Decimal BASIC Open Source Project

[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