Faulhaber's triangle: Difference between revisions
Content added Content deleted
m (Updated description and link for Fōrmulæ solution) |
|||
Line 1,680: | Line 1,680: | ||
faulhaber(4, 1000) |
faulhaber(4, 1000) |
||
200500333333300</pre> |
200500333333300</pre> |
||
=={{header|jq}}== |
|||
'''Works with [[jq]]''' (*) |
|||
'''Works with gojq, the Go implementation of jq''' |
|||
The solution presented here requires a "rational arithmetic" package, such |
|||
as the "Rational" module at [[Arithmetic/Rational#jq]]. |
|||
As a reminder, the jq directive for including such a package appears as the first line of |
|||
the program below. |
|||
Note also that the function `bernoulli` is defined here in a way that ensures B(1) is `1 // 2`. |
|||
(*) The C implementation of jq does not have sufficient numeric precision for the "extra credit" task. |
|||
<lang jq>include "Rational"; |
|||
# Preliminaries |
|||
def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .; |
|||
# for gojq |
|||
def idivide($j): |
|||
. as $i |
|||
| ($i % $j) as $mod |
|||
| ($i - $mod) / $j ; |
|||
# use idivide for precision |
|||
def binomial(n; k): |
|||
if k > n / 2 then binomial(n; n-k) |
|||
else reduce range(1; k+1) as $i (1; . * (n - $i + 1) | idivide($i)) |
|||
end; |
|||
# Here we conform to the modern view that B(1) is 1 // 2 |
|||
def bernoulli: |
|||
if type != "number" or . < 0 then "bernoulli must be given a non-negative number vs \(.)" | error |
|||
else . as $n |
|||
| reduce range(0; $n+1) as $i ([]; |
|||
.[$i] = r(1; $i + 1) |
|||
| reduce range($i; 0; -1) as $j (.; |
|||
.[$j-1] = rmult($j; rminus(.[$j-1]; .[$j])) ) ) |
|||
| .[0] # the modern view |
|||
end; |
|||
# Input: a non-negative integer, $p |
|||
# Output: an array of Rationals corresponding to the |
|||
# Faulhaber coefficients for row ($p + 1) (counting the first row as row 1). |
|||
def faulhabercoeffs: |
|||
def altBernoulli: # adjust B(1) for this task |
|||
bernoulli as $b |
|||
| if . == 1 then rmult(-1; $b) else $b end; |
|||
. as $p |
|||
| r(1; $p + 1) as $q |
|||
| { coeffs: [], sign: -1 } |
|||
| reduce range(0; $p+1) as $j (.; |
|||
.sign *= -1 |
|||
| binomial($p + 1; $j) as $b |
|||
| .coeffs[$p - $j] = ([ .sign, $q, $b, ($j|altBernoulli) ] | rmult)) |
|||
| .coeffs |
|||
; |
|||
# Calculate the sum for ($k|faulhabercoeffs) |
|||
def faulhabersum($n; $k): |
|||
($k|faulhabercoeffs) as $coe |
|||
| reduce range(0;$k+1) as $i ({sum: 0, power: 1}; |
|||
.power *= $n |
|||
| .sum = radd(.sum; rmult(.power; $coe[$i])) |
|||
) |
|||
| .sum; |
|||
# pretty print a Rational assumed to have the {n,d} form |
|||
def rpp: |
|||
if .n == 0 then "0" |
|||
elif .d == 1 then .n | tostring |
|||
else "\(.n)/\(.d)" |
|||
end; |
|||
def testfaulhaber: |
|||
(range(0; 10) as $i |
|||
| ($i | faulhabercoeffs | map(rpp | lpad(6)) | join(" "))), |
|||
"\nfaulhabersum(1000; 17):", |
|||
(faulhabersum(1000; 17) | rpp) ; |
|||
testfaulhaber</lang> |
|||
{{out}} |
|||
<pre> |
|||
1 |
|||
1/2 1/2 |
|||
1/6 1/2 1/3 |
|||
0 1/4 1/2 1/4 |
|||
-1/30 0 1/3 1/2 1/5 |
|||
0 -1/12 0 5/12 1/2 1/6 |
|||
1/42 0 -1/6 0 1/2 1/2 1/7 |
|||
0 1/12 0 -7/24 0 7/12 1/2 1/8 |
|||
-1/30 0 2/9 0 -7/15 0 2/3 1/2 1/9 |
|||
0 -3/20 0 1/2 0 -7/10 0 3/4 1/2 1/10 |
|||
faulhabersum(1000; 17): |
|||
56056972216555580111030077961944183400198333273050000 |
|||
</pre> |
|||
=={{header|Julia}}== |
=={{header|Julia}}== |