/* This file is in the public domain. */

#include <stdio.h>
#include <unistd.h>
#include <strings.h>
#include <sys/wait.h>

/*

if table[e] = { e2, mul }, then N * 10^e = (N * mul) * 2^e2
mul is represented in decimal, with 1<=mul<2; the decimal point is
assumed to be just after the leading 1, which is present explicitly.

NRD must be at least log10(NDB) less than ND; roundoff error in bitdigs[][]
	can cause the result to be in error by up to NDB in the LSD, so to
	ensure an error bound of 1/2 LSD in the resulting decimal expansions,
	at least log10(NDB) digits must be dropped in rounding.  This is
	not checked for.

NDB must be high enough that the least bit falls right off the end of ND.
	If this is not true, resulting inaccuracies will be on the order of
	2^-NDB instead of 10^-ND.  A warning is printed if this is not true.

NB must be large enough to compute
	2^MAXEXP without overflow
	2^MINEXP with at least NDB significant bits left
	If either is found to be untrue, the program will print a message
	and die partway through the run.

NRD must be no smaller than as11's MAXDDIG define.  If this is not true,
	as11 built with the output will attempt to access outside the
	generated string constants.  This obviously cannot be checked for
	by this program.

If MINEXP is > 0, or if MAXEXP is < 0, this program will run fine, but
	as11 built with the output will attempt to do illegal things
	with array addresses.  This obviously cannot be checked for by
	this program.

For VAX-11/PDP-11 floats, which range from approximately .587747e-38 to
	.170141e+39, a MINEXP of -38 and a MAXEXP of 39 are appropriate,
	since any number requiring a more extreme exponent is guaranteed
	to overflow/underflow.  (Numbers with exponents in range may still
	overflow/underflow; as11 must check for that when converting from
	decimal to binary.)

*/

#define MINEXP  (-38) /* most negative exponent to compute */
#define MAXEXP   (39) /* most positive exponent to compute */
#define NB     (1024) /* number of bits to carry computations to */
#define NDB      (97) /* number of bits to use when computing decimal values */
#define ND       (30) /* number of digits to work out decimal values to */
#define NRD      (26) /* number of digits to round to */

static char bitdigs[NDB][ND];

static void initbitdigs(void)
{
 int b;
 int d;
 int t;
 char digs[ND];

 for (d=0;d<ND;d++) digs[d] = 0;
 digs[0] = 1;
 for (b=0;b<NDB;b++)
  { bcopy(&digs[0],&bitdigs[b][0],ND);
    t = 0;
    for (d=0;d<ND;d++)
     { t = (10 * t) + digs[d];
       digs[d] = t / 2;
       t &= 1;
     }
  }
#ifdef DUMP
 for (b=0;b<NDB;b++)
  { printf("%3d ",b);
    for (d=0;d<ND;d++) putchar('0'+bitdigs[b][d]);
    printf("\n");
  }
 exit(0);
#endif
 for (d=0;d<ND-1;d++)
  { if (bitdigs[NDB-1][d])
     { fprintf(stderr,"Warning: NDB needs to be increased\n");
       break;
     }
  }
 if (d == ND-1)
  { switch (bitdigs[NDB-1][ND-1])
     { case 0:
	  fprintf(stderr,"Warning: NDB is wastefully high\n");
	  break;
       case 1:
	  break;
       default:
	  fprintf(stderr,"Warning: NDB needs to be increased\n");
	  break;
     }
  }
}

static void printdigits(const char *bp)
{
 int b;
 int d;
 int c;
 char digs[ND];

 c = 0;
 for (d=ND-1;d>=0;d--)
  { for (b=0;b<NDB;b++) if (bp[b]) c += bitdigs[b][d];
    digs[d] = c % 10;
    c /= 10;
  }
 if (c > 0) printf("<<%d>>",c);
 c = 5;
 d = NRD;
 while ((d >= 0) && (c > 0))
  { c += digs[d];
    digs[d] = c % 10;
    c /= 10;
    d --;
  }
 if (c > 0) printf("<<%d>>",c);
 for (d=0;d<NRD;d++) putchar('0'+digs[d]);
}

static void outfilter(void)
{
 int p[2];

 pipe(p);
 if (fork() == 0)
  { if (p[0] != 0)
     { dup2(p[0],0);
       close(p[0]);
     }
    close(p[1]);
    execlp("sort","sort","-n","+1",(char *)0);
    perror("sort");
    exit(1);
  }
 if (p[1] != 1)
  { dup2(p[1],1);
    close(p[1]);
  }
 close(p[0]);
}

int main(void);
int main(void)
{
 int i;
 int j;
 int e2;
 int n;
 int tz;
 char bits[NB];
 char dbits[NDB];

 initbitdigs();
 printf("#define EXPTABMAX (%d)\n",MAXEXP);
 printf("#define EXPTABMIN (%d)\n",MINEXP);
 printf("#define EXPTABOFF (-EXPTABMIN)\n");
 fflush(stdout);
 outfilter();
 for (i=0;i<NB;i++) bits[i] = 0;
 bits[0] = 1;
 for (i=0;i<39;i++)
  { for (e2=0;(e2<NB)&&!bits[e2];e2++) ;
    if (e2+NDB > NB)
     { fprintf(stderr,"Fatal: NB needs to be increased (-ve powers)\n");
       exit(1);
     }
    printf("{ %4d, \""/*}*/,-e2);
    printdigits(&bits[e2]);
    printf(/*{*/"\" }, /* %d */\n",-i);
    n = 0;
    for (j=0;j<NB;j++)
     { n = n + n + bits[j];
       if (n >= 10)
	{ bits[j] = 1;
	  n -= 10;
	}
       else
	{ bits[j] = 0;
	}
     }
  }
 for (i=0;i<NB;i++) bits[i] = 0;
 bits[0] = 1;
 e2 = 0;
 tz = 0;
 for (i=1;i<40;i++)
  { n = 0;
    for (j=tz;(j<=e2)||n;j++)
     { if (j >= NB)
	{ fprintf(stderr,"Fatal: NB needs to be increased (+ve powers)\n");
	  exit(1);
	}
       if (bits[j]) n += 10;
       bits[j] = n & 1;
       n >>= 1;
     }
    for (e2=NB-1;!bits[e2];e2--) ;
    for (tz=0;!bits[tz];tz++) ;
    printf("{ %4d, \""/*}*/,e2);
    for (j=0;j<NDB;j++) dbits[j] = (j>e2) ? 0 : bits[e2-j];
    printdigits(&dbits[0]);
    printf(/*{*/"\" }, /* %d */\n",i);
  }
 fflush(stdout);
 close(1);
 wait(0);
 exit(0);
}