Revision history [back]

In such cases it is better to use the complementary error function. First of all note that pari/gp does a pretty good job to compute the integral $$ J = \int_{255.5}^{256.5}\frac{1}{\sqrt{2\pi}}\;\exp -\frac{x^2}2\, dx $$ "as is", without contorsions:

? intnum( x= 255.5, 256.5, exp( -x^2/2 ) / sqrt(2*Pi) )
%18 = 5.8524320023056441353545578190269084186749602193913 E-14179

Same can be used in sage (as a poor man solution), e.g. as

sage: pari( "intnum( x=255.5, 256.5, exp( -x^2/2 ) / sqrt(2*Pi) )" )
5.85243200230564 E-14179

or

sage: a, b = 255.5, 256.5
sage: pari( "intnum( x=%s, %s, exp( -x^2/2 ) / sqrt(2*Pi) )" % (a,b) )
5.85243200230564 E-14179

Recall that the complementary error function is $$ \textrm{erfc}(z) = 1-\textrm{erf}(z) = \frac 2\pi\int_z^\infty\exp(-t^2)\; dt$$ (change of variables, $t=u/sqrt2$...) $$=2\cdot \frac 1{\sqrt{2\pi}}\int_{z/\sqrt 2}^\infty\exp-\frac {u^2}2\; du\ .$$

The better solution is then:

def MyPhi( a,b ):
    return ( erfc(a/sqrt(2)) - erfc(b/sqrt(2)) ) / 2

a, b = 255.5, 256.5
print MyPhi( a,b ).n()

This gives:

5.85243200228551e-14179

For more...

sage: print MyPhi( a,b ).n( digits=100 )
5.852432002305644135354557819026908418674960218659348608849721844978854485293380389274750960850111602e-14179