i5r.c

calc大数库

/* i5R.c */
/* programs for doing arithmetic with MPR's. */

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

#ifdef DEBUG
extern long int nettbytes;
#endif

MPR *BUILDMPR( )
/*
 * mallocs space for an MPR.
 */
{
	MPR *Mptr;

	Mptr = (MPR *)mmalloc((USL)sizeof(MPR));
	return (Mptr);
}

void FREEMPR(MPR *Mptr)
/*
 * deallocates the space alloted for the MPR Mptr.
 */
{
	if (Mptr == NULL)
		return;
	FREEMPI(Mptr->N);
	FREEMPI(Mptr->D);
	ffree ((char *)Mptr, sizeof(MPR));
	return;
}

MPR *COPYR(MPR *Aptr)
/*
 * returns a pointer to a copy of *Aptr.
 */
{
	MPR *Bptr;
	Bptr = BUILDMPR();
	Bptr->N = COPYI(Aptr->N);
	Bptr->D = COPYI(Aptr->D);
	return (Bptr);
}

unsigned int EQUALR(MPR *Aptr, MPR *Bptr)
/*
 * returns 1 if *Aptr = *Bptr, 0 otherwise.
 */
{
	if (EQUALI(Aptr->N, Bptr->N) * EQUALI(Aptr->D, Bptr->D))
		return (1);
	else
		return (0);
}

MPR *MINUSR(MPR *Aptr)
/*
 * *Bptr = -(*Aptr).
 */
{
	MPR *Bptr;

	Bptr = BUILDMPR();
	Bptr->N = MINUSI(Aptr->N);
	Bptr->D = COPYI(Aptr->D);
	return (Bptr);
}

MPR *RATIOI(MPI *Aptr, MPI *Bptr)
/*
 * *Cptr = (*Aptr) / (*Bptr).
 */
{
	MPI *G;
	MPR *Cptr;

	Cptr = BUILDMPR();
	G = GCD(Aptr, Bptr);
	Cptr->N = INT(Aptr, G);
	Cptr->D = INT(Bptr, G);
	if ((Cptr->D)->S == -1)
	{
		(Cptr->D)->S = 1;
		(Cptr->N)->S = -((Cptr->N)->S);
	}
	FREEMPI(G);
	return (Cptr);
}

unsigned long LENGTHR(MPR *Mptr)
/*
 * returns the sum of the lengths of numerator and denominator + 1,
 * if Mptr is not an integer, otherwise returns the length of *Mptr.
 */
{
	if ((Mptr->D)->D || (Mptr->D)->V[0] > 1)
		return (LENGTHI(Mptr->N) + LENGTHI(Mptr->D) + 1);
	else
		return (LENGTHI(Mptr->N));
}

MPR *FINPUTR(FILE *f, unsigned int *uptr)
/*
 * Converts the ratio of two decimal inputs from stream into an MPR. 
 * *uptr =  0 if input fails, 1 if successful.
 */
{
	int t,  c;
	MPI *G, *H;
	MPR *Aptr;

	G = FINPUTI(f, uptr);
	if (*uptr == 0)
	{
		FREEMPI(G);
		return ZEROR();
	}
	t = c = fgetc(f);
	while (c == ' ')
		c = fgetc(f);
	if (c != '\n' && c != '/' && c != 'i'  && c != '+' && c != '-' && (c < '0' || c > '9'))
	{
		printf("inside INPUTSR: illegal character %c encountered:\n", c);
		FREEMPI(G);
		*uptr = 0;
		return ZEROR();
	}
	if (c == '/')
        {
	 	H = FINPUTI(f, uptr);
		if (*uptr == 0)
		{
			FREEMPI(G);
			FREEMPI(H);
			return ZEROR();
		}
		else if (H->S == 0)
		{
			printf("zero denominator entered\n");
			*uptr = 0;
			FREEMPI(G);
			FREEMPI(H);
			return ZEROR();
		}
		else
		{
			Aptr = RATIOI(G, H);
			FREEMPI(G);
			FREEMPI(H);
			c = fgetc(f);
			if (c != '\n')
				ungetc(c, f);
		}
	}
        else 
	{
		if (c != '\n' || t == ' ')
			ungetc(c, f);
		Aptr = BUILDMPR();
		Aptr->N = G;
		Aptr->D = ONEI();
	}
	return (Aptr);
}

MPR *INPUTR(unsigned int *uptr)
{
	return FINPUTR(stdin, uptr);
}

MPR *INPUTSR(char **ptr, unsigned int *uptr)
/*
 * Converts the ratio of two decimal inputs from stream into an MPR. 
 * *uptr = 0 if input fails, 1 if successful. For use with complex rationals.
 */
{
	char t;
	MPI *G, *H;
	MPR *Aptr;

	G = INPUTSI(ptr, uptr);
	if (*uptr == 0) 
	{
		FREEMPI(G);
		return ZEROR();
	}
	t = **ptr;
	while (**ptr == ' ')
		(*ptr)++;
	if (**ptr != '\0' && **ptr != '/' && **ptr != 'i'  && **ptr != '+' && **ptr != '-' && (**ptr < '0' || **ptr > '9'))
	{
		printf("inside INPUTSR: illegal character %c entered:\n", **ptr);
		*uptr = 0;
		FREEMPI(G);
		return ZEROR();
	}
	else if (**ptr == '/')
	{
		(*ptr)++;
		H = INPUTSI(ptr, uptr);
		if (*uptr == 0)
		{
			FREEMPI(G);
			FREEMPI(H);
			return ZEROR();
		}
		else if (H->S == 0)
		{
			printf("zero denominator entered\n");
			*uptr = 0;
			FREEMPI(G);
			FREEMPI(H);
			return ZEROR();
		}
		else
		{
			Aptr = RATIOI(G, H);
			FREEMPI(G);
			FREEMPI(H);
		}
	}
	else
	{
		if (t == ' ')
			(*ptr)--;
		Aptr = BUILDMPR();
		Aptr->N = G;
		Aptr->D = ONEI();
	}
	return (Aptr);
}

void FPRINTR(FILE *outfile, MPR *Aptr)
/*
 * prints the MPR *Aptr as (Aptr->N)/(Aptr->D). 
 */
{
	FPRINTI(outfile, Aptr->N);
	if ((Aptr->D)->D || (Aptr->D)->V[0] > 1)
	{
		fprintf(outfile, "/");
		FPRINTI(outfile, Aptr->D);
	}
	return;
}

void PRINTR(MPR *Aptr)
{
	FPRINTR(stdout, Aptr);
	return;
}

MPR *ZEROR()
/*
 * Returns the MPR ZERO.
 */
{
	MPR *Eptr;

	Eptr = BUILDMPR();
	Eptr->N = ZEROI();
	Eptr->D = ONEI();
	return Eptr;
}

MPR *ONER()
/*
 * Returns the MPR ONE.
 */
{
	MPR *Eptr;

	Eptr = BUILDMPR();
	Eptr->N = ONEI();
	Eptr->D = ONEI();
	return Eptr;
}

MPR *TWOR()
     /* Returns the MPR TWO */
{
  MPR *Eptr;
  Eptr = BUILDMPR();
  Eptr->N = TWOI();
  Eptr->D = ONEI();
  return Eptr;
}
MPR *MINUS_ONER()
/*
 * Returns the MPR MINUS_ONE.
 */
{
	MPR *Eptr;

	Eptr = BUILDMPR();
	Eptr->N = MINUS_ONEI();
	Eptr->D = ONEI();
	return (Eptr);
}

void SPRINTR(char *buffer, MPR *Aptr)
/*
 * prints the MPR *Aptr as (Aptr->N)/(Aptr->D). 
 */
{
	SPRINTI(buffer, Aptr->N);
	if ((Aptr->D)->D || (Aptr->D)->V[0] > 1)
	{
		buffer += strlen(buffer);
		sprintf(buffer, "/");
		SPRINTI(buffer + 1, Aptr->D);
	}
	return;
}

void FPRINTDR(FILE *outfile, unsigned int d, MPR *Mptr)
/*
 * prints *Mptr truncated to d (>= 1) decimal places to outfile.
 */
{
	int s;
	unsigned int i;
	unsigned long l;
	MPI *F, *G, *TempI;
	MPR *X, *TempR;

	s = (Mptr->N)->S;
	X = COPYR(Mptr);
	if (s == -1)
	{
		TempR = X;
		X = MINUSR(X);	
		FREEMPR(TempR);
	}
	G = INT0(X->N, X->D);
	if (s == -1)
		fprintf(outfile, "-");
	FPRINTI(outfile, G);
	fprintf(outfile, ".");
	FREEMPI(G);
	G = MOD0(X->N, X->D);
	F = POWER_I(10L, d);
	TempI = F;
	F = MULTI(F, G);
	FREEMPI(TempI);
	TempI = G;
	G = INT0(F, X->D);
	FREEMPI(F);
	FREEMPI(TempI);
	l = LENGTHI(G);
	for (i = 1; i <= d - l; i++)
		fprintf(outfile, "0");
	FPRINTI(outfile, G);
	FREEMPI(G);
	FREEMPR(X);
	return;
}

void PRINTDR(unsigned int d, MPR *Mptr)
/*
 * prints *Mptr truncated to d (>= 1) decimal places to stdout.
 */
{
	FPRINTDR(stdout, d, Mptr);
	return;
}

void SPRINTDR(char *buffer, unsigned int d, MPR *Mptr)
/*
 * prints *Mptr truncated to d (>= 1) decimal places to buffer.
 */
{
	int s;
	unsigned int i;
	unsigned long l;
	MPI *F, *G, *TempI;
	MPR *X, *TempR;

	s = (Mptr->N)->S;
	X = COPYR(Mptr);
	if (s == -1)
	{
		TempR = X;
		X = MINUSR(X);	
		FREEMPR(TempR);
	}
	G = INT0(X->N, X->D);
	if (s == -1)
	{
		sprintf(buffer, "-");
		buffer++;
	}
	SPRINTI(buffer, G);
	buffer += strlen(buffer);
	sprintf(buffer, ".");
	buffer++;
	FREEMPI(G);
	G = MOD0(X->N, X->D);
	F = POWER_I(10L, d);
	TempI = F;
	F = MULTI(F, G);
	FREEMPI(TempI);
	TempI = G;
	G = INT0(F, X->D);
	FREEMPI(TempI);
	FREEMPI(F);
	l = LENGTHI(G);
	for (i = 1; i <= d - l; i++)
	{
		sprintf(buffer, "0");
		buffer++;
	}
	SPRINTI(buffer, G);
	FREEMPI(G);
	FREEMPR(X);
}

unsigned long LENGTHDR(unsigned int d, MPR *Mptr)
/*
 * returns LENGTHI(int(*Mptr)) + d + 1, if *Mptr >= 0.
 * returns LENGTHI(int(abs(*Mptr))) + d + 2, if *Mptr < 0.
 */
{
	int s;
	unsigned long l;
	MPI *G;
	MPR *X, *TempR;

	s = (Mptr->N)->S;
	X  = COPYR(Mptr);
	if (s == -1)
	{
		TempR = X;
		X = MINUSR(X);	
		FREEMPR(TempR);
	}
	G = INT0(X->N, X->D);
	l = LENGTHI(G);
	FREEMPR(X);
	FREEMPI(G);
	if (s == -1)
		return (l + d + (USL)2);
	else
		return (l + d + (USL)1);
}

MPR *RECIPROCAL(unsigned long n)
/*
 * input: unsigned int n, 0 < n < R0, output:  MPR *Aptr = 1 / n. 
 */
{
	MPR *Aptr;

	Aptr = BUILDMPR();
	Aptr->N = ONEI();
	Aptr->D = CHANGE(n);
	return (Aptr);
}

MPR *ADDR(MPR *Aptr, MPR *Bptr)
/*
 * *Eptr = *Aptr + *Bptr.
 * P. Henrici's method: see Computer Algebra Symbolic and Algebraic
 * Computation, Springer 1982,p. 200.
 */
{
	MPI *D, *E, *R, *S, *T1, *T2, *X, *Y;
	MPR *Eptr;

	if ((Aptr->N)->S == 0)