📄 i5i.c
字号:
fprintf(outfile, "P[%lu]=", i);
FPRINTI(outfile, (*P)->A[i]);
fprintf(outfile, ";Q[%lu]=", i);
FPRINTI(outfile, (*Q)->A[i]);
fprintf(outfile, "\n");
}
fclose(outfile);
FREEMPI(ONE);
return;
}
MPI *SUMI(MPI *A[], USL n)
/* Returns A[0]+...+A[n] */
{
MPI *S, *temp;
USL i;
S = ZEROI();
for (i = 0; i <= n; i++)
{
temp = S;
S = ADDI(S, A[i]);
FREEMPI(temp);
}
return (S);
}
void LAGRANGE(POLYI P, MPIA *AA, MPI *M)
/*
* P(x)=A[n]x^n+...+A[0], A[n]>0, is a polynomial
* with integer coefficients, having no rational roots
* and having exactly one real positive root x, this being > 1.
* The method of Lagrange (1797) is used to find the
* the first m+1 partial quotients AA[0],...,AA[m] of x.
* (See Knuth, Art of computer programming, volume2,
* problem 13, 4.5.3. Also S. Lang and H. Trotter,
* 'Continued fractions for some algebraic numbers'
* J. fuer Math. 255 (1972) 112-134; Addendum 267 (1974) ibid. 219-220.)
* Also page 261, Number Theory with Applications, by R. Kumanduri and
* C. Romero.
*/
{
int i,j,m,d;
MPR *HIGH;
MPR *ONE = ONER();
MPI *ZERO = ZEROI();
POLYI Q = COPYPI(P), CURSOR;
MPI *TMP, *HI, *LO, *ON, *TEMP, *TEMP1, *ROOT, *MID;
MPIA COEF;
ON=ONEI();
m = CONVERTI(M);
d = DEGREEPI(P);
/* I should a check in here for a single positive root. I will get around
* to it */
*AA = BUILDMPIA();
for(i = 0;i <= m; i++){
LO=ONEI();
HIGH = CAUCHY(Q); /* this is a power of 2 */
HI=COPYI(HIGH->N);
FREEMPR(HIGH);
TEMP=ADD0I(LO,ON);
while(RSV(TEMP,HI)){
FREEMPI(TEMP);
TEMP1=ADD0I(LO,HI);
MID=INT0_(TEMP1,2);
FREEMPI(TEMP1);
TEMP1=VALPI(Q, MID);
if(TEMP1->S>0){
TEMP=HI;
HI=MID;
FREEMPI(TEMP);
}else{
TEMP=LO;
LO=MID;
FREEMPI(TEMP);
}
FREEMPI(TEMP1);
TEMP=ADD0I(LO,ON);
}
FREEMPI(TEMP);
ROOT=LO;
FREEMPI(HI); /* HI=LOW+1 here */
ADD_TO_MPIA(*AA, ROOT , i);
P_OF_X_PLUS(&Q, ROOT);
FREEMPI(ROOT);
/* This performs transformation p(x) = -x^d*p(1/x +a) */
COEF = BUILDMPIA();
for (j=0; j<=d; j++)
ADD_TO_MPIA(COEF, ZERO, j);
CURSOR = Q;
while(CURSOR) {
TMP = MINUSI(CURSOR->COEF);
ADD_TO_MPIA(COEF, TMP, DEGREEPI(CURSOR));
FREEMPI(TMP);
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 */
if(DEGREEPI(Q)<d){/* noticed by KRM mid Sept 2002 */
printf("rational root\n");
break;
}
}
DELETEPI(Q);
FREEMPI(ON);
FREEMPR(ONE);
FREEMPI(ZERO);
}
void EUCLID(MPI *Aptr, MPI *Bptr, MPIA *Q, MPIA *R, MPIA *S, MPIA *T, MPI **Dptr)
/*
* Returns Q[0]=NULL,Q[1],...Q[n],Q[n+1]=NULL,
* R[0],...R[n + 1],
* S[0],...S[n + 1], T[0],...T[n + 1]
* for Euclid's algorithm on R[0]=Aptr, R[1]=Bptr.
* R[0]=R[1]*Q[1]+R[2]
* R[1]=R[1]*Q[2]+R[3]
* .....
* R[n-2]=R[n-1]*Q[n-1]+R[n]
* R[n-1]=R[n]*Q[n], R[n+1]=0.
* S[0]=1,S[1]=0, S[j]=s[j-2]-Q[j-1]*S[j-1],
* T[0]=0,T[1]=1, T[j]=T[j-2]-Q[j-1]*T[j-1], j=2,...,n+1
* Here *Dptr = n.
*/
{
MPI *QQ;
USL n, i, j;
char buff[20];
FILE *outfile;
MPI *ONE= ONEI();
MPI *ZERO = ZEROI();
MPI *TMPI;
strcpy(buff, "euclid.out");
if(access(buff, R_OK) == 0)
unlink(buff);
outfile = fopen(buff, "a");
n = 2 * BINARYB(Bptr) + 1;
(*Q) = BUILDMPIA();
(*R) = BUILDMPIA();
(*S) = BUILDMPIA();
(*T) = BUILDMPIA();
ADD_TO_MPIA((*R), Aptr, 0);
ADD_TO_MPIA((*R), Bptr, 1);
ADD_TO_MPIA((*S), ONE, 0);
ADD_TO_MPIA((*S), ZERO, 1);
ADD_TO_MPIA((*T), ZERO, 0);
ADD_TO_MPIA((*T), ONE, 1);
ADD_TO_MPIA(*Q, NULL, 0);
TMPI = INT0((*R)->A[0], (*R)->A[1]);
ADD_TO_MPIA((*Q), TMPI, 1);
FREEMPI(TMPI);
TMPI = MOD0((*R)->A[0], (*R)->A[1]);
ADD_TO_MPIA((*R), TMPI, 2);
FREEMPI(TMPI);
i = 2;
while ((*R)->A[i]->S == 1)
{
QQ = MINUSI((*Q)->A[i - 1]);
TMPI = MULTABC((*S)->A[i - 2], QQ, (*S)->A[i - 1]);
ADD_TO_MPIA(*S, TMPI, i);
FREEMPI(TMPI);
TMPI = MULTABC((*T)->A[i - 2], QQ, (*T)->A[i - 1]);
ADD_TO_MPIA(*T, TMPI, i);
FREEMPI(TMPI);
FREEMPI(QQ);
TMPI = INT0((*R)->A[i - 1], (*R)->A[i]);
ADD_TO_MPIA(*Q, TMPI, i);
FREEMPI(TMPI);
TMPI = MOD0((*R)->A[i - 1], (*R)->A[i]);
ADD_TO_MPIA(*R, TMPI, i+1);
FREEMPI(TMPI);
i++;
}
QQ = MINUSI((*Q)->A[i - 1]);
TMPI = MULTABC((*S)->A[i - 2], QQ, (*S)->A[i - 1]);
ADD_TO_MPIA(*S, TMPI, i);
FREEMPI(TMPI);
TMPI = MULTABC((*T)->A[i - 2], QQ, (*T)->A[i - 1]);
ADD_TO_MPIA(*T, TMPI, i);
FREEMPI(TMPI);
FREEMPI(QQ);
ADD_TO_MPIA(*Q, NULL, i);
i++;
*Dptr = CHANGE(i - 1);
for (j = 0; j < i; j++){
if (j == 0 || j == i - 1)
fprintf(outfile, "Q[%lu]=-", j);
else
{
fprintf(outfile, "Q[%lu]=", j);
FPRINTI(outfile, (*Q)->A[j]);
}
fprintf(outfile, ", R[%lu]=", j); FPRINTI(outfile, (*R)->A[j]);
fprintf(outfile, ", S[%lu]=", j); FPRINTI(outfile, (*S)->A[j]);
fprintf(outfile, ", T[%lu]=", j); FPRINTI(outfile, (*T)->A[j]);
fprintf(outfile, " \n");
}
fclose(outfile);
FREEMPI(ZERO);
FREEMPI(ONE);
return;
}
/* Allocates space for an array initially of size 11 (enough to hold a[0] to
* a[10] and sets these slots to contain the zero MPI*/
MPIA BUILDMPIA()
{
MPIA MA;
int i;
if (!(MA = mmalloc(sizeof(struct _MPIA)))) {
printf("Not enough memory to allocate MPIA\n");
exit(1);
}
MA->A = (MPI **)mmalloc((USL)((MIN_ARRAY_SIZE)*sizeof(MPI *)));
MA->size = 0;
MA->slots = MIN_ARRAY_SIZE;
for (i=0; i < MIN_ARRAY_SIZE; i++)
MA->A[i] = ZEROI();
return MA;
}
/* ADD_TO_MPIA
* Adds the supplied MPI at the subscript n
* If slot already exists MPI at that slot is freed and the new one is added.
* If the n is greater than the number of slots then the array is correctly
* resized and the new MPI is added. Slots between the previous last slots
* and the new subscript n are initialised to zero.
*/
void ADD_TO_MPIA(MPIA MA, MPI *V, USL n)
{
int i;
if (n > MAX_ARRAY_SIZE-1) {
printf("Error: Array index must be in range 0 < n < %d\n",
MAX_ARRAY_SIZE-1 );
return;
}
if (n > MA->slots-1) {
MA->A = (MPI **)rrealloc(MA->A,(n+1)*sizeof(MPI *),
(n + 1 - (MA->slots)) * sizeof(MPI*));
for (i=MA->slots; i < n+1; i++)
(MA->A)[i] = ZEROI();
MA->slots = n+1;
}
FREEMPI(MA->A[n]);
if (V != NULL)
MA->A[n] = COPYI(V);
else MA->A[n] = NULL;
if (n+1 > MA->size)
MA->size = n+1;
}
/* Frees an MPIA previously returned by BUILDMPIA. Undefined behaviour occurs
* if this is not the case.
* It will free all the slots of the MPI's they contain.
*/
void FREEMPIA(MPIA MA)
{
int i;
if (!MA)
return;
for (i=0; i < MA->slots; i++)
if (MA->A != NULL)
FREEMPI(MA->A[i]);
ffree(MA->A, (MA->slots)*sizeof(MPI *));
ffree(MA, sizeof(struct _MPIA));
}
/* Inserts the term V at position n. If n > size of array - 1 then this
* function behaves the same as ADD_TO_MPIA
* By insert I mean that the object is placed in at that subscript and the
* rest are shuffled up.
*/
void MPIA_INSERT(MPIA MA, MPI *V, USL n)
{
/* add an extra slot */
if (n > MA->size -1)
ADD_TO_MPIA(MA, V, n);
else { /* n < MA->size */
if (MA->slots == MA->size) { /* then realloc space */
if (!(MA->A = rrealloc(MA->A, (++MA->slots) * sizeof(MPI *), \
sizeof (MPI *)))) {
printf("Not enough memory to increase size of MPIA\n");
exit(1);
}
memmove(MA->A+n+1, MA->A+n, sizeof(MPI *) * (MA->slots-n-1));
} else { /* no need to realloc */
FREEMPI(MA->A[MA->size]);
memmove(MA->A+n+1, MA->A+n, sizeof(MPI *) * (MA->size-n));
}
MA->size++;
MA->A[n] = COPYI(V);
}
}
MPI *EUCLIDI1(MPI *M, MPI *N)
/*
Returns the length of Euclid's algorithm
*/
{
MPI *I, *A, *B, *R, *TEMP, *TEMP1, *T1, *T2, *T3;
USL i, k;
k = 0;
i = 1;
A = COPYI(M);
B = COPYI(N);
R = MOD0(A, B);
while (R->S > 0){
TEMP = A;
A = B;
FREEMPI(TEMP);
B = R;
R = MOD0(A, B);
i++;
if(i%1000==0){
printf("i=%lu\n", i);
}
if(i == 2147483647){
i = 0;
k++;
if(k == 2147483647){
printf("exiting prematurely, k = 2147483647\n");
return(NULL);
}
}
}
FREEMPI(A);
FREEMPI(B);
FREEMPI(R);
T1 = CHANGE(k);
T2 = CHANGE(2147483648UL);
T3 = CHANGE(i);
TEMP1 = MULTI(T1,T2);
I = ADD0I(TEMP1,T3); /*I = k*2147483648 + i */
FREEMPI(T1);
FREEMPI(T2);
FREEMPI(T3);
FREEMPI(TEMP1);
return (I);
}
MPI *CFRACN(MPI *N){
USL n;
MPI *L, *Z, *ONE, *TWO, *AA, *BB, *TEMP;
system("date");
n = CONVERTI(N);
Z = POWER_I(2, n);
ONE = ONEI();
TWO = TWOI();
printf("reached here\n");
POWERD(ONE,ONE,TWO,Z,&AA,&BB);
printf("completed exponentation\n");
TEMP = BB;
BB = INT0(BB, Z);
FREEMPI(TEMP);
L = EUCLIDI1(AA, BB);
FREEMPI(Z);
FREEMPI(ONE);
FREEMPI(TWO);
FREEMPI(AA);
FREEMPI(BB);
if((L->V[0]) % 2){
TEMP = L;
L = ADD0_I(L, 1);
FREEMPI(TEMP);
}
system("date");
return(L);
}
MPI *CFRAC_PERIOD(MPI *D)
/*
* This is a program for finding the period of the cfrac of sqrt(D).
* The algorithm is based on K. Rosen, Elementary number theory
* and its applications, p382, B.A. Venkov, Elementary Number theory, p.62
* and D. Knuth, Art of computer programming, Vol.2, p359, with Pohst's trick
* of using half the period.
*/
{
MPI *P, *Q, *E, *F, *G, *X, *Y, *TEMP, *TEMP1, *H, *PERIOD;
MPI *T2, *T3;
USL h, k;
USI t;
system("date");
X = BIG_MTHROOT(D, 2);
G = MULTI(X, X);
TEMP = ADD0_I(G, (USL)1);
FREEMPI(G);
t = EQUALI(D, TEMP);
FREEMPI(TEMP);
if(t) /* period 1, exceptional case */
{
FREEMPI(X);
return(ONEI());
}
P = ZEROI();
Q = ONEI();
h = 0;
k = 0;
while(1)
{
TEMP = ADD0I(X, P);
Y = INT0(TEMP, Q);
FREEMPI(TEMP);
F = P;
TEMP = MULTI(Y, Q);
FREEMPI(Y);
P = SUBI(TEMP, P);
FREEMPI(TEMP);
t = EQUALI(P, F);
FREEMPI(F);
if(t){ /*P_h=P_{h+1}, even period 2h */
printf("period=");
T2 = CHANGE(2147483648UL);
/*T2 = CHANGE(32768);*/
T3 = CHANGE(h);
TEMP1 = MULT_I(T2,k);
H = ADD0I(TEMP1,T3); /*H = k*2147483648 + h */
FREEMPI(T2);
FREEMPI(T3);
FREEMPI(TEMP1);
PERIOD = MULT_I(H, 2);
FREEMPI(H);
FREEMPI(X);
FREEMPI(P);
FREEMPI(Q);
PRINTI(PERIOD);
printf("\n");
system("date");
return(PERIOD);
}
E = Q;
TEMP = MULTI(P, P);
TEMP1 = SUBI(D, TEMP);
FREEMPI(TEMP);
Q = INT0(TEMP1, Q);
t = EQUALI(Q, E);
FREEMPI(TEMP1);
FREEMPI(E);
if(t){ /*Q_h=Q_{h-1}, odd period 2h+1 */
printf("period=");
T2 = CHANGE(2147483648UL);
/*T2 = CHANGE(32768);*/
T3 = CHANGE(h);
TEMP1 = MULT_I(T2,k);
H = ADD0I(TEMP1,T3); /*H = k*2147483648 + h */
FREEMPI(T2);
FREEMPI(T3);
FREEMPI(TEMP1);
TEMP = MULT_I(H, 2);
FREEMPI(H);
PERIOD = ADD0_I(TEMP, 1);
FREEMPI(TEMP);
FREEMPI(X);
FREEMPI(P);
FREEMPI(Q);
PRINTI(PERIOD);
printf("\n");
system("date");
return(PERIOD);
}
h++;
if(h == 2147483648UL){
h = 0;
k++;
if(k == 2147483648UL){
printf("exiting prematurely, k = 2147483648\n");
return(NULL);
}
}
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -