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