/************************************************************
 **
 ** cmul.c : generates JAL code to multiply by a constant
 **          fixed point number
 **
 ** Copyright (c) 2008, Kyle A. York
 ** All rights reserved
 **
 ************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <errno.h>
#include <ctype.h>

#define CGI
#ifdef CGI
#include "form.h"
#include "cgi_fns.h"
#include "html_fns.h"
#include "misc.h"
#endif

#define BIT_TOP 0x80000000
#define UNKNOWN (-1UL)
#define UNUSED(x) (void) (x)

static unsigned detail;

/* interesting numbers:
 *   1.234
 *   3.578
 *   0.1
 */

/* The bits[] array below holds the binary representation of
 * the floating value. Initially it simply contains a '1'
 * for each bit set in the value. Later it will contain a '1'
 * where addition is required, and a '-1' where subtraction
 * is required. Finally it will contain shift distances. */
typedef struct cmul_ {
  signed char bits[64];  /* which bits are used */
  size_t      bit_first; /* first bit used      */
  size_t      bit_last;  /* last bit used       */
} cmul_t;

/* Dump the bit array in 32.32 format. This starts looking
 * really odd if the shift distance is greater than 9, but
 * such is life. */
static unsigned cmul_bits_dump(char sign, const cmul_t *cmul)
{
  size_t   ii;
  unsigned sub;

  fputc(sign, stdout);
  for (ii = 0, sub = 0; ii < 64; ii++) {
    if (32 == ii) {
      fputc('.', stdout);
    } else if (ii && !(ii & 3)) {
      fputc(' ', stdout);
    }
    if (('+' == sign) && (cmul->bits[ii] > 0)) {
      fputc('0' + cmul->bits[ii], stdout);
    } else if (cmul->bits[ii] < 0) {
      sub++;
      if ('-' == sign) {
        fputc('0' - cmul->bits[ii], stdout);
      } else {
        fputc('0', stdout);
      }
    } else {
      fputc('0', stdout);
    }
  }
  fputc('\n', stdout);
  return sub;
}

/* Dump the bit array. For readability, the positive shift values
 * are dumped on one line, and the negative on a diffent one. */
static void cmul_dump(const char *hdr, const cmul_t *cmul)
{
  if (detail) {
    if (hdr) {
      printf("%s\n", hdr);
    }
    if (cmul_bits_dump('+', cmul)) {
      cmul_bits_dump('-', cmul);
    }
    fputc('\n', stdout);
  }
}

/* Convert a double into 32.32 fixed. Note that this just sets
 * the used bits to '1' */
static void cmul_init(cmul_t *cmul, double val)
{
  unsigned long intr;
  unsigned long frac;
  size_t        ii;

  intr    = (unsigned long) val;
  frac    = (unsigned long) ((val - intr) * 4294967296.0);

  for (ii = 0; ii < 32; ii++) {
    cmul->bits[ii]      = (intr & BIT_TOP) != 0;
    cmul->bits[ii + 32] = (frac & BIT_TOP) != 0;
    intr <<= 1;
    frac <<= 1;
  }
  cmul->bit_last  = UNKNOWN;
  cmul->bit_first = UNKNOWN;
  cmul_dump("-- initial bits", cmul);
}

/* When a string of contiguous set bits is found, sometimes more
 * optimal code will be generated by doing an addition and a
 * subtraction instead of multiple additions.
 *   Let's say there's a string of '1's representing 2^n...2^m
 * (n > m). Instead of:
 *    2^n + 2^(n - 1) + ... + 2^m
 * we can do:
 *    2^(n + 1) - 2^m
 * Note this doesn't always produce superiour results!
 */
static void cmul_sets_optimize(cmul_t *cmul)
{
  size_t ii;

  for (ii = 0; ii < 64; ii++) {
    if (cmul->bits[ii]) {
      size_t ct;

      for (ct = 1; (ii + ct < 64) && cmul->bits[ii + ct]; ct++)
        ;
      if (ct >= 3) {
        if (ii) {
          cmul->bits[ii - 1] = 1;
        }
        ct--;
        cmul->bits[ii + ct] = -1;
        while (ct--) {
          cmul->bits[ii + ct] = 0;
        }
      }
    }
  }
  cmul_dump("-- split add/sub", cmul);
}

/* look for identities (defined in comments below) */
static void cmul_identities_optimize(cmul_t *cmul)
{
  size_t ii;

  for (ii = 0; ii < 63; ii++) {
    if ((cmul->bits[ii] > 0) && (cmul->bits[ii + 1] < 0)) {
      /* + x/2^n - x/2^(n+1) == x/2^(n+1)  */
      cmul->bits[ii]     = 0;
      cmul->bits[ii + 1] = 1;
    } else if ((cmul->bits[ii] < 0) && (cmul->bits[ii + 1] > 0)) {
      /* - x/2^n + x/2^(n+1) == -x/2^(n+1) */
      cmul->bits[ii]     = 0;
      cmul->bits[ii + 1] = -1;
    }
  }
  cmul_dump("-- fixed identities", cmul);
}

/* setup shift distances */
static void cmul_shift_distance_calc(cmul_t *cmul)
{
  size_t ii;
  size_t bit_last;

  /* setup the fractional shift distances */
  for (ii = 32, bit_last = 31; ii < 64; ii++) {
    if (cmul->bits[ii]) {
      cmul->bits[ii] = cmul->bits[ii] / abs(cmul->bits[ii]) * (ii - bit_last);
      bit_last = ii;
    }
  }
  /* setup the integer shift distances */
  /* nb: these values will be one *greater* than needed
   *     to account for a possible zero shift */
  for (ii = 32, bit_last = 31; ii > 0; ii--) {
    if (cmul->bits[ii - 1]) {
      cmul->bits[ii - 1] = cmul->bits[ii - 1] / abs(cmul->bits[ii - 1])
        * (bit_last - ii + 2);
      bit_last = ii - 1;
    }
  }
  cmul_dump("-- set shift distances", cmul);
}

/* setup number of bits used and last bit used
 * (calculated according to the required err) */
static void cmul_precision_calc(cmul_t *cmul, double val, double err_max)
{
  double err;
  double vbit;
  double calc;
  size_t ct;

  /* determine how many bits are required to reach the desired
   * margin of error */
  err_max = err_max / 100.0;
  err     = err_max + 1.0;
  calc    = 0.0;
  for (ct = 0, vbit = BIT_TOP; 
       (ct < 63) && (err > err_max); 
       ct++, vbit = vbit / 2.0) {
    if (cmul->bits[ct]) {
      if (UNKNOWN == cmul->bit_first) {
        cmul->bit_first = ct;
      }
      if (cmul->bits[ct] < 0) {
        calc -= vbit;
      } else {
        calc += vbit;
      }
    }
    err = fabs(calc - val) / val;
  }
  cmul->bit_last = ct - 1;
  printf("Requested(%f) attained(%f) error(%f%%)\n\n", 
    val, calc, err * 100.0);
}

/* if the bit *past* the last bit used is '1', propogate it back */
static void cmul_round(cmul_t *cmul, double val, double err_max)
{
  if ((cmul->bit_last < 64) && cmul->bits[cmul->bit_last + 1]) {
    size_t pos;

    pos = cmul->bit_last;
    while (pos && (cmul->bits[pos])) {
      cmul->bits[pos] = 0;
      pos--;
    }
    cmul->bits[pos] = 1;
    cmul_shift_distance_calc(cmul);
    cmul_precision_calc(cmul, val, err_max);
  }
}

/* Dump the basic equation. This is human readable but would
 * generate awful code! */
static void cmul_basic_equation_dump(const cmul_t *cmul, double val)
{
  size_t        ii;
  unsigned long bit;

  fputs("-- basic equation\n", stdout);
  if (val < 0.0) {
    fputs("- ( ", stdout);
  }
  for (ii = 0, bit = BIT_TOP; ii <= cmul->bit_last; ii++) {
    if (cmul->bits[ii] > 0) {
      printf("%s%s%lux", (ii == cmul->bit_first) ? "" : " + ",
        (ii >= 32) ? "1/" : "", bit);
    } else if (cmul->bits[ii] < 0) {
      printf(" - %s%lux", (ii >= 32) ? "1/" : "", bit);
    }
    if (ii < 31) {
      bit >>= 1;
    } else if (ii == 31) {
      bit = 2;
    } else {
      bit <<= 1;
    }
  }
  if (val < 0.0){
    fputs(" )", stdout);
  }
  fputs("\n\n", stdout);
}

/* Dump the factored equation. This can be used directly to
 * generate good code. */
static void cmul_factored_equation_dump(const cmul_t *cmul, double val)
{
  size_t ii;
  size_t paren_ct;
  size_t pos;
  int    is_first;

  /* first the integer part */
  fputs("-- factored equation\n", stdout);
  if (val < 0.0) {
    fputs("- ( ", stdout);
  }
  if (cmul->bit_first < 32) {
    is_first = 1;
    pos = 31;
    if (cmul->bit_last < pos) {
      pos = cmul->bit_last;
    }
    fputs("( ", stdout);
    for (ii = cmul->bit_first; ii <= pos; ii++) {
      if (cmul->bits[ii]) {
        if (!is_first) {
          fputc(' ', stdout);
        }
        fputc('(', stdout);
        is_first = 0;
      }
    }
    is_first = 1;
    pos = 31;
    if (cmul->bit_last < pos) {
      pos = cmul->bit_last;
    }
    for (ii = cmul->bit_first; ii <= pos; ii++) {
      if (cmul->bits[ii]) {
        if (cmul->bits[ii] > 0) {
          printf(" %sx ) ", (is_first) ? "" : "+ ");
          if (cmul->bits[ii] != 1) {
            printf("* %lu", 1ul << (cmul->bits[ii] - 1));
          }
        } else {
          printf(" - x ) * %lu", 1ul << (-cmul->bits[ii] - 1));
        }
        is_first = 0;
      }
    }
    fputs(" )", stdout);
  }
  /* now the fractional part */
  for (ii = cmul->bit_last, paren_ct = 0; ii >= 32; ii--) {
    if (cmul->bits[ii]) {
      paren_ct++;
    }
  }
  if (paren_ct && (cmul->bit_first < 32)) {
    fputs(" + ", stdout);
  }

  for (ii = 0; ii < paren_ct; ii++) {
    fputs("( ", stdout);
  }
  for (ii = cmul->bit_last; ii >= 32; ii--) {
    if (cmul->bits[ii] > 0) {
      printf("%sx ) / %lu", 
        (ii == cmul->bit_last) ? "" : " + ", 1ul << cmul->bits[ii]);
    } else if (cmul->bits[ii] < 0) {
      printf(" - x ) / %lu", 
        1ul << (-cmul->bits[ii]));
    }
  }
  if (val < 0.0) {
    fputs(" )", stdout);
  }
  fputs("\n\n", stdout);
}

/* Dump one representation of the code. This uses the following
 * operators, translate into individual operators in the JAL
 * language. */
static void cmul_code_dump(const cmul_t *cmul, double val)
{
  size_t ii;
  size_t pos;

  /* first the fractional part part */
  fputs("--jal code\n", stdout);
  if (cmul->bit_first != cmul->bit_last) {
    fputs("tmp = x\n", stdout);
  }
  if (cmul->bit_last >= 32) {
    for (ii = cmul->bit_last; ii >= 32; ii--) {
      if (cmul->bits[ii] > 0) {
        if (ii != cmul->bit_last) {
          fputs("x = x + tmp\n", stdout);
        }
        printf("x = x >> %u\n", cmul->bits[ii]);
      } else if (cmul->bits[ii] < 0) {
        if (ii == cmul->bit_last) {
          fputs("x = -x\n", stdout);
        } else {
          fputs("x = x - tmp\n", stdout);
        }
        printf("x = x >> %u\n", -cmul->bits[ii]);
      }
    }
    for (ii = ii + 1; ii; ii--) {
      if (cmul->bits[ii - 1]) {
        unsigned char shift;

        shift = abs(cmul->bits[ii - 1]) - 1;

        if (shift) {
          printf("tmp = tmp << %u\n", shift);
        }
        printf("x = x %c tmp\n", (cmul->bits[ii - 1] < 0) ? '-' : '+');
      }
    }
  } else {
    /* now the integer part */
    if (cmul->bit_first < 32) {
      pos = 31;
      if (cmul->bit_last < pos) {
        pos = cmul->bit_last;
      }
      for (ii = cmul->bit_first; ii <= pos; ii++) {
        if (cmul->bits[ii]) {
          if (cmul->bits[ii] > 0) {
            if (ii != cmul->bit_first) {
              fputs("x = x + tmp\n", stdout);
            }
            if (cmul->bits[ii] != 1) {
              printf("x = x << %u\n", cmul->bits[ii] - 1);
            }
          } else {
            fputs("x = x - tmp\n", stdout);
            if (cmul->bits[ii] != -1) {
              printf("x = x << %u\n", -cmul->bits[ii] - 1);
            }
          }
        }
      }
    }
  }
  if (val < 0.0) {
    fputs("x = -x\n", stdout);
  }
  fputs("\n\n", stdout);
}

#ifdef CGI
enum {
  FORM_ELEMENT_VALUE,
  FORM_ELEMENT_ERROR,
  FORM_ELEMENT_DETAIL,
  FORM_ELEMENT_ROUND
};

static void exec(void *arg, const char *opt)
{
  form_handle frm;
  const char *str;

  UNUSED(opt);

  frm = *(form_handle *) arg;
  if (0 == form_handle_str_get(frm, FORM_ELEMENT_VALUE, &str)) {
    double val;
    char  *eptr;

    val = strtod(str, &eptr);
    while (isspace((unsigned char) *eptr)) {
      eptr++;
    }
    if (*eptr || (val == 0.0)) {
      fputs("<span class=\"error\">Invalid value."
             " Value must not be equal to 0</span>\n", stdout);
    } else {
      double err;
      cmul_t cmul;

      err = 0.0;
      if (0 == form_handle_str_get(frm, FORM_ELEMENT_ERROR, &str)) {
        err = strtod(str, &eptr);
        if (*eptr || (err < 0.0) || (err >= 100.0)) {
          fputs("<span class=\"error\">Invalid maximum error value."
                 " This must be a % &gt;= 0</span>\n", stdout);
          err = 0.0;
        }
      }

      if (EEXIST == form_handle_checkbox_test(frm, FORM_ELEMENT_DETAIL, 0)) {
        detail = 1;
      }

      cmul_init(&cmul, val);
      cmul_sets_optimize(&cmul);
      cmul_identities_optimize(&cmul);
      cmul_shift_distance_calc(&cmul);
      cmul_precision_calc(&cmul, val, err);
      if (EEXIST == form_handle_checkbox_test(frm, FORM_ELEMENT_ROUND, 0)) {
        cmul_round(&cmul, val, err);
      }
      cmul_basic_equation_dump(&cmul, val);
      cmul_factored_equation_dump(&cmul, val);
      cmul_code_dump(&cmul, val);
    }
  }
}

void cgiMain()
{
  static const form_element_t elements[] = {
    { "value",  FORM_ELEMENT_TYPE_TEXT,     0 },
    { "error",  FORM_ELEMENT_TYPE_TEXT,     0 },
    { "detail", FORM_ELEMENT_TYPE_CHECKBOX, 0 },
    { "round",  FORM_ELEMENT_TYPE_CHECKBOX, 0 }
  };
  static html_form_key_value_t kv[] = {
    { HTML_FORM_KV_VALUE,    "form", 0, 0 },
    { HTML_FORM_KV_CALLBACK, "exec", 0, exec }
  };
  form_handle frm;

  form_handle_alloc(&frm, COUNT(elements), elements);
  cgi_form_read(frm);
  kv[0].value = &frm;
  kv[1].value = &frm;
  html_header_contenttype("text/html");
  html_form_show("cmul.frm", COUNT(kv), kv);
  form_handle_free(frm);
}

#else

int main(int argc, char **argv)
{
  const char *prog;
  char       *eptr;
  const char *estr;
  double      val;
  double      err_max;
  int         round;

  UNUSED(argc);

  estr    = 0;
  err_max = 0.0;
  round   = 0;
  prog    = *argv;
  argv++;

  if (*argv && (0 == strcmp(*argv, "-detail"))) {
    detail = 1;
    argv++;
  }
  if (*argv) {
    errno = 0;
    val = strtod(*argv, &eptr);
    if ((ERANGE == errno) || *eptr || (val == 0.0)) {
      estr = "Invalid value";
    } else {
      argv++;
    }
  } else {
    estr = "No value supplied"; /* show help */
  }
  if (!estr && *argv) {
    errno = 0;
    err_max = strtod(*argv, &eptr);
    if ((ERANGE == errno) || *eptr || (err_max < 0.0) || (err_max >= 100.0)) {
      estr = "Invalid maximum error";
    } else {
      argv++;
      if (*argv && (0 == strcmp(*argv, "-round"))) {
        round = 1;
        argv++;
      }
    }
  }
  if (!estr && *argv) {
    estr = "Too many or invalid parameters";
  }

  if (estr) {
    printf("%s: %s\n"
           "Format: %s [-detail] val [ err [ -round ] ]\n"
           "  Create the output necessary to perform a constant\n"
           "multiplication. Parameters:\n"
           "  -detail : display the internal state at each transform\n"
           "  val     : the constant value for multiplication\n"
           "  err     : allowable error, as a percent (eg 0.5 = 0.5%%)\n"
           "            if missing all bits are used\n"
           "  -round  : If the first bit after the last bit used is one,\n"
           "            it it propogated up through the result. This can\n"
           "            result in a smaller error.\n",
           prog, estr, prog);
  } else {
    cmul_t cmul;

    cmul_init(&cmul, fabs(val));
    cmul_sets_optimize(&cmul);
    cmul_identities_optimize(&cmul);
    cmul_shift_distance_calc(&cmul);
    cmul_precision_calc(&cmul, fabs(val), err_max);
    if (round) {
      cmul_round(&cmul, fabs(val), err_max);
    }
    cmul_basic_equation_dump(&cmul, val);
    cmul_factored_equation_dump(&cmul, val);
    cmul_code_dump(&cmul, val);
  }
  return 0;
}

#endif