Faulhaber's formula: Difference between revisions

=={{header|C}}==

{{trans|Modula-2}}

#include <stdio.h>

#include <stdlib.h>

 

return 0;

{{out}}

<pre>0 : n

=={{header|C sharp|C#}}==

{{trans|Java}}

 

namespace FaulhabersFormula {

}

}

{{out}}

<pre>0 : n

{{trans|D}}

Uses C++17

#include <numeric>

#include <sstream>

 

return 0;

{{out}}

<pre>0 : n

=={{header|D}}==

{{trans|Kotlin}}

import std.exception : enforce;

import std.format : formattedWrite;

faulhaber(i);

}

{{out}}

<pre>0 : n

 

=={{header|EchoLisp}}==

(lib 'math) ;; for bernoulli numbers

(string-delimiter "")

(define (Faulcomp n p)

(printf "Σ(1..%d) n^%d = %d" n p (/ (poly n (Faulhaber p)) (1+ p) )))

{{out}}

<pre>

 

=={{header|Factor}}==

math.functions regexp sequences ;

 

: poly>str ( seq -- str ) (poly>str) clean-up ;

 

{{out}}

<pre>

9: 1/10n^10 + 1/2n^9 + 3/4n^8 - 7/10n^6 + 1/2n^4 - 3/20n^2

</pre>

 

=={{header|Fōrmulæ}}==

 

 

 

 

=={{header|GAP}}==

Straightforward implementation using GAP polynomials, and two different formulas: one based on Stirling numbers of the second kind (sum1, see Python implementation below in this page), and the usual Faulhaber formula (sum2). No optimization is made (one could compute Stirling numbers row by row, or the product in sum1 may be kept from one call to the other). Notice the Bernoulli term in the first formula is here only to correct the value of sum1(0), which is off by one because sum1 computes sums from 0 to n.

 

sum1 := p -> Sum([0 .. p], k -> Stirling2(p, k) * Product([0 .. k], j -> n + 1 - j) / (k + 1)) + 2 * Bernoulli(2 * p + 1);

sum2 := p -> Sum([0 .. p], j -> (-1)^j * Binomial(p + 1, j) * Bernoulli(j) * n^(p + 1 - j)) / (p + 1);

1/8*n^8+1/2*n^7+7/12*n^6-7/24*n^4+1/12*n^2

1/9*n^9+1/2*n^8+2/3*n^7-7/15*n^5+2/9*n^3-1/30*n

 

=={{header|Go}}==

 

import (

fmt.Println()

}

{{out}}

<pre>

=={{header|Groovy}}==

{{trans|Java}}

 

class FaulhabersFormula {

}

}

{{out}}

<pre>0 : n

=={{header|Haskell}}==

====Bernouilli polynomials====

import Data.List (intercalate, transpose)

import Data.Bifunctor (bimap)

unsignedLength xs =

let l = length xs

{{Out}}

<pre>0 -> n

Implementation:

 

+/,(<:*(_1^[)*!*(y^~1+[)%1+])"0/~i.1x+y

)

Faul=:adverb define

(0,|.(%m+1x) * (_1x&^ * !&(m+1) * Bfirst) i.1+m)&p.

 

Task example:

 

0 1x&p.

1 Faul

0 _1r30 0 2r9 0 _7r15 0 2r3 1r2 1r9&p.

9 Faul

 

Double checking our work:

 

9 Faul i.5

0 1 513 20196 282340

0 1 5 14 30

2 Fcheck i.5

 

=={{header|Java}}==

{{trans|Kotlin}}

{{works with|Java|8}}

import java.util.stream.IntStream;

 

}

}

{{out}}

<pre>0 : n

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|Julia}}==

 

'''Module''':

 

function bernoulli(n::Integer)

end

 

 

'''Main''':

 

{{out}}

=={{header|Kotlin}}==

As Kotlin doesn't have support for rational numbers built in, a cut-down version of the Frac class in the Arithmetic/Rational task has been used in order to express the polynomial coefficients as fractions.

 

fun gcd(a: Long, b: Long): Long = if (b == 0L) a else gcd(b, a % b)

fun main(args: Array<String>) {

for (i in 0..9) faulhaber(i)

 

{{out}}

=={{header|Lua}}==

{{trans|C}}

if n<0 or k<0 or n<k then return -1 end

if n==0 or k==0 then return 1 end

for i=0,9 do

faulhaber(i)

{{out}}

<pre>0 : n

 

=={{header|Mathematica}} / {{header|Wolfram Language}}==

Faulhaber[n_, 0] := n

Faulhaber[n_, p_] := n^(p + 1)/(p + 1) + 1/2 n^p + Sum[BernoulliB[k]/k! p!/(p - k + 1)! n^(p - k + 1), {k, 2, p}]

{{out}}

<pre>0 n

 

=={{header|Maxima}}==

sum2(p):=sum((-1)^j*binomial(p+1,j)*bern(j)*n^(p-j+1),j,0,p)/(p+1)$

 

[0,0,0,0,0,0,0,0,0,0]

 

{{out}}

<pre>

=={{header|Modula-2}}==

{{trans|C#}}

FROM EXCEPTIONS IMPORT AllocateSource,ExceptionSource,GetMessage,RAISE;

FROM FormatString IMPORT FormatString;

END;

ReadChar

{{out}}

<pre>0 : n

=={{header|Nim}}==

{{trans|Kotlin}}

 

type

 

for n in 0..9:

 

{{out}}

It's not worth using Bernoulli numbers in PARI/GP, because too much cleaning if you are avoiding "dirty" (but correct) result.<br>

Note: Find ssubstr() function here on RC.

\\ Using "Faulhaber's" formula based on Bernoulli numbers. aev 2/7/17

\\ In str string replace all occurrences of the search string ssrch with the replacement string srepl. aev 3/8/16

{\\ Testing:

for(i=0,9, Faulhaber2(i))}

{{Output}}

<pre>

This version is using the sums of pth powers formula from [[wp:Bernoulli_polynomials| Bernoulli polynomials]].

It has small, simple and clear code, and produces instant result.

\\ Using a formula based on Bernoulli polynomials. aev 2/5/17

Faulhaber1(m)={

{\\ Testing:

for(i=0,9, Faulhaber1(i))}

{{Output}}

<pre>

> ##

*** last result computed in 0 ms

</pre>

 

=={{header|Perl}}==

use Math::Algebra::Symbols;

 

foreach my $i (0 .. 9) {

say "$i: ", faulhaber_s_formula($i);

{{out}}

<pre>

=={{header|Phix}}==

{{trans|C#}}

{{out}}

<pre>

 

 

 

def nextu(a):

 

for i, p in enumerate(sumpol(10)):

 

{{out}}

Racket will simplify rational numbers; if this code simplifies the expressions too much for your tastes (e.g. you like <code>1/1 * (n)</code>) then tweak the simplify... clauses to taste.

 

 

(require racket/match

(for ((p (in-range 0 (add1 9))))

(printf "f(~a) = ~a~%" p (expression->infix-string (faulhaber p))))

 

{{out}}

(formerly Perl 6)

{{works with|Rakudo|2018.04.01}}

 

return 1/2 if $n == 1;

for 0..9 -> $p {

say "f($p) = ", faulhaber_s_formula($p);

{{out}}

<pre>

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|Ruby}}==

{{trans|C}}

if n < 0 or k < 0 or n < k then

return -1

end

 

{{out}}

<pre>0 : n

=={{header|Scala}}==

{{trans|Java}}

 

abstract class Frac extends Comparable[Frac] {

}

}

{{out}}

<pre>0 : n

 

=={{header|Sidef}}==

(p+1).of { |j|

}

 

for p in (^10) {

printf("%2d: %s\n", p, faulhaber_formula(p))

{{out}}

<pre>

</pre>

 

=={{header|Visual Basic .NET}}==

{{trans|C#}}

Function Gcd(a As Long, b As Long)

If b = 0 Then

Next

End Sub

{{out}}

<pre>0 : n

{{libheader|Wren-math}}

{{libheader|Wren-rat}}

 

var bernoulli = Fn.new { |n|

}

 

 

{{out}}

{{libheader|GMP}} GNU Multiple Precision Arithmetic Library

Uses code from the Bernoulli numbers task (copied here).

 

fcn faulhaberFormula(p){ //-->(Rational,Rational...)

[p..0,-1].pump(List(),'wrap(k){ B(k)*BN(p+1).binomial(k) })

.apply('*(Rational(1,p+1)))

println("F(%d) --> %s".fmt(p,polyRatString(faulhaberFormula(p))))

fcn init(_a,_b){ var a=_a, b=_b; normalize(); }

fcn toString{

}

fcn __opDiv(n){ self(a*n.b,b*n.a) } // Rat / Rat

var A=List.createLong(100,0); // aka static aka not thread safe

foreach m in (N+1){

if(str[0]=="+") str[1,*]; // leave leading space

else String("-",str[2,*]);

{{out}}

<pre>