pi.c

calc大数库

/* pI.c */

#include <stdio.h>
#include <stdlib.h>
#include "integer.h"
#include "pI.h"
#include "fun.h"


#ifndef DEBUG

POLYI BANKI = NULL;

#endif

#ifndef DEBUG
void SAVEI(Ptr)
POLYI Ptr;
/*
This routine takes a POLYI (which is a linked list of TERMI's) that is
no longer needed, and links it onto the start of the free list BANKI,
which itself is just a linked list opf TERMI's.
*/
{
	POLYI tmp;

	tmp=Ptr;
	while(tmp->NEXT)
	  tmp=tmp->NEXT;
	tmp->NEXT = BANKI;
	BANKI = Ptr;
	Ptr = NULL;
}
#endif

POLYI CREATEPI()
/*
 * If there is anything in BANKI (the memory bank of TERMI's) then use that,
 * otherwise allocate a new one.
 * No initialisation is performed
*/
{
	POLYI Ptr;
	#ifndef DEBUG
	if (BANKI)
	{
		Ptr = BANKI;
		BANKI = BANKI->NEXT;
	}
	else
	#endif
	Ptr = (POLYI)mmalloc(sizeof(struct TERMI));
	return (Ptr);	
}

int DEGREEPI(POLYI Pptr)
/*
 * returns the degree of the polynomial Pptr.
 */
{
	if (Pptr == NULL)
		return (-1);
	else
		return (Pptr->DEG);
}

MPI *LEADPI(POLYI Pptr)
/*
 * *Aptr is the  leading coefficient of the polynomial Pptr.
 */
{
	if (Pptr == NULL)
		return ZEROI();
	else
		return COPYI(Pptr->COEF);
}

void DELETEPI(POLYI Pptr)
/*
 * returns back to the system storage dynamically allocated to the linked list
 * which Pptr points to.
 */
{
	POLYI Qtmp, orig;
	
	if(!Pptr) return;
	orig = Pptr;
	while (Pptr != NULL)
	{
 		Qtmp = Pptr->NEXT;
		FREEMPI(Pptr->COEF);
#ifdef DEBUG /* This fix was due to Peter Adams on 30/6/94 */
		SAVEPI(Pptr);
#endif
		Pptr= Qtmp;
	}
#ifndef DEBUG
	SAVEPI(orig);
#endif
	return;
}

void PINSERTPI(int DEG, MPI *COEF, POLYI *Pptr, USI OPT)
/*
 * inserts a copy of a single term *Qptr into polynomial *Pptr.
 * It has two options that control the insertion in case a term of *Pptr already
 * has the degree of *Qptr, If OPT = 0, then the term is replaced; if OPT = 1,
 * then the coefficient of *Qptr is added to that of the corresponding term of
 *  *Pptr.
 */
{
	POLYI CURSOR, Ttmp;
	MPI *TempI;

	if (*Pptr == NULL)
	{
		*Pptr = CREATEPI();
		(*Pptr)->DEG = DEG;
		(*Pptr)->COEF = COPYI(COEF);
		(*Pptr)->NEXT = NULL;
		return;
	}
	else
	{
		CURSOR = *Pptr;
		while ((CURSOR->DEG > DEG) && (CURSOR->NEXT != NULL))
			CURSOR = CURSOR->NEXT;
		if (CURSOR->DEG < DEG)
		{
			/* Qptr inserted before CURSOR */
			Ttmp = CREATEPI();
			Ttmp->DEG = CURSOR->DEG;
			Ttmp->COEF = CURSOR->COEF;
			Ttmp->NEXT = CURSOR->NEXT;
			CURSOR->DEG = DEG;
			CURSOR->COEF = COPYI(COEF);
			CURSOR->NEXT = Ttmp;
			return;
		}
		else if (CURSOR->DEG == DEG)
		{
			if (OPT == 0)
				CURSOR->COEF = COPYI(COEF);
			else
			{
				TempI = CURSOR->COEF;
				CURSOR->COEF = ADDI(CURSOR->COEF, COEF);
				FREEMPI(TempI);
			}
			return;
		}
		else
		{
			/* CURSOR->DEG > Qptr->DEG and CURSOR->NEXT === NULL */
			/* Qptr inserted at end */
			Ttmp = CREATEPI();
			Ttmp->DEG = DEG;
			Ttmp->COEF = COPYI(COEF);
			Ttmp->NEXT = NULL;
			CURSOR->NEXT = Ttmp;
			return;
		}
	}
}

void PURGEPI(POLYI *Pptr)
/*
 * gets rid of any zero coefficients in the polynomial Pptr.
 */
{
	POLYI R1tmp, R2tmp;

	if (*Pptr != NULL)
	{
		R1tmp = *Pptr;
		R2tmp = (*Pptr)->NEXT;
		while (R2tmp != NULL)
		{
			/* purge non-leading terms */
			if ((R2tmp->COEF)->S == 0)
			{
				R1tmp->NEXT = R2tmp->NEXT;
				R2tmp->NEXT = NULL;
				DELETEPI(R2tmp);
				R2tmp = R1tmp->NEXT;
			}
			else
			{
				R1tmp = R2tmp;
				R2tmp = R2tmp->NEXT;
			}
		}
		if (((*Pptr)->COEF)->S == 0)
		{
			/* purge leading term */
			R1tmp = *Pptr;
			*Pptr = (*Pptr)->NEXT;
			R1tmp->NEXT = NULL;
			DELETEPI(R1tmp);
		}
	}
	return;
}

POLYI COPYPI(POLYI Pptr)
/*
 * makes a copy of the polynomial Pptr.
 */
{
	POLYI Rtmp, Sptr;

	Sptr = NULL;
	Rtmp = Pptr;
	while (Rtmp != NULL)
	{
		PINSERTPI(Rtmp->DEG, Rtmp->COEF, &Sptr, 0);
		Rtmp = Rtmp->NEXT;
	}
	return (Sptr);
}
 
void PRINTPI(POLYI Pptr)
{
  FPRINTPI(stdout, Pptr);
}

void FPRINTPI(FILE *outfile, POLYI Pptr)
/*
 * outputs the terms of the polynomial Pptr.
 */
{
	POLYI CURSOR;
	int k, l;

	if (Pptr == NULL)
	{
		fprintf(outfile, "0");
		return;
	}
	else
	{
		CURSOR = Pptr;
		k = DEGREEPI(Pptr);
		while (CURSOR != NULL)
		{
			l = (CURSOR->COEF)->S;
			if (l == -1)
			{
				(CURSOR->COEF)->S = 1;
				fprintf(outfile, "- ");
			}
			else
			{
				if (CURSOR->DEG != k)
					fprintf(outfile, "+ ");
			}
			if (CURSOR->DEG == 0 || !EQONEI(CURSOR->COEF))
				FPRINTI(outfile, CURSOR->COEF);
			if (l == -1)
				(CURSOR->COEF)->S = -1;
			if (CURSOR->DEG == 1)
				fprintf(outfile, "X");
			if (CURSOR->DEG > 1)
				fprintf(outfile, "X^%d", CURSOR->DEG);
			if (CURSOR->NEXT != NULL)
				fprintf(outfile, " ");
			CURSOR = CURSOR->NEXT;
		}
		return;
	}
}



POLYI SCALARPI(MPI *Cptr, POLYI Pptr)
/*
 * multiplying the polynomial Pptr by the MPI *Cptr.
 */
{
	POLYI Sptr, Rtmp;
	MPI *TempI;

	if (Pptr == NULL || Cptr->S == 0)
		Sptr = NULL;
	else
	{
		Sptr = COPYPI(Pptr);
		Rtmp = Sptr;
		while (Rtmp != NULL)
		{
			TempI = Rtmp->COEF;
			Rtmp->COEF = MULTI(Rtmp->COEF, Cptr);
			FREEMPI(TempI);
			Rtmp = Rtmp->NEXT;
		}
	}
	return (Sptr);
}

POLYI ADDPI(POLYI Pptr, POLYI Qptr)
/*
 * returns the sum of two polynomials.
 */
{
	POLYI CURSOR, Sptr;

		Sptr = COPYPI(Pptr);
		CURSOR = Qptr;
		while (CURSOR != NULL)
		{
			PINSERTPI(CURSOR->DEG, CURSOR->COEF, &Sptr, 1);
			CURSOR = CURSOR->NEXT;
		}
		PURGEPI(&Sptr);
		return (Sptr);
}

POLYI SUBPI(POLYI Pptr, POLYI Qptr)
/*
 * returns the difference Pptr - Qptr of two polynomials.
 */
{
	POLYI CURSOR, Sptr;
	MPI *TempI;

		Sptr = COPYPI(Pptr);
		if (Qptr == NULL)
			return (Sptr);
		CURSOR = Qptr;
		while (CURSOR != NULL)
		{
			TempI = MINUSI(CURSOR->COEF);
			PINSERTPI(CURSOR->DEG, TempI, &Sptr, 1);
			FREEMPI(TempI);
			CURSOR = CURSOR->NEXT;
		}
		PURGEPI(&Sptr);
		return (Sptr);
}

POLYI MULTPI(POLYI Pptr, POLYI Qptr)
/*
 * returns the product of two polynomials.
 */
{
	POLYI Sptr, CURSOR_P, CURSOR_Q;
	MPI *TempI;
	int k;

	if (Pptr == NULL || Qptr == NULL)
		return (NULL);
	Sptr = NULL;
	CURSOR_P = Pptr;
	while (CURSOR_P != NULL)
	{
		CURSOR_Q = Qptr;
		while (CURSOR_Q != NULL)
		{
			k = CURSOR_P->DEG + CURSOR_Q->DEG;
			TempI = MULTI(CURSOR_P->COEF, CURSOR_Q->COEF);
			PINSERTPI(k, TempI, &Sptr, 1);
			FREEMPI(TempI);
			CURSOR_Q = CURSOR_Q->NEXT;
		}
		CURSOR_P = CURSOR_P->NEXT;
	}
	PURGEPI(&Sptr);
	return (Sptr);
}

MPI *VALPI(POLYI Pptr, MPI *Cptr)
/*
 * Evaluating the polynomial Pptr at the MPI *Cptr.
 */
{
	int k;
	MPI *B, *Aptr, *TempI;
	POLYI CURSOR;

	if (Pptr == NULL)
		return ZEROI();
	else
		Aptr = COPYI(Pptr->COEF);
	CURSOR = Pptr;
	while (CURSOR != NULL)
	{
		if (CURSOR->NEXT == NULL)
			k = DEGREEPI(CURSOR);
		else
			k = DEGREEPI(CURSOR) - DEGREEPI(CURSOR->NEXT);
		B = POWERI(Cptr, (unsigned int)k);
		TempI = Aptr;
		Aptr = MULTI(Aptr, B);
		FREEMPI(B);
		FREEMPI(TempI);
		CURSOR = CURSOR->NEXT;
		if (CURSOR != NULL)
		{
			TempI = Aptr;
			Aptr = ADDI(Aptr, CURSOR->COEF);
			FREEMPI(TempI);
		}
	}
	return (Aptr);
}

POLYI DIVPI(POLYI Fptr, POLYI Gptr)
/*
 * Pseudo-division of polynomials: see Knuth, Vol 2, p.407 
 * and H. Flanders, Scientific Pascal, p.510.
 */
{
	POLYI F, Tmp, Htmp, CURSOR, Qptr;
	int j, d;
	MPI *A, *G, *TempI;

	Qptr = NULL;
	d = DEGREEPI(Gptr);
	j = DEGREEPI(Fptr) - d;
	if (j < 0)
		return (Qptr);
	else
	{
		F = COPYPI(Fptr);
		G = LEADPI(Gptr);
		A = POWERI(G, (unsigned int)(j + 1));
		Tmp = F;
		F = SCALARPI(A, Tmp);
		FREEMPI(A);
		DELETEPI(Tmp);
		while ((j = DEGREEPI(F)) >= d)
		{
			j = j - d;
			A = LEADPI(F);
			TempI = A;
			A = INTI(A, G);
			FREEMPI(TempI);
			Htmp = COPYPI(Gptr);
			CURSOR = Htmp;
			do
			{
				CURSOR->DEG = j + CURSOR->DEG;
				TempI = CURSOR->COEF;
				CURSOR->COEF = MINUSI(CURSOR->COEF);
				FREEMPI(TempI);
				TempI = CURSOR->COEF;
				CURSOR->COEF = MULTI(CURSOR->COEF, A);
				FREEMPI(TempI);
				CURSOR = CURSOR->NEXT;
			}
			while (CURSOR != NULL);
			Tmp = F;
			F = ADDPI(Tmp, Htmp);
			DELETEPI(Tmp);
			DELETEPI(Htmp);
			PINSERTPI(j, A, &Qptr, 0);
			FREEMPI(A);
		}
		DELETEPI(F);
		FREEMPI(G);
		return (Qptr);
	}
}

POLYI MODPI(POLYI Fptr, POLYI Gptr)
/*
 * Pseudo-division of polynomials: see Knuth, Vol 2, p.407 
 * and H. Flanders, Scientific Pascal, p.510.
 * NOTE: This returns a polynomial that is a POSITIVE multiple of the true remainder.
 * i.e. The sign of the remainder is correct. 
 */
{
	POLYI F, Tmp, Htmp, CURSOR;
	int j, d, degDiff;
	MPI *A, *G, *Temp1I, *Temp2I;

	d = DEGREEPI(Gptr); /* n */
	j = DEGREEPI(Fptr) - d; /* m - n */
	degDiff=j;
	F = COPYPI(Fptr);	
	if (j < 0)
		return (F);
	else
	{
		G = LEADPI(Gptr);
		A = POWERI(G, (unsigned int)(j + 1));
		Tmp = F;
		F= SCALARPI(A, Tmp);
		FREEMPI(A);
		DELETEPI(Tmp);
		while ((j = DEGREEPI(F)) >= d)
		{
			j = j - d;
			A = LEADPI(F);
			if (!EQONEI(G))
			{
				Temp2I = A;
				A = INTI(A, G);
				FREEMPI(Temp2I);
			}
			Htmp = COPYPI(Gptr);
			CURSOR = Htmp;
			do
			{
				CURSOR->DEG = j + CURSOR->DEG;
				Temp1I = CURSOR->COEF;
				CURSOR->COEF = MINUSI(CURSOR->COEF);
				FREEMPI(Temp1I);
				Temp1I = CURSOR->COEF;
				CURSOR->COEF = MULTI(CURSOR->COEF, A);
				FREEMPI(Temp1I);
				CURSOR = CURSOR->NEXT;
			}
			while (CURSOR != NULL);
			Tmp = F;
			F= ADDPI(Tmp, Htmp);
			DELETEPI(Tmp);
			DELETEPI(Htmp);
			FREEMPI(A);
		}
		
		if ((degDiff % 2 == 0) && G->S == -1) {
		  POLYI TMPI;
		  TMPI = F;
		  F = MINUSPI(F);
		  DELETEPI(TMPI);
		  }	
		FREEMPI(G);
		return (F);
	}
}

MPI *CONTENTPI(POLYI Pptr)
/*
 * *Cptr is the (positive) content of the polynomial Pptr.
 */
{
	POLYI CURSOR;
	MPI *TempI, *Cptr;

	if (Pptr == NULL)
		return (ZEROI());
	else if (Pptr->DEG == 0)
	{
		Cptr = COPYI(Pptr->COEF);
		if (Cptr->S == -1)
			Cptr->S = 1;
		return (Cptr);
	}
	else
	{
		CURSOR = Pptr;
		Cptr = COPYI(CURSOR->COEF);
		if (CURSOR->NEXT != NULL)
		{
			do
			{
				TempI = Cptr;