View Full Version : Riemann zeta function in O2
also works for fractional arguments.
<edit> just timed the zeta function in the gsl library and it's almost three times faster than this implementation, bummer.
some timings
O2 858 nSec
gsl 343 nSec
FreeBASIC 378 nSec
my routine compiled with FreeBASIC is only 10% slower than gsl, that's not bad.</edit>
Uses "Console" , "oxygen"
Dim i , j, t As Long
Dim scr As String
Dim As Long pZeta
Dim As Long pFinish
Dim x, y As Double
scr = "function zeta(byval n as double) as double link #pZeta
'Riemann zeta function evaluated using the Euler-Maclaurin summation formula.
'first we sum the first ten terms of the Zeta series,
'then use the Euler-Maclaurin summation formula.
dim as double c(7)
dim as integer i, k, k1, nc
dim as double s, dx, n1
c(1) = 8.3333333333333330e-02 'B(2)/2!
c(2) = -1.3888888888888890e-03 'B(4)/4!
c(3) = 3.3068783068783070e-05 'B(6)/6!
c(4) = -8.2671957671957670e-07 'B(8)/8!
c(5) = 2.0876756987868100e-08 'B(10)/10!
c(6) = -5.2841901386874930e-10 'B(12)/12!
c(7) = 1.3382536530684680e-11 'B(14)/14!
nc=10
s=0
for i=1 to nc
s=s+1/i^n
next
'f(k)=(k+10)^-n
s=s+nc^(1-n)/(n-1) 'int(f(k)dk,k=0..infinity)=10^(1-n)/(n-1)
dx=-n*nc^(-n-1) 'first derivative of f(k)
s=s-.5*nc^(-n) 'f(0)*B(1), first term of the Euler-Maclaurin summation formula.
n1=nc^(-2)
k=1
k1=2
for i=1 to 7
s=s-dx*c(i) 'next term of the series
dx=dx*n1*(n+k1)*(n+k) 'next derivative of f(k)
k=k+2
k1=k1+2
next
function = s
end function
sub finish() link #pFinish
terminate
end sub
"
O2_Asm scr
If O2_Error Then
MsgBox 0,O2_Error
Stop
Else
O2_Exec
End If
Declare Function Zeta(ByVal Double ) As Double At pZeta
Declare Sub Finish() At pFinish
For j=2 To 20
x=Zeta(j)
Console_WriteLine j, x
Next
t=GetTickCount
For i=1 To 10000
For j=2 To 50
x=Zeta(j)
'Console_WriteLine x
Next
Next
t=GetTickCount-t
Console_WriteLine "10000 iterations of computing the Zeta function from 2 to 50 in "&t&" nSec."
Console_WriteLine "All done. Press any key to finish"
Console_WaitKey
Finish()
ErosOlmi
30-12-2013, 11:02
I think 1/4 of the time is taken by thinBasic just for interpreting and calling O2 function in the loop:
For i=1 To 10000
For j=2 To 50
x=Zeta(j)
Next
Next
Below code is 100% O2 code and it takes about 600 ms instead of 800 ms
I moved constant declarations outside zeta function (they are evaluated just once during O2 just in time compilation.
Created a test lRun function to see how much time was taken by thinBasic for interpretation in the main loop.
Uses "Console" , "oxygen"
Dim i , j, t As Long
Dim scr As String
Dim As Long pZeta
Dim As Long pRun
Dim As Long pFinish
Dim x, y As Double
'---O2 code
scr = "
Dim As Double c(7)
c(1) = 8.3333333333333330e-02 'B(2)/2!
c(2) = -1.3888888888888890e-03 'B(4)/4!
c(3) = 3.3068783068783070e-05 'B(6)/6!
c(4) = -8.2671957671957670e-07 'B(8)/8!
c(5) = 2.0876756987868100e-08 'B(10)/10!
c(6) = -5.2841901386874930e-10 'B(12)/12!
c(7) = 1.3382536530684680e-11 'B(14)/14!
function zeta(byval n as double) as double link #pZeta
'Riemann zeta function evaluated using the Euler-Maclaurin summation formula.
'first we sum the first ten terms of the Zeta series,
'then use the Euler-Maclaurin summation formula.
Dim As Integer i, k, k1, nc
Dim As Double s, dx, n1
nc=10
s=0
For i=1 To nc
s=s+1/i^n
Next
'f(k)=(k+10)^-n
s=s+nc^(1-n)/(n-1) 'int(f(k)dk,k=0..infinity)=10^(1-n)/(n-1)
dx=-n*nc^(-n-1) 'first derivative of f(k)
s=s-.5*nc^(-n) 'f(0)*B(1), first term of the Euler-Maclaurin summation formula.
n1=nc^(-2)
k=1
k1=2
For i=1 To 7
s=s-dx*c(i) 'next term of the series
dx=dx*n1*(n+k1)*(n+k) 'next derivative of f(k)
k=k+2
k1=k1+2
Next
Function = s
End Function
function lRun() as long link #pRun
dim i as long
dim j as long
dim x as double
For i=1 To 10000
For j=2 To 50
x=Zeta(j)
Next
Next
end function
Sub finish() link #pFinish
terminate
End Sub
"
'---JIT compile O2 code
O2_Asm scr
If O2_Error Then
MsgBox 0,O2_Error
Stop
Else
O2_Exec
End If
'---Start thinBasic code
Declare Function Zeta(ByVal Double ) As Double At pZeta
Declare Function lRun() As Long At pRun
Declare Sub Finish() At pFinish
For j=2 To 20
x=Zeta(j)
Console_WriteLine j, x
Next
t=GetTickCount
lRun()
'For i=1 To 10000
' For j=2 To 50
' x=Zeta(j)
' Next
'Next
t=GetTickCount-t
Console_WriteLine "10000 iterations of computing the Zeta function from 2 to 50 in "&t&" mSec."
Console_WriteLine "All done. Press any key to finish"
Console_WaitKey
Finish()
Do you know where I can get a pre-compiled version (DLL) of GSL library?
Thanks a lot
Eros
mike lobanovsky
30-12-2013, 18:11
1. Smart optimizing static C compilers like GCC used to compile GSL may unroll user loops transparently in full or in batches. This is done in order to reduce the overhead and resources involved in reproducing such high-level constructs in low-level machine code they are resolved to by the compiler. This may apply to both the inner loop of Zeta function proper and the outer loop of test script.
2. Such compilers may also transparently optimize the low-level code to use a faster SSE floating-point instruction set if such a capability is detected automatically. Modern Pentium and Athlon CPU's do support this instruction set in some form or other.
These two points may account for why original GSL performs faster than FreeBasic whose compiler is somewhat "dumber". JIT compilers like O2 (and DynC, for that matter ;) ) won't perform any such optimizations at all for the sake of much, much faster compilation speed.
Do you know where I can get a pre-compiled version (DLL) of GSL library?
Here's everything (http://gnuwin32.sourceforge.net/packages/gsl.htm) you may need to get a 32-bit Windows binary of GSL, the GNU scientific library for numerical computing, working for you.
http://www.fbsl.net/phpbb2/images/smilies/icon_ml_noel.gif
thank you Eros and mike for your comments, Eros if you have MinGW installed on your system then you can compile the GSL sources no problem but you may want to run the strip --unneeded command afterwards to reduce the size of the dll.
<off topic> one of my favorite tricks is to convert a static library to a dll by first running ar x libname.a and then gcc -shared -static -o libname.dll *.o </off topic>
I have attached the latest version of the GSL dll.
forgot post the gsl test code
'gsl test
Uses "Math", "Console"
Type gsl_sf_result
Vl As Double
Er As Double
End Type
Declare Function gamma CDECL Lib "gsl.dll" Alias "gsl_sf_gamma" ( ByVal x As Double) As Double
Declare Function gamma_e CDECL Lib "gsl.dll" Alias "gsl_sf_gamma_e" ( ByVal x As Double, ByVal r As gsl_sf_result) As Long
Declare Function zeta CDECL Lib "gsl.dll" Alias "gsl_sf_zeta" ( ByVal n As Double) As Double
Declare Function zeta_e CDECL Lib "gsl.dll" Alias "gsl_sf_zeta_e" ( ByVal x As Double, ByVal r As gsl_sf_result) As Long
Dim i, j, t, status As Long
Dim x, y As Double
Dim gsl_result As gsl_sf_result
'two ways to call the functions depending on your need
'y=gamma(0.5)
'Console_WriteLine "status = "&status, "Gamma 0.5 = "&y
'or the following
status=gamma_e(0.5, gsl_result)
Console_WriteLine "status = "&status, "Gamma 0.5 = "&gsl_result.vl, "Error = "&gsl_result.er
status=zeta_e(3, gsl_result)
Console_WriteLine "status = "&status, "Zeta 3 = "&gsl_result.vl, "Error = "&gsl_result.er
Console_WriteLine
t=GetTickCount
For i=1 To 10000
For j=2 To 50
x=zeta(j)
'Console_WriteLine j,x
Next
Next
t=GetTickCount-t
Console_WriteLine "10000 iterations of computing the Zeta function from 2 to 50 in "&t&" mSec."
Console_WriteLine "All done. Press any key to finish"
Console_WaitKey
John Spikowski
30-12-2013, 23:44
Hi Jack,
Back in the summer of 2011, I started a thread on how to build a ScriptBasic extension module. I thought the GSL library would have practical use after the tutorial was completed. I forgot about it until your thread here about GSL. FWIW this is where I left off with that project.
ScriptBasic GSL extension mudule (http://www.allbasic.info/forum/index.php?topic=115.msg1650#msg1650)
I did a quick hypot test to see if the extension module was still working and I have the GSL dependencies installed.
The hypotenuse of a right triangle with sides the length of 3.25 and 4.5
DECLARE SUB hypot ALIAS "_hypot" LIB "gsl"
x = 3.25
y = 4.5
PRINT "SB hypot: ",FORMAT("%2.8f",SQR(x^2 + y^2)),"\n"
PRINT "GSL hypot: ",FORMAT("%2.8f",hypot(x, y)),"\n"
jrs@laptop:~/sb/sb22/gsl$ scriba testgsl.sb
SB hypot: 5.55090083
GSL hypot: 5.55090083
jrs@laptop:~/sb/sb22/gsl$
Good to see SB is maintaining accuracy out to 8 decimal places.
mike lobanovsky
30-12-2013, 23:45
@Jack:
Yep, your DLL seems to be somewhat newer than the one which can be downloaded via the link in my message. Thanks!
@All:
Gentlemen, please be warned that upon closer investigation the both DLL's are not 100% perfect from Windows' standpoint. Matter is, they've been compiled with MinGW/GCC which stands for Minimum GNU Compiler Collection for Windows. Some versions of this GNU software generate an imperfect section table in the Windows PE header of their binaries (both EXE and DLL). Sometimes they also append some trash at the end of trailing .rsrc section in the file image.
I don't know what the real reason for such negligeance is but I'm inclined to attribute it to general blind hatred among the GNU GPL Linux community towards its glorious and far more fortunate MS Windows rival. It often leads to trivial ignorance of basic concepts of their principal enemy they should've known inside out by now.
For instance, FBSL is compiled with MinGW/GCC v4.3.3. So far it's been the best version to generate speed-critical Windows binaries with. All the later versions generate Windows binaries that are noticeably fatter and/or run up to 40% slower despite all the advanced optimizations of GCC's -O3 option. Regretfully, this GCC version generates faulty section tables and wrong PE checksums, and it also appends trash at the end of the file. I have to use my own FBSL script at the end of the toolchain to patch the resultant binaries to be fully conformant with what Windows expects from a binary that it launches for execution.
Another example of such a careless GNU C compiler is Fabrice Bellard's Tiny C Compiler that's used as a prototype for FBSL's Dynamic C layer.
Summing up the aforesaid, the both libraries will run on your PC's but they can hardly be accepted as high-quality Windows product. Redistributing them in this form on the net is evil and detrimental to the reputation of MS Windows. In case professionals do confirm the importance of this particular piece of software for public at large, it's source code should be adapted to the MSVC family of compilers and recompiled accordingly.
thank you John and mike.
mike, this is the first time I read of this, but you are right about the linux community disdain for windows os.
mike, will your script work on any exe or dll created with the GNU tool chain?
if so, would you share it with me?
btw, here's vc build of gsl-1.15 https://code.google.com/p/gsl-w32/downloads/list
hello Eros, here's a version that's about 30% faster than the one you posted.
Uses "Console" , "oxygen"
Dim i , j, t As Long
Dim scr As String
Dim As Long pZeta
Dim As Long pRun
Dim As Long pFinish
Dim x, y As Double
'---O2 code
scr = "
dim As Double c(12)
c(1) = 8.3333333333333330e-02 'B(2)/2!
c(2) = -1.3888888888888890e-03 'B(4)/4!
c(3) = 3.3068783068783070e-05 'B(6)/6!
c(4) = -8.2671957671957670e-07 'B(8)/8!
c(5) = 2.0876756987868100e-08 'B(10)/10!
c(6) = -5.2841901386874930e-10 'B(12)/12!
c(7) = 1.3382536530684680e-11 'B(14)/14!
c(8) = -3.3896802963225829e-13 'B(16)/16!
c(9) = 8.5860620562778450e-15 'B(16)/16!
c(10)= -2.1748686985580620e-16 'B(16)/16!
c(11)= 5.5090028283602300e-18 'B(16)/16!
c(12)= -1.3954464685812520e-19 'B(18)/18!
Function zeta(ByVal n As Double) As Double link #pZeta
'Riemann zeta function evaluated using the Euler-Maclaurin summation formula.
'first we sum the first seven terms of the Zeta series,
'then use the Euler-Maclaurin summation formula.
Dim As Integer i, k, k1, nc
Dim As Double s, dx, n1
nc=7
' f(k)=(k+nc)^-n
' inf nc inf inf
' ==== ==== / ==== B
' \ 1 \ 1 [ \ (2*k) (2*k-1)
' > ----- = > ----- + I f(k) dk + B * f(0) - > ------ f (0)
' / n / n ] 1 / (2*k)!
' ==== k ==== k / ====
' k = 1 k = 1 0 k = 1
s=0
For i=2 To nc
s=s+1/i^n
Next
' inf
' /
' [
s=s+nc^(1-n)/(n-1) ' I f(k) dk = nc^(1-n)/(n-1)
' ]
' /
' 0
dx=-n*nc^(-n-1) 'first derivative of f(0)
s=s-.5*nc^(-n) 'f(0)*B(1), first term of the Euler-Maclaurin summation formula.
n1=nc^(-2)
k=1
k1=2
For i=1 To 12
s=s-dx*c(i) 'next term of the series
' (2*i+1)
dx=dx*n1*(n+k1)*(n+k) 'f (0)
k=k+2
k1=k1+2
Next
Function = 1+s
End Function
Function lRun() As Long link #pRun
Dim i As Long
Dim j As Long
Dim x As Double
For i=1 To 10000
For j=2 To 50
x=Zeta(j)
Next
Next
End Function
Sub finish() link #pFinish
terminate
End Sub
"
'---JIT compile O2 code
O2_ASM scr
If O2_ERROR Then
MsgBox 0,O2_ERROR
Stop
Else
O2_EXEC
End If
'---Start thinBasic code
Declare Function Zeta(ByVal Double ) As Double At pZeta
Declare Function lRun() As Long At pRun
Declare Sub Finish() At pFinish
For j=2 To 20
x=Zeta(j)
Console_WriteLine j, x
Next
t=GetTickCount
lRun()
'For i=1 To 10000
' For j=2 To 50
' x=Zeta(j)
' Next
'Next
t=GetTickCount-t
Console_WriteLine "10000 iterations of computing the Zeta function from 2 to 50 in "&t&" mSec."
Console_WriteLine "All done. Press any key to finish"
Console_WaitKey
Finish()
mike lobanovsky
31-12-2013, 02:55
He-he Jack,
Everything in this world has its beginning, including my public GNU GPL criticism. But there's only one thing in this world that has no end, and that's unmotivated and incompetent linuxoid disdain. http://www.fbsl.net/phpbb2/images/smilies/icon_ml_winkwink.gif
And no, Jack, my script is only universal for the FBSL toolchain purposes in the sense that it doesn't depend on exactly which built-in layers this or that FBSL compilation contains; it'll do its job to fix Fbsl.exe, Fbsl.dll and Fbsl_Tiny.exe regardless. But I can illustrate, and elaborate on, what actually takes place in a lopsided binary generated by MinGW/GCC.
First, please have a look at what's there in a raw Fbsl.dll using PE Explorer (http://www.heaventools.com/) and any hex editor:
http://i1240.photobucket.com/albums/gg490/FbslGeek/Fbsldll_BAD.png
As you may see, there's a bad NULL pointer (i.e. offset within the file image on disk) to an empty .bss section in the DLL's section header. For a zero raw data length .bss section, its pointer in the table should in fact be the same as the pointer to the next section in the table that follows, whichever section that might be. There's always at least one section that follows because the trailing one is always section number 15 which is reserved by Windows presumably for such cases when the linker proscribes a zero-length section into the last permissible entry of section header table - section number 14.
You may also deduce that the end of file would be expected at a file offset equal to the pointer to the last section actually proscribed in the section header table plus that section's raw data size, which in this particular case would be at 0009DC00h. But as seen in the hex editor snapshot, there are some extra 29 bytes of bla bla rubbish written in there. These bytes are unexpected by Windows whose PE loader is governed exclusively by the PE header information where these extra bytes are not accounted for in any way.
Now have a look at the same DLL after having been fixed by my script:
http://i1240.photobucket.com/albums/gg490/FbslGeek/Fbsldll_GOOD.png
Do you feel the difference? http://www.fbsl.net/phpbb2/images/smilies/icon_ml_shine.gif
And finally, have a look at Gsl.dll that you uploaded to this thread:
http://i1240.photobucket.com/albums/gg490/FbslGeek/Gsldll_BAD.png
Seems like you've already seen some such trash somewhere before, ainchya? http://www.fbsl.net/phpbb2/images/smilies/icon_ml_winkwink.gif How am I supposed to regard the GCC devs professionally after all that? And this is going on for years! http://www.fbsl.net/phpbb2/images/smilies/icon_m_tomato.gif
A Happy New Year to you and all the best in tidying up your GNU software archives even if only manually with a Windows calculator and a hex editor in hand!
http://www.fbsl.net/phpbb2/images/smilies/icon_ml_noel.gif
P.S. Thanks for this extra MSVC link but its version number looks somewhat alarming...
thank you mike for the info.
Happy new year :)
ErosOlmi
31-12-2013, 11:43
Jack,
I've developed Zeta function inside thinBasic Math module using you latest O2 code.
It seems it is a little bit faster even if interpreted: 430 ms here on my PC and results seem OK.
To test the script, leave thinBasic_Math.dll in the same directory of the script: thinBasic will search for needed modules in script directory first.
If OK for you, I will release it in next thinBasic update.
Ciao
Eros
Eros
yes of course it's OK to use my code, but perhaps some boundary checks should be included, maybe you did that already, for double x>52 the result is simply 1 also my function does not handle negative arguments, perhaps it should.
on second thought, if it's OK with you I will write a Zeta function using Chebyshev approximation, it should be about 8 times faster and I will make it more bullet proof.
hello Eros.
here's the zeta function using Chebyshev approximations, it's 3 times faster for positive arguments than the one you used in the Math module, also works for negative arguments.
code is in PowerBasic Console compiler.
#COMPILE EXE
#DIM ALL
FUNCTION gamma(BYVAL y AS DOUBLE) AS DOUBLE
DIM Pi AS STATIC DOUBLE
DIM sq2pi AS STATIC DOUBLE
DIM g AS STATIC DOUBLE
DIM t AS DOUBLE, w AS DOUBLE, gam AS DOUBLE, x AS DOUBLE
DIM i AS LONG, cg AS LONG
DIM c(15) AS STATIC DOUBLE
Pi = 3.1415926535897932385
sq2pi = 2.50662827463100050241577 'sqrt(2Pi)
g = 607/128 ' best resu'ts when 4<=g<=5
x = y-1
cg = 14
'Lanczos approximation for the complex plane
'calculated using vpa digits(256)
'the best set of coeffs was selected from
'a solution space of g=0 to 32 with 1 to 32 terms
'these coeffs really give superb performance
'of 15 sig. digits for 0<=real(z)<=171
'coeffs should sum to about g*g/2+23/24
'http://www.numericana.com/answer/info/godfrey.htm
c( 1) = 0.99999999999999709182
c( 2) = 57.156235665862923517
c( 3) = -59.597960355475491248
c( 4) = 14.136097974741747174
c( 5) = -0.49191381609762019978
c( 6) = 0.000033994649984811888699
c( 7) = 0.000046523628927048575665
c( 8) = -0.000098374475304879564677
c( 9) = 0.00015808870322491248884
c(10) = -0.00021026444172410488319
c(11) = 0.00021743961811521264320
c(12) = -0.00016431810653676389022
c(13) = 0.000084418223983852743293
c(14) = -0.000026190838401581408670
c(15) = 0.0000036899182659531622704
IF ( x < 0 ) THEN
x=-x
IF FRAC(x)=0 THEN
gam=10^4932
ELSE
t = c(1)
FOR i=1 TO cg
t = t + c(i+1)/(x+i)
NEXT
w = x + g + 0.5
gam=w^(x+0.5)*EXP(-w)*sq2pi*t
gam = Pi*x/(gam*SIN(Pi*x))
END IF
ELSE
t = c(1)
FOR i=1 TO cg
t = t + c(i+1)/(x+i)
NEXT
w = x + g + 0.5
gam=w^(x+0.5)*EXP(-w)*sq2pi*t
END IF
FUNCTION = gam
END FUNCTION
FUNCTION zeta(BYVAL n AS DOUBLE) AS DOUBLE
'Riemann zeta function evaluated using tables from
'Chebyshev Approximations for the Riemann Zeta Function
'By W. J. Cody,* K. E. Hillstrom,* and Henry C. Thacher, Jr.**
'except for the second set which was computed using Maple.
DIM p(8) AS STATIC DOUBLE
p(0) = 12871681214.82446392809
p(1) = 13753969320.37025111825
p(2) = 5106655918.364406103683
p(3) = 856147100.2433314862469
p(4) = 74836181.24380232984824
p(5) = 4860106.585461882511535
p(6) = 273957.4990221406087728
p(7) = 4631.710843183427123061
p(8) = 57.87581004096660659109
DIM q(8) AS STATIC DOUBLE
q(0) = 25743362429.64846244667
q(1) = 5938165648.679590160003
q(2) = 900633037.3261233439089
q(3) = 80425366.34283289888587
q(4) = 5609711.759541920062814
q(5) = 224743.1202899137523543
q(6) = 7574.578909341537560115
q(7) =-23.73835781373772623086
q(8) = 1.000000000000000000000
DIM p1(8) AS STATIC DOUBLE
p1(0) = 0.1811021757722157193018e-1
p1(1) =-0.1781077292672229847925e-2
p1(2) = 0.3271270368846556844317e-2
p1(3) =-0.1183785392295491855720e-2
p1(4) = 0.2811461476445431284369e-3
p1(5) =-0.4540699985455522705413e-4
p1(6) = 0.4921410048947204581619e-5
p1(7) =-3.3060968398926623965732e-7
p1(8) = 1.0649769599519517944407e-8
DIM q1(8) AS STATIC DOUBLE
q1(0) =-.95084345744468539815114
q1(1) = .23724094500966523346054
q1(2) = .12719149298560634692396
q1(3) =-0.3742718269184690837860e-1
q1(4) =-0.2745688015667984911198e-1
q1(5) =-0.7479966479870980001960e-2
q1(6) =-0.1221861376076105986461e-2
q1(7) =-0.1189441500804843551022e-3
q1(8) =-0.6171014074208758609628e-5
DIM p2(7) AS STATIC DOUBLE
p2(0) = 1.66156480515774675916e-11
p2(1) =-4.68068827660654526862e-09
p2(2) = 5.83519727319147047318e-07
p2(3) =-4.17644012643145602124e-05
p2(4) = 1.85468422843597959483e-03
p2(5) =-5.11288800220490240591e-02
p2(6) = 8.10450231751100353193e-01
p2(7) =-5.69951948768478922618e+00
DIM q2(9) AS STATIC DOUBLE
q2(0) =-6.99562633519191654964e-10
q2(1) =-1.77757961895149256941e-08
q2(2) =-9.82231825734078036442e-07
q2(3) =-2.84927282759096487594e-05
q2(4) =-5.81727909388048093531e-04
q2(5) =-1.15848749169766585807e-02
q2(6) =-1.28149124051978195742e-01
q2(7) =-1.11913057349097709324e+00
q2(8) =-7.67928761604628812537e-01
q2(9) = 1.00000000000000000000e+00
DIM p3(8) AS STATIC DOUBLE
p3(0) = 1.0314487718885971168e-15
p3(1) =-5.1258461396468824062e-13
p3(2) = 1.1294879419487354786e-10
p3(3) =-1.4423466537313095228e-08
p3(4) = 1.1682467698445809766e-06
p3(5) =-6.1497516799031480614e-05
p3(6) = 2.0559467798883032750e-03
p3(7) =-3.9933942939466886853e-02
p3(8) = 3.4523497673617845708e-01
DIM q3(8) AS STATIC DOUBLE
q3(0) = 5.9395941728841905020e-11
q3(1) =-6.0475535907999180572e-09
q3(2) = 3.6468020866838856275e-07
q3(3) =-1.2945690556801181241e-05
q3(4) = 3.2018949847022925001e-04
q3(5) =-5.0780155709999407748e-03
q3(6) = 5.4962890788158726560e-02
q3(7) =-3.2451761115597241852e-01
q3(8) = 1.0000000000000000000e+00
DIM pi AS STATIC DOUBLE
DIM pn AS DOUBLE, qn AS DOUBLE, n2 AS DOUBLE, y AS DOUBLE
DIM j AS LONG
pi=3.141592653589793
IF n<0 THEN
n2=ABS(n)
IF FRAC(n2)=0 THEN
j=n2
IF j MOD 2=0 THEN
FUNCTION=0
EXIT FUNCTION
END IF
END IF
n2=1+n2
FUNCTION=2^n*pi^(-n2)*COS(0.5*pi*n2)*gamma(n2)*zeta(n2)
ELSEIF n>55 THEN
FUNCTION=1
ELSE
SELECT CASE n
CASE n=0.5 TO 5
pn=p(0)
qn=q(0)
n2=n
FOR j=1 TO 8
pn+=n2*p(j)
qn+=n2*q(j)
n2*=n
NEXT
qn*=(n-1)
FUNCTION=pn/qn
CASE n=5 TO 11
y=n/3-8/3
n2=y*y
pn=p1(0)+p1(1)*n
qn=q1(0)+q1(1)*n
FOR j=2 TO 8
pn+=n2*p1(j)
qn+=n2*q1(j)
n2*=y
NEXT
FUNCTION=pn/qn+1
CASE n=11 TO 25
y=1/n
n2=y
pn=p2(0)+p2(1)*n2
qn=q2(0)+q2(1)*n2
FOR j=2 TO 7
n2*=y
pn+=n2*p2(j)
qn+=n2*q2(j)
NEXT
n2*=y
qn+=n2*q2(8)
n2*=y
qn+=n2*q2(9)
FUNCTION=1+2^(y*pn/qn-n)
CASE n=25 TO 55
y=1/n
n2=y
pn=p3(0)+p3(1)*n2
qn=q3(0)+q3(1)*n2
FOR j=2 TO 8
n2*=y
pn+=n2*p3(j)
qn+=n2*q3(j)
NEXT
FUNCTION=1+2^(y*pn/qn-n)
END SELECT
END IF
END FUNCTION
FUNCTION PBMAIN () AS LONG
DIM s AS DOUBLE, t AS DOUBLE
DIM i AS LONG, j AS LONG
t=TIMER
FOR i=1 TO 10000
FOR j=2 TO 50
s=zeta(j)
' ?j,s
NEXT
NEXT
t=TIMER-t
PRINT t*1000;" nSec"
INPUT t
END FUNCTION
fixed copy/paste mistake.