integral competition HP50g vs. DM42

[Albert Chan]

(08-23-2020 08:56 AM)peacecalc Wrote:  @Albert Chan, do you know how this is implemented in XCAS? So, I don't have to invent the wheel a second time and maybe it is possible for me to translate it in code for the HP50g.

I just downloaded giac source, usual.cc, function erf0():

Quote:if (absx<=3){
      // numerical computation of int(exp(-t^2),t=0..x)
      // by series expansion at x=0
      // x*sum( (-1)^n*(x^2)^n/n!/(2*n+1),n=0..inf)

Code:

def hyp1f1(a,b,x):
    s = t = c = 1
    while 1:
        t *= a/(b*c)*x
        if s == s+t: return s
        a, b, c, s = a+1, b+1, c+1, s+t

>>> k = 2/pi**0.5
>>> erf = lambda x: k*x * hyp1f1(1/2, 3/2, -x*x)
>>> CiS = lambda x, erf=erf: (1+1j)/2 * erf((1-1j)/k * x)
>>> CiS(3.9)
(0.42233270272270274 + 0.47520241152278453j)

CiS(3.9) called erf(z = (1-1j) * 3.4562850092657564), |z| ≈ 5 > 3
As expected, we lose many digits precision due to massive cancellation.

For big |z|, XCas use continued fraction form:

Quote:// 2*exp(z^2)*int(exp(-t^2),t=z..inf) = 1/(z+1/2/(z+1/(z+3/2/(z+...))))

int(exp(-t^2),t=z..inf) = sqrt(pi)/2 * erfc(z), we got this:

Code:

def erfc(x, n=20, res=0):   # continued fraction, bottom up
    while n>0: res = n/(x+res); n -= 0.5
    return 0.5*k * exp(-x*x) / (x+res)

>>> CiS(3.9, erf = lambda x: 1 - erfc(x))
(0.42233271026093344 + 0.47520240235068834j)

Or, we could calculate generalized continued fraction top down, without hard-coded maximum n.

Code:

def erfc2(x, eps=1/16):     # continued fraction, top down
    a1, b1, a2, b2 = 1, 0, x, 1
    f1, f2, n = 0, 1/x, 0.5
    while f2 != f2 + (f1-f2)*eps:
        a1, b1, a2, b2 = a2, b2, n*a1+x*a2, n*b1+x*b2
        f1, f2, n = f2, b2/a2, n+0.5
    return 0.5*k * exp(-x*x) * f2

>>> CiS(3.9, erf = lambda x: 1 - erfc2(x))
(0.42233271026093344 + 0.47520240235068834j)