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