Bernoulli numbers: Difference between revisions

m
→‎{{header|REXX}}: added/changed comments and whitespace, simplified some statements, optimized some functions.
m (Phix/mpfr)
m (→‎{{header|REXX}}: added/changed comments and whitespace, simplified some statements, optimized some functions.)
Line 3,284:
<lang rexx>/*REXX program calculates N number of Bernoulli numbers expressed as fractions. */
parse arg N .; if N=='' | N=="," then N= 60 /*Not specified? Then use the default.*/
d=numeric digits max(9, n*2;) if d>digits() then numeric digits d /*increase the decimal digits if needed*/
!.=0; w= max(length(N), 4); Nw= N + w + N % 4 /*used for aligning (output) fractions.*/
say 'B(n)' center("Bernoulli number expressed as a fraction", max(78-w, Nw)) /*title.*/
say copies('─',w) copies("─", max(78-w,Nw+2*w)) /*display 2nd line of title, separators*/
!.= 0; do #=0 to N /*process the numbers from 0 ──► N. */
b= bern(#); if b==0 then iterate /*calculate Bernoulli number, skip if 0*/
indent= max(0, nW - pos('/', b) ) /*calculate the alignment (indentation)*/
say right(#, w) left('', indent) b /*display the indented Bernoulli number*/
end /*#*/ /* [↑] align the Bernoulli fractions. */
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
Line 3,298:
if x==1 then return '-1/2' /* " " " " " one. */
if x//2 then return 0 /* " " " " " odds > 1.*/
do j=2 to x by 2; jp= j+1; d= j+j /*process the positive integers up to X*/
sn= 1 - j /*define the numerator. */
sd= 2 /* " " denominator. */
Line 3,304:
parse var @.k bn '/' ad /*get a previously calculated fraction.*/
an= comb(jp, k) * bn /*use COMBination for the next term. */
$LCMlcm= LCM(sd, ad) /*use Least Common Denominator function*/
sn= $LCMlcm % sd * sn; sd= $LCMlcm /*calculate the current numerator. */
an= $LCMlcm % ad * an; ad= $LCMlcm /* " " next " */
sn= sn + an /* " " current " */
end /*k*/ /* [↑] calculate the SN/SD sequence.*/
sn= -sn /*adjustflip the sign for the numerator. */
sd= sd * jp /*calculate the denominator. */
if sn\==1 then do; _= GCD(sn, sd) /*get the Greatest Common Denominator.*/
sn= sn%_; sd= sd%_ /*reduce the numerator and denominator.*/
end /* [↑] done with the reduction(s). */
@.j= sn'/'sd /*save the result for the next round. */
end /*j*/ /*j*/ /* [↑] done calculating Bernoulli #'s.*/
return sn'/'sd
/*──────────────────────────────────────────────────────────────────────────────────────*/
comb: procedure expose !.; parse arg x,y; if x==y then return 1
if !.c.x.y \== 0 then return !.c.x.y /*combination computed before?*/
if x-y < y then y= x-y; z= perm(x, y); do j=2 to y; z= z % j; end /*Jx-y < y? Then use a new Y.*/
z= perm(x, y); do j=2 for y-1; z= z % j
end /*j*/
!.c.x.y= z; return z /*assign memoization & return.*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
GCD: procedure; parse arg x,y; x= abs(x)
do until y==0; parse value x//y y with y x; end; return x
/*──────────────────────────────────────────────────────────────────────────────────────*/
LCM: procedure; parse arg x,y /*xX=absABS(xX); yY=absABS(yY) not ¬ needed for Bernoulli numbers#s.*/
if y /*IF Y==0 THEN then returnRETURN 0 " " " /*if zero," then LCM is also zero." */
d$= x * y /*calculate part of the LCM here. */
do until y==0; parse value x//y y with y x
end /*until*/ /* [↑] this is a short & fast GCD*/
return d$ % x /*divide the pre─calculated value.*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
perm: procedure expose !.; parse arg x,y; if !.p.x.y \== 0 then return !.p.x.y
z= 1; do j=x-y+1 to x; z= z*j; end; !.p.x.y= z; return z</lang>
</lang>
{{out|output|text=&nbsp; when using the default input:}}
<pre>