accuracy of gamma function algorithm (2017-04-11 16:21 by bspinoza #79690)
I am using the following program to calculate the gamma function, which uses the Stirling approximation with a Gaussian correction.
With the value n in line 330 you can set the accuracy of the calculation.
n is set to 55, when I use n = 56 then you will get the error message:
EXTYPE 1002
^
numerical operation overflow
For comparison here are some results with other gamma calculations:
calculation the gammafunction with the value 1.25:
Keisan , CASIO online calculator:
gamma (1.25) =
0.90640247705547707798267128896691800074879192072
QB64 (BASIC compiler) with _FLOAT - Type precision:
gamma(1.25) = 0.906402477055477
MS Excel:
GAMMA(1.25) = 0.906402477055477
Decimal BASIC
gamma (x) = 0.906402477054182
I think with n = 55 (see line 330) I reached the maximum for the accuracy in Decimal BASIC or has somebody an idea to increase the accuracy?
100 ! PROGRAM: GAMMA.BAS
110 ! This Programm calculates the gammafunction for a specified value.
120 ! created April 04, 2017 by Bruno Schaefer, Losheim am See, Germany
130 ! last review: April 11, 2017
140 DECLARE EXTERNAL FUNCTION gamma
150 INPUT PROMPT "value input: ": VALUE
160 LET result = gamma(value)
170 PRINT result
180 END
190 ! -------------- gamma function --------------------------
200 ! GAMMA returns the value of the Gamma function (Euler-Mascheroni)
210 ! by using the Stirling's approximation (or Stirling's formula)
220 ! with correction factor
320 EXTERNAL FUNCTION gamma(x)
330 LET n = 55
340 LET a = (x +n)
350 LET b= x
360 FOR I = 1 TO (n-1)
370 LET b = b *(x+I)
380 NEXT I
390 LET gamma = SQR((2*PI)/a)*a^a*EXP(a*(-1)+(1/(12*a))-(1/(360*a^3))) /b
400 END FUNCTION
Re: accuracy of gamma function algorithm (2017-04-13 20:26 by SHIRAISHI Kazuo #79705)
When a=57.25, a^a = 4.2918729809699E100 .
In decimal mode, MAXNUM=1E100.
That's a limit of decimal mode.
Use native arithmetic.
135 OPTION ARITHMETIC NATIVE
140 DECLARE EXTERNAL FUNCTION gamma
150 INPUT PROMPT "value input: ": VALUE
160 LET result = gamma(value)
170 PRINT result
180 END
190 ! -------------- gamma function --------------------------
200 ! GAMMA returns the value of the Gamma function (Euler-Mascheroni)
210 ! by using the Stirling's approximation (or Stirling's formula)
220 ! with correction factor
320 EXTERNAL FUNCTION gamma(x)
325 OPTION ARITHMETIC NATIVE
330 LET n = 56
340 LET a = (x +n)
350 LET b= x
360 FOR I = 1 TO (n-1)
370 LET b = b *(x+I)
380 NEXT I
382 PRINT a
385 PRINT a^a
390 LET gamma = SQR((2*PI)/a)*a^a*EXP(a*(-1)+(1/(12*a))-(1/(360*a^3))) /b
400 END FUNCTION
Re: accuracy of gamma function algorithm (2017-04-17 00:57 by jack #79711)
[Reply To Message #79690]
> I am using the following program to calculate the gamma function, which uses the Stirling approximation with a Gaussian correction.
> With the value n in line 330 you can set the accuracy of the calculation.
> n is set to 55, when I use n = 56 then you will get the error message:
> EXTYPE 1002
> ^
> numerical operation overflow
>
you can use option arithmetic decimal high by using mp_exp instead of exp
External Function mp_exp(x)
OPTION ARITHMETIC DECIMAL_HIGH
LET euler=2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952605956307381323286279434907632338298807531952510190115738341879307021540891499348841675092447614606680822648001684774118537423454424371075390777449920695517027618386062613313845830007520449338265602976067371132007093287091274437470472306969772093101416928368190255151086574637721112523897844250569536967707854499699679468644549059879316368892300987931277361782154249992295763514822082698951936680331825288693984964651058209392398294887933203625094431173012381970684161403970198376793206832823764648042953118023287825098194558153017567173613320698112509961818815930416903515988885193458072738667385894228792284998920868058257492796104841984443634632449684875602336248270419786232090021609902353043699418491463140934317381436405462531520961836908887070167683964243781405927145635490613031072085103837505101157477041718986106873969655212671546889570350354021234078
!taylor series
DECLARE NUMERIC fac,inv,x_,temp1,temp2,accum,strc,one,p,term
DECLARE NUMERIC c
LET fac=1
LET x_=fp(x)/1125899906842624
LET temp1=1
LET accum=1
LET one=1
LET p=1
Do
LET c=c+1
LET temp2=accum
LET strC=c
LET fac=fac*strc
LET p=p*x_
LET term=p/fac
LET Accum=temp2+Term
Loop Until abs(accum-temp2)<1e-1000
LET temp2=euler^IP(x)
LET temp1=Accum^1125899906842624
LET mp_exp= temp1*temp2
end Function
Re: accuracy of gamma function algorithm (2017-04-17 15:25 by bspinoza #79713)
Hi jack,
good idea with the taylor serie, but using it with my gamma calculation it seem to make a conflict with the ^sign:
"
EXTYPE 1000
power must be an integer of range from - 9223372036854775807 to...
^
"
[Reply To Message #79711]
> [Reply To Message #79690]
> > I am using the following program to calculate the gamma function, which uses the Stirling approximation with a Gaussian correction.
> > With the value n in line 330 you can set the accuracy of the calculation.
> > n is set to 55, when I use n = 56 then you will get the error message:
> > EXTYPE 1002
> > ^
> > numerical operation overflow
> >
> you can use option arithmetic decimal high by using mp_exp instead of exp
>
Re: accuracy of gamma function algorithm (2017-04-18 00:20 by jack #79715)
[Reply To Message #79713]
>
> Hi jack,
> good idea with the taylor serie, but using it with my gamma calculation it seem to make a conflict with the ^sign:
> "
> EXTYPE 1000
> power must be an integer of range from - 9223372036854775807 to...
> ^
> "
>
> [Reply To Message #79711]
> > [Reply To Message #79690]
> > > I am using the following program to calculate the gamma function, which uses the Stirling approximation with a Gaussian correction.
> > > With the value n in line 330 you can set the accuracy of the calculation.
> > > n is set to 55, when I use n = 56 then you will get the error message:
> > > EXTYPE 1002
> > > ^
> > > numerical operation overflow
> > >
> > you can use option arithmetic decimal high by using mp_exp instead of exp
> >
I made a small change to your code, now it should work.
10 OPTION ARITHMETIC DECIMAL_HIGH
100 ! PROGRAM: GAMMA.BAS
110 ! This Programm calculates the gammafunction for a specified value.
120 ! created April 04, 2017 by Bruno Schaefer, Losheim am See, Germany
130 ! last review: April 11, 2017
140 DECLARE EXTERNAL FUNCTION gamma
150 INPUT PROMPT "value input: ": VALUE
160 LET result = gamma(value)
170 PRINT result
180 END
190 ! -------------- gamma function --------------------------
200 ! GAMMA returns the value of the Gamma function (Euler-Mascheroni)
210 ! by using the Stirling's approximation (or Stirling's formula)
220 ! with correction factor
320 EXTERNAL FUNCTION gamma(x)
325 OPTION ARITHMETIC DECIMAL_HIGH
330 LET n = 55
340 LET a = (x +n)
350 LET b= x
360 FOR I = 1 TO (n-1)
370 LET b = b *(x+I)
380 NEXT I
384 LET temp=mp_log(a)*a
386 LET temp=mp_exp(temp)
390 LET gamma = SQR((2*PI)/a)*temp*mp_exp(a*(-1)+(1/(12*a))-(1/(360*a^3))+(1/(1260*a^5)-(1/(1680*a^7)))) /b
400 END FUNCTION
EXTERNAL FUNCTION mp_exp(x)
OPTION ARITHMETIC DECIMAL_HIGH
LET euler=2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952605956307381323286279434907632338298807531952510190115738341879307021540891499348841675092447614606680822648001684774118537423454424371075390777449920695517027618386062613313845830007520449338265602976067371132007093287091274437470472306969772093101416928368190255151086574637721112523897844250569536967707854499699679468644549059879316368892300987931277361782154249992295763514822082698951936680331825288693984964651058209392398294887933203625094431173012381970684161403970198376793206832823764648042953118023287825098194558153017567173613320698112509961818815930416903515988885193458072738667385894228792284998920868058257492796104841984443634632449684875602336248270419786232090021609902353043699418491463140934317381436405462531520961836908887070167683964243781405927145635490613031072085103837505101157477041718986106873969655212671546889570350354021234078
!taylor series
DECLARE NUMERIC fac,inv,x_,temp1,temp2,accum,strc,one,p,term
DECLARE NUMERIC c
LET fac=1
LET x_=fp(x)/1125899906842624
LET temp1=1
LET accum=1
LET one=1
LET p=1
Do
LET c=c+1
LET temp2=accum
LET strC=c
LET fac=fac*strc
LET p=p*x_
LET term=p/fac
LET Accum=temp2+Term
Loop Until abs(accum-temp2)<1e-1000
LET temp2=euler^IP(x)
LET temp1=Accum^1125899906842624
LET mp_exp= temp1*temp2
end Function
External Function mp_logTaylor(x)
OPTION ARITHMETIC DECIMAL_HIGH
!'taylor series
!'====================Log Guard==================
DECLARE NUMERIC g,zero,INV,XX,Term,Accum,strC,x_,tmp,tmp2
DECLARE NUMERIC T,B,one,Q,two,minus1
DECLARE NUMERIC invflag, c
LET g=x
LET zero=0
LET x_=x
LET one=1
LET two=2
LET minus1=-1
LET c=1
LET g=abs(g)
if g=x and x>0 then
If x<one Then
LET invflag=1
LET x_=one/x_
End If
LET T=x_-one
LET B=x_+one
LET accum=T/B
LET Q=T/B
LET tmp=Q
LET XX=Q*Q
Do
LET c=c+2
LET tmp2=tmp
LET Q=Q*XX
LET term=Q/c
LET Accum=tmp+Term
swap tmp,Accum
Loop Until abs(tmp-tmp2)<1e-1000
LET accum=accum*two
If invflag<>0 Then
LET accum=minus1*accum
End If
end if
LET mp_logTaylor= accum
End Function
External Function mp_log(x)
OPTION ARITHMETIC DECIMAL_HIGH
!====================Log Guard==================
DECLARE NUMERIC g,zero,approx,ans,logx,factor
LET g=x
LET zero=0
LET approx=x
LET factor=8192
LET g=abs(g)
if g=x and x>zero then
LET ry=EXP(LOG(approx)/8192)
FOR j=0 TO 5
LET r0=ry^8191
LET ry=ry-(r0*ry-approx)/(8192*r0)
NEXT j
LET r0=ry^8191
LET approx=ry-(r0*ry-approx)/(8192*r0)
LET logx=mp_logTaylor(approx)
LET ans=factor*logx
end if
LET mp_log= ans
end Function
Re: accuracy of gamma function algorithm (2017-04-18 12:40 by bspinoza #79716)
> I made a small change to your code, now it should work.
> 10 OPTION ARITHMETIC DECIMAL_HIGH
> 100 ! PROGRAM: GAMMA.BAS
> 110 ! This Programm calculates the gammafunction for a specified value.
> 120 ! created April 04, 2017 by Bruno Schaefer, Losheim am See, Germany
> 130 ! last review: April 11, 2017
> 140 DECLARE EXTERNAL FUNCTION gamma
> 150 INPUT PROMPT "value input: ": VALUE
> 160 LET result = gamma(value)
> 170 PRINT result
> 180 END
> 190 ! -------------- gamma function --------------------------
> 200 ! GAMMA returns the value of the Gamma function (Euler-Mascheroni)
> 210 ! by using the Stirling's approximation (or Stirling's formula)
> 220 ! with correction factor
> 320 EXTERNAL FUNCTION gamma(x)
> 325 OPTION ARITHMETIC DECIMAL_HIGH
> 330 LET n = 55
> 340 LET a = (x +n)
> 350 LET b= x
> 360 FOR I = 1 TO (n-1)
> 370 LET b = b *(x+I)
> 380 NEXT I
> 384 LET temp=mp_log(a)*a
> 386 LET temp=mp_exp(temp)
> 390 LET gamma = SQR((2*PI)/a)*temp*mp_exp(a*(-1)+(1/(12*a))-(1/(360*a^3))+(1/(1260*a^5)-(1/(1680*a^7)))) /b
> 400 END FUNCTION
> .....
Thanks again, Jack!
Now the program works without error but there seems to be a small mistake in the code, because the result isn't correct:
e.g. value input: 1.25 results in: 0.5907613718183358... , but it should be: 0.906402477055477...
I will check the code in the next days ...
Re: accuracy of gamma function algorithm (2017-04-18 20:27 by jack #79720)
[Reply To Message #79719]
> Sorry, Jack,
> your version works fine with Decimal BASIC Version 7.7.8.1!
>
> By coincidendeI found a strange behaviour of Decimal BASIC with this program:
>
> Decimal BASIC Version 7.7.8.1:
> value input: 1.25; result: .906402477055477077...
>
> Decimal BASIC Version 7.8.0 also:
> value input: 1.25; result: .906402477055477077...
>
> Decimal BASIC Version 6.6.2.0 however:
> value input: 1.25; result: .590761371818335894... ???
>
> and Decimal BASIC Version 6.6.2.2 too:
> value input: 1.25; result: .590761371818335894... ???
>
> Can somebody confirm this?
> Is it an error in my logic or is it a bug?
>
> The Decimal BASIC Version 6.6.2.2 was downloaded from https://osdn.net/projects/decimalbasic/releases/p8178. It is the English version BASIC6622En_Win.zip from 2017-04-10.
there's a bug somewhere, but while that gets resolved here's a temporary solution
in the function mp_exp replace the constant 1125899906842624 to 1073741824
there are two occurrences of the constant.
Re: accuracy of gamma function algorithm (2017-04-18 22:00 by jack #79721)
there's problem with the ^ operator, if the exponent is >=2147483648 it fails.
with exponent=1073741824 it works ok, have not tested for values in-between.
Re: accuracy of gamma function algorithm (2017-04-19 08:40 by SHIRAISHI Kazuo #79725)
[Reply To Message #79721]
> there's problem with the ^ operator, if the exponent is >=2147483648 it fails.
> with exponent=1073741824 it works ok, have not tested for values in-between.
I found FPC 3.0.0 had a bug on the INT function for COMP type.
FPC 3.0.2 might fix it, compiling with Lazarus 1.6.4 +FPC 3.0.2 , the above bug on Decimal BASIC vanished.
Re: accuracy of gamma function algorithm (2017-04-19 02:08 by bspinoza #79724)
[Reply To Message #79723]
> [Reply To Message #79720]
> > [Reply To Message #79719]
> > > Sorry, Jack,
> > > your version works fine with Decimal BASIC Version 7.7.8.1!
>
> Just curious. Why would anyone use 6.x when 7.x is available?
>
> Tom L
Question: "Why would anyone use 6.x when 7.x is available?"
Answer: The Linux version !