DataMuseum.dk

top - download

    Length: 5701 (0x1645)
    Types: TextFile
    Names: »prodlist.c«

└─⟦a0efdde77⟧ Bits:30001252 EUUGD11 Tape, 1987 Spring Conference Helsinki
    └─ ⟦this⟧ »EUUGD11/stat-5.3/eu/stat/src/prodlist.c« 

TextFile

/*  Copyright 1986 Gary Perlman */

/*LATTICE BUG
	prodlist did not compile under that Lattice large model.
	It was necessary to rewrite the routines to use a simple
	array, not a structure, and keep the list size in a global.
*/

/*LINTLIBRARY*/

#include "prodlist.h"
#include <math.h>

#ifdef	TRACE
#include <stdio.h>
#endif	TRACE

#ifndef	NULL
#define	NULL 0
#endif	NULL

#ifndef MSDOS
static char sccsfid[] = "@(#) prodlist.c 5.1 (|stat) 12/22/86";
#endif

\f

/*FUNCTION prod_new: create a new product list */
PLIST *
prod_new (n)
int 	n;
	{
	PLIST	*new = (PLIST *) malloc (sizeof (PLIST));
	if (new)
		{
		new->n = n;
		new->power = (int  *) malloc (((unsigned) n + 1) * sizeof (int ));
		if (new->power == NULL)
			new = NULL;
		else
			prod_init (new);
		}
	return (new);
	}

\f

/*FUNCTION prod_rel: release a product list */
void
prod_rel (list)
PLIST	*list;
	{
	if (list)
		{
		if (list->power)
			free ((char *) list->power);
		free ((char *) list);
		}
	}

\f

/*FUNCTION prod_init: set all exponents in a list to 0 */
void
prod_init (list)
PLIST	*list;
	{
	register	int 	i;
	for (i = prod_n(list); i >= 2; i--)
		prod_set (list, i, 0);
	}

\f

/*FUNCTION prod_fact: add an expanded factorial to a list */
void
prod_fact (list, fact, sign)
register PLIST	*list;
register	int 	fact;    /* factorial */
int 	sign;                /* if < 0, then divide, else multiply */
	{
	if (sign < 0)
		while (fact >= 2)
			{
			prod_div (list, fact);
			fact--;
			}
	else
		while (fact >= 2)
			{
			prod_mult (list, fact);
			fact--;
			}
	}

\f

/*FUNCTION prod_compute: return the computed value of a product list */
double
prod_compute (list)
PLIST 	*list;
	{
	static	double	loghuge = 0.0;    /* HUGE is max double number */
	int 	p;                        /* a specific power exponent */
	register	int 	power;        /* loop through all the powers */
	double	logresult = 0.0;          /* = sum of logs of factors */
	double	log ();                   /* standard log function */
	double	logprime ();              /* returns logs of primes */
	
	if (loghuge == 0.0)               /* set static over/underflow value */
		loghuge = log (HUGE);
	
	/*ALGORITHM
		First we reduce the product list to a product of powers of primes.
		This is done to reduce the number of terms.
		Then we
			add the (product of the exponents times the logs of the primes)
		to obtain the
			log of the (product of the powers of the primes).
		We then compare the logsum to underflow and overflow limits,
		and then return the exponential of the logsum.
	*/
	primefactor (list);
	for (power = prod_n (list); power >= 2; power--)
		if (p = prod_get (list, power))
			logresult += logprime (power) * p;
	if (logresult <= -loghuge)
		return (0.0);
	if (logresult >= loghuge)
		return (HUGE); /* overflow */
	return (exp (logresult));
	}

\f

/*FUNCTION primefactor: remove composite factors of a product list */
/*
	shift the representation of each power count
	to smaller numbers, so that, for example,
	the powers of 6 turn into powers of 2 and 3.
	
	Note that the list is reduced in reverse order, so that if we
	find that an index has factors, and one of them is composite,
	the that composite will be factored later.  We only need to check
	down to the minimum composite number: 4 because 3 and 2
	will not be factorable.
	
	After the procedure is done, the list is a product of powers of primes.
*/
static
primefactor (list)
register	PLIST	*list;
	{
	static	int 	divn = 0;            /* number of divisors set */
	static	int 	smalldiv[MAXN+1];    /* smallest divisor of numbers */
	static	int 	largediv[MAXN+1];    /* largediv[i] = i / smalldiv[i] */
	register	int 	i;
	int 	factor1, factor2;
	register	int 	p;
	void	dofactor ();
	
#ifdef	TRACE
	printlist ("shift  in: ", list);
#endif
	
	/* factor down to minimum composite number: 4 */
	if (prod_n (list) > divn)
		{
		dofactor (smalldiv, largediv, divn+1, prod_n (list));
		divn = prod_n (list);
		}
	for (i = prod_n (list); i >= 4; i--)       /* reverse order to retry */
		{
		if ((p = prod_get (list, i)) && (factor1 = smalldiv[i]) > 1)
			{
			factor2 = largediv[i];
			prod_pow (list, factor1, p);
			prod_pow (list, factor2, p);
			prod_set (list, i, 0);
			}
		}
#ifdef	TRACE
	printlist ("shift out: ", list);
#endif
	}

\f

/*FUNCTION printlist: print a product list */
#ifdef	TRACE
printlist (msg, list)
char	*msg;
PLIST	*list;
	{
	int 	i;
	int 	power;
	int 	count = 0;

	printf ("%s", msg);
	for (i = 2; i <= prod_n (list); i++)
		{
		count += abs (prod_get (list, i));
		if (prod_get (list, i))
			printf ("%d|%d ", i, prod_get (list, i));
		}
	printf ("(%d)\n", count);
	}
#endif	TRACE

\f


#ifdef	STANDALONE

#include <stdio.h>
#define	N 1000
main ()
	{
	char	line[BUFSIZ];
	PLIST	*list = prod_new (N+1);
	int 	small[N+1];
	int 	big[N+1];
	int 	i;
	int 	n;

	prod_init (list);
	while (gets (line))
		{
		i = atoi (line+1);
		switch (*line)
			{
			default:
				puts ("[*/!?nr] N   or =");
				break;
			case '*': prod_mult (list, i);        break;
			case '/': prod_div (list, i);         break;
			case '!': prod_fact (list, i, 1);     break;
			case '?': prod_fact (list, i, -1);    break;
			case '=':
				printf ("	%f\n", prod_compute (list));
				break;
			case '0': prod_init (list);           break;
			case 'n': printf ("n = %d\n", n = i); break;
			case 'r':
				printf ("n = %d\n", n);
				printf ("r = %d\n", i);
				prod_init (list);
				prod_fact (list, n, 1);
				prod_fact (list, n-i, -1);
				printf ("permutations = %f\n", prod_compute (list));
				prod_fact (list, i, -1);
				printf ("combinations = %f\n", prod_compute (list));
				prod_pow (list, 2, -n);
				printf ("binomial     = %f\n", prod_compute (list));
				break;
			}
		}
	}

#endif	STANDALONE