Faulhaber's formula: Difference between revisions
→{{header|jq}}
(→{{header|Sidef}}: switched to built-in polynomials) Tag: Made through Tor |
|||
Line 1,411:
8 : 1/9n^9 + 1/2n^8 + 2/3n^7 - 7/15n^5 + 2/9n^3 - 1/30n
9 : 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2</pre>
=={{header|jq}}==
{{works with|jq}}
'''Works with gojq, the Go implementation of jq, and with fq'''
The following assumes the following rational arithmetic functions
as provided by the jq "rational" module at [[Arithmetic/Rational#jq]]:
* r/2 (constructor)
* requal/2 (equality)
* rlessthan/1
* rminus/1
* rminus/2 (subtraction)
* rmult/0 and rmult/2 (multiplication)
In addition, the pretty-printing function defined here (rpp) assumes
rationals are JSON objects of the form {n,d}
<syntaxhighlight lang=jq>
include "rational" {search: "."}; # see [[Arithmetic/Rational#jq]]:
# Preliminaries
# 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;
# 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;
# The following definition reflects 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;
# The task
def faulhaber($p):
# The traditional view:
def bernouilli($n):
$n | bernouilli | if $n==1 then rminus else . end;
r(1; $p+1) as $q
| { sign: -1 }
| reduce range(0; 1+$p) as $j (.;
.sign *= -1
| r(binomial($p+1; $j); 1) as $b
| ([$q, .sign, $b, ($j|bernoulli)] | rmult) as $coeff
| if requal($coeff; r(0;1))|not
then .emit +=
(if $j == 0
then (if requal($coeff; r( 1;1)) then ""
elif requal($coeff; r(-1;1)) then "-"
else "\($coeff|rpp) " end)
else (if requal($coeff; r(1;1)) then " + "
elif requal($coeff; r(-1;1)) then " - "
elif r(0;1)|rlessthan($coeff) then " + \($coeff|rpp) "
else " - \($coeff|rminus|rpp) "
end)
end)
| ($p + 1 - $j) as $pwr
| .emit += (if 1 < $pwr then "n^\($pwr)" else "n" end)
else .
end )
| .emit ;
range(0;10) | "\(.) : \(faulhaber(.))"
</syntaxhighlight>
{{output}}
<pre>
0 : n
1 : 1/2 n^2 - 1/2 n
2 : 1/3 n^3 - 1/2 n^2 + 1/6 n
3 : 1/4 n^4 - 1/2 n^3 + 1/4 n^2
4 : 1/5 n^5 - 1/2 n^4 + 1/3 n^3 - 1/30 n
5 : 1/6 n^6 - 1/2 n^5 + 5/12 n^4 - 1/12 n^2
6 : 1/7 n^7 - 1/2 n^6 + 1/2 n^5 - 1/6 n^3 + 1/42 n
7 : 1/8 n^8 - 1/2 n^7 + 7/12 n^6 - 7/24 n^4 + 1/12 n^2
8 : 1/9 n^9 - 1/2 n^8 + 2/3 n^7 - 7/15 n^5 + 2/9 n^3 - 1/30 n
9 : 1/10 n^10 - 1/2 n^9 + 3/4 n^8 - 7/10 n^6 + 1/2 n^4 - 3/20 n^2
</pre>
=={{header|Julia}}==
| |||