identification division.
program-id. zeta.
* Riemann zeta function of two arguments
* from Cephes-Library, zetaf.c
author. "JD".
date-written. 27.8.2006.
date-compiled.
data division.
working-storage section.
78 MAXNUMF value 1e8.
78 MACHEPF value 1e8.
78 DOMAIN value 7.
* Expansion coefficients
* for Euler-Maclaurin summation formula
* (2k)! / B2k
* where B2k are Bernoulli numbers
77 AR pic s9(4)v9(4) occurs 12 values are
12.0,
-720.0,
30240.0,
-1209600.0,
47900160.0,
-1.892437580318379e9,
7.47242496e10,
-2.9501307279181e12,
1.1646782814350e14,
-4.5979787224074e15,
1.8152105401943e17,
-7.1661652561756e18.
77 i pic 9(8).
77 x pic s9(4)v9(4).
77 q pic s9(4)v9(4).
77 a pic s9(4)v9(4).
77 b pic s9(4)v9(4).
77 k pic s9(4)v9(4).
77 s pic s9(4)v9(4).
77 w pic s9(4)v9(4).
77 t pic s9(4)v9(4).
77 te pic s9(4)v9(4).
linkage section.
77 xx pic s9(4)v9(4).
77 qq pic s9(4)v9(4).
77 res pic s9(4)v9(4).
procedure division using xx, qq returning res.
move xx to x
move qq to q
if( x = 1.0 )
move MAXNUMF to res
stop run.
if( x < 1.0 )
call mtherr using "zetaf", DOMAIN
move 0.0 to res stop run
* Euler-Maclaurin summation formula */
if( x < 25.0 )
move 9.0 to w
move function powf( q, -x ) to s
move q to a
perform varying i from 0 by 1 until i>=9
add 1.0 to a
move function powf( a, -x ) to b
add b to s
if b/s < MACHEPF
goto done
end-perform.
move a to w
compute s=s + b*w/(x-1.0);
compute s=s - 0.5 * b;
move 1.0 to a
move 0.0 to k
perform varying i from 0 by 1 until i>=12
compute a=a * x + k
compute b=b / w;
compute t = a*b/AR(i)
add t to s
compute te=t/s
move function fabsf(te) to t
if t < MACHEPF
goto done
add 1.0 to k
compute a = a* x + k
compute b = b/ w
add 1.0 to k
end-perform.
done.
move s to res
stop run
*
* Basic sum of inverse powers */
pseres.
move function powf(q,-x) to s
move q to a
perform with test after
until b/s <= MACHEPF
add 2.0 to a
move function powf( a, -x ) to b
add b to s
end-perform
move function powf( 2.0, -x ) to b
compute s = (s + b)/(1.0-b);
move s to res
stop run
end-program saegezahn.
¤ Dauer der Verarbeitung: 0.27 Sekunden
(vorverarbeitet)
¤
|
schauen Sie vor die Tür
Fenster
Die Firma ist wie angegeben erreichbar.
Die farbliche Syntaxdarstellung ist noch experimentell.
|