Faulhaber's formula: Difference between revisions
Added FreeBASIC
(→{{header|Sidef}}: switched to built-in polynomials) Tag: Made through Tor |
(Added FreeBASIC) |
||
| (9 intermediate revisions by 5 users not shown) | |||
Line 826:
9: 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2
</pre>
=={{header|FreeBASIC}}==
{{trans|C}}
<syntaxhighlight lang="vbnet">Type Fraction
As Integer num
As Integer den
End Type
Function Binomial(n As Integer, k As Integer) As Integer
If n < 0 Or k < 0 Then Print "Arguments must be non-negative integers": Exit Function
If n < k Then Print "The second argument cannot be more than the first.": Exit Function
If n = 0 Or k = 0 Then Return 1
Dim As Integer i, num, den
num = 1
For i = k + 1 To n
num *= i
Next i
den = 1
For i = 2 To n - k
den *= i
Next i
Return num / den
End Function
Function GCD(n As Integer, d As Integer) As Integer
Return Iif(d = 0, n, GCD(d, n Mod d))
End Function
Function makeFrac(n As Integer, d As Integer) As Fraction
Dim As Fraction result
Dim As Integer g
If d = 0 Then
result.num = 0
result.den = 0
Return result
End If
If n = 0 Then
d = 1
Elseif d < 0 Then
n = -n
d = -d
End If
g = Abs(gcd(n, d))
If g > 1 Then
n /= g
d /= g
End If
result.num = n
result.den = d
Return result
End Function
Function negateFrac(f As Fraction) As Fraction
Return makeFrac(-f.num, f.den)
End Function
Function subFrac(lhs As Fraction, rhs As Fraction) As Fraction
Return makeFrac(lhs.num * rhs.den - lhs.den * rhs.num, rhs.den * lhs.den)
End Function
Function multFrac(lhs As Fraction, rhs As Fraction) As Fraction
Return makeFrac(lhs.num * rhs.num, lhs.den * rhs.den)
End Function
Function equalFrac(lhs As Fraction, rhs As Fraction) As Integer
Return (lhs.num = rhs.num) And (lhs.den = rhs.den)
End Function
Function lessFrac(lhs As Fraction, rhs As Fraction) As Integer
Return (lhs.num * rhs.den) < (rhs.num * lhs.den)
End Function
Sub printFrac(f As Fraction)
Print Str(f.num);
If f.den <> 1 Then Print "/" & f.den;
End Sub
Function Bernoulli(n As Integer) As Fraction
If n < 0 Then Print "Argument must be non-negative": Exit Function
Dim As Fraction a(16)
Dim As Integer j, m
If (n < 0) Then
a(0).num = 0
a(0).den = 0
Return a(0)
End If
For m = 0 To n
a(m) = makeFrac(1, m + 1)
For j = m To 1 Step -1
a(j - 1) = multFrac(subFrac(a(j - 1), a(j)), makeFrac(j, 1))
Next j
Next m
If n <> 1 Then Return a(0)
Return negateFrac(a(0))
End Function
Sub Faulhaber(p As Integer)
Dim As Fraction coeff, q
Dim As Integer j, pwr, sign
Print p & " : ";
q = makeFrac(1, p + 1)
sign = -1
For j = 0 To p
sign = -1 * sign
coeff = multFrac(multFrac(multFrac(q, makeFrac(sign, 1)), makeFrac(Binomial(p + 1, j), 1)), Bernoulli(j))
If (equalFrac(coeff, makeFrac(0, 1))) Then Continue For
If j = 0 Then
If Not equalFrac(coeff, makeFrac(1, 1)) Then
If equalFrac(coeff, makeFrac(-1, 1)) Then
Print "-";
Else
printFrac(coeff)
End If
End If
Else
If equalFrac(coeff, makeFrac(1, 1)) Then
Print " + ";
Elseif equalFrac(coeff, makeFrac(-1, 1)) Then
Print " - ";
Elseif lessFrac(makeFrac(0, 1), coeff) Then
Print " + ";
printFrac(coeff)
Else
Print " - ";
printFrac(negateFrac(coeff))
End If
End If
pwr = p + 1 - j
Print Iif(pwr > 1, "n^" & pwr, "n");
Next j
Print
End Sub
For i As Integer = 0 To 9
Faulhaber(i)
Next i
Sleep</syntaxhighlight>
{{out}}
<pre>Same as C entry.</pre>
=={{header|Fōrmulæ}}==
{{FormulaeEntry|page=https://formulae.org/?script=examples/Faulhaber}}
'''Solution'''
The following function creates the Faulhaber's coefficients up to a given number of rows, according to the [http://www.ww.ingeniousmathstat.org/sites/default/files/Torabi-Dashti-CMJ-2011.pdf paper] of of Mohammad Torabi Dashti:
(This is exactly the as than the task [[Faulhaber%27s_triangle#F%C5%8Drmul%C3%A6|Faulhaber's triangle]])
[[File:Fōrmulæ - Faulhaber 01.png]]
The following function creates the sum of the p-th powers of the first n positive integers as a (p + 1)th-degree polynomial function of n:
Notes. The -1 index means the last element (-2 is the penultimate element, and so on). So it retrieves the last row of the triangle. |x| is the cardinality (number of elements) of x.
(This is exactly the as than the task [[Faulhaber%27s_triangle#F%C5%8Drmul%C3%A6|Faulhaber's triangle]])
This function can be used for both symbolic or numeric computation of the polynomial:
[[File:Fōrmulæ - Faulhaber 06.png]]
To generate the first 10 closed-form expressions, starting with p = 0:
[[File:Fōrmulæ - Faulhaber 07.png]]
[[File:Fōrmulæ - Faulhaber 08.png]]
=={{header|GAP}}==
Line 1,411 ⟶ 1,581:
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}}==
Line 2,177 ⟶ 2,444:
> ##
*** last result computed in 0 ms
</pre>
=={{header|Pascal}}==
A console program that runs under Lazarus or Delphi. Does not make use of Bernoulli numbers.
<syntaxhighlight lang="pascal">
program Faulhaber;
{$IFDEF FPC} // Lazarus
{$MODE Delphi} // ensure Lazarus accepts Delphi-style code
{$ASSERTIONS+} // by default, Lazarus does not compile 'Assert' statements
{$ELSE} // Delphi
{$APPTYPE CONSOLE}
{$ENDIF}
uses SysUtils;
type TRational = record
Num, Den : integer; // where Den > 0 and Num, Den are coprime
end;
const
ZERO : TRational = ( Num: 0; Den : 1);
HALF : TRational = ( Num: 1; Den : 2);
// Construct rational a/b, assuming b > 0.
function Rational( const a, b : integer) : TRational;
var
t, x, y : integer;
begin
if b <= 0 then raise SysUtils.Exception.Create( 'Denominator must be > 0');
// Find HCF of a and b (Euclid's algorithm) and cancel it out.
x := Abs(a);
y := b;
while y <> 0 do begin
t := x mod y;
x := y;
y := t;
end;
result.Num := a div x;
result.Den := b div x
end;
function Prod( r, s : TRational) : TRational; // result := r*s
begin
result := Rational( r.Num*s.Num, r.Den*s.Den);
end;
procedure DecRat( var r : TRational;
const s : TRational); // r := r - s
begin
r := Rational( r.Num*s.Den - s.Num*r.Den, r.Den * s.Den);
end;
// Write a term such as ' - (7/10)n^6' to the console.
procedure WriteTerm( coeff : TRational;
index : integer;
printPlus : boolean);
begin
if Coeff.Num = 0 then exit;
with coeff do begin
if Num < 0 then Write(' - ')
else if printPlus then Write(' + ');
// Put brackets round a fractional coefficient
if (Den > 1) then Write('(');
// If coefficient is 1, don't write it
if (Den > 1) or (Abs(Num) > 1) then Write( Abs(Num));
// Write denominator if it's not 1
if (Den > 1) then Write('/', Den, ')');
end;
Write('n');
if index > 1 then Write('^', index);
end;
{-------------------------------------------------------------------------------
Main routine. Calculation of Faulhaber polynomials
F_p(n) = 1^p + 2^p + ... + n^p, p = 0, 1, ..., p_max
}
var
p_max : integer;
c : array of array of TRational;
i, j, p : integer;
coeff_of_n : TRational;
begin
// User types program name, optionally followed by maximum power p (defaults to 9)
if ParamCount = 0 then p_max := 9
else p_max := SysUtils.StrToInt( ParamStr(1));
// c[p, i] is coefficient of n^i in the polynomial F_p(n).
// Initialize all coefficients to 0.
SetLength( c, p_max + 1, p_max + 2);
for i := 0 to p_max do
for j := 0 to p_max + 1 do
c[i, j] := ZERO;
c[0, 1] := Rational(1, 1); // F_0(n) = n, special case
for p := 1 to p_max do begin
// Initialize calculation of coefficient of n, needed if p is even.
// If p is odd, still calculate it as a check on the working (should be 0).
// Calculation uses the fact that F_p(1) = 1.
coeff_of_n := Rational(1, 1);
c[p, p+1] := Rational(1, p + 1);
DecRat( coeff_of_n, c[p, p + 1]);
c[p, p] := HALF;
DecRat( coeff_of_n, c[p, p]);
i := p - 1;
while (i >= 2) do begin
c[p, i] := Prod( Rational(p, i), c[p - 1, i - 1]);
DecRat( coeff_of_n, c[p, i]);
dec(i, 2);
end;
if i = 1 then // p is even
c[p, 1] := coeff_of_n // = the Bernoulli number B_p
else // p is odd
Assert( coeff_of_n.Num = 0); // just checking
end; // for p
// Print the result
for p := 0 to p_max do begin
Write( 'F_', p, '(n) = ');
for j := p + 1 downto 1 do WriteTerm( c[p, j], j, j <= p);
WriteLn;
end;
end.
</syntaxhighlight>
{{out}}
<pre>
F_0(n) = n
F_1(n) = (1/2)n^2 + (1/2)n
F_2(n) = (1/3)n^3 + (1/2)n^2 + (1/6)n
F_3(n) = (1/4)n^4 + (1/2)n^3 + (1/4)n^2
F_4(n) = (1/5)n^5 + (1/2)n^4 + (1/3)n^3 - (1/30)n
F_5(n) = (1/6)n^6 + (1/2)n^5 + (5/12)n^4 - (1/12)n^2
F_6(n) = (1/7)n^7 + (1/2)n^6 + (1/2)n^5 - (1/6)n^3 + (1/42)n
F_7(n) = (1/8)n^8 + (1/2)n^7 + (7/12)n^6 - (7/24)n^4 + (1/12)n^2
F_8(n) = (1/9)n^9 + (1/2)n^8 + (2/3)n^7 - (7/15)n^5 + (2/9)n^3 - (1/30)n
F_9(n) = (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>
Line 2,534 ⟶ 2,938:
f(8) = 1/9 * (n^9 + (9/2)*n^8 + (6/1)*n^7 + (-21/5)*n^5 + (2/1)*n^3 + (-3/10)*n)
f(9) = 1/10 * (n^10 + (5/1)*n^9 + (15/2)*n^8 + (-7/1)*n^6 + (5/1)*n^4 + (-3/2)*n^2)
</pre>
=={{header|RPL}}==
{{works with|HP|49}}
≪ → p
≪ 0
p 0 '''FOR''' m
p 1 + m 1 + COMB
p m - IBERNOULLI
'''IF''' LASTARG 1 == '''THEN''' NEG '''END'''
* EVAL + 'n' *
-1 '''STEP'''
p 1 + / EXPAN
≫ ≫ '<span style="color:blue>FAULH</span>' STO
≪ P <span style="color:blue>FAULH</span> ≫ 'P' 0 9 1 SEQ
{{out}}
<pre>
1: {'n' '1/2n^2+1/2n' '1/3n^3+1/2n^2+1/6n' '1/4n^4+1/2n^3+1/4n^2' '1/5n^5+1/2n^4+1/3n^3-1/30n' '1/6n^6+1/2n^5+5/12n^4-1/12n^2' '1/7n^7+1/2n^6+1/2n^5-1/6n^3+1/42n' '1/8n^8+1/2n^7+7/12n^6-7/24n^4+1/12n^2' '1/9n^9+1/2n^8+2/3n^7-7/15n^5+2/9n^3-1/30n' '1/10n^10+1/2n^9+3/4n^8-7/10n^6+1/2n^4-3/20n^2'}
</pre>
Line 3,028 ⟶ 3,451:
{{libheader|Wren-math}}
{{libheader|Wren-rat}}
<syntaxhighlight lang="
import "./rat" for Rat
var bernoulli = Fn.new { |n|
| |||