1/*
2 * bernoulli - calculate the Nth Bernoulli number B(n)
3 *
4 * Copyright (C) 2000,2021  David I. Bell and Landon Curt Noll
5 *
6 * Calc is open software; you can redistribute it and/or modify it under
7 * the terms of the version 2.1 of the GNU Lesser General Public License
8 * as published by the Free Software Foundation.
9 *
10 * Calc is distributed in the hope that it will be useful, but WITHOUT
11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
13 * Public License for more details.
14 *
15 * A copy of version 2.1 of the GNU Lesser General Public License is
16 * distributed with calc under the filename COPYING-LGPL.  You should have
17 * received a copy with calc; if not, write to Free Software Foundation, Inc.
18 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
19 *
20 * Under source code control:	1991/09/30 11:18:41
21 * File existed as early as:	1991
22 *
23 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
24 */
25
26/*
27 * Calculate the Nth Bernoulli number B(n).
28 *
29 * NOTE: This is now a builtin function.
30 *
31 * The non-builtin code used the following symbolic formula to calculate B(n):
32 *
33 *	(b+1)^(n+1) - b^(n+1) = 0
34 *
35 * where b is a dummy value, and each power b^i gets replaced by B(i).
36 * For example, for n = 3:
37 *
38 *	(b+1)^4 - b^4 = 0
39 *	b^4 + 4*b^3 + 6*b^2 + 4*b + 1 - b^4 = 0
40 *	4*b^3 + 6*b^2 + 4*b + 1 = 0
41 *	4*B(3) + 6*B(2) + 4*B(1) + 1 = 0
42 *	B(3) = -(6*B(2) + 4*B(1) + 1) / 4
43 *
44 * The combinatorial factors in the expansion of the above formula are
45 * calculated interactively, and we use the fact that B(2i+1) = 0 if i > 0.
46 * Since all previous B(n)'s are needed to calculate a particular B(n), all
47 * values obtained are saved in an array for ease in repeated calculations.
48 */
49
50
51/*
52static Bnmax;
53static mat Bn[1001];
54*/
55
56define B(n)
57{
58/*
59	local	nn, np1, i, sum, mulval, divval, combval;
60
61	if (!isint(n) || (n < 0))
62		quit "Non-negative integer required for Bernoulli";
63
64	if (n == 0)
65		return 1;
66	if (n == 1)
67		return -1/2;
68	if (isodd(n))
69		return 0;
70	if (n > 1000)
71		quit "Very large Bernoulli";
72
73	if (n <= Bnmax)
74		return Bn[n];
75
76	for (nn = Bnmax + 2; nn <= n; nn+=2) {
77		np1 = nn + 1;
78		mulval = np1;
79		divval = 1;
80		combval = 1;
81		sum = 1 - np1 / 2;
82		for (i = 2; i < np1; i+=2) {
83			combval = combval * mulval-- / divval++;
84			combval = combval * mulval-- / divval++;
85			sum += combval * Bn[i];
86		}
87		Bn[nn] = -sum / np1;
88	}
89	Bnmax = n;
90	return Bn[n];
91*/
92	return bernoulli(n);
93}
94