www.pudn.com > calc.rar > cyclo.c

/* cyclo.c */ 
 
#include  
#include  
#include "integer.h" 
#include "fun.h" 
#define S0 65536 
 
void NEXTPHI(USL p, USL *gptr, MPIA *b, MPIA *coef) 
{ 
	USL cycle, grad, t, lastgrad; 
        int i, j, next; 
	MPI *l, *ZERO, *TEMP, *w; 
 
	ZERO = ZEROI(); 
	cycle=p; 
	grad = *gptr; 
	next=grad; 
	 
	for(i=grad; i>=0; i--){ 
		  if(cycle==p) 
		  { 
		    cycle=1; 
		    TEMP = COPYI(((*coef)->A)[next]); 
		    ADD_TO_MPIA(*b, TEMP, i); 
		    FREEMPI(TEMP); 
		    next=next-1; 
		  }else{ 
		    cycle++; 
		    ADD_TO_MPIA(*b, ZERO, i); 
		  } 
	} 
	lastgrad=grad; 
	grad=grad*(p-1); 
 
	t=grad/2; 
	for(i=grad; i>=t; i--){ 
		  /* PHI = PHI(x^p)/PHI.  First half of the coefficients. */ 
		  l = COPYI(((*b)->A)[lastgrad]); 
		  ADD_TO_MPIA(*coef, l, i); 
		  for(j=lastgrad-1; j>=0; j--) 
		  { 
		    w = MULTI(l, ((*coef)->A)[j]); 
		    TEMP = SUBI(((*b)->A)[j], w); 
		    FREEMPI(w); 
		    ADD_TO_MPIA(*b, TEMP, j+1); 
		    FREEMPI(TEMP); 
		  } 
		  FREEMPI(l); 
		  if(cycle==p) 
		  { 
		    cycle=1; 
		    TEMP = COPYI(((*coef)->A)[next]); 
		    ADD_TO_MPIA(*b, TEMP, 0); 
		    FREEMPI(TEMP); 
	/*	    next--;*/ /*commented out by KRM on 1h26/9/02 */ 
		  } 
		  else 
		  { 
		    cycle++; 
		    ADD_TO_MPIA(*b, ZERO, 0); 
		  } 
	} 
	FREEMPI(ZERO); 
 
	t=grad/2 - 1; 
	for(j = 0; j <= t; j++){ 
            TEMP=COPYI(((*coef)->A)[grad-j]); 
            ADD_TO_MPIA(*coef, TEMP, j); 
            FREEMPI(TEMP); 
	} 
	*gptr = grad; 
} 
 
int NOTYET(USL n, USL i, USL *qptr) 
{ 
	USL q; 
 
	q = n / i; 
	*qptr = q; 
	if ((n > q * i) && (i <= q)) 
		return 1; 
	else 
		return 0; 
} 
 
int LPF(USL n) 
{ 
	USL c, i, q; 
 
	if (n % 2 == 0) return (2); 
	if (n % 3 == 0) return (3); 
	else 
	{ 
		i = 5; 
		c = 0; 
		while (NOTYET(n, i, &q)) 
		{ 
			i = (c == 0) ? i + 2 : i + 4; 
			c = 1 - c; 
		} 
		if (q < i) 
			return (n); 
		else 
			return (i);	 
	} 
} 
 
POLYI CYCLOTOMIC(MPI *N) 
/* This algorithm is from Heinz Luneburg's book Galoisfelder, Kreisteilungs 
 * korper and Schieberegisterfolgen. It was originally coded in CMAT by perhaps Peter Adams. 
 * I imported the code to CALC on 25-26 Seotember 2002, tightening up the code. 
 */ 
 
{ 
	USL n, n1, socn, p, grad; 
	int i; 
	POLYI P; 
	MPI *H, *TEMP, *MINUSONE, *ONE, *ZERO; 
	MPIA b, coef; 
 
	if(N->D){ 
		printf("N >= 65536\n"); 
		return (POLYI)NULL; 
	} 
	MINUSONE = MINUS_ONEI(); 
	ONE = ONEI(); 
	ZERO = ZEROI(); 
 
	n = CONVERTI(N); 
	b = BUILDMPIA(); 
        coef = BUILDMPIA(); 
	n1 = n; 
	if (n1 == 1) 
	{ 
		grad = 1; 
		ADD_TO_MPIA(coef, MINUSONE, 0); 
		ADD_TO_MPIA(coef, ONE, 1); 
	} 
	else 
	{ 
		p = LPF(n1); 
		socn = p; 
		while (n1 % p == 0) 
			n1 = n1 / p; 
		grad = p - 1; 
		for (i = 0; i <= grad; i++){ 
			ADD_TO_MPIA(coef, ONE, i); 
		} 
		/* PHI is the pth cyclotomic polynomial. */ 
			 
		while (n1 > 1) 
		{ 
			p = LPF(n1); 
			socn = socn * p; 
			while (n1 % p == 0) 
				n1 = n1 / p; 
			NEXTPHI(p, &grad, &b, &coef); 
		/* PHI is the socnth cyclotomic polynomial. */		 
		} 
		n1 = n / socn; 
		if (n1 > 1) 
		{ 
		/* PHI = PHI(x^n). */ 
			grad = n1 * grad; 
			for (i = grad; i >= 0; i--) 
			{ 
				if (i % n1 == 0){ 
					TEMP = COPYI((coef->A)[i/n1]); 
					ADD_TO_MPIA(coef, TEMP, i); 
					FREEMPI(TEMP); 
				} 
				else{ 
					ADD_TO_MPIA(coef, ZERO, i); 
				} 
			} 
		}	 
	} 
	FREEMPIA(b); 
	P = NULL; 
	for (i = 0; i <= grad; i++) 
	{ 
		if (((coef->A)[i])->S) 
		{ 
			H = COPYI((coef->A)[i]); 
			PINSERTPI((int)i, H, &P, 0); 
			FREEMPI(H); 
		} 
	} 
	FREEMPIA(coef); 
	FREEMPI(MINUSONE); 
	FREEMPI(ONE); 
	FREEMPI(ZERO); 
	return (P); 
}