Revision history [back]

This is not a satisfactory answer... but anyway: The integrals you are trying to compute are expressible in terms of the error function. For some reason that I do not know (and should be investigated), SageMath is unable to correctly compute the error function for large values, when this error function is very close to $1$. Even increasing the precision does not work (that is working with RealField(1000) or more, or with RealIntervalField(1000), etc. Note that the expression of your integral in terms of erf is the best (in some sense) exact answer one can give, since the error function is not elementary.

Actually, the library Arb of Fredrik Johansson is able to perform the computation. The simplest way (to my taste, you may disagree) to use this C library is through the CAS Nemo, written for the language Julia. I give you below an example of use in Nemo (but you can of course write an equivalent C code directly using Arb if you prefer).

The computation I perform is a simplified version of your need:

$$\sum_{c=1}^{d}\int_{256c-0.5}^{256c+0.5}\frac{2}{\sqrt{\pi}}e^{-x^2}dx$$

The code below shows that for $d=50$, the value of the sum is something like $3.1\cdot 10^{-28354}$. The closeness to zero explains why one needs so much precision to be able to distinguish from $0$!

julia> using Nemo

Welcome to Nemo version 0.5.0

Nemo comes with absolutely no warranty whatsoever


julia> function s(d, prec = 1000)       # using default precision 1000 bits
           F = ComplexField(prec)        # For some reason, erf does not work over RealField
           sum = F(0)
           for c = 1:d
               u = F(256*c+.5)
               v = F(256*c-.5)
               sum += erf(u) - erf(v)
           end
           return sum
       end
s (generic function with 2 methods)

julia> s(50) # Not enough precision: get 0 +/- something
[+/- 8.40e-300] + i*0

julia> s(50, 10000) # Still not enough!
[+/- 4.52e-3009] + i*0

julia> s(50, 100000) # Finally!
[3.102262678278053076398652035217490023670965971936356980763793132685150297353942594965955722762805573392153330406480953855731387398023115984622786616977425540305272338247087011441401529284717218316248016664526853584800590168918124704184766482548557105353651652368296375023569553822793725443102971750875855897752226561077609876866162061711562923329377439218548528916354400384468966317197680008340200698258455463251434793207006670012945890581001148258637826042374760859045345143168726986415658022488421604607834339264038625902663546478278814980043653958176479672213608838686103342175667992448852731310595727573672549161581892009023641358983676574601083376637614888507061187244777656151985731856257941976475105038054800870251422661227024669221253196375384625938232368807030347647045465382519932395272845244305245675307763521925808514638809267924250263503220341147139844099847262900841324568603016291331669298541013247881551863448642400104609478241062200239333495268610840375558431182176503376839668955839352645002793527349614716665957443338984442304676833795005909267696213605538509355001746049287714536707616120292584911365481512840133877579846868859068920299888550789717502048145517061355005325552635169752304689718047493686303050935112466379110899771012977020553628991752738579948440664336445330854288254372713613062555003022198940042814571068769124858706877287484803032872382430528486856936500533970245684004943194205532248702178101002288809615625387993450186953306376601302771382662537102238863693773280261194027940865405478664244791943840999854120140399491991893785064480156556621092433860406215704328279792938871986918708868424499764959378290014074093197179867944667683870991035370327938743608633506239623813944138924732150370855135636489875222e-28354 +/- 5.52e-30101] + i*0

This is not a satisfactory answer... but anyway: The integrals you are trying to compute are expressible in terms of the error function. For some reason that I do not know (and should be investigated), SageMath is unable to correctly compute the error function for large values, when this error function is very close to $1$. Even increasing the precision does not work (that is working with RealField(1000) or more, or with RealIntervalField(1000), etc. Note that the expression of your integral in terms of erf is the best (in some sense) exact answer one can give, since the error function is not elementary.

SageMath is able to perform the computation if you allow for enough precision (thanks mforest for testing!), but it takes a very long time. An example of use in SageMath explains why. Below, the result is around $3.1\cdot 10^{-28354}$, very close to zero, so hard to distinguish from zero.

sage: F = RealField(100000)
sage: u = F(256.5)
sage: v = F(255.5)
sage: d = erf(u) - erf(v) 
sage: d
3.10226267827805307639865203[... removed loooots of digits! ...]95e-28354

My previous answer, as I thought SageMath was not able to perform the computation. Yet, it is still interesting since you can have the same result much faster using directly (or though Nemo) the library Arb.

Actually, the library Arb of Fredrik Johansson is able to perform the computation. The simplest way (to my taste, you may disagree) to use this C library is through the CAS Nemo, written for the language Julia. I give you below an example of use in Nemo (but you can of course write an equivalent C code directly using Arb if you prefer).

The computation I perform is a simplified version of your need:

$$\sum_{c=1}^{d}\int_{256c-0.5}^{256c+0.5}\frac{2}{\sqrt{\pi}}e^{-x^2}dx$$

The code below shows that for $d=50$, the value of the sum is something like $3.1\cdot 10^{-28354}$. The closeness to zero explains why one needs so much precision to be able to distinguish from $0$!

julia> using Nemo

Welcome to Nemo version 0.5.0

Nemo comes with absolutely no warranty whatsoever


julia> function s(d, prec = 1000)       # using default precision 1000 bits
           F = ComplexField(prec)        # For some reason, erf does not work over RealField
           sum = F(0)
           for c = 1:d
               u = F(256*c+.5)
               v = F(256*c-.5)
               sum += erf(u) - erf(v)
           end
           return sum
       end
s (generic function with 2 methods)

julia> s(50) # Not enough precision: get 0 +/- something
[+/- 8.40e-300] + i*0

julia> s(50, 10000) # Still not enough!
[+/- 4.52e-3009] + i*0

julia> s(50, 100000) # Finally!
[3.102262678278053076398652035217490023670965971936356980763793132685150297353942594965955722762805573392153330406480953855731387398023115984622786616977425540305272338247087011441401529284717218316248016664526853584800590168918124704184766482548557105353651652368296375023569553822793725443102971750875855897752226561077609876866162061711562923329377439218548528916354400384468966317197680008340200698258455463251434793207006670012945890581001148258637826042374760859045345143168726986415658022488421604607834339264038625902663546478278814980043653958176479672213608838686103342175667992448852731310595727573672549161581892009023641358983676574601083376637614888507061187244777656151985731856257941976475105038054800870251422661227024669221253196375384625938232368807030347647045465382519932395272845244305245675307763521925808514638809267924250263503220341147139844099847262900841324568603016291331669298541013247881551863448642400104609478241062200239333495268610840375558431182176503376839668955839352645002793527349614716665957443338984442304676833795005909267696213605538509355001746049287714536707616120292584911365481512840133877579846868859068920299888550789717502048145517061355005325552635169752304689718047493686303050935112466379110899771012977020553628991752738579948440664336445330854288254372713613062555003022198940042814571068769124858706877287484803032872382430528486856936500533970245684004943194205532248702178101002288809615625387993450186953306376601302771382662537102238863693773280261194027940865405478664244791943840999854120140399491991893785064480156556621092433860406215704328279792938871986918708868424499764959378290014074093197179867944667683870991035370327938743608633506239623813944138924732150370855135636489875222e-28354 +/- 5.52e-30101] + i*0