## Numerical integration

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

### Numerical integration

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

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 string(25, "-")
next
``````
jdebord
Posts: 550
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

### Re: Numerical integration

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

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

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

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

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

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
-------------------------
n = 103
n! =   9.902900716486181e+163
-------------------------
n = 104
n! =   1.029901674514563e+166
-------------------------``````
hhr
Posts: 220
Joined: Nov 29, 2019 10:41

### Re: Numerical integration

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

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

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 ``````