roots.c

calc大数库

/* A wrapper function that takes a single polynomial as its argument.
 * (Handled by parser).  Returns open rational intervals that are
 * guaranteed to enclose the real roots of the supplied polynomial.
 */
void STURM_W(Stack s)
{
  POLYI P=stackPop(s);
  Stack Results=STURM(P);
  FILE *outfile;

  if (!(outfile = fopen("sturm.out", "w"))) {
    printf("Could not open file sturm.out for writing\n");
    exit(1);
  }

  fprintf(outfile, "The intervals containing only one root for the polynomial:\n");
  FPRINTPI(outfile, P);
  fprintf(outfile, " are:\n");

  DELETEPI(P);
  while (!stackEmpty(Results)) {
    Interval intvl;
    intvl = stackPop(Results);
    
    printf("(");
    PRINTR(intvl->left);
    printf(", ");
    PRINTR(intvl->right);
    printf(")\n");

    fprintf(outfile, "(");
    FPRINTR(outfile, intvl->left);
    fprintf(outfile, ", ");
    FPRINTR(outfile, intvl->right);
    fprintf(outfile, ")\n");

    freeInterval(intvl);
  } 
  fclose(outfile);
  stackFree(&Results);
}


/* Using a bisection method this finds the integer part root of the
 * supplied polynomial in the supplied _OPEN_ interval.  This function is
 * only guaranteed to return the correct value if the interval
 * supplied contains exactly one root, and this root does not fall at the
 * end points. */
MPI *rootInIntervalI(POLYI P, MPI *LEFT, MPI *RIGHT)
{
  MPI *TMP, *LOW=COPYI(LEFT), *HIGH=COPYI(RIGHT), *MID;
  MPI *ONE = ONEI();
  MPI *TWO = TWOI();
  TMP=ADDI(LOW, ONE);

  while (COMPAREI(TMP, HIGH) != 0) {
    MPI *HIGHVAL, *MIDVAL;
    FREEMPI(TMP);
    MID=ADDI(LOW,HIGH);
    
    TMP = MID;
    MID=INTI(MID, TWO);
    FREEMPI(TMP);
    MIDVAL = VALPI(P, MID);
    HIGHVAL = VALPI(P, HIGH);

    if (SIGNI(MIDVAL) != SIGNI(HIGHVAL)) {
      TMP = LOW;
      LOW = MID;
      FREEMPI(TMP);
    } else {
      TMP = HIGH;
      HIGH = MID;
      FREEMPI(TMP);
    }

    FREEMPI(MIDVAL);
    FREEMPI(HIGHVAL);
    TMP = ADDI(LOW, ONE);
  }
  
  FREEMPI(ONE);
  FREEMPI(TWO);
  FREEMPI(TMP);
  FREEMPI(HIGH);
  return LOW;
}

/* Using a bisection method this finds the integer part root of the
 * supplied polynomial in the supplied _OPEN_ interval.  This function is
 * only guaranteed to return the correct value if the interval
 * supplied contains exactly one root, and this root does not fall at the
 * end points. 
 * If RIGHT == NULL this signifies that we are searching for a root in
 * the open interval (LEFT, INFINITY).
*/
MPI *rootInIntervalR(POLYI P, MPR *LEFT, MPR *RIGHT)
{
  MPR *TMP, *LOW=COPYR(LEFT), *HIGH, *MID;
  MPR *TWO = TWOR();
  MPI *INT_LOW, *INT_HIGH; /* stores the integer part of LOW and HIGH */

  if (RIGHT != NULL)
    HIGH = COPYR(RIGHT);
  else 
    HIGH = CAUCHY(P);
 
  INT_LOW = INTI(LOW->N, LOW->D);
  INT_HIGH= INTI(HIGH->N, HIGH->D);

  while (COMPAREI(INT_LOW, INT_HIGH) != 0) {
    FREEMPI(INT_LOW);
    FREEMPI(INT_HIGH);
    MID=ADDR(LOW,HIGH);
    TMP = MID;
    MID=RATIOR(MID, TWO);
    FREEMPR(TMP);
    if (SIGN_AT_R_PI(P, MID) != SIGN_AT_R_PI(P, HIGH)) {
      TMP = LOW;
      LOW = MID;
      FREEMPR(TMP);
    } else {
      TMP = HIGH;
      HIGH = MID;
      FREEMPR(TMP);
    }
    INT_LOW = INTI(LOW->N, LOW->D);
    INT_HIGH= INTI(HIGH->N, HIGH->D);
  }
  FREEMPR(TWO);
  FREEMPI(INT_HIGH);
  FREEMPR(HIGH);
  FREEMPR(LOW);
  return INT_LOW;
}


/* Finds the continued fraction expansion of all real roots of the supplied
 * polynomial using Lagrange's method and methods first presented in a paper
 * by Cantor, Galyean and Zimmer called A Continued Fraction Algorithm for
 * Real Algebraic Number 
 */
void ROOTEXPANSION(Stack s)
{
  POLYI P = stackPop(s);
  MPI *M  = stackPop(s);
  POLYI Q, CURSOR;
  /* Making P square free */
  POLYI D=DERIVPI(P);
  POLYI G=GCDPI(P,D);
  POLYI TMPP;

  Stack intvlStack;
  Interval intvl;
  FILE *outfile;
  MPR *ONE = ONER();
  MPI *NUMQUO, *TMPI;
  MPI *ZERO =ZEROI();
  MPR *TMP;
  MPIA LANGFRAC; /* the lagrange partial quotients */
  int i=0, j, d;
  /* main loop */

  TMPP = P;
  P = DIVPI(P,G);
  DELETEPI(TMPP);
  DELETEPI(G);
  DELETEPI(D);

  intvlStack = STURM(P);


  if (!(outfile = fopen("rootexp.out", "w"))) {
    printf("Error: Could not open file rootexp.out for writing\n");
    exit(1);
  }

  while (!stackEmpty(intvlStack)) {
    MPI *A=NULL; /* next partial quotient */
    MPR *AR=NULL;
    MPR *RN, *RP=NULL; /* next and previous left bound */
    MPR *SN, *SP=NULL; /* next and previous right bound */
    MPIA CONTFRAC=BUILDMPIA(); /* the zimmer partial quotients */

    MPIA COEF; /* coefficients for the transformation */
    /* SN will equal NULL to represent INFINITY */
    intvl = stackPop(intvlStack);

    Q = COPYPI(P);
    RN = COPYR(intvl->left);
    SN = COPYR(intvl->right);
    i = 0;
    while(COMPARER(RN, ONE) != 0 || SN != NULL) {
      A = rootInIntervalR(Q, RN, SN);
      AR = MPI_TO_MPR(A);

      ADD_TO_MPIA(CONTFRAC, A, i++);

      FREEMPR(RP); /* ok even if null */
      FREEMPR(SP);

      RP = RN;
      SP = SN;

      TMP = ADDR(AR, ONE);
      if (SP != NULL && COMPARER(SP, TMP) < 0) {
	MPR *TMP2;
	RN = SUBR(SP, AR);
	TMP2 = RN;
	RN = INVERSER(RN);
	FREEMPR(TMP2);
      } else 
	RN = ONER();
      FREEMPR(TMP);

      if (COMPARER(RP, AR) > 0) {
	MPR *TMP2;
	SN = SUBR(RP, AR);
	TMP2 = SN;
	SN = INVERSER(SN);
	FREEMPR(TMP2);
      } else
	SN = NULL;


      P_OF_X_PLUS(&Q, A);

      FREEMPI(A);
      FREEMPR(AR);

      /* This performs transformation p(x) = x^d*p(1/x +a) */
      COEF = BUILDMPIA();
      d = DEGREEPI(Q);
      for (j=0; j<= d; j++)
	ADD_TO_MPIA(COEF, ZERO, j);
      
      CURSOR = Q;
      while(CURSOR) {
	ADD_TO_MPIA(COEF, CURSOR->COEF, DEGREEPI(CURSOR));
	CURSOR=CURSOR->NEXT;
      } 
      
      DELETEPI(Q);
      Q = NULL;
      for (j=0; j <= d; j++) 
	PINSERTPI(d-j, COEF->A[j], &Q, 0);
      PURGEPI(&Q);
      FREEMPIA(COEF);

      /* end of transformation */

    } /* the zimmer method is finished. We can use Lagrange's method now */
    FREEMPR(RN);
    FREEMPR(RP);
    FREEMPR(SN);
    FREEMPR(SP);

    if (SIGNI(Q->COEF) == -1) {
      POLYI TMPPOLY;
      TMPPOLY = Q;
      Q = MINUSPI(Q);
      DELETEPI(TMPPOLY);
    }
    
    /* this bit runs lagrange */
    NUMQUO = CHANGEI(i);
    TMPI = NUMQUO;
    NUMQUO = SUBI(M, NUMQUO);
    FREEMPI(TMPI);
    LAGRANGE(Q, &LANGFRAC, NUMQUO);
    FREEMPI(NUMQUO);

    /* now to output the continued fractions */
    fprintf(outfile, "Continued fraction expansion of roots for the polynomial: \n");
    FPRINTPI(outfile, P);
    fprintf(outfile, "\n");
    printf("The root in the interval (");
    PRINTR(intvl->left); 
    printf(", ");
    PRINTR(intvl->right);
    printf(") has the continued fraction expansion: \n");
    fprintf(outfile, "The root in the interval (");
    FPRINTR(outfile, intvl->left); 
    fprintf(outfile, ", ");
    FPRINTR(outfile, intvl->right);
    fprintf(outfile, ") has the continued fraction expansion: \n");
    freeInterval(intvl);
    for (j=0; j < CONTFRAC->size; j++) {
      printf("[%d]: ", j); 
      PRINTI(CONTFRAC->A[j]); 
      printf("\n");

      fprintf(outfile, "[%d]: ", j); 
      FPRINTI(outfile, CONTFRAC->A[j]); 
      fprintf(outfile, "\n");
    }
    for (j=0; j< LANGFRAC->size; j++) {
      printf("[%d]: ", j+CONTFRAC->size); 
      PRINTI(LANGFRAC->A[j]);
      printf("\n");

      fprintf(outfile, "[%d]: ", j); 
      FPRINTI(outfile, LANGFRAC->A[j]);
      fprintf(outfile, "\n");

    }
    printf("--------------------------------------------------\n");
    fprintf(outfile, "--------------------------------------------------\n");
    
    
    FREEMPIA(CONTFRAC);
    FREEMPIA(LANGFRAC);
    DELETEPI(Q);
  }  /* end of main loop */
  fclose(outfile);
  stackFree(&intvlStack);
  FREEMPR(ONE);
  FREEMPI(ZERO);
  DELETEPI(P); 
  FREEMPI(M);
}

USI SIGN_COUNT(MPIA A){
/* This counts the number of sign changes in the array A */
	USI i,j,n, changes;
	MPIA S;

        j=0;
	n=A->size;
	S=BUILDMPIA();
        for(i=0;i<n;i++){
                if((A->A[i]->S)==0){
                        continue;
                }else{
                        ADD_TO_MPIA(S, A->A[i], j);
                        j=j+1;
                }
        }
        changes=0;
        for(i=0;i<j;i++){
                if((S->A[i]->S)*(S->A[i+1]->S)<0){
                        changes=changes+1;
                }
        }
	FREEMPIA(S);
        return(changes);

}

MPI *STURM_SEQUENCE(POLYI P, MPI *B, MPI *E){
	POLYI TEMPPOLYI, PA, PB;
	MPIA VALUES;
	int n;
	MPI *TEMP, *MINUSONE;
	USI count, i, changes;

	MINUSONE = MINUS_ONEI();
	VALUES = BUILDMPIA();
	PA = PRIMITIVEPI_(P);
	if(E->S){
		printf("STURM[0]=");
		PRINTPI(PA);
		printf("\n");
	}
	TEMP = VALPI(PA, B);
	ADD_TO_MPIA(VALUES, TEMP, 0); 
	FREEMPI(TEMP);

	TEMPPOLYI = DERIVPI(PA);
	PB = PRIMITIVEPI_(TEMPPOLYI);
	DELETEPI(TEMPPOLYI);

	if(E->S){
		printf("STURM[1]=");
		PRINTPI(PB);
		printf("\n");
	}

	TEMP = VALPI(PB, B);
	ADD_TO_MPIA(VALUES, TEMP, 1); 
	FREEMPI(TEMP);

	n = DEGREEPI(PB);
	count = 2;
	while(n > 0){
		TEMPPOLYI = PB;
		PB = MODPI(PA, PB);
		n = DEGREEPI(PB);
		DELETEPI(PA);
		PA = TEMPPOLYI;

		TEMPPOLYI = PB;
		PB = PRIMITIVEPI_(PB);
		DELETEPI(TEMPPOLYI);
		TEMPPOLYI = PB;
		PB = SCALARPI(MINUSONE, PB);
		DELETEPI(TEMPPOLYI);
		if(E->S){
			printf("STURM[%u]=",count);
			PRINTPI(PB);
			printf("\n");
		}
		TEMP = VALPI(PB, B);
		ADD_TO_MPIA(VALUES, TEMP, count); 
		FREEMPI(TEMP);
		count++;
	}
	if(E->S){
		for(i = 0;i < count;i++){
			PRINTI(VALUES->A[i]);
			printf(", ");
		}
		printf("\n");
	}
	changes = SIGN_COUNT(VALUES);
	FREEMPIA(VALUES);
	FREEMPI(MINUSONE);
	DELETEPI(PA);
	DELETEPI(PB);
	return(CHANGE(changes));
}