239 lines
6.2 KiB
C
239 lines
6.2 KiB
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;
|
|
}
|