i5r.c

calc大数库

		return COPYR(Bptr);
	else if ((Bptr->N)->S == 0)
		return COPYR(Aptr);
	else
	{
		D = GCD(Aptr->D, Bptr->D);	
		if (EQONEI(D))
		{
			FREEMPI(D);
			X = MULTI(Aptr->N, Bptr->D);
			Y = MULTI(Aptr->D, Bptr->N);
			Eptr = BUILDMPR();	
			Eptr->N = ADDI(X, Y);
			Eptr->D = MULTI(Aptr->D, Bptr->D);
			FREEMPI(X);
			FREEMPI(Y);
			return (Eptr);
		}	
		else
		{
			R = INT0(Aptr->D, D);
			S = INT0(Bptr->D, D);
			X = MULTI(Aptr->N, S);
			Y = MULTI(Bptr->N, R);
			T1 = ADDI(X, Y);
			T2 = MULTI(Aptr->D, S);
			FREEMPI(R);
			FREEMPI(S);
			FREEMPI(X);
			FREEMPI(Y);
			if (T1->S == 0)
			{
				FREEMPI(D);
				FREEMPI(T1);
				FREEMPI(T2);
				return ZEROR();
			}
			else
			{
				Eptr = BUILDMPR();	
				E = GCD(T1, D);				
				FREEMPI(D);
				if (EQONEI(E))
				{
					FREEMPI(E);
					Eptr->N = T1;
					Eptr->D = T2;
					return (Eptr);
				}
				else
				{
					Eptr->N = INT(T1, E);
					Eptr->D = INT(T2, E);
					FREEMPI(E);
					FREEMPI(T1);
					FREEMPI(T2);
					return (Eptr);
				}
			}
		}
	}
}

MPR *MULTR(MPR *Aptr, MPR *Bptr)
/*
 * *Eptr = *Aptr * *Bptr.
 * P. Henrici's method: see Computer Algebra Symbolic and Algebraic
 * Computation, Springer 1982,p. 202.
 */
{
	MPI *D1=NULL, *D2=NULL, *R1=NULL, *R2=NULL, *S1=NULL, *S2=NULL;
	MPR *Eptr;
	if ((Aptr->N)->S == 0 || (Bptr->N)->S == 0)
		return ZEROR();
	else
	{
	  D1 = GCD(Aptr->N, Bptr->D);
	  D2 = GCD(Bptr->N, Aptr->D);
	  if (EQONEI(D1))
	    {
	      R1 = COPYI(Aptr->N);
	      S2 = COPYI(Bptr->D);
	    }
	  else
	    {
	      R1 = INT(Aptr->N, D1);
	      S2 = INT0(Bptr->D, D1);
	    }
	  if (EQONEI(D2))
	    {

	      S1 = COPYI(Bptr->N);
	      R2 = COPYI(Aptr->D);
	    }
	  else
	    {
	      S1 = INT(Bptr->N, D2);

	      R2 = INT0(Aptr->D, D2);


	    }
	  Eptr = BUILDMPR();
	  Eptr->N = MULTI(R1, S1);
	  Eptr->D = MULTI(R2, S2);
	  FREEMPI(D1); 
	  FREEMPI(D2); 
	  FREEMPI(R1); 
	  FREEMPI(R2); 
	  FREEMPI(S1); 
	  FREEMPI(S2); 


	  return (Eptr);
	}
}

MPR *POWERR(MPR *Aptr, unsigned int n)
/*
 * *Eptr = (*Aptr) ^ n, where 0 <= n < R0 * R0.
 */
{
  MPR *B, *Temp, *Eptr;
  
  if (n == 0)
    return ONER();
  B = COPYR(Aptr);
  Eptr = ONER();
  while (1)
    {
      if (n & 1)
	{
	  
	  Temp = Eptr;
	  Eptr = MULTR(Eptr, B);
	  FREEMPR(Temp);
	  if (n == 1)
	    {
	      FREEMPR(B);
	      return (Eptr);
	    }
	}
      Temp = B;
      B = MULTR(B, B);
      FREEMPR(Temp);
      n = n >> 1;
    }
}

/* As above but can find powers of negative numbers */
MPR *POWERR_2(MPR *R, int n)
{
  MPR *RETVAL;
  if (n < 0) {
    MPR *TMPR;
    RETVAL= POWERR(R, -n);   
    TMPR=RETVAL;
    RETVAL=INVERSER(RETVAL);
    FREEMPR(TMPR);
  } else {
    RETVAL = POWERR(R, n);
  }
  return RETVAL;
}

MPR *SUBR(MPR *Aptr, MPR *Bptr)
/*
 * *Eptr = *Aptr - *Bptr.
 */
{
	MPI *X, *Y, *U, *V, *G;
	MPR *Eptr;

	X = MULTI(Aptr->N, Bptr->D);
	Y = MULTI(Aptr->D, Bptr->N);
	V = MULTI(Aptr->D, Bptr->D);
	U = SUBI(X, Y);
	G = GCD(U, V);
	Eptr = BUILDMPR();
	Eptr->N = INT(U, G);
	Eptr->D = INT(V, G);
	FREEMPI(X);
	FREEMPI(Y);
	FREEMPI(U);
	FREEMPI(V);
	FREEMPI(G);
	return (Eptr);
}

MPR *INVERSER(MPR *Aptr)
/*
 * *Eptr = 1 / *Aptr.
 */
{
	MPR *Eptr;
	int s;

	Eptr = BUILDMPR();
	Eptr->D = ABSI(Aptr->N); 
	Eptr->N = COPYI(Aptr->D);
	s = (Aptr->N)->S;
	(Eptr->N)->S = ((Eptr->N)->S) * s;
	return (Eptr);
}
 
MPR *RATIOR(MPR *Aptr, MPR *Bptr)
/*
 * *Eptr = *Aptr / *Bptr.
 */
{
	MPR *X, *Eptr;
 
	X = INVERSER(Bptr);
	Eptr = MULTR(Aptr, X);
	FREEMPR(X);
	return (Eptr);
}

MPR *FRAC_PARTR(MPR *Aptr)
/*
 * *Bptr = fractional part of *Aptr.
 */
{
	MPI *X;
	MPR *Bptr;

	if (EQONEI(Aptr->D))
		Bptr = ZEROR();
	else
	{
		X = MOD(Aptr->N, Aptr->D);
		Bptr = RATIOI(X, Aptr->D);
		FREEMPI(X);
	}
	return (Bptr);
}

MPR *FRAC_PARTI(MPI *Aptr, MPI *Bptr)
/*
 * *Cptr = fractional part of *Aptr/(*Bptr).
 */
{
	MPR *C, *Cptr;

	C = RATIOI(Aptr, Bptr);
	Cptr = FRAC_PARTR(C);
	FREEMPR(C);
	return (Cptr);
}

MPI *NEAREST_INTR(MPR *Aptr)
/*
 * *Bptr is the nearest integer to *Aptr.
 */
{
	MPI *Y, *Z, *Temp, *Bptr;
	MPR *X;

	Y = INT(Aptr->N, Aptr->D);
	X = FRAC_PARTR(Aptr);
	Z = MULT_I(X->N, 2L);
	if (RSV(Z, X->D) <= 0)
		Bptr = Y;
	else
	{
		Temp = ONEI();
		Bptr = ADDI(Y, Temp);
		FREEMPI(Y);
		FREEMPI(Temp);
	}
	FREEMPR(X);
	FREEMPI(Z);
	return (Bptr);
} 

MPI *ABS_NEAREST_INTR(MPR *Aptr)
/*
 * *Bptr is the nearest integer to *Aptr, taking the integer closer to 0
 * if half an odd integer.
 */
{
	MPI *Y, *Z, *Bptr, *O;
	MPR *X;
	int t;

	Y = INT(Aptr->N, Aptr->D);
	X = FRAC_PARTR(Aptr);
	Z = MULT_I(X->N, 2L);
	t = RSV(Z, X->D);
	FREEMPR(X);
	FREEMPI(Z);
	if (t < 0)
		Bptr = Y;
	else if (t == 0)
	{
		if (Y->S >= 0)	
			Bptr = Y;
		else
		{
			O = ONEI();
			Bptr = ADDI(Y, O);
			FREEMPI(O);
			FREEMPI(Y);
		}
	}
	else
	{
		O = ONEI();
		Bptr = ADDI(Y, O);
		FREEMPI(Y);
		FREEMPI(O);
	}
	return (Bptr);
} 

MPI *INPUTSID(char **ptr, unsigned int *uptr)
/*
 * For use in converting the decimals after the decimal point to an MPI.
 * *uptr is the number of decimals met after the decimal point.
 */
{
	unsigned int long n;
	char *t;
	MPI *Mptr, *Temp;

	Mptr = ZEROI();
	t = *ptr;
	while ((**ptr >= '0' && **ptr <= '9')) 
	{
		Temp = Mptr;
		Mptr = MULT_I(Temp, 10L);
		FREEMPI(Temp);
		n = (unsigned long)(**ptr - '0');
		Temp = Mptr;
		Mptr = ADD0_I(Temp, n);
		FREEMPI(Temp);
		(*ptr)++;
	}
	*uptr = *ptr - t;
	return (Mptr);
}

MPR *INPUTSRD(char **ptr)
/*
 * Converts a floating point decimal to an MPI. Used in inputting a file of
 * matrices.
 */
{
	unsigned int n, length;
	MPI *G, *H, *K, *L, *M, *N;
	MPR *Aptr;

	G = INPUTSI(ptr, &n);
	if (n == 0) 
	{
		FREEMPI(G);/* default - junk at start of number */
		return ZEROR();
	}
	else if (**ptr == '.')
	{
		(*ptr)++;
		H = INPUTSID(ptr, &length);
		K = POWER_I(10L, length);
		M = ABSI(G);
		L = MULTI(M, K);
		N = ADD0I(L, H);
		if (G->S == -1)
			N->S = -(N->S);
		Aptr = RATIOI(N, K);
		FREEMPI(H);
		FREEMPI(K);
		FREEMPI(M);
		FREEMPI(L);
		FREEMPI(N);
		FREEMPI(G);
	}
	else 
	{
		Aptr = BUILDMPR();
		Aptr->N = G;
		Aptr->D = ONEI();
	}
	return (Aptr);
}

void FPRINTMATR(FILE *outfile, USI i1, USI i2, USI j1, USI j2, MPMATR *Mptr)
/*
 * Modified 8/1/93 using Peter Adams' improvement from, August 1992.
 */
{
	char **str, *buff;
	unsigned int tmp, i, j, nrow, ncol, len, nstr = 0;
	int *colwidth;
	unsigned int ct;
	
	nrow = i2 - i1 + 1;
	ncol = j2 - j1 + 1;
	str = (char **)mmalloc(nrow * ncol * sizeof(char *));
	colwidth=(int *)mmalloc(ncol*sizeof(int));

	for(i=0; i<ncol; i++)
	  colwidth[i]=0;
	for (i = i1; i <= i2; i++)
	{
		for (j = j1; j <= j2; j++)
		{
			tmp = 1 + LENGTHR(elt(Mptr, i, j));
			buff = (char *)mmalloc(tmp * sizeof(char));
			SPRINTR(buff, elt(Mptr, i, j));
		    	if ((len = strlen(buff)) > colwidth[j-j1])
				colwidth[j-j1] = len;
			str[nstr++] = strcpy((char *)mmalloc(len + 1), buff);
			ffree(buff, tmp * sizeof(char));
		}
	}
	if (outfile != stdout)
		fprintf(outfile, "%u %u\n", nrow, ncol);
	ct=0;
	for (i = i1; i <= i2; i++)
	{
		for (j = j1; j <= j2; j++)
		{
	    		fprintf(outfile, "%*s ", colwidth[j-j1], str[ct]);
	    		ffree(str[ct], (unsigned int) (strlen(str[ct]) + 1) * sizeof(char));
	    		if ((ct % ncol) == (ncol - 1))
				fprintf(outfile, "\n");
	    		ct++;
		}
	}
	ffree((char *) str, nrow * ncol * sizeof(char *));
	ffree((char *) colwidth, ncol*sizeof(int));
	return;
}

void PRINTMATR(USI i1, USI i2, USI j1, USI j2, MPMATR *Mptr)
/*
 * prints *Mptr from rows i1 to i2 and cols j1 to j2.
 */
{
	unsigned int i, j;

	i = i2 - i1 + 1;
	j = j2 - j1 + 1;
	if (MAX((int)i, (int)j) > 20)
	{
		printf("(The number of rows or columns to be printed exceeds 20;\n");
		printf("there is no point in printing this matrix on the screen.)\n");
		return;
	}
	FPRINTMATR(stdout, i1, i2, j1, j2, Mptr);
	return;
}


MPR *INTR(MPR *Aptr)
/*
 *  Returns the integer part of *Aptr.
 */
{
	MPI *X;
	MPR *Bptr, *Temp;

	if (EQONEI(Aptr->D))
		Bptr = COPYR(Aptr);
	else
	{
		X = MOD(Aptr->N, Aptr->D);
		Temp = RATIOI(X, Aptr->D);
		FREEMPI(X);
		Bptr = SUBR(Aptr, Temp);
		FREEMPR(Temp);
	}
	return (Bptr);
}