/*
 * bernoulli - clculate the Nth Bernoulli number B(n)
 *
 * Copyright (C) 1999  David I. Bell
 *
 * Calc is open software; you can redistribute it and/or modify it under
 * the terms of the version 2.1 of the GNU Lesser General Public License
 * as published by the Free Software Foundation.
 *
 * Calc is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
 * Public License for more details.
 *
 * A copy of version 2.1 of the GNU Lesser General Public License is
 * distributed with calc under the filename COPYING-LGPL.  You should have
 * received a copy with calc; if not, write to Free Software Foundation, Inc.
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 *
 * @(#) $Revision: 29.1 $
 * @(#) $Id: bernoulli.cal,v 29.1 1999/12/14 09:15:30 chongo Exp $
 * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/bernoulli.cal,v $
 *
 * Under source code control:	1991/09/30 11:18:41
 * File existed as early as:	1991
 *
 * Share and enjoy!  :-)	http://reality.sgi.com/chongo/tech/comp/calc/
 */

/*
 * Calculate the Nth Bernoulli number B(n).
 * This uses the following symbolic formula to calculate B(n):
 *
 *	(b+1)^(n+1) - b^(n+1) = 0
 *
 * where b is a dummy value, and each power b^i gets replaced by B(i).
 * For example, for n = 3:
 *	(b+1)^4 - b^4 = 0
 *	b^4 + 4*b^3 + 6*b^2 + 4*b + 1 - b^4 = 0
 *	4*b^3 + 6*b^2 + 4*b + 1 = 0
 *	4*B(3) + 6*B(2) + 4*B(1) + 1 = 0
 *	B(3) = -(6*B(2) + 4*B(1) + 1) / 4
 *
 * The combinatorial factors in the expansion of the above formula are
 * calculated interatively, and we use the fact that B(2i+1) = 0 if i > 0.
 * Since all previous B(n)'s are needed to calculate a particular B(n), all
 * values obtained are saved in an array for ease in repeated calculations.
 */


static Bnmax;
static mat Bn[1001];

define B(n)
{
	local	nn, np1, i, sum, mulval, divval, combval;

	if (!isint(n) || (n < 0))
		quit "Non-negative integer required for Bernoulli";

	if (n == 0)
		return 1;
	if (n == 1)
		return -1/2;
	if (isodd(n))
		return 0;
	if (n > 1000)
		quit "Very large Bernoulli";

	if (n <= Bnmax)
		return Bn[n];

	for (nn = Bnmax + 2; nn <= n; nn+=2) {
		np1 = nn + 1;
		mulval = np1;
		divval = 1;
		combval = 1;
		sum = 1 - np1 / 2;
		for (i = 2; i < np1; i+=2) {
			combval = combval * mulval-- / divval++;
			combval = combval * mulval-- / divval++;
			sum += combval * Bn[i];
		}
		Bn[nn] = -sum / np1;
	}
	Bnmax = n;
	return Bn[n];
}