Numerical integration

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
Post Reply
jdebord
Posts: 550
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: 220
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: 550
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: 7996
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: 220
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: 220
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: 220
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: 220
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: 7996
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 
Post Reply