Faulhaber's triangle: Difference between revisions
→{{header|jq}}
m (Updated description and link for Fōrmulæ solution) |
|||
Line 1,680:
faulhaber(4, 1000)
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}}==
| |||