// 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 <stdio.h>
#include <math.h>

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<s;k++)
    {
      quasi[k]=0.0;
      r=rqs;
      for (j=0;j<=hisum;j++)
        {
          ztemp=0;
          for (i=j;i<=hisum;i++)
            ztemp+=coef[i][j]*ytemp[i];
// no apparent alternative one-dimensional coefficient array except via
// subscript address computations and equivalencing

// new ytemp[j] is the sum over i from j to hisum
// of (old ytemp[i]*binom(i,j))mod qs
          ytemp[j]=ztemp%qs;
          quasi[k]+=ytemp[j]*r;
          r*=rqs;
        }
     }

// update nextn and, if needed, testn and hisum

  nextn++;
  if (nextn==testn)
    {
      testn*=qs;
      hisum++;
    }

  return(0);
}

void main()
{
  double quasi[40];
  int dimen, n;
  int i,j;

  dimen=2;
  n=50000;
  infaur(dimen,n);
  
  for (i=0;i<n;i++)
    {
      gofaur(quasi);
      for (j=0;j<dimen;j++)
        printf("%f ",quasi[j]);
      printf("\n");
    }
}