roots.c

calc大数库

#include "roots.h"
#include "integer.h"
#include <stdio.h>
#include <stdlib.h>
#include "fun.h"
#include "stack.h"

/* This code was written by Sean Seefried when he was a vacation student
 * with me at UQ Maths. Dept.
 */
/* Simple linked list implementation of polynomial sequence */
typedef struct _PolySeq {
  POLYI poly;
  struct _PolySeq *next;
} *PolySeq;

/* 
 * Adds a POLYI to the end of a PolySeq.  If supplied PolySeq is NULL then a 
 *  new PolySeq is created.
 */
PolySeq addPolySeqTerm(PolySeq ps, POLYI P)
{
  PolySeq CURSOR=ps;
  if (ps == NULL) {
    if ((ps=mmalloc(sizeof (struct _PolySeq))) == NULL)
      return NULL;
    else {
      ps->poly=COPYPI(P);
      ps->next=NULL;
    }
  } else {
    while (CURSOR->next) 
      CURSOR=CURSOR->next; 
    /* CURSOR now points to last term */
    CURSOR->next=mmalloc(sizeof (struct _PolySeq));
    CURSOR->next->poly=COPYPI(P);
    CURSOR->next->next=NULL;
  }
  return ps;
}

/* 
 * Frees a polynomial sequence previously created with createPolySeq 
 */
void freePolySeq(PolySeq ps)
{
  PolySeq CURSOR=ps;
  PolySeq prev;
  while(CURSOR) {
    prev=CURSOR;
    DELETEPI(CURSOR->poly);
    CURSOR=CURSOR->next;
    ffree(prev, sizeof (struct _PolySeq)); 
  }
}
void printPolySeq(PolySeq ps) 
{
  PolySeq CURSOR=ps;
  printf("{");
  if (ps) {
    PRINTPI(CURSOR->poly);
    while (CURSOR->next) {
      printf(", ");
      CURSOR=CURSOR->next;
      PRINTPI(CURSOR->poly);
    }
  }
  printf("}\n");
}

void printPolySeqEval(PolySeq ps, MPR *R) 
{
  PolySeq CURSOR=ps;
  int i;
  printf("{");
  if (ps) {
    i = SIGN_AT_R_PI(CURSOR->poly, R);
    printf("%d", i);
    while (CURSOR->next) {
      printf(", ");
      CURSOR=CURSOR->next;
      i = SIGN_AT_R_PI(CURSOR->poly, R);
      printf("%d", i);
    }
  }
  printf("}\n");
}
/*
 * Returns the number of sign variations in the Polynomial Sequence at the 
 * evualated at the supplied value.
 * A sign variation exists between two non-zero numbers c[p] and c[q] (p<q)
 * of a finite or infinite sequence of numbers, c[1], c[2], c[3],... if the
 * following hold:
 * i. for q = p + 1, c[p] and c[q] have opposite signs.
 * ii. for q >= p +2, the numbers c[p+1], ..., c[q-1] are all zero and c[p]
 * and c[q] have opposite signs.
 */


int signVar(PolySeq PS, MPR *R) 
{
  int prev=SIGN_AT_R_PI(PS->poly, R);
  /* previous polynomial in sequence's value at I */
  int vars=0; /* number of variations */
  MPR *Z=ZEROR();
  int sign;
  PolySeq CURSOR=PS;
  while (CURSOR) {
    /* skip polynomials that evaluate to zero */
    sign = SIGN_AT_R_PI(CURSOR->poly, R);
    while (sign == 0 && CURSOR->next) {
      CURSOR=CURSOR->next;
      sign = SIGN_AT_R_PI(CURSOR->poly, R);
    }
    if (prev != 0 && (prev*sign < 0))
      vars++;
    prev=sign;
    CURSOR=CURSOR->next;
  }
  FREEMPR(Z);
  return vars;
}

/* 
 * Returns a polynomial sequence that is the sturm sequence of the supplied 
 * polynomial.  See Akritas, Elements of Computer Algebra, p341 
 */ 
PolySeq sturmSeq(POLYI Pptr)
{
  POLYI R1=COPYPI(Pptr);
  POLYI R2=DERIVPI(R1);
  POLYI TR, R; /* temp remainder, remainder, modified remainder */
  MPI *CONT;
  PolySeq ps=NULL;
  
  ps = addPolySeqTerm(ps, R1);
  /* It is possible that Pptr was a constant polynomial. Hence its derivative
   * is zero or NULL.  The Sturm sequence should consists of one term - Pptr.
   */
  if (R2 != NULL) 
    ps = addPolySeqTerm(ps, R2);
  else {
    DELETEPI(R1);
    return ps;
  }
  while ((DEGREEPI(TR=MODPI(R1,R2)) != 0) && TR != NULL) {
    CONT=CONTENTPI2(TR);
    if (SIGNI(CONT) > 0) {
      POLYI TMPI;
      R = PRIMITIVEPI(TR);
      DELETEPI(TR); /* free TR. not needed */

      TMPI = R;
      R = MINUSPI(R);
      DELETEPI(TMPI);
    } else {
      R = PRIMITIVEPI(TR);
      DELETEPI(TR);
    }
    ps = addPolySeqTerm(ps, R);
    TR=R1;
    R1=R2;
    R2=R;
    DELETEPI(TR); /* deletes R1 */
    FREEMPI(CONT);
  }

  if (TR != NULL) { /* can guarantee that this is a constant polynomial */
    MPI *TMP;
    POLYI ONE_P=ONEPI();
    TMP = LEADPI(TR);
    if (SIGNI(TMP) > 0) {
      POLYI NEGONE_P;
      MPI *CONST=MINUS_ONEI();
      NEGONE_P=CONSTANTPI(CONST);
      addPolySeqTerm(ps, NEGONE_P);
      DELETEPI(NEGONE_P);
      FREEMPI(CONST);
    } else {
      addPolySeqTerm(ps, ONE_P);
    }
    DELETEPI(ONE_P);
    FREEMPI(TMP);
  }
  DELETEPI(TR);
  DELETEPI(R1);
  DELETEPI(R2);
  return ps;
}

/*
 * Calculates an upperbound on the positve roots of the supplied polynomial
 * using Cauchy's Rule. See Akritas, Elements of Computer Algebra.
 */
MPR *CAUCHY(POLYI P)
{
  POLYI Q=COPYPI(P);
  POLYI CURSOR;
  MPR *RETVAL=NULL;
  MPI *LEAD;
  MPR *TWO=TWOR();
  int t, deg=DEGREEPI(P);
  /* t a flag, deg is degree of polynomial */
  int  k, k_, k__, lambda=0, i, j, p, q, r;
  MPI *CI_=NULL;
  MPR *CN_=NULL;
  
  /* if the leading coefficient is less than zero then multiply polynomial by
   * -1 */

  LEAD=LEADPI(Q);
  if (SIGNI(LEAD) < 0) {
    MPI *NEGONE=MINUS_ONEI();
    POLYI TMPI;
    TMPI=Q;
    Q=SCALARPI(NEGONE, Q);
    DELETEPI(TMPI); 
    FREEMPI(NEGONE);
  }
  CURSOR=Q;
  while (CURSOR){ 
    if (SIGNI(CURSOR->COEF) < 0) /* count negative terms in poly */
      lambda++;
    CURSOR=CURSOR->NEXT;
  }

  k__ = 0 ; 
  j = BINARYB(LEAD)-1; /* j = log base 2 of LEAD */
  t = 0;
  if (!(lambda==0 || deg==0)) {
    CURSOR=Q;
    while (CURSOR) {
      if (SIGNI(CURSOR->COEF) < 0) {
	k=Q->DEG - CURSOR->DEG;
	CI_ = MULT_I(CURSOR->COEF, -lambda);
	i = BINARYB(CI_)-1; /* i = log base 2 of CI_ */
	p=i-j-1;
	q=p/k;
	r=p-k*q;
	if (r < 0) {
	  r+=k;
	  q--;
	}
	k_ = q + 1;

	if (r == (k-1)) {
	  MPR *TMPR;
	  MPR *POW;
	  POW=POWERR_2(TWO, k*k_);
	  TMPR=MPI_TO_MPR(Q->COEF);
	  CN_ = MULTR(TMPR, POW);
	  FREEMPR(POW);
	  FREEMPR(TMPR);
	  TMPR=MPI_TO_MPR(CI_);
	  if (COMPARER(TMPR, CN_) > 0) {
	    k_ = k_ + 1;
	    } 
	  FREEMPR(CN_);
	  FREEMPR(TMPR); 
	}
	if (t==0 || k_ > k__) {
	  k__ = k_;
	  t = 1;
	}
	FREEMPI(CI_); 
      }
      CURSOR=CURSOR->NEXT;
    }
  }
  RETVAL = POWERR_2(TWO, k__);
  
  FREEMPR(TWO);
  FREEMPI(LEAD);
  DELETEPI(Q);
  return RETVAL;
  
}


/* Returns a newly created Interval if possible.  
 * Returns NULL on error.
 */
Interval createInterval(MPR* left, MPR *right) 
{
  Interval intvl;
  if (!(intvl = malloc (sizeof (struct _Interval))))
    return NULL;
  intvl->left = COPYR(left);
  intvl->right = COPYR(right);
  return intvl;
}
/*-------------------------------------------------------------------*/
/* Frees an interval previously created with createInterval 
 * Undefined behaviour occurs if 'intvl' does not point to a structure of 
 * this type. */
void freeInterval(Interval intvl) 
{
  FREEMPR(intvl->left);
  FREEMPR(intvl->right);
  free(intvl);
}
/*-------------------------------------------------------------------*/

/* Returns a stack of intervals that contain only one root each respectively 
 * of the supplied polynomial.  If there are no roots, it returns an empty
 * stack */
Stack STURM(POLYI P)
{
  
  Stack Bounds=stackNew();
  Stack Results=stackNew();
  PolySeq sseq; /* sturm sequence */
  MPI *ZI, *V;
  MPR *BR, *ZERO=ZEROR(), *TWO = TWOR();
  unsigned debugCounter=0;
  unsigned numRoots, pn_flag=0;
  
  /* must find square free version of polynomial or else 
   * we cannot use the theorm of Sturm.  This is done simply by dividing 
   * the supplied polynomial by the gcd of it and it's derivative */
  POLYI D=DERIVPI(P);
  POLYI G=GCDPI(P,D);
  POLYI Pw=DIVPI(P, G);
  DELETEPI(D);
  DELETEPI(G);

  /* if there is a zero at zero then divide polynomial by X */
  ZI = ZEROI();
  V = VALPI(Pw, ZI);
  if (COMPAREI(V, ZI) == 0) {
    MPI *ONE=ONEI();
    POLYI X=NULL, TMPI;
    PINSERTPI(1, ONE, &X, 0);
    TMPI = Pw;
    Pw = DIVPI(Pw, X);
    DELETEPI(TMPI);
    DELETEPI(X);
    FREEMPI(ONE); 
  }
    FREEMPI(V);
    FREEMPI(ZI); 
  
    while (pn_flag < 2) {
      POLYI TMPPOLY;
      sseq = sturmSeq(Pw);

      BR = CAUCHY(Pw); 
      stackPush(Bounds, createInterval(ZERO, BR));
      FREEMPR(BR);
      /* main loop */
      
      while (!stackEmpty(Bounds) && debugCounter < 20) {
	Interval intvl;
	intvl = stackPop(Bounds);
	numRoots = signVar(sseq, intvl->left) - signVar(sseq, intvl->right);
	debugCounter++;
	/* DEBUG INFO */

	/*	printf("The number of roots in the interval (");
	 PRINTR(intvl->left);
	 printf(", ");
	 PRINTR(intvl->right);
	 printf(") is %d\n", numRoots);     
	 printPolySeqEval(sseq, intvl->left);
	 printPolySeqEval(sseq, intvl->right);         */
	 /* DEBUG END */

	 if (numRoots == 1) {
	   
	   MPR *L, *R;
	   L = intvl->left;
	   R = intvl->right;
	   if (pn_flag == 1) {
	     MPR *TMPR;
	     TMPR = MINUSR(L); /* swap left and right of interval because  
				* they are now negative values */
	     L = MINUSR(R);
	     R = TMPR; 
	}
	   stackPush(Results, createInterval(L, R));
	   if (pn_flag == 1) {
	     FREEMPR(L);
	     FREEMPR(R);
	   } 
	 }
	 else if (numRoots > 1) { /* numroots > 1 */
	   MPR *T, *CENTRE;
	   T = ADDR(intvl->left,intvl->right);
	   CENTRE = RATIOR(T,TWO);
	   FREEMPR(T);  /* New intervals are (L, (L+R)/2) and ((L+R)/2, R) */
	   stackPush(Bounds, createInterval(intvl->left, CENTRE));
	   stackPush(Bounds, createInterval(CENTRE, intvl->right));
	   
	/* debug code */
	   /*	   printf("This interval has been subdivided into: (");
		   PRINTR(intvl->left);
	     printf(", ");
	     PRINTR(CENTRE);
	     printf(") , (");
	     PRINTR(CENTRE);
	     printf(", ");
	     PRINTR(intvl->right);
	     printf(")\n");    */
	   /* debug code end */
	   FREEMPR(CENTRE); 
	 }
	 freeInterval(intvl); 
      }
      TMPPOLY=Pw;
      Pw=P_OF_NEG_X(Pw);
      DELETEPI(TMPPOLY);   
      freePolySeq(sseq);
      pn_flag++;
    }
    stackFree(&Bounds);

    FREEMPR(ZERO);
    FREEMPR(TWO);
    DELETEPI(Pw);
    return Results;
}
/*-------------------------------------------------------------------*/