/* copyright 1987,2003 David G. Long. All rights reserved. */ /* *** LAST REVISED ON 9-FEB-1987 12:26:00.06 */ /* *** SOURCE FILE: [DL.TXT.STURM]STURM.C */ #include #include #include #define TRUE 1 #define FALSE 0 #define HUGE 1.e33 /* largest positive float */ /* test program for modified McGowan's phase unwrapping */ #define KIMAX 100000000 /* max value of one word in K representation */ /* note KIMAX < 1/10 of largest integer in type long */ #define KDIGPW 8 /* digits per word */ #define KISQR 10000 /* sqrt of KIMAX */ #define KMDTR 20 /* minimum number of words for poly coef trunc */ #define KXDTR 10 /* maximum number of words to truncate poly coef */ #define MAXPTS 30 /* maximum number of points time sequence */ #define debug TRUE /* enable printout */ struct k_long { /* extended integer k representation */ long k_num, *k_p; }; typedef struct k_long* INTK ; /* global variables for storing modified Chebychev polynomial values */ INTK v_n[64], v_nn[MAXPTS]; main() /* time domain sequence is McGowan and Kuc's example */ { long x[] = {10000,-36000,48600,-29160,6561}; /* scaled by 10000 */ long n = 4; /* number of time samples < MAXPTS */ long m = 8; /* number of phase samples in 0 to Pi*/ unwrap(x,n,m); /* m is limited to powers of 2 <= 32 */ } unwrap(x,n,m) /* generate unwrapped phase */ long x[],n,m; { INTK p[MAXPTS * MAXPTS],c[MAXPTS]; long i,ns,n0,l,test[MAXPTS]; double chebs(),w,phi,phi1,r,atan(),sin(),pi = 3.141592654; sineinit(m,n); /* initialize sine value table */ sturm(x,n,p); /* generate sturm sequence */ for (i = 0; i < m; i++) { w = chebs(n,c,i,m,test); ns = vsign(p,n,c,&r,test); if (i == 0) n0 = ns; l = ns - n0; /* multiples of two pi for unwrapped phase */ if (i != 0) { phi1 = atan(2.0 * r * sin(w)); /* principle value */ phi = phi1 + l * pi - n * w; /* unwrapped phase */ } else { phi1 = 0.0; phi = l * pi; } printf("\n%d w = %g phi1 = %g l = %d phi = %g\n",i,w,phi1,l,phi,l); } } sturm(x,n,p) /* generate sturm sequence from time series */ long x[],n; INTK p[]; { INTK kint(); long i,k0,k1,k2,nt,nw,sizepoly(); /* first polynomial element of sturm sequence */ p[0] = kint(2 * x[n] - x[n - 2]); for(i = 1; i < n - 1; i++) p[i] = kint(x[n - i] - x[n - i - 2]); p[n-1] = kint(x[1]); p[n] = kint(x[0]); if (debug) {printf("\np0 %d\n",n); typepoly(&p[0],n);} /* second polynomial element of sturm sequence */ k1 = ipoff(1,n); for (i = 0; i < n; i++) p[k1 + i] = kint(x[n - i - 1]); if (debug) {printf("\np1 %d\n",n); typepoly(&p[k1],n - 1);} /* generate remainder of sturm sequence */ for (i = 2; i <= n; i++) { k0 = ipoff(i - 2,n); k1 = ipoff(i - 1,n); k2 = ipoff(i,n); nt = n - i + 2; divide(&p[k0],nt,&p[k1],&p[k2]); /* remove common factors in polynomial coefficients */ simple(&p[k2],nt - 2); if (debug) {printf("\np%d %d\n",i,n); typepoly(&p[k2],nt-2);} /* truncate off least significant words of polynomial if all coefficients are of minimum size -- not used in example */ if ((nw = sizepoly(&p[k2],nt - 2)) > KMDTR) { truncpoly(&p[k2],nt - 2,nw - KXDTR); if (debug) printf("Truncating Polynomial\n"); } } } divide(p0,n,p1,p2) /* chebychev division algorithm for polynomials */ INTK p0[],p1[],p2[]; { INTK q,r,t1,t2,t3,t4,kmul(),ksub(),kadd(); long j; /* compute q = p0[n]*p1[n-1], r = p0[n-1]*p1[n-1]-p0[n]*p1[n-2], p2[0] = r*p1[0]+q*p1[1]-p0[0]*p1[n-1]*p1[n-1] */ q = kmul(p0[n],p1[n - 1]); t1 = kmul(p0[n - 1],p1[n - 1]); t2 = kmul(p0[n],p1[n - 2]); r = ksub(t1,t2); kdeal(t1); kdeal(t2); if (debug) { printf("\n q = "); kout(q); printf(" r = "); kout(r); } t1 = kmul(r,p1[0]); t2 = kmul(q,p1[1]); t3 = kadd(t1,t2); kdeal(t1); kdeal(t2); t4 = kmul(p1[n - 1],p1[n - 1]); t1 = kmul(t4,p0[0]); p2[0] = ksub(t3,t1); kdeal(t1); kdeal(t3); if (n == 0) return; /* compute p2[j]=q*(p1[j-1]+p1[j+1])+r*p1[j]-p0[j]*p1[n-1]*p1[n-1] */ for (j = 1;j < n - 1;j++) { t1 = kadd(p1[j - 1],p1[j + 1]); t2 = kmul(q,t1); kdeal(t1); t1 = kmul(r,p1[j]); t3 = kadd(t2,t1); kdeal(t2); kdeal(t1); t1 = kmul(p0[j],t4); p2[j] = ksub(t3,t1); kdeal(t3); kdeal(t1); } kdeal(t4); } ipoff(i,n) /* comput offset of ith polynomial */ long i,n; { long k, j = 0; for(k = 1; k <= i; k++) j = j + (n - k + 2); return(j); } double chebs(n,c,j,m,test) /* returns list modified chebychev polynomials */ long n,j,m,test[]; INTK c[]; { double w,pi = 3.14159265386; long i,k; w = j * pi / m; if(debug)printf("\nModified Chebychev Polynomials at %d of %d w=%g:\n" ,j,m,w); for(i = 0;i <= n; i++) { k = ((i + 1) * j) % (2 * m); if (j == 0) { /* w=0 */ test[i] = 0; /* computation exact, accuracy test unneeded */ c[i] = v_nn[i]; } else { /* w<>0 */ if (k > m) { k = 2 * m - k; if (k > m / 2) k = m - k; k = k + m / 2; } else if (k > m / 2) k = m - k; if (k == m / 2 || k == 0 || k == m) test[i] = 0; else test[i] = 1; /* accuracy test if v <> 1 or 0 */ c[i] = v_n[k]; } if (debug) {printf(" c[%d],%d : ",i,k); kout(c[i]);} } return(w); } vsign(p,n,c,r,test) /* compute number of sign changes */ INTK p[],c[]; /* in sturm sequence */ long n,test[]; double *r; { long i,is0,is1,k,ns; INTK evalpoly(),pw,p0,p1; double krmul(),krratio(); p0 = evalpoly(&p[0],c,n,test); k = ipoff(1,n); p1 = evalpoly(&p[k],c,n - 1,test); if (debug) { printf("\nP0(w) = "); kout(p0); printf("P1(w) = "); kout(p1); } /* signs of first two values in sturm sequence */ is0 = ksign(p0); is1 = ksign(p1); *r = krratio(p1,p0); /* ratio p1/p0 for principle value */ kdeal(p1); kdeal(p0); /* count sign changes */ ns = 0; if (is0 * is1 < 0) ns = 1; for (i = 2;i <= n; i++) { k = ipoff(i,n); pw = evalpoly(&p[k],c,n - i,test); if (debug) { printf("P%d(w) = ",i); kout(pw); } is0 = ksign(pw); kdeal(pw); if (is0 * is1 < 0) ns = ns + 1; is1 = is0; } return(ns); } INTK evalpoly(p,c,n,test) /* evaluate polynomial */ INTK p[],c[]; long n,test[]; { INTK v,u,w,kmul(),kadd(),a,ktok(),kneg(); long i,nt; nt = 0; for(i = 0;i <= n; i++) { u = kmul(p[i],c[i]); if (i == 0) v = u; else{ w = kadd(v,u); kdeal(u); kdeal(v); v = w; } if (test[i]) { /* compute maximum error */ if (nt == 0) { /* test[i]=0 indicates cheby eval exact */ a = ktok(p[i]); if (ksign(a) < 0) a = kneg(a); } else { if (ksign(p[i]) < 0) u = kneg(kadd(p[i],kneg(a))); else u = kadd(p[i],a); kdeal(a); a = u; } nt++; } } /* check maximum error against sum */ if (nt != 0) { u = ktok(v); if (ksign(u) < 0) u = kneg(u); if (klt(u,a)) { printf("*** error threshold check failed in evalpoly ***\n"); printf("max error = "); kout(a); printf("sum = "); kout(v); } kdeal(u); kdeal(a); } return(v); } long sizepoly(p,n) /* determine minimum size of poly coefs */ INTK p[]; long n; { long i,cnt; cnt = knum(p[0]); for(i = 1;i <= n; i++) cnt = min(knum(p[i]),cnt); return(cnt); } truncpoly(p,n,nw) /* truncate of nw least significant words */ INTK p[]; long n; /* from each coefficient in polynomial */ { long i,cnt,sizepoly(); INTK ktrunc(); if (sizepoly(p,n) < nw) return; for(i = 0;i <= n; i++) p[i] = ktrunc(p[i],nw); } typepoly(p,n) /* type polynomial coefficients out */ INTK p[]; long n; { long i; for(i = 0;i <= n; i++) { printf("%d: ",i); kout(p[i]); } } simple(p,n) /* factor out common factors from polynomial */ INTK p[]; long n; /* coefficients */ { long i,j,ok; long v[] = {2,3,5,7,11,0}; /* consider only some prime factors */ INTK t, kidiv(), kint(); /* special case of only one coefficient -- return sign of coefficient */ if (n == 0) { i = ksign(p[0]); t = kint(i); kdeal(p[0]); p[0] = t; return; } /* check for common factors */ j = 0; while (v[j] != 0) { ok = TRUE; for(i = 0;i <= n;i++) if (kimod(p[i],v[j])) ok = FALSE; if (ok) { /* remove common factor v[j] */ for(i = 0;i <= n;i++) { t = kidiv(p[i],v[j]); if (t != p[i]) kdeal(p[i]); p[i] = t; } } else j++; } } min(i,j) /* returns the minimum of i and j */ long i,j; { if (i < j) return(i); else return(j); } sineinit(m,n) /* generate sine value table */ long m,n; /* m = number of phase values from 0 to pi */ { long i,j; /* n = number of points */ INTK kint(),kstr(),ktok(),kneg(),kimul(); char *s[17]; /* integerized sine table for D=28 */ s[1] = "009801714032956060199419556389"; /* sin(pi/32) */ s[2] = "019509032201612826784828486848"; /* sin(pi/16) */ s[3] = "029028467725446236763619237582"; /* sin(3pi/32)*/ s[4] = "038268343236508977172845998403"; /* sin(pi/8) */ s[5] = "047139673682599764855638762591"; /* sin(5pi/32) */ s[6] = "055557023301960222474283081395"; /* sin(3pi/16 */ s[7] = "063439328416364549821517161323"; /* sin(7pi/32) */ s[8] = "070710678118654752440084436210"; /* sin(pi/4) */ s[9] = "077301045336273696081090660976"; /* sin(9pi/32) */ s[10] = "083146961230254523707878837762"; /* sin(5pi/16) */ s[11] = "088192126434835502971275686366"; /* sin(11pi/32) */ s[12] = "092387953251128675612818318940"; /* sin(3pi/8) */ s[13] = "095694033573220886493579788698"; /* sin(13pi/32) */ s[14] = "098078528040323044912618223613"; /* sin(7pi/16) */ s[15] = "099518472667219688624483695311"; /* sin(15pi/32) */ s[16] = "100000000000000000000000000000"; /* sin(pi/2) */ v_n[0] = kint(0); /* sin(0) */ for (i = 1; i <= m/2; i++) { j = i * (32 / m); v_n[i] = kstr(s[j]); /* positive 1/4 sine wave */ v_n[i+m/2] = kneg(kstr(s[j])); /* negative 1/4 sine wave */ } for (i = 0; i <= n; i++) v_nn[i] = kimul(v_n[m/2],i+1); /* w=0 cases */ } /* computational routines used for extended integer (k integer) */ double krratio(n,d) /* provides n/d as a real number */ INTK n,d; /* destroys contents of n and d integers in the process */ { long nl,dl; INTK ktrunc(); double r,nr,dr,rp,krmul(); if (n == NULL || d == NULL) { kerr(0,"krratio"); return(0.0); } else if (kzero(d)) { kerr(3,"krratio"); return(HUGE); /* largest real value */ } /* reduce length of extended integers to managable size */ nl = knum(n) - 2; dl = knum(d) - 2; if (nl > 0) n = ktrunc(n,nl); if (dl > 0) d = ktrunc(d,dl); nr = krmul(n,1.0); dr = krmul(d,1.0); rp = (nl - dl) * KDIGPW; r = nr/dr * pow(10.0,rp); /* compute ratio */ return(r); } kerr(err,routine) /* k integer error handling routine */ long err; char *routine; { printf("*** Error %d has occured in routine %s\n",err,routine); if (err == 0) printf("*** NULL K ***\n"); else if (err == 1) printf("*** Out of K space ***\n"); else if (err == 2) printf("*** Out of string space ***\n"); else if (err == 4) printf("*** Division by zero I ***\n"); else if (err == 5) printf("*** Internal Error ***\n"); else if (err == 6) printf("*** NULL string input ***\n"); else printf("*** UNKNOWN ERROR CODE ***\n"); return(-99); } isign(i) long i; /* integer sign function */ { if (i < 0) return(-1); else if (i == 0) return(0); else return(1);} kout(k) /* print out K integer */ INTK k; /* output form is "<#words * sign> sign words" */ { long *i; if (k == NULL) return(kerr(0,"kout")); else if (k->k_p == NULL) return(kerr(5,"kout")); else { printf("<%d> ",k->k_num); if (k->k_num < 0) printf("-"); else printf("+"); for (i = k->k_p + abs(k->k_num)-1; i >= k->k_p;i--) printf(" %08d",*i); /* show leading zeros */ printf(".\n"); } } ksign(k) /* return sign of k */ INTK k; { if (k == NULL) return(kerr(0,"ksign")); else if (abs(k->k_num) == 1 && *k->k_p == 0) return(0); else return(isign(k->k_num)); } knum(k) /* return number of words in k */ INTK k; { if (k == NULL) return(kerr(0,"knum")); else return(abs(k->k_num)); } INTK kneg(k) /* returns -k */ INTK k; /* done in place */ { if (k == NULL) { kerr(0,"kneg"); return(NULL); } else { k->k_num = -k->k_num; return(k);} } INTK ktok(k) /* returns a copy of k */ INTK k; { INTK tk,kalloc(); long *j,*l; if (k == NULL) return(NULL); else { tk = kalloc(knum(k)); tk->k_num = k->k_num; for (j = tk->k_p, l = k->k_p ; j < tk->k_p + knum(k) ;j++,l++) *j = *l; } return(tk); } INTK kalloc(n) /* allocate space for an n word k integer */ long n; { INTK k; char *calloc(); k = (struct k_long *) calloc(1,sizeof(struct k_long)); if (k == NULL) { kerr(1,"kalloc 1"); return(NULL);} k->k_num = n; k->k_p = (long *) calloc(n,sizeof(long)); if (k->k_p == NULL) { kerr(1,"kalloc 2"); return(NULL);} return(k); } kdeal(k) /* deallocate storage for k */ INTK k; { if (k == NULL) kerr(0,"kdeal"); else { cfree(k->k_p); cfree(k); } } INTK kint(i) /* set k integer to value i */ long i; { INTK tk,kalloc(); tk = kalloc(1); if (tk == NULL) { kerr(1,"kint"); return(NULL);} tk->k_num = isign(i); if (tk->k_num == 0) tk->k_num = 1; *tk->k_p = abs(i); return(tk); } INTK kstr(s) /* convert a string into a k integer */ char *s; /* no leading or trailing blanks permitted */ { char *p; long i,j,*l,len,is,n; INTK tk, kalloc(); if (s == NULL) { kerr(6,"kstr"); return(kint(0)); } len = strlen(s); p = s + len - 2; is = 1; if (*s == '-') { len--; is = -1; } else if (*s == '+') len--; n = len / KDIGPW; while (n * KDIGPW < len) n++; tk = kalloc(n); if (tk == NULL) { kerr(1,"kstr"); return(NULL); } tk->k_num = is * n; l = tk->k_p; i = 0; j = 1; while(p >= s) { if (isdigit(*p)) { i = i + j * (*p - 48); j = j * 10; if (j == KIMAX) { j = 1; *(l++) = i; i = 0; }}; p--; } if (i != 0) *l = i; return(tk); } INTK ksub(k1,k2) /* returns k1-k2 */ INTK k1,k2; { INTK kneg(),kadd(),tk,ktok(); if (kzero(k2)) return(ktok(k1)); tk=kadd(k1,kneg(k2)); k2=kneg(k2); return(tk); } keq(k1,k2) /* logical test k1=k2 */ INTK k1,k2; { long *i1,*i2; if (k1 == NULL) return(kerr(0,"keq arg 1")); if (k2 == NULL) return(kerr(0,"keq arg 2")); if (knum(k1) != knum(k2)) return(FALSE); if (ksign(k1) != ksign(k2)) return(FALSE); for (i1 = k1->k_p, i2 = k2->k_p; i1 <= k1->k_p + knum(k1); i1++, i2++) if (*i1 != *i2) return(FALSE); return(TRUE); } klt(k1,k2) /* logical test k1 < k2 */ INTK k1,k2; { long *i1,*i2; if (k1 == NULL) return(kerr(0,"keq arg 1")); if (k2 == NULL) return(kerr(0,"keq arg 2")); if (knum(k1) < knum(k2)) return(TRUE); else if (knum(k1) > knum(k2)) return (FALSE); if (ksign(k1) < ksign(k2)) return(TRUE); if (keq(k1,k2)) return(FALSE); for (i1 = k1->k_p + knum(k1) - 1, i2 = k2->k_p + knum(k2) - 1; i1 >= k1->k_p; i1--, i2--) if (*i1 > *i2) return(FALSE); else if (*i1 < *i2) return (TRUE); return(TRUE); } kzero(k) /* logical test k=0 */ INTK k; { if (k == NULL) return(kerr(0,"kzero")); else if (ksign(k) == 0) return(TRUE); else return(FALSE); } INTK kadd(k1,k2) /* returns k=k1+k2 */ INTK k1,k2; { long l1,l2,n,icarry; long ik1,ik2,ik3,*p1,*p2,*p3,i; INTK tk,kalloc(),ktok(); if (k1 == NULL) if (k2 != NULL) return(ktok(k2)); else return(NULL); else if (k2 == NULL) return(ktok(k1)); if (kzero(k1)) return(ktok(k2)); if (kzero(k2)) return(ktok(k1)); l1 = knum(k1); l2 = knum(k2); p1 = k1->k_p; p2 = k2->k_p; n = (l1 > l2 ? l1 : l2) + 1; if (ksign(k1) == ksign(k2)) { /* add same signs */ tk = kalloc(n); if (tk == NULL){ kerr(1,"kadd"); return(NULL); } p3 = tk->k_p; icarry = 0; for (i = 1;i <= n;i++) { ik1 = 0; ik2 = 0; if (i <= l1) ik1 = *(p1++); if (i <= l2) ik2 = *(p2++); ik3 = ik1 + ik2 + icarry; icarry = ik3 / KIMAX; ik3 = ik3 - icarry * KIMAX; *(p3++) = ik3; } if (icarry != 0) kerr(5,"kadd"); while (*(--p3) == 0 && p3 > tk->k_p) n--; tk->k_num = n * ksign(k1); } else { /* different signs */ n--; tk = kalloc(n); if (tk == NULL){ kerr(1,"kadd"); return(NULL); } p3 = tk->k_p; icarry = 1; for (i = 1;i <= n;i++) { ik1 = 0; ik2 = KIMAX - 1; if (i <= l1) ik1 = *(p1++); if (i <= l2) ik2 = KIMAX - 1 - *(p2++); ik3 = ik1 + ik2 + icarry; icarry = ik3 / KIMAX; ik3 = ik3 - icarry * KIMAX; *(p3++) = ik3; } tk->k_num = n * ksign(k1); if (icarry == 0) { tk->k_num = - tk->k_num; icarry = 1; p3 = tk->k_p; for (i = 1; i <= n; i++) { ik1 = KIMAX - 1 - *p3 + icarry; icarry = ik1 / KIMAX; *(p3++) = ik1 % KIMAX; } } p3 = tk->k_p + n; while (*(--p3) == 0 && p3 > tk->k_p) n--; tk->k_num = n * isign(tk->k_num); } return(tk); } double krmul(k,r) /* return k*r */ INTK k; double r; { long *i; double sum,pow; sum = 0.0; pow = 1.0; if (k == NULL) { kerr(0,"krmul"); return(sum); } for (i = k->k_p; i < k->k_p + knum(k); i++) { sum = sum + *i * pow * r; pow = pow * KIMAX; } if (ksign(k) < 0) sum = -sum; return(sum); } INTK kmul(k1,k2) /*returns k1*k2 */ INTK k1,k2; { long n,*p1,*p2,*p3,i1,i2,i3; INTK tk,kalloc(),kint(); if (k1 == NULL) kerr(0,"kmul 1"); if (k2 == NULL) kerr(0,"kmul 2"); if (k1 == NULL || k2 == NULL) return(kint(0)); if (kzero(k1) || kzero(k2)) return(kint(0)); n = knum(k1) + knum(k2); tk = kalloc(n); if (tk == NULL) { kerr(2,"kmul"); return(NULL); } p3 = tk->k_p; for (p3 = tk->k_p,i3 = 1; i3 <= n; i3++,p3++) *p3 = 0; kkkmul(tk->k_p,k1->k_p,knum(k1),k2->k_p,knum(k2)); p3 = tk->k_p + n - 1; while(n > 1 && *(p3--) == 0) n--; tk->k_num = n * ksign(k1) * ksign(k2); return(tk); } kkkmul(r,i1,n1,i2,n2) /* multiply 2 arrays of integers */ long r[],i1[],n1,i2[],n2; { long i,j,t,ic,a,b,c,d,rr; for (i = 0; i < n1 ; i++) { ic = 0; for (j = 0; j <= n2; j++) { t = 0; if (j < n2) t = i2[j]; ic = iiimul(&rr,t,i1[i],ic); rr = r[i + j] + rr; ic = ic + rr / KIMAX; r[i + j] = rr % KIMAX; } } } iiimul(r,i1,i2,ic) /* multiplies two large int without overflow */ long *r,i1,i2,ic; { long a,b,c,d,t,rr; a = i1 / KISQR; b = i1 % KISQR; c = i2 / KISQR; d = i2 % KISQR; t = (a * d + c * b); rr = b * d + ic + (t % KISQR) * KISQR; ic = a * c + t / KISQR + rr / KIMAX; *r = rr % KIMAX; return(ic); } INTK kidiv(k,i) /* returns k/i */ INTK k; /* done in place */ long i; { INTK t,kneg(),kidiv1(); long n,*m,l,ll,ic; if (k == NULL) { kerr(0,"kidiv"); return(NULL); } if (i == 0) { kerr(4,"kidiv"); return(NULL); } if (kzero(k)) return(k); if (abs(i) == 1) if (i > 0) return(k); else return(kneg(k)); /* for larger i call specialized routine */ if (abs(i) > 10) return(kidiv1(k,i)); n = knum(k); l = ksign(k); ic = 0; for (m = k->k_p + n - 1; m >= k->k_p; m--) { ll = *m; *m = (ll + ic * KIMAX) / abs(i); ic = (ll + ic * KIMAX) % abs(i); } m = k->k_p + n - 1; while (*m == 0 && n > 1) n--; k->k_num = n * isign(i) * l; return(k); } iiidiv(r,ic,i1,i2) /* r = i1/i2, returns i1 % i2 for larger ints */ long *r,ic,i1,i2; { long a,b,c,d,e,qd,qr,t; qd = KIMAX / i2; qr = KIMAX % i2; a = i1 / i2 + qd * ic; b = i1 % i2; e = ic; while ((c = iiimul(&d,e,qr,0)) > 0 ) { b = b + d % i2; a = a + d / i2 + qd * c + b / i2; b = b % i2; e = c; } b = b + (qr * e ) % i2; a = a + (qr * e ) / i2 + b / i2; *r = a; b = b % i2; return(b); } INTK kidiv1(k,i) /* returns k/i for larger i */ INTK k; /* done in place */ long i; { INTK t,kneg(),ktok(),kalloc(); long n,*l,m,ic; if (k == NULL) { kerr(0,"kidiv1"); return(NULL); } if (i == 0) { kerr(4,"kidiv1"); return(NULL); } if (kzero(k)) return(kint(0)); if (abs(i) == 1) if (i > 0) return(ktok(k)); else return(kneg(ktok(k))); ic = 0; n = knum(k); m = ksign(k); for (l = k->k_p + n - 1; l >= k->k_p; l--) ic = iiidiv(l,ic,*l,abs(i)); l = k->k_p + n -1; while (*l == 0 && n > 1) n--; k->k_num = n * isign(i) * m; return(k); } INTK ktrunc(k,nw) /* divides k by KIMAX**nw */ INTK k; long nw; /* done in-place */ { long *i,*j,n; if (k == NULL) { kerr(0,"ktrunc"); return(k); } if (nw <= 1) return(k); if ((n = knum(k)) < nw) { *k->k_p = 0; k->k_num = isign(k->k_num); } for(i = k->k_p + nw,j = k->k_p; i <= k->k_p + n - 1;i++, j++) *j = *i; k->k_num = (n - nw) * isign(k->k_num); return(k); } kimod(k,i) /* returns k mod i */ INTK k; long i; { long *l,ic; if (k == NULL) return(kerr(0,"kimod")); if (i == 0) return(kerr(4,"kimod")); if (kzero(k)) return(0); if (abs(i) == 1) return(ksign(k)*isign(i)); /* for larger i call specialized routine */ if (abs(i) > 10) return(kimod1(k,i)); ic = 0; for (l = k->k_p + knum(k) - 1; l >= k->k_p; l--) ic = (*l + ic * KIMAX) % abs(i); return(ic * ksign(k) * isign(i)); } kimod1(k,i) /* returns k mod i for larger i */ INTK k; long i; { long *l,ic,r; if (k == NULL) return(kerr(0,"kimod1")); if (i == 0) return(kerr(4,"kimod1")); if (kzero(k)) return(0); if (abs(i) == 1) return(ksign(k)*isign(i)); ic = 0; for (l = k->k_p + knum(k) - 1; l >= k->k_p; l--) ic = iiidiv(&r,ic,*l,abs(i)); return(ic * ksign(k) * isign(i)); } INTK kimul(k,i) /* returns k*i */ INTK k; long i; { INTK tk,kalloc(),kint(),ktok(),kneg(); long n,*l,*m,icarry,ik; if (k == NULL) { kerr(0,"kimul"); return(NULL); } if (i == 0 || ksign(k) == 0) return(kint(0)); if (abs(i) == 1) { tk = ktok(k); if (ksign(k) * isign(i) < 0) tk = kneg(tk); return(tk); } n = knum(k) + 1; tk = kalloc(n); if (tk == NULL) { kerr(2,"kimul"); return(NULL); } icarry = 0; for (l = k->k_p, m = tk->k_p; l <= k->k_p + n - 1; l++,m++) { ik = 0; if (l < k->k_p + knum(k)) ik = *l; icarry = iiimul(m,ik,abs(i),icarry); } l = tk->k_p + n - 1; while (*(l--) == 0 && n > 1) n--; tk->k_num = (abs(tk->k_num) - 1) * isign(tk->k_num); tk->k_num = n * ksign(k) * isign(i); return(tk); } /* sample output from program */ /* p0 4 0: <-1> - 00035478. 1: <1> + 00006840. 2: <1> + 00038600. 3: <-1> - 00036000. 4: <1> + 00010000. p1 4 0: <-1> - 00029160. 1: <1> + 00048600. 2: <-1> - 00036000. 3: <1> + 00010000. q = <2> + 00000001 00000000. r = <-1> - 00000000. p2 4 0: <1> + 00042039. 1: <-1> - 00036000. 2: <1> + 00010000. q = <2> + 00000001 00000000. r = <-1> - 00000000. p3 4 0: <-1> - 00006840. 1: <1> + 00003439. q = <1> + 34390000. r = <-1> - 55404000. p4 4 0: <1> + 00000001. Modified Chebychev Polynomials at 0 of 8 w=0: c[0],0 : <4> + 00010000 00000000 00000000 00000000. c[1],0 : <4> + 00020000 00000000 00000000 00000000. c[2],0 : <4> + 00030000 00000000 00000000 00000000. c[3],0 : <4> + 00040000 00000000 00000000 00000000. c[4],0 : <4> + 00050000 00000000 00000000 00000000. P0(w) = <4> + 00020000 00000000 00000000 00000000. P1(w) = <4> + 00400000 00000000 00000000 00000000. P2(w) = <4> + 00390000 00000000 00000000 00000000. P3(w) = <4> + 00380000 00000000 00000000 00000000. P4(w) = <4> + 00010000 00000000 00000000 00000000. 0 w = 0 phi1 = 0 l = 0 phi = 0 Modified Chebychev Polynomials at 1 of 8 w=0.392699: c[0],1 : <4> + 00003826 83432365 08977172 84599840. c[1],2 : <4> + 00007071 06781186 54752440 08443621. c[2],3 : <4> + 00009238 79532511 28675612 81831894. c[3],4 : <4> + 00010000 00000000 00000000 00000000. c[4],3 : <4> + 00009238 79532511 28675612 81831894. P0(w) = <4> + 01603128 49915866 49334917 51292520. P1(w) = <-4> - 00533224 92506131 27833545 19537800. P2(w) = <-4> - 01294199 84406834 40345583 58742240. P3(w) = <-4> - 01858144 56876677 10220816 25292981. P4(w) = <4> + 00003826 83432365 08977172 84599840. 1 w = 0.392699 phi1 = -0.249278 l = 2 phi = 4.46311 Modified Chebychev Polynomials at 2 of 8 w=0.785398: c[0],2 : <4> + 00007071 06781186 54752440 08443621. c[1],4 : <4> + 00010000 00000000 00000000 00000000. c[2],2 : <4> + 00007071 06781186 54752440 08443621. c[3],0 : <1> + 00000000. c[4],6 : <-4> - 00007071 06781186 54752440 08443621. P0(w) = <4> + 19765195 58998926 12717099 24774762. P1(w) = <4> + 25249221 37884563 31004098 13655640. P2(w) = <4> + 07971297 86166746 62229553 97593219. P3(w) = <-4> - 13976103 83315985 06690177 54367640. P4(w) = <4> + 00007071 06781186 54752440 08443621. 2 w = 0.785398 phi1 = 1.06525 l = 2 phi = 4.20684 Modified Chebychev Polynomials at 3 of 8 w=1.1781: c[0],3 : <4> + 00009238 79532511 28675612 81831894. c[1],2 : <4> + 00007071 06781186 54752440 08443621. c[2],5 : <-4> - 00003826 83432365 08977172 84599840. c[3],8 : <-4> - 00010000 00000000 00000000 00000000. c[4],5 : <-4> - 00003826 83432365 08977172 84599840. P0(w) = <-5> - 00000001 05392024 84062809 37301706 29791692. P1(w) = <5> + 00000001 12016659 62780319 65940777 36191560. P2(w) = <4> + 95562932 20875375 34515769 62235866. P3(w) = <-4> - 38875957 81876664 47550226 92542341. P4(w) = <4> + 00009238 79532511 28675612 81831894. 3 w = 1.1781 phi1 = -1.09982 l = 3 phi = 3.61257 Modified Chebychev Polynomials at 4 of 8 w=1.5708: c[0],4 : <4> + 00010000 00000000 00000000 00000000. c[1],0 : <1> + 00000000. c[2],8 : <-4> - 00010000 00000000 00000000 00000000. c[3],0 : <1> + 00000000. c[4],4 : <4> + 00010000 00000000 00000000 00000000. P0(w) = <-5> - 00000006 40780000 00000000 00000000 00000000. P1(w) = <4> + 68400000 00000000 00000000 00000000. P2(w) = <5> + 00000003 20390000 00000000 00000000 00000000. P3(w) = <-4> - 68400000 00000000 00000000 00000000. P4(w) = <4> + 00010000 00000000 00000000 00000000. 4 w = 1.5708 phi1 = -1.5708 l = 3 phi = 1.5708 Modified Chebychev Polynomials at 5 of 8 w=1.9635: c[0],3 : <4> + 00009238 79532511 28675612 81831894. c[1],6 : <-4> - 00007071 06781186 54752440 08443621. c[2],5 : <-4> - 00003826 83432365 08977172 84599840. c[3],4 : <4> + 00010000 00000000 00000000 00000000. c[4],5 : <-4> - 00003826 83432365 08977172 84599840. P0(w) = <-5> - 00000009 22124232 50694779 50682061 38526972. P1(w) = <-5> - 00000003 75291131 68552099 71235429 83769640. P2(w) = <5> + 00000006 04679814 66306797 10201849 02947866. P3(w) = <-4> - 87510762 22877738 34833127 67767579. P4(w) = <4> + 00009238 79532511 28675612 81831894. 5 w = 1.9635 phi1 = 0.644787 l = 3 phi = 2.21558 Modified Chebychev Polynomials at 6 of 8 w=2.35619: c[0],2 : <4> + 00007071 06781186 54752440 08443621. c[1],8 : <-4> - 00010000 00000000 00000000 00000000. c[2],2 : <4> + 00007071 06781186 54752440 08443621. c[3],0 : <1> + 00000000. c[4],6 : <-4> - 00007071 06781186 54752440 08443621. P0(w) = <-5> - 00000001 17034804 41001073 87282900 75225238. P1(w) = <-5> - 00000009 46750778 62115436 68995901 86344360. P2(w) = <5> + 00000007 27971297 86166746 62229553 97593219. P3(w) = <-4> - 82756103 83315985 06690177 54367640. P4(w) = <4> + 00007071 06781186 54752440 08443621. 6 w = 2.35619 phi1 = 1.48361 l = 3 phi = 1.48361 Modified Chebychev Polynomials at 7 of 8 w=2.74889: c[0],1 : <4> + 00003826 83432365 08977172 84599840. c[1],6 : <-4> - 00007071 06781186 54752440 08443621. c[2],3 : <4> + 00009238 79532511 28675612 81831894. c[3],8 : <-4> - 00010000 00000000 00000000 00000000. c[4],3 : <4> + 00009238 79532511 28675612 81831894. P0(w) = <5> + 00000006 24870920 83283896 35954562 42557240. P1(w) = <-5> - 00000008 87841016 23838550 65009752 39499000. P2(w) = <5> + 00000005 07822682 61024587 35340495 81969760. P3(w) = <-4> - 50492948 97877750 97503717 00518219. P4(w) = <4> + 00003826 83432365 08977172 84599840. 7 w = 2.74889 phi1 = -0.827273 l = 4 phi = 0.743523 */