tool/mkfptab.c

/*
** This program generates C code for the tables used to generate powers
** of 10 in the powerOfTen() subroutine in util.c.
**
** The objective of the powerOfTen() subroutine is to provide the most
** significant 96 bits of any power of 10 between -348 and +347.  Rather
** than generate a massive 8K table, three much smaller tables are constructed,
** which can then generate the requested power of 10 using a single
** 160-bit multiple.
**
** This program works by internally generating a table of powers of
** 10 accurate to 256 bits each.  It then that full-sized, high-accuracy
** table to construct the three smaller tables needed by powerOfTen().
**
** LIMITATION:
**
** This program uses the __uint128_t datatype, available in gcc/clang.
** It won't build using other compilers.
*/
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <assert.h>

typedef unsigned __int128 u128;          /* 128-bit unsigned integer */
typedef unsigned long long int u64;      /* 64-bit unsigned integer */

/* There is no native 256-bit unsigned integer type, so synthesize one
** using four 64-bit unsigned integers.  Must significant first. */
struct u256 {
  u64 a[4];  /* big-endian */
};
typedef struct u256 u256;

/*
** Return a u64 with the N-th bit set.
*/
#define U64_BIT(N)  (((u64)1)<<(N))

/* Multiple *pX by 10, in-place */
static void u256_times_10(u256 *pX){
  u64 carry = 0;
  int i;
  for(i=3; i>=0; i--){
    u128 y = (10 * (u128)pX->a[i]) + carry;
    pX->a[i] = (u64)y;
    carry = (u64)(y>>64);
  }
}

/* Multiple *pX by 2, in-place.  AKA, left-shift */
static void u256_times_2(u256 *pX){
  u64 carry = 0;
  int i;
  for(i=3; i>=0; i--){
    u64 y = pX->a[i];
    pX->a[i] = (y<<1) | carry;
    carry = y>>63;
  }
}

/* Divide *pX by 10, in-place */
static void u256_div_10(u256 *pX){
  u64 rem = 0;
  int i;
  for(i=0; i<4; i++){
    u128 acc = (((u128)rem)<<64) | pX->a[i];
    pX->a[i] = acc/10;
    rem = (u64)(acc%10);
  }
}

/* Divide *pX by 2, in-place,  AKA, right-shift */
static void u256_div_2(u256 *pX){
  u64 rem = 0;
  int i;
  for(i=0; i<4; i++){
    u64 y = pX->a[i];
    pX->a[i] = (y>>1) | (rem<<63);
    rem = y&1;
  }
}

/* Note: The main table is a little larger on the low end than the required
** range of -348..+347, since we need the -351 value for the reduced tables.
*/
#define SCALE_FIRST  (-351)                          /* Least power-of-10 */
#define SCALE_LAST   (+347)                          /* Largest power-of-10 */
#define SCALE_COUNT  (SCALE_LAST - SCALE_FIRST + 1)  /* Size of main table */
#define SCALE_ZERO   (351)

int main(int argc, char **argv){
  int i, j, iNext;
  int e;
  int bRound = 0;
  int bTruth = 0;
  const u64 top = ((u64)1)<<63;
  u256 v;
  u64 aHi[SCALE_COUNT];
  u64 aLo[SCALE_COUNT];
  int aE[SCALE_COUNT];

  for(i=1; i<argc; i++){
    const char *z = argv[i];
    if( z[0]=='-' && z[1]=='-' && z[2]!=0 ) z++;
    if( strcmp(z,"-round")==0 ){
      bRound = 1;
    }else
    if( strcmp(z,"-truth")==0 ){
      bTruth = 1;
    }else
    {
      fprintf(stderr, "unknown option: \"%s\"\n", argv[i]);
      exit(1);
    }
  }

  /* Generate the master 256-bit power-of-10 table */
  v.a[0] = top;
  v.a[1] = 0;
  v.a[2] = 0;
  v.a[3] = 0;
  for(i=SCALE_ZERO, e=63; i>=0; i--){
    aHi[i] = v.a[0];
    aLo[i] = v.a[1];
    aE[i] = e;
    u256_div_10(&v);
    while( v.a[0] < top ){
      e++;
      u256_times_2(&v);
    }
  }
  v.a[0] = 0;
  v.a[1] = top;
  v.a[2] = 0;
  v.a[3] = 0;
  for(i=SCALE_ZERO+1, e=63; i<SCALE_COUNT; i++){
    u256_times_10(&v);
    while( v.a[0]>0 ){
      e--;
      u256_div_2(&v);
    }
    aHi[i] = v.a[1];
    aLo[i] = v.a[2];
    aE[i] = e;
  }

  if( bTruth ){
    /* With the --truth flag, also output the aTruth[] table that
    ** contains 128 bits of every power-of-two in the range */
    printf("  /* Powers of ten, accurate to 128 bits each */\n");
    printf("  static const struct {u64 hi; u64 lo;} aTruth[] = {\n");
    for(i=0; i<SCALE_COUNT; i++){
      u64 x = aHi[i];
      u64 y = aLo[i];
      int e = aE[i];
      char *zOp;
      if( e>0 ){
        zOp = "<<";
      }else{
        e = -e;
        zOp = ">>";
      }
      printf("   {0x%016llx, 0x%016llx}, /* %2d: 1.0e%+d %s %d */\n",
           x, y, i, i+SCALE_FIRST, zOp, e);
    }
    printf("  };\n");
  }

  /* The aBase[] table contains powers of 10 between 0 and 26.  These
  ** all fit in a single 64-bit integer.
  */
  printf("  static const u64 aBase[] = {\n");
  for(i=SCALE_ZERO, j=0; i<SCALE_ZERO+27; i++, j++){
    u64 x = aHi[i];
    int e = aE[i];
    char *zOp;
    if( e>0 ){
      zOp = "<<";
    }else{
      e = -e;
      zOp = ">>";
    }
    printf("    UINT64_C(0x%016llx), /* %2d: 1.0e%+d %s %d */\n",
           x, j, i+SCALE_FIRST, zOp, e);
  }
  printf("  };\n");

  /* For powers of 10 outside the range [0..26], we have to multiple
  ** on of the aBase[] entries by a scaling factor to get the true
  ** power of ten.  The scaling factors are all approximates accurate
  ** to 96 bytes, represented by a 64-bit integer in aScale[] for the
  ** most significant bits and a 32-bit integer in aScaleLo[] for the
  ** next 32 bites.
  **
  ** The scale factors are at increments of 27.  Except, the entry for 0
  ** is replaced by the -1 value as a special case.
  */
  printf("  static const u64 aScale[] = {\n");
  for(i=j=0; i<SCALE_COUNT; i=iNext, j++){
    const char *zExtra = "";
    iNext = i+27;
    if( i==SCALE_ZERO ){ i--; zExtra = " (special case)"; }
    u64 x = aHi[i];
    int e = aE[i];
    char *zOp;
    if( e>0 ){
      zOp = "<<";
    }else{
      e = -e;
      zOp = ">>";
    }
    printf("    UINT64_C(0x%016llx), /* %2d: 1.0e%+d %s %d%s */\n",
           x, j, i+SCALE_FIRST, zOp, e, zExtra);
  }
  printf("  };\n");
  printf("  static const unsigned int aScaleLo[] = {\n");
  for(i=j=0; i<SCALE_COUNT; i=iNext, j++){
    const char *zExtra = "";
    iNext = i+27;
    if( i==SCALE_ZERO ){ i--; zExtra = " (special case)"; }
    u64 x = aLo[i];
    int e = aE[i];
    char *zOp;
    if( bRound && (x & U64_BIT(31))!=0 && i!=SCALE_ZERO-1 ) x += U64_BIT(32);
    if( e>0 ){
      zOp = "<<";
    }else{
      e = -e;
      zOp = ">>";
    }
    printf("    0x%08llx, /* %2d: 1.0e%+d %s %d%s */\n",
           x>>32, j, i+SCALE_FIRST, zOp, e, zExtra);
  }
  printf("  };\n");
  return 0;
}