// Alogrithm 659, Collected Algorithm from ACM // This is the C version of program // By Yaohang Li // Dept. of Comp. Sci. // Florida State Univ. // 7/2000 #include #include int s, qs, coef[20][20], nextn, testn, hisum; double rqs; int primes[]={1,2,3,5,5,7,7,11,11,11,11,13,13,17,17,17,17,19,19,23,23,23,23,29,29,29,29,29,29,31,31,37,37,37,37,37,37,41,41,41}; int infaur(int dimen, int atmost) { int i,j; // Check s s=dimen; if (s<1||s>40) return(-1); qs=primes[s-1]; testn=qs*qs*qs*qs; // Compute log(atmost+testn) in base qs using a ratio of natural logs to get // an upper bound on (the number of digits in the base qs representation of // atmost+testn) minus one. hisum=(int)(log((double)(atmost+testn))/log((double)(qs))); if (hisum>=20) return(-2); // now find binomial coefficients mod qs in a lower-triangular matrix "coef" // using recursion binom(i,j)=binom(i-1,j)+binom(i-1,j-1) and a=b+c implies // mod(a,d)=mod(mod(b,d)+mod(c,d),d) coef[0][0]=1; for (j=1;j<=hisum;j++) { coef[j][0]=1; coef[j][j]=1; } for (j=1;j<=hisum;j++) for (i=j+1;i<=hisum;i++) coef[i][j]=(coef[i-1][j]+coef[i-1][j-1])%qs; // calculating these coefficients mod qs avoids possible overflow problems // with raw binomial coefficients // now complete initialization as described in section 2. nextn has 4 digits // in base qs, so hisum equals 3. nextn=testn-1; hisum=3; rqs=1.0/((double)qs); return(0); } int gofaur(double *quasi) { int i, j, k, ytemp[20], ztemp, ktemp, ltemp, mtemp; double r; // find quasi[1] using faure (section 3.3) // nextn has a representation in base // qs of the form: sum over j from zero to hisum of ytemp(j)*(qs**j) // we now compute the ytemp(j)'s. ktemp=testn; ltemp=nextn; for (i=hisum;i>=0;i--) { ktemp=ktemp/qs; mtemp=ltemp%ktemp; ytemp[i]=(ltemp-mtemp)/ktemp; ltemp=mtemp; } // quasi[k] has the form sum over j // from zero to hisum of ytemp[j]*(qs**(-(j+1))) // ready to compute quasi[0] // using nested multiplication r=(double)(ytemp[hisum]); for (i=hisum-1;i>=0;i--) r=ytemp[i]+rqs*r; quasi[0]=r*rqs; // find the other s-1 components of quasi using "faure" (section 3.2 and 3.3) for (k=1;k