pi.c

calc大数库

				Cptr = GCD(Cptr, (CURSOR->NEXT)->COEF);
				FREEMPI(TempI);
				CURSOR = CURSOR->NEXT;
			}
			while (CURSOR->NEXT != NULL);
		}
		return (Cptr);
	}
}

MPI *CONTENTPI2(POLYI Pptr)
/* Returns content of Pptr such that content(Pptr)*primitive(Pptr) = Pptr */
{
  MPI *C, *C2, *L, *Z=ZEROI();
  C = CONTENTPI(Pptr);
  L=LEADPI(Pptr);
  if (COMPAREI(L,Z) < 0 && COMPAREI(C,Z) > 0) {
    C2 = MINUSI(C);
    FREEMPI(C);
    FREEMPI(L);
    FREEMPI(Z);
    return C2;
  } else {
    FREEMPI(L);
    FREEMPI(Z);
    return C;
  }
}
    
    
POLYI PRIMITIVEPI(POLYI Pptr)
/*
 * returns the primitive part of the polynomial Pptr.
 * leading coefficient is positive.
 */
{
	POLYI Tmp, Sptr, CURSOR;
	MPI *C, *TempI;

	if (Pptr == NULL)
		return (NULL);
	else
	{
		Sptr = COPYPI(Pptr);
		C = CONTENTPI(Sptr);
		if ((Sptr->COEF)->S == -1)
		{
			Tmp = Sptr;
			TempI = MINUS_ONEI();
			Sptr = SCALARPI(TempI, Tmp);
			FREEMPI(TempI);
			DELETEPI(Tmp);
		}
		CURSOR = Sptr;
		while (CURSOR != NULL)
		{
			TempI = CURSOR->COEF;
			CURSOR->COEF = INT(CURSOR->COEF, C);
			FREEMPI(TempI);
			CURSOR = CURSOR->NEXT;
		}
		FREEMPI(C);
		return (Sptr);
	}
}

POLYI GCDPI(POLYI Pptr, POLYI Qptr)
/*
 * returns the gcd of polynomials Pptr and Qptr.
 * see Knuth, Vol 2, p.403-408 and H. Flanders, Scientific Pascal, p.510.
 */
{
	POLYI Ptmp, Qtmp, Rtmp, Sptr;
	MPI *D, *P, *Q, *TempI;
	int j;
	
	if (Qptr == NULL)
	{
		if (Pptr == NULL)
			return (NULL);
		else
			return (COPYPI(Pptr));
	}
	if (Pptr == NULL)
	{
		if (Qptr == NULL)
			return (NULL);
		else
			return (COPYPI(Qptr));
	}
	P = CONTENTPI(Pptr);
	Q = CONTENTPI(Qptr);
	D = GCD(P, Q);
	FREEMPI(P);
	FREEMPI(Q);
	Ptmp = PRIMITIVEPI(Pptr);
	Qtmp = PRIMITIVEPI(Qptr);
        j = DEGREEPI(Qtmp);
	while ((Qtmp != NULL) && ((j = DEGREEPI(Qtmp)) != 0))
	{
		Rtmp =	MODPI(Ptmp, Qtmp);
		DELETEPI(Ptmp);
		Ptmp = Qtmp;
		Qtmp = PRIMITIVEPI(Rtmp);
		DELETEPI(Rtmp);
	}
	if (j == 0)
	{
		TempI = Qtmp->COEF;
		Qtmp->COEF = ONEI();
		FREEMPI(TempI);
	}
	else/* Qtmp = NULL*/
		Qtmp = COPYPI(Ptmp);
	Sptr = SCALARPI(D, Qtmp);
	FREEMPI(D);
	DELETEPI(Ptmp);
	DELETEPI(Qtmp);
	return (Sptr);
}

POLYI DERIVPI(POLYI Pptr)
/*
 * returns the derivative of a polynomial.
 */
{
	POLYI CURSOR, Sptr;
	MPI *TempI;

	if (DEGREEPI(Pptr) == 0)
		return (NULL);
	else
	{
		Sptr = COPYPI(Pptr);
  		CURSOR = Sptr;
		while (CURSOR != NULL)
		{
			TempI = CURSOR->COEF;
			CURSOR->COEF = MULT_I(CURSOR->COEF, CURSOR->DEG);
			FREEMPI(TempI);
			CURSOR->DEG = CURSOR->DEG - 1;
			if (CURSOR->NEXT == NULL)
				break;
			else if ((CURSOR->NEXT)->DEG == 0)
			{
				DELETEPI(CURSOR->NEXT);
				CURSOR->NEXT = NULL;
				break;
			}
			else
				CURSOR = CURSOR->NEXT;
		}
		return (Sptr);
	}
}

POLYI *SQUAREFREEPI(POLYI Pptr, USI **D, USI *tptr)
/*
 * returns the "squarefree decomposition" of Pptr, a non-constant polynomial:
 * Pptr = G[0]^D[0]...G[*tpr - 1]^D[*tptr - 1],
 * where G[i] is squarefree and is the product of the irreducible factors of
 * multiplicity D[i]. See L. Childs, Introduction to Higher Algebra, p.159.
 */
{
	POLYI Ttmp, Ttmp1, Dtmp, Ftmp, Gtmp, Htmp, Ptmp, *G;
	unsigned int i = 1, j = 0, n;

	n = (unsigned int)(Pptr->DEG);
	Ptmp = PRIMITIVEPI(Pptr);
	Ftmp = DERIVPI(Ptmp);
	Dtmp = GCDPI(Ptmp, Ftmp);
	G = (POLYI *)mmalloc(n * sizeof(POLYI));
	*D = (unsigned int *)mmalloc(n * sizeof(unsigned int));
	if (DEGREEPI(Dtmp) == 0)
	{
		G[j] = Ptmp;
		(*D)[j] = i;
		j++;
		DELETEPI(Ftmp);
		DELETEPI(Dtmp);
		*tptr = j;
		return (G);
	}
	else
	{
		Gtmp = DIVPI(Ptmp, Dtmp);
		while (DEGREEPI(Dtmp))
		{
			DELETEPI(Ptmp);
			DELETEPI(Ftmp);
			Htmp = Gtmp;
			Ptmp = Dtmp;
			Ftmp = DERIVPI(Ptmp);
			Dtmp = GCDPI(Ptmp, Ftmp);
			Gtmp = DIVPI(Ptmp, Dtmp);
			Ttmp = DIVPI(Htmp, Gtmp);
			DELETEPI(Htmp);
			Ttmp1 = PRIMITIVEPI(Ttmp);
			if (Ttmp1->DEG > 0)
			{
				G[j] = COPYPI(Ttmp1);
				(*D)[j] = i;
				j++;
			}
			DELETEPI(Ttmp);
			DELETEPI(Ttmp1);
			i++;
		}
		DELETEPI(Ptmp);
		DELETEPI(Ftmp);
		G[j] = PRIMITIVEPI(Gtmp);
		(*D)[j] = i;
		j++;
		DELETEPI(Gtmp);
		DELETEPI(Dtmp);
		*tptr = j;
		return (G);
	}
}

POLYI POWERPI(POLYI Pptr, USI n)
/*
 * returns Pptr ^ n, where 0 <= n < R0 * R0.
 */
{
	POLYI B, Eptr, Tmp;

	Eptr = ONEPI();
	if (n == 0)
		return (Eptr);
	B = COPYPI(Pptr);
	while (1)
        {
		if (n & 1) 
		{
			Tmp = Eptr;
			Eptr = MULTPI(Tmp, B);
			DELETEPI(Tmp);
			if (n == 1)
			{
				DELETEPI(B);
				return (Eptr);
			}
		}
		Tmp = B;
		B = MULTPI(Tmp, Tmp);
		DELETEPI(Tmp);
		n = n >> 1;
	}
}

POLYI ONEPI()
/*
 * returns the constant polynomial 1.
 */
{
	POLYI Pptr;

	Pptr = CREATEPI();
	Pptr->DEG = 0;
	Pptr->COEF = ONEI();
	Pptr->NEXT = NULL;
	return (Pptr);
}
POLYI ZEROPI()
     /* returns the zero polynomial */
{
  POLYI Pptr;
  Pptr= CREATEPI();
  Pptr->DEG=0;
  Pptr->COEF=ZEROI();
  Pptr->NEXT=NULL;
  return (Pptr);
}
POLYI CONSTANTPI(MPI *Aptr)
/*
 * returns the constant polynomial whose constant term is *Aptr.
 */
{
	POLYI Pptr;	

	if (Aptr->S == 0)
		return (NULL);
	else
	{
		Pptr = CREATEPI();
		Pptr->DEG = 0;
		Pptr->COEF = COPYI(Aptr);
		Pptr->NEXT = NULL;
		return (Pptr);
	}
}

MPI *SUBRESULTANT(POLYI Pptr, POLYI Qptr)
/*
 * *Aptr is the resultant of Pptr and Qptr, found using the subresultant
 * algorithm, p. 130. E. Kaltofen, G. E. Collins etc, Algorithm 4.5.
 * similar to Knuth, Algorithm C, p. 410.
 */
{
	unsigned delta;
	MPI *G, *G1, *H1, *H2, *Aptr, *TempI;
	POLYI Utmp, Vtmp, Rtmp, CURSOR;

	G = ONEI();
	Aptr = ONEI();
	Utmp = COPYPI(Pptr);
	Vtmp = COPYPI(Qptr);
	while (Vtmp != NULL)
	{
		Rtmp = MODPI(Utmp, Vtmp);
		if (Rtmp == NULL && DEGREEPI(Vtmp) > 0)
		{
			FREEMPI(G);
			FREEMPI(Aptr);
			DELETEPI(Utmp);
			DELETEPI(Vtmp);
			return (ZEROI());
		}
		delta = DEGREEPI(Utmp) - DEGREEPI(Vtmp);
		DELETEPI(Utmp);
		Utmp = Vtmp;
		G1 = MINUSI(G);
		H1 = MINUSI(Aptr);
		TempI = H1;
		H1 = POWERI(H1, delta);
		FREEMPI(TempI);
		TempI = H1;
		H1 = MULTI(H1, G1);
		FREEMPI(TempI);
		CURSOR = Rtmp;
		while (CURSOR != NULL)
		{
			TempI = CURSOR->COEF;
			CURSOR->COEF = INTI(CURSOR->COEF, H1);
			FREEMPI(TempI);
			CURSOR = CURSOR->NEXT;
		}
		Vtmp = Rtmp;
		TempI = G;
		G = LEADPI(Utmp);
		FREEMPI(TempI);
		TempI = G1;
		G1 = POWERI(G, delta);
		FREEMPI(TempI);
		TempI = H1;
		H1 = MULTI(Aptr, G1);
		FREEMPI(TempI);
		H2 = POWERI(Aptr, delta);
		TempI = Aptr;
		Aptr = INTI(H1, H2);
		FREEMPI(TempI);
		FREEMPI(H1);
		FREEMPI(H2);
		FREEMPI(G1);
	}
	DELETEPI(Utmp);
	FREEMPI(G);
	return (Aptr);
}

MPI *DISCRIMINANTPI(POLYI Pptr)
/*
 * *Dptr = Discriminant of Pptr = (1/a[n])*(-1)^{n*(n-1)/2 * Res(Pptr, Pptr').
 * See O. Perron, Algebra, Vol 1, p.212.
 */
{
	POLYI P;
	MPI *A, *Dptr, *TempI;
	int n;

	P = DERIVPI(Pptr);
	Dptr = SUBRESULTANT(Pptr, P);
	A = LEADPI(Pptr);
	TempI = Dptr;
	Dptr = INTI(Dptr, A);
	FREEMPI(TempI);
	FREEMPI(A);
	DELETEPI(P);
	n = (DEGREEPI(Pptr)) % 4;
	if (n == 2 || n == 3)
		Dptr->S = -(Dptr->S);
	return (Dptr);
}

MPI *LENGTHSQPI(POLYI Pptr)
/*
 * *Cptr is the sum of the squares of the coefficients of Pptr.
 */
{
	POLYI CURSOR;
	MPI *C, *Cptr, *TempI;

	Cptr = ZEROI();
	if (Pptr == NULL)
		return (Cptr);
	CURSOR = Pptr;
	while (CURSOR != NULL)
	{
		C = MULTI(CURSOR->COEF, CURSOR->COEF);
		TempI =  Cptr;
		Cptr = ADDI(C, Cptr);
		FREEMPI(C);
		FREEMPI(TempI);
		CURSOR = CURSOR->NEXT;
	}
	return (Cptr);
}

MPI *VAL0PI(POLYI Pptr)
/*
 * *Mptr = Pptr(0)
 */
{
	POLYI CURSOR;

	if (Pptr == NULL)
		return ZEROI();
	CURSOR = Pptr;
	while (CURSOR->NEXT != NULL)
		CURSOR = CURSOR->NEXT;
	if (DEGREEPI(CURSOR) != 0)
		return ZEROI();
	else
		return COPYI(CURSOR->COEF);	
}

/* Returns respectively positive, or negative if the sign of the polynomial 
 * evaluated at the supplied rational number is positive or negative.
 * Returns zero if polynomial has a root at the supplied rational */
int SIGN_AT_R_PI(POLYI P, MPR *R)
{
  MPI *TMPI, *RESULT=ZEROI(), *X_TO_I, *TERM;
  MPI *A=COPYI(R->N), *D=COPYI(R->D);
  POLYI CURSOR=P;
  unsigned int n=DEGREEPI(P), i;

  if (P == NULL) 
    i=0;
  else while (CURSOR) {
    i = DEGREEPI(CURSOR);
    TERM=POWERI(A, i);
    
    TMPI=TERM;
    TERM=MULTI(CURSOR->COEF, TERM);
    FREEMPI(TMPI);
    
    X_TO_I=POWERI(D, n-i);
    
    TMPI=TERM;
    TERM=MULTI(TERM, X_TO_I);
    FREEMPI(TMPI);
    FREEMPI(X_TO_I);

    TMPI=RESULT;
    RESULT=ADDI(RESULT, TERM);
    FREEMPI(TMPI);
    FREEMPI(TERM);

    CURSOR=CURSOR->NEXT;
  }
  i = SIGNI(RESULT); /* shouldn't really use this as return val but it works */
  FREEMPI(RESULT);
  FREEMPI(A);
  FREEMPI(D);
  return i;
}

/* If P = p(x) then this function returns p(-x) */
POLYI P_OF_NEG_X(POLYI P)
{
  POLYI Q=COPYPI(P);
  POLYI CURSOR=Q;
  while (CURSOR) {
    if (DEGREEPI(CURSOR)%2 == 1) /* if odd degree */
      CURSOR->COEF->S = - CURSOR->COEF->S;
    CURSOR=CURSOR->NEXT;
  }
  return Q;
}

/* Returns -1*P(x) */
POLYI MINUSPI(POLYI P)
{
  POLYI Q=COPYPI(P);
  POLYI CURSOR=Q;
  while (CURSOR) {
    CURSOR->COEF->S = -CURSOR->COEF->S;
    CURSOR=CURSOR->NEXT;
  }
  return Q;
}

/* Modifies the supplied polynomial so that p(x) = p(x+A) */

void P_OF_X_PLUS(POLYI *P, MPI *A)
{
  MPI *ZERO = ZEROI();
  MPIA COEF = BUILDMPIA();
  POLYI CURSOR = *P;
  int i, j,  n = DEGREEPI(*P);

  for (i=n; i >=0; i--) {
	ADD_TO_MPIA(COEF, ZERO, i);
  }
  while (CURSOR) {
    ADD_TO_MPIA(COEF, CURSOR->COEF, DEGREEPI(CURSOR));
    CURSOR=CURSOR->NEXT;
  }

  for (i=0; i <= n-1; i++)
    for (j = n-1;j >= i; j--) {
      MPI *TMP, *MUL;
      MUL = MULTI(A, COEF->A[j+1]);
      TMP = COEF->A[j];
      COEF->A[j] = ADDI(MUL, COEF->A[j]);
      FREEMPI(TMP);
      FREEMPI(MUL);
    }
  
  DELETEPI(*P);
  *P = NULL;
  for (i=0; i <= n; i++)
    PINSERTPI(i, COEF->A[i], P, 0);
  PURGEPI(P);
  FREEMPIA(COEF);
  FREEMPI(ZERO);
}
    
POLYI PRIMITIVEPI_(POLYI Pptr)
/*
 * returns the primitive part of the polynomial Pptr.
 * retaining the sign of the leading coefficient.
 */
{
	POLYI Sptr, CURSOR;
	MPI *C, *TempI;

	if (Pptr == NULL)
		return (NULL);
	else
	{
		Sptr = COPYPI(Pptr);
		C = CONTENTPI(Sptr);
		CURSOR = Sptr;
		while (CURSOR != NULL)
		{
			TempI = CURSOR->COEF;
			CURSOR->COEF = INT(CURSOR->COEF, C);
			FREEMPI(TempI);
			CURSOR = CURSOR->NEXT;
		}
		FREEMPI(C);
		return (Sptr);
	}
}