Numerical integration

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
Post Reply
jdebord
Posts: 551
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Numerical integration

Post by jdebord »

Here is my FreeBASIC version of an integration program in C by Robert van Engelen. It is posted here with permission of the original author.

Code: Select all

' ************************************************************************
' Numerical integration by Tanh-Sinh, Sinh-Sinh and Exp-Sinh formulas
' ************************************************************************
' Ref. : R. A. van Engelen, https://www.genivia.com/files/qthsh.pdf
' ************************************************************************

#include "crt.bi"

#define isfinite(x) (abs(x) <> INFINITY)

function exp_sinh_opt_d(func as function (x as double) as double, _
                        a as double, eps as double, d as double) as double
  
  dim as double fl, fr, h, h2, lfl, lfr, lr, r, s

  dim as long ev, i, j
  
  ev = 2 : i = 1 : j = 32  ' j=32 is optimal to search for r
  
  h2 = func(a + d / 2) - func(a + d * 2) * 4
  
  if (isfinite(h2) and abs(h2) > 1e-5) then  ' if |h2| > 2^-16
    s = 0
    lr = 2
    
    do
      ' find max j such that fl and fr are finite
      j /= 2
      r = 1 shl (i + j)
      fl = func(a + d / r)
      fr = func(a + d * r) * r * r
      ev += 2
      h = fl - fr
    loop while (j > 1 and not(isfinite(h)))
    
    if (j > 1 and isfinite(h) and sgn(h) <> sgn(h2)) then
      lfl = fl          ' last fl=func(a+d/r)
      lfr = fr          ' last fr=func(a+d*r)*r*r
    
      do                ' bisect in 4 iterations
        j /= 2
        r = 1 shl (i + j)
        fl = func(a + d / r)
        fr = func(a + d * r) * r * r
        ev += 2
        h = fl - fr
        
        if (isfinite(h)) then
          s += abs(h)   ' sum |h| to remove noisy cases
          if (sgn(h) = sgn(h2)) then
            i += j      ' search right half
          else          ' search left half
            lfl = fl    ' record last fl=f(a+d/r)
            lfr = fr    ' record last fr=f(a+d*r)*r*r
            lr = r      ' record last r
          end if
        end if
      loop while (j > 1)
      
      if (s > eps) then ' if sum of |h| > eps
        h = lfl - lfr   ' use last fl and fr before the sign change
        r = lr          ' use last r before the sign change
        
        ' if last difference was nonzero, back up r by one step
        if (h <> 0) then r /= 2   
          
        if (abs(lfl) < abs(lfr)) then
          d /= r        ' move d closer to the finite endpoint
        else
          d *= r        ' move d closer to the infinite endpoint
        end if  
      end if
    end if
  end if
  
  return d
end function

   
function quad (func as function(x as double) as double, a as double, b as double, _
                n as long = 6, eps as double = 1.0e-9, byref err_est as double = 0) as double

  dim as double c, d, eh, fp, fm, h, p, q, r, s, sign, t, tol, u, v, w, x, y 

  dim as long k
  dim as long mode  ' Tanh-Sinh = 0, Exp-Sinh = 1, Sinh-Sinh = 2
  
  d = 1 : sign = 1 : h = 2 : tol = 10 * eps
  if (b < a) then swap a, b : sign = -1
  
  if (isfinite(a) and isfinite(b)) then
    c = (a + b) / 2
    d = (b - a) / 2
    v = c
  elseif (isfinite(a)) then
    mode = 1                             ' Exp-Sinh
    d = exp_sinh_opt_d(func, a, eps, d)
    c = a
    v = a + d
  elseif (isfinite(b)) then
    mode = 1                             ' Exp-Sinh
    d = exp_sinh_opt_d(func, b, eps, -d)
    sign = -sign
    c = b
    v = b + d
  else
    mode = 2                             ' Sinh-Sinh
    v = 0
  end if
  
  s = func(v)
  
  do 
    p = 0 : fp = 0 : fm = 0
    h /= 2
    eh = exp(h)
    t = eh
    if (k > 0) then eh *= eh
    
    if (mode = 0) then                   ' Tanh-Sinh
        
      do 
        u = exp(1 / t - t)               ' = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2
        r = 2 * u / (1 + u)              ' = 1 - tanh(sinh(j*h))
        w = (t + 1 / t) * r / (1 + u)    ' = cosh(j*h)/cosh(sinh(j*h))^2
        x = d * r
        
        if (a + x > a) then              ' if too close to a then reuse previous fp
          y = func(a + x)
          if (isfinite(y)) then fp = y   ' if f(x) is finite, add to local sum
        end if
        
        if (b - x < b) then              ' if too close to b then reuse previous fm
          y = func(b - x)
          if (isfinite(y)) then fm = y   ' if f(x) is finite, add to local sum
        end if
        
        q = w * (fp + fm)
        p += q
        t *= eh
      loop while (abs(q) > eps * abs(p))

    else 
        
      t /= 2
      do 
        r = exp(t - 0.25 / t)            ' = exp(sinh(j*h))
        w = r
        q = 0
        
        if (mode = 1) then               ' Exp-Sinh
          x = c + d / r
          if (x = c) then exit do        ' if x hit the finite endpoint then break
          y = func(x)
          if isfinite(y) then q += y / w ' if f(x) is finite, add to local sum
        else                             ' Sinh-Sinh
          r = (r - 1 / r) / 2            ' = sinh(sinh(j*h))
          w = (w + 1 / w) / 2            ' = cosh(sinh(j*h))
          x = c - d * r
          y = func(x) 
          if isfinite(y) then q += y * w ' if f(x) is finite, add to local sum
        end if
        
        x = c + d * r
        y = func(x)
        if isfinite(y) then q += y * w   ' if f(x) is finite, add to local sum
        q *= t + 0.25 / t                ' q *= cosh(j*h)
        p += q
        t *= eh
      loop while (abs(q) > eps * abs(p))
      
    end if

    v = s - p
    s += p
    k += 1
  
  loop while (abs(v) > tol * abs(s) and k <= n)

  ' if the estimated relative error is desired, then return it
  if err_est <> 0 then err_est = abs(v) / (abs(s) + eps)
  
  ' result with estimated relative error err_est
  return sign * d * s * h
end function

' ************************************************************************

' example functions to integrate, exact integrals are 1, 5 and 2

function f1(x as double) as double
  return acos(x)
end function

function f2(x as double) as double
  return exp(- x / 5)
end function

function f3(x as double) as double
  return 1 / cosh(x)^2
end function

printf(!"integrate(acos(x), x=0..1)           = %.15g\n", quad(@f1, 0, 1))
printf(!"integrate(exp(-x/5), x=0..+inf)      = %.15g\n", quad(@f2, 0, INFINITY))
printf(!"integrate(1/cosh(x)^2, x=-inf..+inf) = %.15g\n", quad(@f3, -INFINITY, INFINITY))

sleep
hhr
Posts: 222
Joined: Nov 29, 2019 10:41

Re: Numerical integration

Post by hhr »

Very nice, is there a way to estimate the error?

Code: Select all

function f4(x as double) as double
  return 1-sqr(x*x-1)/x
end function

printf(!"integrate(1-sqr(x*x-1)/x, x=1..+inf) = %.15g\n", quad(@f4, 1, INFINITY))
printf(!"should be                   (pi-2)/2 = %.15g\n", M_PI/2 - 1)

Code: Select all

dim shared as double n

function f5(x as double) as double
  return (x^n)*exp(-x)
end function

for n = 0 to 4
   print "n =";n
   print "n! =  ";tgamma(n+1)
   print "quad =";quad(@f5, 0, INFINITY)
   print string(25, "-")
next
jdebord
Posts: 551
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Re: Numerical integration

Post by jdebord »

The quad function has a sixth parameter err_est with a default value of 0. To return it, initialize it to a different value.

Code: Select all

dim as double err_est = 1

function f4(x as double) as double
  return 1-sqr(x*x-1)/x
end function

printf(!"integrate(1-sqr(x*x-1)/x, x=1..+inf) = %.15g\n", quad(@f4, 1, INFINITY, , , err_est))
printf(!"should be                   (pi-2)/2 = %.15g\n", M_PI/2 - 1)

print "estimated error = ", err_est      
Result:

Code: Select all

integrate(1-sqr(x*x-1)/x, x=1..+inf) = 0.570796326697611
should be                   (pi-2)/2 = 0.570796326794897
estimated error =            2.962455091046144e-011
dodicat
Posts: 7998
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Numerical integration

Post by dodicat »

Using Simpson's rules (as in shipyards for waterplane areas) gets quite close.
...
COMPARE:

Code: Select all




' ************************************************************************
' Numerical integration by Tanh-Sinh, Sinh-Sinh and Exp-Sinh formulas
' ************************************************************************
' Ref. : R. A. van Engelen, https://www.genivia.com/files/qthsh.pdf
' ************************************************************************

#include "crt.bi"

#define isfinite(x) (abs(x) <> INFINITY)

function exp_sinh_opt_d(func as function (x as double) as double, _
                        a as double, eps as double, d as double) as double
  
  dim as double fl, fr, h, h2, lfl, lfr, lr, r, s

  dim as long ev, i, j
  
  ev = 2 : i = 1 : j = 32  ' j=32 is optimal to search for r
  
  h2 = func(a + d / 2) - func(a + d * 2) * 4
  
  if (isfinite(h2) and abs(h2) > 1e-5) then  ' if |h2| > 2^-16
    s = 0
    lr = 2
    
    do
      ' find max j such that fl and fr are finite
      j /= 2
      r = 1 shl (i + j)
      fl = func(a + d / r)
      fr = func(a + d * r) * r * r
      ev += 2
      h = fl - fr
    loop while (j > 1 and not(isfinite(h)))
    
    if (j > 1 and isfinite(h) and sgn(h) <> sgn(h2)) then
      lfl = fl          ' last fl=func(a+d/r)
      lfr = fr          ' last fr=func(a+d*r)*r*r
    
      do                ' bisect in 4 iterations
        j /= 2
        r = 1 shl (i + j)
        fl = func(a + d / r)
        fr = func(a + d * r) * r * r
        ev += 2
        h = fl - fr
        
        if (isfinite(h)) then
          s += abs(h)   ' sum |h| to remove noisy cases
          if (sgn(h) = sgn(h2)) then
            i += j      ' search right half
          else          ' search left half
            lfl = fl    ' record last fl=f(a+d/r)
            lfr = fr    ' record last fr=f(a+d*r)*r*r
            lr = r      ' record last r
          end if
        end if
      loop while (j > 1)
      
      if (s > eps) then ' if sum of |h| > eps
        h = lfl - lfr   ' use last fl and fr before the sign change
        r = lr          ' use last r before the sign change
        
        ' if last difference was nonzero, back up r by one step
        if (h <> 0) then r /= 2   
          
        if (abs(lfl) < abs(lfr)) then
          d /= r        ' move d closer to the finite endpoint
        else
          d *= r        ' move d closer to the infinite endpoint
        end if  
      end if
    end if
  end if
  
  return d
end function

   
function quad (func as function(x as double) as double, a as double, b as double, _
                n as long = 6, eps as double = 1.0e-9, byref err_est as double = 0) as double

  dim as double c, d, eh, fp, fm, h, p, q, r, s, sign, t, tol, u, v, w, x, y 

  dim as long k
  dim as long mode  ' Tanh-Sinh = 0, Exp-Sinh = 1, Sinh-Sinh = 2
  
  d = 1 : sign = 1 : h = 2 : tol = 10 * eps
  if (b < a) then swap a, b : sign = -1
  
  if (isfinite(a) and isfinite(b)) then
    c = (a + b) / 2
    d = (b - a) / 2
    v = c
  elseif (isfinite(a)) then
    mode = 1                             ' Exp-Sinh
    d = exp_sinh_opt_d(func, a, eps, d)
    c = a
    v = a + d
  elseif (isfinite(b)) then
    mode = 1                             ' Exp-Sinh
    d = exp_sinh_opt_d(func, b, eps, -d)
    sign = -sign
    c = b
    v = b + d
  else
    mode = 2                             ' Sinh-Sinh
    v = 0
  end if
  
  s = func(v)
  
  do 
    p = 0 : fp = 0 : fm = 0
    h /= 2
    eh = exp(h)
    t = eh
    if (k > 0) then eh *= eh
    
    if (mode = 0) then                   ' Tanh-Sinh
        
      do 
        u = exp(1 / t - t)               ' = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2
        r = 2 * u / (1 + u)              ' = 1 - tanh(sinh(j*h))
        w = (t + 1 / t) * r / (1 + u)    ' = cosh(j*h)/cosh(sinh(j*h))^2
        x = d * r
        
        if (a + x > a) then              ' if too close to a then reuse previous fp
          y = func(a + x)
          if (isfinite(y)) then fp = y   ' if f(x) is finite, add to local sum
        end if
        
        if (b - x < b) then              ' if too close to b then reuse previous fm
          y = func(b - x)
          if (isfinite(y)) then fm = y   ' if f(x) is finite, add to local sum
        end if
        
        q = w * (fp + fm)
        p += q
        t *= eh
      loop while (abs(q) > eps * abs(p))

    else 
        
      t /= 2
      do 
        r = exp(t - 0.25 / t)            ' = exp(sinh(j*h))
        w = r
        q = 0
        
        if (mode = 1) then               ' Exp-Sinh
          x = c + d / r
          if (x = c) then exit do        ' if x hit the finite endpoint then break
          y = func(x)
          if isfinite(y) then q += y / w ' if f(x) is finite, add to local sum
        else                             ' Sinh-Sinh
          r = (r - 1 / r) / 2            ' = sinh(sinh(j*h))
          w = (w + 1 / w) / 2            ' = cosh(sinh(j*h))
          x = c - d * r
          y = func(x) 
          if isfinite(y) then q += y * w ' if f(x) is finite, add to local sum
        end if
        
        x = c + d * r
        y = func(x)
        if isfinite(y) then q += y * w   ' if f(x) is finite, add to local sum
        q *= t + 0.25 / t                ' q *= cosh(j*h)
        p += q
        t *= eh
      loop while (abs(q) > eps * abs(p))
      
    end if

    v = s - p
    s += p
    k += 1
  
  loop while (abs(v) > tol * abs(s) and k <= n)

  ' if the estimated relative error is desired, then return it
  if err_est <> 0 then err_est = abs(v) / (abs(s) + eps)
  
  ' result with estimated relative error err_est
  return sign * d * s * h
end function

' ************************************************************************

' example functions to integrate, exact integrals are 1, 5 and 2

function f1(x as double) as double
  return acos(x)
end function

function f2(x as double) as double
  return exp(- x / 5)
end function

function f3(x as double) as double
  return 1 / cosh(x)^2
end function

function f4(x as double) as double
  return 1-sqr(x*x-1)/x
end function

function f5(x as double) as double
  return x^2
end function

function f6(x as double) as double
    return 1/(1+exp(x))
    end function


'=================================================
#macro SimpsonArea(fn,Lo,Up,ordinates,sum)
scope
    dim as double x,part1,part2
    dim as double l=lo,u=up,n=iif(ordinates mod 2=0,ordinates+1,ordinates)
    var z=(len(str(n))-1)/2
     if abs(l)=infinity then l=sgn(l)*10^4
     if abs(u)=infinity then u=sgn(u)*10^4
    Var h=(U-L)/n
    x=L+.5*h:Part1=fn
    For i As long=1 To n-1
        x=L+h*i+h/2
        Part1+=fn
        x=L+h*i
        Part2+=fn
    Next i
    x=L:var fn1=fn
    x=U:var fn2=fn
    sum= (h/6)*(fn1+fn2+Part1*4+Part2*2) 
    end scope
#endmacro
'==================================================




printf(!"integrate(acos(x), x=0..1)           = %.15g\n", quad(@f1, 0, 1))
printf(!"integrate(exp(-x/5), x=0..+inf)      = %.15g\n", quad(@f2, 0, INFINITY))
printf(!"integrate(1/cosh(x)^2, x=-inf..+inf) = %.15g\n", quad(@f3, -INFINITY, INFINITY))
printf(!"integrate(1-sqr(x*x-1)/x, x=1..+inf) = %.15g\n", quad(@f4, 1, INFINITY))
printf(!"integrate(1/(1+exp(x)), x=0..+inf)   = %.15g\n", quad(@f6, 0, INFINITY))
printf(!"integrate(x^2, x=-100.005..12.556)   = %.15g\n",quad(@f5, -100.005, 120.556))

print
'Simpson macro
dim as double Result
const ordinates=1000000 ' no less anyway.
SimpsonArea(acos(x),0,1,ordinates,Result)
printf(!"integrate(acos(x), x=0..1)           = %.15g\n", Result)

SimpsonArea(exp(-x/5),0,infinity,ordinates,Result)
printf(!"integrate(exp(-x/5), x=0..+inf)      = %.15g\n", Result)

SimpsonArea(1/cosh(x)^2,-infinity,infinity,ordinates,Result)
printf(!"integrate(1/cosh(x)^2, x=-inf..+inf) = %.15g\n", Result)

SimpsonArea(1-sqr(x*x-1)/x,1,infinity,ordinates,Result)
printf(!"integrate(1-sqr(x*x-1)/x, x=1..+inf) = %.15g\n", Result)

SimpsonArea(1/(1+exp(x)),0,infinity,ordinates,Result)
printf(!"integrate(1/(1+exp(x)), x=0..+inf)   = %.15g\n", Result)

SimpsonArea(x^2,-100.005,120.556,ordinates,Result)
printf(!"integrate(x^2, x=-100.005..120.556)  = %.15g\n", Result)
print
print "actual value for last integral (x^2)  ";120.556^3/3- -100.005^3/3
print



sleep
 
hhr
Posts: 222
Joined: Nov 29, 2019 10:41

Re: Numerical integration

Post by hhr »

First of all, thank you for the beautiful program.

The following example works for n = 0 to 97 and then again for n = 168 to 170.
For n > 170, FreeBASIC-Double can no longer display the result.

Code: Select all

dim shared as double n
dim as double q, g, err_est

function f5(x as double) as double
   return (x^n)*exp(-x)
end function

for n = 0 to 180
   print "n ="; n
   g = tgamma(n+1)
   print "n! =  "; g
   err_est = 1
   q = quad(@f5, 0, INFINITY, , , err_est)
   print "quad ="; q; " "; csng(err_est); " "; csng((g - q) / g)
   print string(25, "-")
   getkey
next
If INFINITY is replaced by 1000, it works better:

Code: Select all

q = quad(@f5, 0, 1000, , , err_est)
I am surprised that the values for INFINITY and n = 98 to 167 are too small. Is there a bug with INFINITY?

I have experimented with gas32.
With gcc and gas64 it behaves differently and with gas64 it crashes with b = 1000 at n = 103.

err_est is very useful.
SARG
Posts: 1790
Joined: May 27, 2005 7:15
Location: FRANCE

Re: Numerical integration

Post by SARG »

hhr wrote: Jun 20, 2024 19:30 with gas64 it crashes with b = 1000 at n = 103.
I'll fix the problem in gas64 but in your example what is tgamma ?
hhr
Posts: 222
Joined: Nov 29, 2019 10:41

Re: Numerical integration

Post by hhr »

tgamma is the gamma function in crt, which is declared in crt\math.bi.

https://en.wikipedia.org/wiki/Gamma_function

In the hope that I can be of some help, I suspect that the program is stuck in the loop that ends with line 141:

Code: Select all

loop while (abs(q) > eps * abs(p))
If I switch to any other data type, it works, for example:

Code: Select all

loop while cbyte(abs(q) > eps * abs(p))
or

Code: Select all

loop while cdbl(abs(q) > eps * abs(p))
SARG
Posts: 1790
Joined: May 27, 2005 7:15
Location: FRANCE

Re: Numerical integration

Post by SARG »

ok, thks.

I don't get a crash with b=1000. Using very last fbc 1.20 on Windows.

Code: Select all

n = 102
n! =   9.614466715035127e+161
quad = 9.614466715023597e+161  1.321906e-012  1.199234e-012
-------------------------
n = 103
n! =   9.902900716486181e+163
quad =-1.#IND  1.#QNAN -1.#IND
-------------------------
n = 104
n! =   1.029901674514563e+166
quad =-1.#IND  1.#QNAN -1.#IND
-------------------------
hhr
Posts: 222
Joined: Nov 29, 2019 10:41

Re: Numerical integration

Post by hhr »

Indeed, it works with fbc 1.20.0, I had used fbc 1.10.1.

Many thanks for gas64 and the effort.
hhr
Posts: 222
Joined: Nov 29, 2019 10:41

Re: Numerical integration

Post by hhr »

With the latest fbc 1.20.0, gcc32, gcc64 and gas64 work exactly the same except for slight differences in the calculated numbers.
gas32 has a different behavior.

What distinguishes gas32 from the other compilers?
dodicat
Posts: 7998
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Numerical integration

Post by dodicat »

32 bit -gen gas works out better for these areas.
I have added simpson's shipyard area (with 500000 iterations -which is more than any ship would get) and a string factorial. (for fun)
infinity at 10^4 for simpson
Choose your compiler option at line 1.

Code: Select all

'#cmdline "-gen gas64"
#include "crt.bi"
#define isfinite(x) (abs(x) <> INFINITY)
dim shared as double n
dim as double q, g, err_est


function exp_sinh_opt_d(func as function (x as double) as double, _
                        a as double, eps as double, d as double) as double
  
  dim as double fl, fr, h, h2, lfl, lfr, lr, r, s

  dim as long ev, i, j
  
  ev = 2 : i = 1 : j = 32  ' j=32 is optimal to search for r
  
  h2 = func(a + d / 2) - func(a + d * 2) * 4
  
  if (isfinite(h2) and abs(h2) > 1e-5) then  ' if |h2| > 2^-16
    s = 0
    lr = 2
    
    do
      ' find max j such that fl and fr are finite
      j /= 2
      r = 1 shl (i + j)
      fl = func(a + d / r)
      fr = func(a + d * r) * r * r
      ev += 2
      h = fl - fr
    loop while (j > 1 and not(isfinite(h)))
    
    if (j > 1 and isfinite(h) and sgn(h) <> sgn(h2)) then
      lfl = fl          ' last fl=func(a+d/r)
      lfr = fr          ' last fr=func(a+d*r)*r*r
    
      do                ' bisect in 4 iterations
        j /= 2
        r = 1 shl (i + j)
        fl = func(a + d / r)
        fr = func(a + d * r) * r * r
        ev += 2
        h = fl - fr
        
        if (isfinite(h)) then
          s += abs(h)   ' sum |h| to remove noisy cases
          if (sgn(h) = sgn(h2)) then
            i += j      ' search right half
          else          ' search left half
            lfl = fl    ' record last fl=f(a+d/r)
            lfr = fr    ' record last fr=f(a+d*r)*r*r
            lr = r      ' record last r
          end if
        end if
      loop while (j > 1)
      
      if (s > eps) then ' if sum of |h| > eps
        h = lfl - lfr   ' use last fl and fr before the sign change
        r = lr          ' use last r before the sign change
        
        ' if last difference was nonzero, back up r by one step
        if (h <> 0) then r /= 2   
          
        if (abs(lfl) < abs(lfr)) then
          d /= r        ' move d closer to the finite endpoint
        else
          d *= r        ' move d closer to the infinite endpoint
        end if  
      end if
    end if
  end if
  
  return d
end function

   

function quad (func as function(x as double) as double, a as double, b as double, _
                n as long = 6, eps as double = 1.0e-9, byref err_est as double = 0) as double

  dim as double c, d, eh, fp, fm, h, p, q, r, s, sign, t, tol, u, v, w, x, y 

  dim as long k
  dim as long mode  ' Tanh-Sinh = 0, Exp-Sinh = 1, Sinh-Sinh = 2
  
  d = 1 : sign = 1 : h = 2 : tol = 10 * eps
  if (b < a) then swap a, b : sign = -1
  
  if (isfinite(a) and isfinite(b)) then
    c = (a + b) / 2
    d = (b - a) / 2
    v = c
  elseif (isfinite(a)) then
    mode = 1                             ' Exp-Sinh
    d = exp_sinh_opt_d(func, a, eps, d)
    c = a
    v = a + d
  elseif (isfinite(b)) then
    mode = 1                             ' Exp-Sinh
    d = exp_sinh_opt_d(func, b, eps, -d)
    sign = -sign
    c = b
    v = b + d
  else
    mode = 2                             ' Sinh-Sinh
    v = 0
  end if
  
  s = func(v)
  
  do 
    p = 0 : fp = 0 : fm = 0
    h /= 2
    eh = exp(h)
    t = eh
    if (k > 0) then eh *= eh
    
    if (mode = 0) then                   ' Tanh-Sinh
        
      do 
        u = exp(1 / t - t)               ' = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2
        r = 2 * u / (1 + u)              ' = 1 - tanh(sinh(j*h))
        w = (t + 1 / t) * r / (1 + u)    ' = cosh(j*h)/cosh(sinh(j*h))^2
        x = d * r
        
        if (a + x > a) then              ' if too close to a then reuse previous fp
          y = func(a + x)
          if (isfinite(y)) then fp = y   ' if f(x) is finite, add to local sum
        end if
        
        if (b - x < b) then              ' if too close to b then reuse previous fm
          y = func(b - x)
          if (isfinite(y)) then fm = y   ' if f(x) is finite, add to local sum
        end if
        
        q = w * (fp + fm)
        p += q
        t *= eh
      loop while (abs(q) > eps * abs(p))

    else 
        
      t /= 2
      do 
        r = exp(t - 0.25 / t)            ' = exp(sinh(j*h))
        w = r
        q = 0
        
        if (mode = 1) then               ' Exp-Sinh
          x = c + d / r
          if (x = c) then exit do        ' if x hit the finite endpoint then break
          y = func(x)
          if isfinite(y) then q += y / w ' if f(x) is finite, add to local sum
        else                             ' Sinh-Sinh
          r = (r - 1 / r) / 2            ' = sinh(sinh(j*h))
          w = (w + 1 / w) / 2            ' = cosh(sinh(j*h))
          x = c - d * r
          y = func(x) 
          if isfinite(y) then q += y * w ' if f(x) is finite, add to local sum
        end if
        
        x = c + d * r
        y = func(x)
        if isfinite(y) then q += y * w   ' if f(x) is finite, add to local sum
        q *= t + 0.25 / t                ' q *= cosh(j*h)
        p += q
        t *= eh
      loop while (abs(q) > eps * abs(p))
      
    end if

    v = s - p
    s += p
    k += 1
  
  loop while (abs(v) > tol * abs(s) and k <= n)

  ' if the estimated relative error is desired, then return it
  if err_est <> 0 then err_est = abs(v) / (abs(s) + eps)
  
  ' result with estimated relative error err_est
  return sign * d * s * h
end function




Function factorial(num As String) As String 
    static  As Ubyte _Mod(0 To 99),_Div(0 To 99),flag=0
    if flag=0 then
    For z As Integer=0 To 99:_Mod(z)=(z Mod 10+48):_Div(z)=z\10:next
    flag=1
    end if
    Var fact="1",a="",b="",c=""
    Dim As Ubyte n,carry,ai
    For z As Integer=1 To Valint(num)
        a=fact:b=Str(z):Var la =Len(a),lb =Len(b)
        c =String(la+lb,"0")
        For i As Integer =la-1 To 0 Step -1
            carry=0:ai=a[i]-48
            For j As Integer =lb-1 To 0 Step -1
                n =ai*(b[j]-48)+(c[i+j+1]-48)+carry
                carry =_Div(n):c[i+j+1]=_Mod(n)
            Next j
            c[i]+=carry
        Next i
        fact=Ltrim(c,"0")
    Next z
    Return fact
End Function


Function Simpson(fn As Function(x As Double) As Double,L As Double,U As Double,Ordinates As Integer) As Double
    if abs(l)=infinity then l=sgn(l)*10^4
    if abs(u)=infinity then u=sgn(u)*10^4
    var n=Ordinates
    If n Mod 2=0 Then n=n+1 
    Var h=(U-L)/n
    Var Part1=fn(L+.5*h)
    Var Part2=0.0
    For i As Integer=1 To n-1
        Part1+=fn(L+h*i+h/2)
        Part2+=fn(L+h*i)
    Next i
    Function= (h/6)*(fn(L)+fn(U)+Part1*4+Part2*2)  
End Function



function f5(x as double) as double
   return (x^n)*exp(-x)
end function

for n = 0 to 180
   print "n = "; n
   g = tgamma(n+1)
   var f=factorial(str(n))
   print "n! =   "; g;" f = ";f
   err_est = 1
   q = quad(@f5, 0, INFINITY, , , err_est)
   print "quad = "; q; " "; csng(err_est); " "; csng((g - q) / g)
   var area=simpson(@f5,0,infinity,500000)
   print "simpson";area;" "; " err not given"; " "; csng((val(f) - area) / val(f))
   print string(25, "-")
   getkey
next 
hhr
Posts: 222
Joined: Nov 29, 2019 10:41

Re: Numerical integration

Post by hhr »

Simpson works well, it's just a bit slow.
I also find the factorial function interesting. I didn't realize that you can easily convert the string to Single or Double.

Code: Select all

#include "crt.bi"

Function factorial(num As String) As String 'by dodicat
   Static  As Ubyte _Mod(0 To 99),_Div(0 To 99),flag=0
   If flag=0 Then
      For z As Integer=0 To 99
         _Mod(z)=(z Mod 10+48)
         _Div(z)=z\10
      Next
      flag=1
   End If
   Var fact="1",a="",b="",c=""
   Dim As Ubyte n,carry,ai
   For z As Integer=1 To Valint(num)
      a=fact:b=Str(z):Var la =Len(a),lb =Len(b)
      c =String(la+lb,"0")
      For i As Integer =la-1 To 0 Step -1
         carry=0:ai=a[i]-48
         For j As Integer =lb-1 To 0 Step -1
            n =ai*(b[j]-48)+(c[i+j+1]-48)+carry
            carry =_Div(n):c[i+j+1]=_Mod(n)
         Next j
         c[i]+=carry
      Next i
      fact=Ltrim(c,"0")
   Next z
   Return fact
End Function

Dim As Long n
Dim As String f

For n=0 To 200
   Print "n =";n
   f = factorial(Str(n))
   Print "f = ";f
   Print "cdbl(f) ="; Cast(Double,f)
   Print "tgamma = ";tgamma(n+1)
   Print String(20,"-")
   Getkey
Next

Sleep
hhr
Posts: 222
Joined: Nov 29, 2019 10:41

Re: Numerical integration

Post by hhr »

I wanted to replace Double with another data type and while searching I noticed this topic:

https://www.freebasic.net/forum/viewtopic.php?t=27314

I copied the bi files from the first three posts and rebuilt the example from jdebord's first post.
This gives me the same results for all four compilers.
This even works with Windows XP, but unfortunately not with Linux.

However, you have to refrain from including crt.bi.

Code: Select all

' ************************************************************************
' Numerical integration by Tanh-Sinh, Sinh-Sinh and Exp-Sinh formulas
' ************************************************************************
' Ref. : R. A. van Engelen, https://www.genivia.com/files/qthsh.pdf
' ************************************************************************

'#cmdline "-gen gcc"

Print __fb_signature__;" ";
Print __fb_backend__;" ";
#IFDEF __FB_64BIT__
  Print "64Bit"
#ELSE
  Print "32Bit"
#ENDIF
Print string(26,"-")

#include "clongdouble.bi"

#define INFINITY (clongdouble(1) / clongdouble(0))

#define isfinite(x) (abs(x) <> INFINITY)

function exp_sinh_opt_d(func as function (x as clongdouble) as clongdouble, _
                        a as clongdouble, eps as clongdouble, d as clongdouble) as clongdouble
  
  dim as clongdouble fl, fr, h, h2, lfl, lfr, lr, r, s

  dim as long ev, i, j
  
  ev = 2 : i = 1 : j = 32  ' j=32 is optimal to search for r
  
  h2 = func(a + d / 2) - func(a + d * 2) * 4
  
  if (isfinite(h2) and abs(h2) > 1e-5) then  ' if |h2| > 2^-16
    s = 0
    lr = 2
    
    do
      ' find max j such that fl and fr are finite
      j /= 2
      r = 1 shl (i + j)
      fl = func(a + d / r)
      fr = func(a + d * r) * r * r
      ev += 2
      h = fl - fr
    loop while (j > 1 and not(isfinite(h)))
    
    if (j > 1 and isfinite(h) and sgn(h) <> sgn(h2)) then
      lfl = fl          ' last fl=func(a+d/r)
      lfr = fr          ' last fr=func(a+d*r)*r*r
    
      do                ' bisect in 4 iterations
        j /= 2
        r = 1 shl (i + j)
        fl = func(a + d / r)
        fr = func(a + d * r) * r * r
        ev += 2
        h = fl - fr
        
        if (isfinite(h)) then
          s += abs(h)   ' sum |h| to remove noisy cases
          if (sgn(h) = sgn(h2)) then
            i += j      ' search right half
          else          ' search left half
            lfl = fl    ' record last fl=f(a+d/r)
            lfr = fr    ' record last fr=f(a+d*r)*r*r
            lr = r      ' record last r
          end if
        end if
      loop while (j > 1)
      
      if (s > eps) then ' if sum of |h| > eps
        h = lfl - lfr   ' use last fl and fr before the sign change
        r = lr          ' use last r before the sign change
        
        ' if last difference was nonzero, back up r by one step
        if (h <> 0) then r /= 2   
          
        if (abs(lfl) < abs(lfr)) then
          d /= r        ' move d closer to the finite endpoint
        else
          d *= r        ' move d closer to the infinite endpoint
        end if  
      end if
    end if
  end if
  
  return d
end function

   
function quad (func as function(x as clongdouble) as clongdouble, a as clongdouble, b as clongdouble, _
                n as long = 6, eps as clongdouble = 1.0e-9, byref err_est as clongdouble = 0) as clongdouble

  dim as clongdouble c, d, eh, fp, fm, h, p, q, r, s, sign, t, tol, u, v, w, x, y 

  dim as long k
  dim as long mode  ' Tanh-Sinh = 0, Exp-Sinh = 1, Sinh-Sinh = 2
  
  d = 1 : sign = 1 : h = 2 : tol = 10 * eps
  if (b < a) then swap a, b : sign = -1
  
  if (isfinite(a) and isfinite(b)) then
    c = (a + b) / 2
    d = (b - a) / 2
    v = c
  elseif (isfinite(a)) then
    mode = 1                             ' Exp-Sinh
    d = exp_sinh_opt_d(func, a, eps, d)
    c = a
    v = a + d
  elseif (isfinite(b)) then
    mode = 1                             ' Exp-Sinh
    d = exp_sinh_opt_d(func, b, eps, -d)
    sign = -sign
    c = b
    v = b + d
  else
    mode = 2                             ' Sinh-Sinh
    v = 0
  end if
  
  s = func(v)
  
  do 
    p = 0 : fp = 0 : fm = 0
    h /= 2
    eh = exp(h)
    t = eh
    if (k > 0) then eh *= eh
    
    if (mode = 0) then                   ' Tanh-Sinh
        
      do 
        u = exp(1 / t - t)               ' = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2
        r = 2 * u / (1 + u)              ' = 1 - tanh(sinh(j*h))
        w = (t + 1 / t) * r / (1 + u)    ' = cosh(j*h)/cosh(sinh(j*h))^2
        x = d * r
        
        if (a + x > a) then              ' if too close to a then reuse previous fp
          y = func(a + x)
          if (isfinite(y)) then fp = y   ' if f(x) is finite, add to local sum
        end if
        
        if (b - x < b) then              ' if too close to b then reuse previous fm
          y = func(b - x)
          if (isfinite(y)) then fm = y   ' if f(x) is finite, add to local sum
        end if
        
        q = w * (fp + fm)
        p += q
        t *= eh
      loop while (abs(q) > eps * abs(p))

    else 
        
      t /= 2
      do 
        r = exp(t - 0.25 / t)            ' = exp(sinh(j*h))
        w = r
        q = 0
        
        if (mode = 1) then               ' Exp-Sinh
          x = c + d / r
          if (x = c) then exit do        ' if x hit the finite endpoint then break
          y = func(x)
          if isfinite(y) then q += y / w ' if f(x) is finite, add to local sum
        else                             ' Sinh-Sinh
          r = (r - 1 / r) / 2            ' = sinh(sinh(j*h))
          w = (w + 1 / w) / 2            ' = cosh(sinh(j*h))
          x = c - d * r
          y = func(x) 
          if isfinite(y) then q += y * w ' if f(x) is finite, add to local sum
        end if
        
        x = c + d * r
        y = func(x)
        if isfinite(y) then q += y * w   ' if f(x) is finite, add to local sum
        q *= t + 0.25 / t                ' q *= cosh(j*h)
        p += q
        t *= eh
      loop while (abs(q) > eps * abs(p))
      
    end if

    v = s - p
    s += p
    k += 1
  
  loop while (abs(v) > tol * abs(s) and k <= n)

  ' if the estimated relative error is desired, then return it
  if err_est <> 0 then err_est = abs(v) / (abs(s) + eps)
  
  ' result with estimated relative error err_est
  return sign * d * s * h
end function

' ************************************************************************

function factorial_d(n as long) as clongdouble
   dim i as long
   dim as clongdouble p = 1
   for i  =  2 to n
      p  *=  i
   next
   return p
end function

dim shared as long n
dim as clongdouble q, g, err_est

function f5(x as clongdouble) as clongdouble
   return (x^n)*exp(-x)
end function

for n = 0 to 500
   print "n ="; n
   g = factorial_d(n)
   print "n! =  "; g
   err_est = 1
   q = quad(@f5, 0, INFINITY, , , err_est)
   print "quad ="; q; " "; csng(cdbl(err_est)); " "; csng(cdbl((g - q) / g))
   print string(25, "-")
   getkey
next

sleep
Here are the examples from dodicat's first post:

Code: Select all

' example functions to integrate, exact integrals are 1, 5 and 2

function f1(x as clongdouble) as clongdouble
  return acos(x)
end function

function f2(x as clongdouble) as clongdouble
  return exp(- x / 5)
end function

function f3(x as clongdouble) as clongdouble
  return 1 / cosh(x)^2
end function

'=================================================

function f4(x as clongdouble) as clongdouble
  return 1-sqr(x*x-1)/x
end function

function f5(x as clongdouble) as clongdouble
  return x^2
end function

function f6(x as clongdouble) as clongdouble
    return 1/(1+exp(x))
    end function

'=================================================
#macro SimpsonArea(fn,Lo,Up,ordinates,sum)
scope
    dim as clongdouble x,part1,part2
    dim as clongdouble l=lo,u=up,n=iif(ordinates mod 2=0,ordinates+1,ordinates)
    var z=(len(str(n))-1)/2
     if abs(l)=infinity then l=sgn(l)*10^4
     if abs(u)=infinity then u=sgn(u)*10^4
    Var h=(U-L)/n
    x=L+.5*h:Part1=fn
    For i As long=1 To n-1
        x=L+h*i+h/2
        Part1+=fn
        x=L+h*i
        Part2+=fn
    Next i
    x=L:var fn1=fn
    x=U:var fn2=fn
    sum= (h/6)*(fn1+fn2+Part1*4+Part2*2) 
    end scope
#endmacro
'==================================================

print "integrate(acos(x), x=0..1)           = ", quad(@f1, 0, 1)
print "integrate(exp(-x/5), x=0..+inf)      = ", quad(@f2, 0, INFINITY)
print "integrate(1/cosh(x)^2, x=-inf..+inf) = ", quad(@f3, -INFINITY, INFINITY)
print "integrate(1-sqr(x*x-1)/x, x=1..+inf) = ", quad(@f4, 1, INFINITY)
print "integrate(1/(1+exp(x)), x=0..+inf)   = ", quad(@f6, 0, INFINITY)
print "integrate(x^2, x=-100.005..12.556)   = ", quad(@f5, -100.005, 120.556)

print
'Simpson macro
dim as clongdouble Result
const ordinates=1000000 ' no less anyway.
SimpsonArea(acos(x),0,1,ordinates,Result)
print "integrate(acos(x), x=0..1)           = ", Result

SimpsonArea(exp(-x/5),0,infinity,ordinates,Result)
print "integrate(exp(-x/5), x=0..+inf)      = ", Result

SimpsonArea(1/cosh(x)^2,-infinity,infinity,ordinates,Result)
print "integrate(1/cosh(x)^2, x=-inf..+inf) = ", Result

SimpsonArea(1-sqr(x*x-1)/x,1,infinity,ordinates,Result)
print "integrate(1-sqr(x*x-1)/x, x=1..+inf) = ", Result

SimpsonArea(1/(1+exp(x)),0,infinity,ordinates,Result)
print "integrate(1/(1+exp(x)), x=0..+inf)   = ", Result

SimpsonArea(x^2,-100.005,120.556,ordinates,Result)
print "integrate(x^2, x=-100.005..120.556)  = ", Result
print
print "actual value for last integral (x^2)  ";120.556^3/3- -100.005^3/3
print

sleep
Arguments must be of type clongdouble:

Code: Select all

#include "clongdouble.bi"

dim as double a
dim as clongdouble b, c, arg = 1

a = 4*atn(1)
b = 4*atn(clongdouble(1))
c = 4*atn(arg)

print a
print b
print c
print clongdouble("3.1415926535897932384626433832795")

sleep
Translated with the help from DeepL.com (free version)
dodicat
Posts: 7998
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Numerical integration

Post by dodicat »

Simpson is a local baker in these parts.
https://www.simpsonsbakery.co.uk/
So here is one of their pies, tied up.(for fun)

Code: Select all


namespace PiBySimpson
const infinity=1.797693134862316e+308
Function Simpson(fn As Function(x As Double) As Double,L As Double,U As Double,Ordinates As Integer) As Double
    if abs(l)=infinity then l=sgn(l)*10^4
    if abs(u)=infinity then u=sgn(u)*10^4
    var n=Ordinates
    If n Mod 2=0 Then n=n+1 
    Var h=(U-L)/n
    Var Part1=fn(L+.5*h)
    Var Part2=0.0
    For i As Integer=1 To n-1
        Part1+=fn(L+h*i+h/2)
        Part2+=fn(L+h*i)
    Next i
    Function= (h/6)*(fn(L)+fn(U)+Part1*4+Part2*2)  
End Function
function PieSlice(x as double) as double
   return exp(-x*x)
end function
dim as const double i=0
sub CreatePi() constructor
    *cptr(double ptr,@i)=simpson(@PieSlice,-infinity,infinity,1000000)^2
end sub
end namespace
'======================================================
'PiBySimpson.i=9 not allowed(const)
dim as const double pi=PiBySimpson.i
const as double fbpi=acos(-1)
for n as long =1 to 5
var t=timer
print pi,
print timer-t,"pi"
t=timer
print fbpi,
print timer-t,"fbpi"
print
next
'fbpi=9 not allowed(const)
'pi=9 not allowed(const)

sleep
 
jdebord
Posts: 551
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Re: Numerical integration

Post by jdebord »

Here is a test using the MPFR library according to srvaldez (https://www.freebasic.net/forum/viewtopic.php?t=24110)

The library has been compiled for 100 digits precision.

The quad routine is called with 10 levels and a tolerance eps = 1E-90

Code: Select all

' ************************************************************************
' Numerical integration by Tanh-Sinh, Sinh-Sinh and Exp-Sinh formulas
' ************************************************************************
' Ref. : R. A. van Engelen, https://www.genivia.com/files/qthsh.pdf
' ************************************************************************

#include "mpfr-fb.bi"

#define INFINITY 1/0

#define isfinite(x) (abs(x) <> INFINITY)

function exp_sinh_opt_d(func as function (x as mpfr) as mpfr, _
                        a as mpfr, eps as mpfr, d as mpfr) as mpfr
  
  dim as mpfr fl, fr, h, h2, lfl, lfr, lr, r, s

  dim as long ev, i, j
  
  ev = 2 : i = 1 : j = 32  ' j=32 is optimal to search for r
  
  h2 = func(a + d / 2) - func(a + d * 2) * 4
  
  if (isfinite(h2) and abs(h2) > 1e-5) then  ' if |h2| > 2^-16
    s = 0
    lr = 2
    
    do
      ' find max j such that fl and fr are finite
      j /= 2
      r = 1 shl (i + j)
      fl = func(a + d / r)
      fr = func(a + d * r) * r * r
      ev += 2
      h = fl - fr
    loop while (j > 1 and not(isfinite(h)))
    
    if (j > 1 and isfinite(h) and sgn(h) <> sgn(h2)) then
      lfl = fl          ' last fl=func(a+d/r)
      lfr = fr          ' last fr=func(a+d*r)*r*r
    
      do                ' bisect in 4 iterations
        j /= 2
        r = 1 shl (i + j)
        fl = func(a + d / r)
        fr = func(a + d * r) * r * r
        ev += 2
        h = fl - fr
        
        if (isfinite(h)) then
          s += abs(h)   ' sum |h| to remove noisy cases
          if (sgn(h) = sgn(h2)) then
            i += j      ' search right half
          else          ' search left half
            lfl = fl    ' record last fl=f(a+d/r)
            lfr = fr    ' record last fr=f(a+d*r)*r*r
            lr = r      ' record last r
          end if
        end if
      loop while (j > 1)
      
      if (s > eps) then ' if sum of |h| > eps
        h = lfl - lfr   ' use last fl and fr before the sign change
        r = lr          ' use last r before the sign change
        
        ' if last difference was nonzero, back up r by one step
        if (h <> 0) then r /= 2   
          
        if (abs(lfl) < abs(lfr)) then
          d /= r        ' move d closer to the finite endpoint
        else
          d *= r        ' move d closer to the infinite endpoint
        end if  
      end if
    end if
  end if
  
  return d
end function

   
function quad (func as function(x as mpfr) as mpfr, a as mpfr, b as mpfr, _
                n as long = 6, eps as mpfr = 1.0e-9, byref err_est as mpfr = 0) as mpfr

  dim as mpfr c, d, eh, fp, fm, h, p, q, r, s, sign, t, tol, u, v, w, x, y 

  dim as long k
  dim as long mode  ' Tanh-Sinh = 0, Exp-Sinh = 1, Sinh-Sinh = 2
  
  d = 1 : sign = 1 : h = 2 : tol = 10 * eps
  if (b < a) then swap a, b : sign = -1
  
  if (isfinite(a) and isfinite(b)) then
    c = (a + b) / 2
    d = (b - a) / 2
    v = c
  elseif (isfinite(a)) then
    mode = 1                             ' Exp-Sinh
    d = exp_sinh_opt_d(func, a, eps, d)
    c = a
    v = a + d
  elseif (isfinite(b)) then
    mode = 1                             ' Exp-Sinh
    d = exp_sinh_opt_d(func, b, eps, -d)
    sign = -sign
    c = b
    v = b + d
  else
    mode = 2                             ' Sinh-Sinh
    v = 0
  end if
  
  s = func(v)
  
  do 
    p = 0 : fp = 0 : fm = 0
    h /= 2
    eh = exp(h)
    t = eh
    if (k > 0) then eh *= eh
    
    if (mode = 0) then                   ' Tanh-Sinh
        
      do 
        u = exp(1 / t - t)               ' = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2
        r = 2 * u / (1 + u)              ' = 1 - tanh(sinh(j*h))
        w = (t + 1 / t) * r / (1 + u)    ' = cosh(j*h)/cosh(sinh(j*h))^2
        x = d * r
        
        if (a + x > a) then              ' if too close to a then reuse previous fp
          y = func(a + x)
          if (isfinite(y)) then fp = y   ' if f(x) is finite, add to local sum
        end if
        
        if (b - x < b) then              ' if too close to b then reuse previous fm
          y = func(b - x)
          if (isfinite(y)) then fm = y   ' if f(x) is finite, add to local sum
        end if
        
        q = w * (fp + fm)
        p += q
        t *= eh
      loop while (abs(q) > eps * abs(p))

    else 
        
      t /= 2
      do 
        r = exp(t - 0.25 / t)            ' = exp(sinh(j*h))
        w = r
        q = 0
        
        if (mode = 1) then               ' Exp-Sinh
          x = c + d / r
          if (x = c) then exit do        ' if x hit the finite endpoint then break
          y = func(x)
          if isfinite(y) then q += y / w ' if f(x) is finite, add to local sum
        else                             ' Sinh-Sinh
          r = (r - 1 / r) / 2            ' = sinh(sinh(j*h))
          w = (w + 1 / w) / 2            ' = cosh(sinh(j*h))
          x = c - d * r
          y = func(x) 
          if isfinite(y) then q += y * w ' if f(x) is finite, add to local sum
        end if
        
        x = c + d * r
        y = func(x)
        if isfinite(y) then q += y * w   ' if f(x) is finite, add to local sum
        q *= t + 0.25 / t                ' q *= cosh(j*h)
        p += q
        t *= eh
      loop while (abs(q) > eps * abs(p))
      
    end if

    v = s - p
    s += p
    k += 1
  
  loop while (abs(v) > tol * abs(s) and k <= n)

  ' result with estimated relative error err_est
  err_est = abs(v) / (abs(s) + eps)
  return sign * d * s * h
end function

' *****************************************************************************

dim shared as long n
dim as mpfr g_ref, g_quad, err_est

function f(x as mpfr) as mpfr
  return x^(n-1) * exp(-x)
end function

n = 200
g_ref = gamma(n)  
g_quad = quad(@f, 0, INFINITY, 10, 1.0E-90, err_est)
? "n        : "; n
? "gamma(n) : "; g_ref
? "quad     : "; g_quad
? "err_est  : "; err_est

sleep
Result :

Code: Select all

n        :  200
gamma(n) :  3.94328933682395251776181606966092531147567988843586631647371266622179724981701671460152142005992312e+372
quad     :  3.94328933682395251776181606966092531147567988843586631647371266622179724981701671460152142005872997e+372
err_est  :  3.0247956024245193836723639732709444567567333141164743994765700329670554374566046179272726029117624e-94
Post Reply