/* makepow5.c (emx+gcc) -- Copyright (c) 1996 by Eberhard Mattes */

#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include <emx/bigint.h>

static const char *out_fname;
static int single_prec;
static int long_mul_threshold;
static int max_exponent;
static int short_factor;
static int long_factor;

struct power
{
  int n;
  int i;
};

static struct power *powers;
static _bi_word *words;
static int words_used;
static int words_alloc;


static void usage (void)
{
  puts ("makepow5 -- Copyright (c) 1996 by Eberhard Mattes\n");
  puts ("Usage:");
  puts ("  makepow5 [-o <output_file>] [-s] [-m <threshold>] <max_exponent>");
  exit (1);
}


static void out_of_memory (void)
{
  fprintf (stderr, "Out of memory\n");
  exit (2);
}


static void compute_all_powers (void)
{
  int i, j, bw;
  _bi_bigint b;

  powers = malloc ((max_exponent + 1) * sizeof (_bi_bigint));
  if (powers == NULL)
    out_of_memory ();

  words_used = 0; words_alloc = max_exponent;
  words = malloc (words_alloc * sizeof (_bi_word));
  if (words == NULL)
    out_of_memory ();

  bw = _BI_WORDS ((max_exponent + 1) * 3);
  b.v = malloc (bw * sizeof (_bi_word));
  if (b.v == NULL)
    out_of_memory ();

  if (_bi_set_w (&b, bw, 1) != 0)
    abort ();
  powers[0].n = b.n;
  powers[0].i = words_used;
  words[words_used++] = b.v[0];

  fprintf (stderr, "Computing powers of 5...\n");
  for (i = 1; i <= max_exponent; ++i)
    {
      if (_bi_mul_bw (&b, bw, &b, 5) != 0)
        abort ();
      if (words_used + b.n > words_alloc)
        {
          words_alloc *= 2;
          words = realloc (words, words_alloc * sizeof (_bi_word));
          if (words == NULL)
            out_of_memory ();
        }
      powers[i].n = b.n;
      powers[i].i = words_used;
      for (j = 0; j < b.n; ++j)
        words[words_used++] = b.v[j];
    }

  fprintf (stderr, "Table size: %d words\n", words_used);

  if (single_prec)
    {
      unsigned long long x = 5;
      short_factor = 1;
      while (x == (_bi_word)x)
        {
          x *= 5;
          ++short_factor;
        }
    }
  else
    short_factor = 1;
}


static void find_size (int factor, int *psize2, int *psize3)
{
  int i, size2, size3;

  size2 = size3 = 0;
  for (i = short_factor; i < factor; i += short_factor)
    size2 += powers[i].n + 1;
  for (i = factor; i <= max_exponent; i += factor)
    size3 += powers[i].n + 1;
  *psize2 = size2; *psize3 = size3;
}


static void optimize (void)
{
  int i, size2, size3, best_factor, best_size2, best_size3;

  fprintf (stderr, "Optimizing...\n");
  best_factor = short_factor * ((max_exponent + short_factor - 1)
                                / short_factor);
  find_size (best_factor, &best_size2, &best_size3);
  for (i = short_factor * ((long_mul_threshold + short_factor - 1)
                           / short_factor);
       i < max_exponent; i += short_factor)
    {
      find_size (i, &size2, &size3);
      if (size2 + size3 < best_size2 + best_size3)
        {
          best_factor = i; best_size2 = size2; best_size3 = size3;
        }
    }
  long_factor = best_factor;
  fprintf (stderr, "Table 1:     %4d entries, %4d words\n",
           short_factor - 1, short_factor - 1);
  fprintf (stderr, "Table 2 & 4: %4d entries, %4d words\n",
           (best_factor - short_factor) / short_factor, best_size2);
  fprintf (stderr, "Table 3 & 4: %4d entries, %4d words\n",
           (max_exponent + 1) / best_factor, best_size3);
  fprintf (stderr, "Total size:                %4d words\n",
           best_size2 + best_size3 + short_factor);
}


static void write_output (void)
{
  FILE *f;
  int i, j, index;
  _bi_word s;

  fprintf (stderr, "Writing output file...\n");
  f = fopen (out_fname, "w");
  if (f == NULL)
    {
      perror (out_fname);
      exit (2);
    }
  fprintf (f, "/* %s */\n\n", out_fname);
  fprintf (f, "/* This file has been generated with\n       makepow5 -o %s ",
           out_fname);
  if (single_prec != 0)
    fprintf (f, "-s ");
  if (long_mul_threshold != 0)
    fprintf (f, "-m%d ", long_mul_threshold);
  fprintf (f, "%d\n   Please do not edit this file! */\n\n", max_exponent);

  if (short_factor != 1)
    fprintf (f, "#define TABLE1_SIZE      %d\n", short_factor - 1);
  fprintf (f, "#define TABLE2_SIZE      %d\n",
           (long_factor - short_factor) / short_factor);
  if (long_factor <= max_exponent)
    fprintf (f, "#define TABLE3_SIZE      %d\n",
             (max_exponent + 1) / long_factor);

  fprintf (f, "\n");
  if (short_factor != 1)
    {
      fprintf (f, "static const _bi_word table1[] =\n{\n");
      s = 5;
      for (i = 1; i < short_factor; ++i)
        {
          fprintf (f, "  %lu, /* 5^%d */\n", s, i);
          s *= 5;
        }
      fprintf (f, "};\n\n");
    }

  fprintf (f, "static const int table2[] =\n{\n");

  index = 0;
  for (i = short_factor; i < long_factor; i += short_factor)
    {
      fprintf (f, "  %d, /* 5^%d */\n", index, i);
      index += powers[i].n;
    }
  fprintf (f, "  %d\n", index);
  fprintf (f, "};\n\n");

  if (long_factor <= max_exponent)
    {
      fprintf (f, "static const int table3[] =\n{\n");
      for (i = long_factor; i <= max_exponent; i += long_factor)
        {
          fprintf (f, "  %d, /* 5^%d */\n", index, i);
          index += powers[i].n;
        }
      fprintf (f, "  %d\n", index);
      fprintf (f, "};\n\n");
    }

  fprintf (f, "static const _bi_word table4[] =\n{\n");
  for (i = short_factor; i < long_factor; i += short_factor)
    {
      fprintf (f, "/* 5^%d */\n", i);
      for (j = 0; j < powers[i].n; ++j)
        fprintf (f, "  0x%.8lx,\n", words[powers[i].i + j]);
    }
  if (long_factor <= max_exponent)
    {
      for (i = long_factor; i <= max_exponent; i += long_factor)
        {
          fprintf (f, "/* 5^%d */\n", i);
          for (j = 0; j < powers[i].n; ++j)
            fprintf (f, "  0x%.8lx,\n", words[powers[i].i + j]);
        }
    }
  fprintf (f, "};\n");

  if (fclose (f) != 0)
    {
      perror (out_fname);
      remove (out_fname);
    }
}


int main (int argc, char *argv[])
{
  int c;
  char *end;

  while ((c = getopt (argc, argv, "m:o:s")) != -1)
    switch (c)
      {
      case 'm':
        long_mul_threshold = strtol (optarg, &end, 10);
        if (long_mul_threshold < 1 || *end != 0)
          usage ();
        break;

      case 'o':
        out_fname = optarg;
        break;

      case 's':
        single_prec = 1;
        break;

      default:
        usage ();
      }

  if (argc - optind != 1)
    usage ();

  max_exponent = strtol (argv[optind], &end, 10);
  if (max_exponent < 1 || *end != 0)
    usage ();

  compute_all_powers ();
  if (long_mul_threshold != 0)
    optimize ();
  else
    long_factor = max_exponent + 1;

  if (out_fname != NULL)
    write_output ();

  return 0;
}