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.
Members who have read this thread: 0
Reply With Quote
Bookmarks