ff2 = ff1 −
¶
0
1
ff1
(
x
)
dx
To find the nth Bernoulli polynomial requires finding all the (n-1) previous polynomials, so this is
time-consuming for higher-order polynomials. For this reason I wrote two different versions of
programs to find the Bernoulli polynomials. One version is a function which returns a single polynomial
of the order specified by the function argument. The second version builds a table of the polynomials
up to a specified order. Individual polynomials can be quickly retrieved from this table by another
function.
bnpolys(n) returns a single Bernoulli polynomial of order n as a function of x:
bnpolys(n)
Func
©(n) return Bernoulli polynomial of order n
©27jun00/dburkett@infinet.com
local k,f,g
x-1/2→f
if n=1:return f
for k,2,n
á(k*f,x)→g
g-á(g,x,0,1)→f
endfor
return f
EndFunc
bnpolys() may not return the polynomial in the form you want. For example, bnpolys(3) returns
x 2x
2
−3x+1
2
Use expand() to convert this result to
x
3
−
3x
2
2
+
x
2
or use comdenom() to get
2x
3
−3x
2
+x
2
bnpolys() slows down for large arguments. For example, bnpolys(20) takes about 22 seconds on my
HW2 92+ with AMS 2.04, and bnpolys(50) takes about 160 seconds. This long execution time occurs
because bnpolys() must calculate all the polynomials of order n-1, to return the polynomial of order n. If
you frequently use higher-order polynomials, it is worthwhile to build a table in advance, then use a
simple routine that just recalls the appropriate polynomial.
These routines perform those functions. bnpoly() builds the list of polynomials, and bpo() is a function
that returns a given polynomial, after bnpoly() is run.
First, use bnpoly() to build the list of polynomials:
bnpoly(nxx)
prgm
©(n) Fill Bernoulli polynomial list bpoly[] up to n.
6 - 41