ext_comp.c

Go to the documentation of this file.

/*
  (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands.
  See the copyright notice in the ACK home directory, in the file "Copyright".
*/

/* $Id: ext_comp.c,v 1.1.1.1 2005/04/21 14:56:05 beng Exp $ */

/* extended precision arithmetic for the strtod() and cvt() routines */

/* This may require some more work when long doubles get bigger than 8
   bytes. In this case, these routines may become obsolete. ???
*/

#include        "ext_fmt.h"
#include        <float.h>
#include        <errno.h>
#include        <ctype.h>

static int b64_add(struct mantissa *e1, struct mantissa *e2);
static b64_sft(struct mantissa *e1, int n);

static
mul_ext(struct EXTEND *e1, struct EXTEND *e2, struct EXTEND *e3)
{
        /*      Multiply the extended numbers e1 and e2, and put the
                result in e3.
        */
        register int    i,j;            /* loop control */
        unsigned short  mp[4];
        unsigned short  mc[4];
        unsigned short  result[8];      /* result */

        register unsigned short *pres;

        /* first save the sign (XOR)                    */
        e3->sign = e1->sign ^ e2->sign;

        /* compute new exponent */
        e3->exp = e1->exp + e2->exp + 1;

        /* check for overflow/underflow ??? */

        /* 128 bit multiply of mantissas        */

        /* assign unknown long formats          */
        /* to known unsigned word formats       */
        mp[0] = e1->m1 >> 16;
        mp[1] = (unsigned short) e1->m1;
        mp[2] = e1->m2 >> 16;
        mp[3] = (unsigned short) e1->m2;
        mc[0] = e2->m1 >> 16;
        mc[1] = (unsigned short) e2->m1;
        mc[2] = e2->m2 >> 16;
        mc[3] = (unsigned short) e2->m2;
        for (i = 8; i--;) {
                result[i] = 0;
        }
        /*
         *      fill registers with their components
         */
        for(i=4, pres = &result[4];i--;pres--) if (mp[i]) {
                unsigned short k = 0;
                unsigned long mpi = mp[i];
                for(j=4;j--;) {
                        unsigned long tmp = (unsigned long)pres[j] + k;
                        if (mc[j]) tmp += mpi * mc[j];
                        pres[j] = tmp;
                        k = tmp >> 16;
                }
                pres[-1] = k;
        }

        if (! (result[0] & 0x8000)) {
                e3->exp--;
                for (i = 0; i <= 3; i++) {
                        result[i] <<= 1;
                        if (result[i+1]&0x8000) result[i] |= 1;
                }
                result[4] <<= 1;
        }       
        /*
         *      combine the registers to a total
         */
        e3->m1 = ((unsigned long)(result[0]) << 16) + result[1];
        e3->m2 = ((unsigned long)(result[2]) << 16) + result[3];
        if (result[4] & 0x8000) {
                if (++e3->m2 == 0) {
                        if (++e3->m1 == 0) {
                                e3->m1 = 0x80000000;
                                e3->exp++;
                        }
                }
        }
}

static
add_ext(struct EXTEND *e1, struct EXTEND *e2, struct EXTEND *e3)
{
        /*      Add two extended numbers e1 and e2, and put the result
                in e3
        */
        struct EXTEND ce2;
        int diff;

        if ((e2->m1 | e2->m2) == 0L) {
                *e3 = *e1;
                return;
        }
        if ((e1->m1 | e1->m2) == 0L) {
                *e3 = *e2;
                return;
        }
        ce2 = *e2;
        *e3 = *e1;
        e1 = &ce2;

        /* adjust mantissas to equal power */
        diff = e3->exp - e1->exp;
        if (diff < 0) {
                diff = -diff;
                e3->exp += diff;
                b64_sft(&(e3->mantissa), diff);
        }
        else if (diff > 0) {
                e1->exp += diff;
                b64_sft(&(e1->mantissa), diff);
        }
        if (e1->sign != e3->sign) {
                /* e3 + e1 = e3 - (-e1) */
                if (e1->m1 > e3->m1 ||
                    (e1->m1 == e3->m1 && e1->m2 > e3->m2)) {
                        /*      abs(e1) > abs(e3) */
                        if (e3->m2 > e1->m2) {
                                e1->m1 -= 1;    /* carry in */
                        }
                        e1->m1 -= e3->m1;
                        e1->m2 -= e3->m2;
                        *e3 = *e1;
                }
                else {
                        if (e1->m2 > e3->m2)
                                e3->m1 -= 1;    /* carry in */
                        e3->m1 -= e1->m1;
                        e3->m2 -= e1->m2;
                }
        }
        else {
                if (b64_add(&e3->mantissa,&e1->mantissa)) {/* addition carry */
                        b64_sft(&e3->mantissa,1);/* shift mantissa one bit RIGHT */
                        e3->m1 |= 0x80000000L;  /* set max bit  */
                        e3->exp++;              /* increase the exponent */
                }
        }
        if ((e3->m2 | e3->m1) != 0L) {
                /* normalize */
                if (e3->m1 == 0L) {
                        e3->m1 = e3->m2; e3->m2 = 0L; e3->exp -= 32;
                }
                if (!(e3->m1 & 0x80000000)) {
                        unsigned long l = 0x40000000;
                        int cnt = -1;

                        while (! (l & e3->m1)) {
                                l >>= 1; cnt--;
                        }
                        e3->exp += cnt;
                        b64_sft(&(e3->mantissa), cnt);
                }
        }
}

static int
cmp_ext(struct EXTEND *e1, struct EXTEND *e2)
{
        struct EXTEND tmp; 
         
        e2->sign = ! e2->sign; 
        add_ext(e1, e2, &tmp);  
        e2->sign = ! e2->sign;
        if (tmp.m1 == 0 && tmp.m2 == 0) return 0; 
        if (tmp.sign) return -1;
        return 1;
}

static
b64_sft(struct mantissa *e1, int n)
{
        if (n > 0) {
                if (n > 63) {
                        e1->l_32 = 0;
                        e1->h_32 = 0;
                        return;
                }
                if (n >= 32) {
                        e1->l_32 = e1->h_32;
                        e1->h_32 = 0;
                        n -= 32;
                }
                if (n > 0) {
                        e1->l_32 >>= n;
                        if (e1->h_32 != 0) {
                                e1->l_32 |= (e1->h_32 << (32 - n));
                                e1->h_32 >>= n;
                        }
                }
                return;
        }
        n = -n;
        if (n > 0) {
                if (n > 63) {
                        e1->l_32 = 0;
                        e1->h_32 = 0;
                        return;
                }
                if (n >= 32) {
                        e1->h_32 = e1->l_32;
                        e1->l_32 = 0;
                        n -= 32;
                }
                if (n > 0) {
                        e1->h_32 <<= n;
                        if (e1->l_32 != 0) {
                                e1->h_32 |= (e1->l_32 >> (32 - n));
                                e1->l_32 <<= n;
                        }
                }
        }
}

static int
b64_add(struct mantissa *e1, struct mantissa *e2)
                /*
                 * pointers to 64 bit 'registers'
                 */
{
        register int    overflow;
        int             carry;

                        /* add higher pair of 32 bits */
        overflow = ((unsigned long) 0xFFFFFFFF - e1->h_32 < e2->h_32);
        e1->h_32 += e2->h_32;

                        /* add lower pair of 32 bits */
        carry = ((unsigned long) 0xFFFFFFFF - e1->l_32 < e2->l_32);
        e1->l_32 += e2->l_32;
        if ((carry) && (++e1->h_32 == 0))
                return(1);              /* had a 64 bit overflow */
        else
                return(overflow);       /* return status from higher add */
}

/* The following tables can be computed with the following bc(1)
   program:

obase=16
scale=0
define t(x){
        auto a, b, c
        a=2;b=1;c=2^32;n=1
        while(a<x) {
                b=a;n+=n;a*=a
        }
        n/=2
        a=b
        while(b<x) {
                a=b;b*=c;n+=32
        }
        n-=32
        b=a
        while(a<x) {
                b=a;a+=a;n+=1
        }
        n-=1
        x*=16^16
        b=x%a
        x/=a
        if(a<=(2*b)) x+=1
        obase=10
        n
        obase=16
        return(x)
}
for (i=1;i<28;i++) {
        t(10^i)
}
0
for (i=1;i<20;i++) {
        t(10^(28*i))
}
0
define r(x){
        auto a, b, c
        a=2;b=1;c=2^32;n=1
        while(a<x) {
                b=a;n+=n;a*=a
        }
        n/=2
        a=b
        while(b<x) {
                a=b;b*=c;n+=32
        }
        n-=32
        b=a
        while(a<x) {
                b=a;a+=a;n+=1
        }
        a=b
        a*=16^16
        b=a%x
        a/=x
        if(x<=(2*b)) a+=1
        obase=10
        -n
        obase=16
        return(a)
}
for (i=1;i<28;i++) {
        r(10^i)
}
0
for (i=1;i<20;i++) {
        r(10^(28*i))
}
0

*/
static struct EXTEND ten_powers[] = {   /* representation of 10 ** i */
        { 0,    0,      0x80000000,     0 },
        { 0,    3,      0xA0000000,     0 },
        { 0,    6,      0xC8000000,     0 },
        { 0,    9,      0xFA000000,     0 },
        { 0,    13,     0x9C400000,     0 },
        { 0,    16,     0xC3500000,     0 },
        { 0,    19,     0xF4240000,     0 },
        { 0,    23,     0x98968000,     0 },
        { 0,    26,     0xBEBC2000,     0 },
        { 0,    29,     0xEE6B2800,     0 },
        { 0,    33,     0x9502F900,     0 },
        { 0,    36,     0xBA43B740,     0 },
        { 0,    39,     0xE8D4A510,     0 },
        { 0,    43,     0x9184E72A,     0 },
        { 0,    46,     0xB5E620F4,     0x80000000 },
        { 0,    49,     0xE35FA931,     0xA0000000 },
        { 0,    53,     0x8E1BC9BF,     0x04000000 },
        { 0,    56,     0xB1A2BC2E,     0xC5000000 },
        { 0,    59,     0xDE0B6B3A,     0x76400000 },
        { 0,    63,     0x8AC72304,     0x89E80000 },
        { 0,    66,     0xAD78EBC5,     0xAC620000 },
        { 0,    69,     0xD8D726B7,     0x177A8000 },
        { 0,    73,     0x87867832,     0x6EAC9000 },
        { 0,    76,     0xA968163F,     0x0A57B400 },
        { 0,    79,     0xD3C21BCE,     0xCCEDA100 },
        { 0,    83,     0x84595161,     0x401484A0 },
        { 0,    86,     0xA56FA5B9,     0x9019A5C8 },
        { 0,    89,     0xCECB8F27,     0xF4200F3A }
};
static struct EXTEND big_ten_powers[] = {  /* representation of 10 ** (28*i) */
        { 0,    0,      0x80000000,     0 },
        { 0,    93,     0x813F3978,     0xF8940984 },
        { 0,    186,    0x82818F12,     0x81ED44A0 },
        { 0,    279,    0x83C7088E,     0x1AAB65DB },
        { 0,    372,    0x850FADC0,     0x9923329E },
        { 0,    465,    0x865B8692,     0x5B9BC5C2 },
        { 0,    558,    0x87AA9AFF,     0x79042287 },
        { 0,    651,    0x88FCF317,     0xF22241E2 },
        { 0,    744,    0x8A5296FF,     0xE33CC930 },
        { 0,    837,    0x8BAB8EEF,     0xB6409C1A },
        { 0,    930,    0x8D07E334,     0x55637EB3 },
        { 0,    1023,   0x8E679C2F,     0x5E44FF8F },
        { 0,    1116,   0x8FCAC257,     0x558EE4E6 },
        { 0,    1209,   0x91315E37,     0xDB165AA9 },
        { 0,    1302,   0x929B7871,     0xDE7F22B9 },
        { 0,    1395,   0x940919BB,     0xD4620B6D },
        { 0,    1488,   0x957A4AE1,     0xEBF7F3D4 },
        { 0,    1581,   0x96EF14C6,     0x454AA840 },
        { 0,    1674,   0x98678061,     0x27ECE4F5 },
        { 0,    1767,   0x99E396C1,     0x3A3ACFF2 }
};

static struct EXTEND r_ten_powers[] = { /* representation of 10 ** -i */
        { 0,    0,      0x80000000,     0 },
        { 0,    -4,     0xCCCCCCCC,     0xCCCCCCCD },
        { 0,    -7,     0xA3D70A3D,     0x70A3D70A },
        { 0,    -10,    0x83126E97,     0x8D4FDF3B },
        { 0,    -14,    0xD1B71758,     0xE219652C },
        { 0,    -17,    0xA7C5AC47,     0x1B478423 },
        { 0,    -20,    0x8637BD05,     0xAF6C69B6 },
        { 0,    -24,    0xD6BF94D5,     0xE57A42BC },
        { 0,    -27,    0xABCC7711,     0x8461CEFD },
        { 0,    -30,    0x89705F41,     0x36B4A597 },
        { 0,    -34,    0xDBE6FECE,     0xBDEDD5BF },
        { 0,    -37,    0xAFEBFF0B,     0xCB24AAFF },
        { 0,    -40,    0x8CBCCC09,     0x6F5088CC },
        { 0,    -44,    0xE12E1342,     0x4BB40E13 },
        { 0,    -47,    0xB424DC35,     0x095CD80F },
        { 0,    -50,    0x901D7CF7,     0x3AB0ACD9 },
        { 0,    -54,    0xE69594BE,     0xC44DE15B },
        { 0,    -57,    0xB877AA32,     0x36A4B449 },
        { 0,    -60,    0x9392EE8E,     0x921D5D07 },
        { 0,    -64,    0xEC1E4A7D,     0xB69561A5 },
        { 0,    -67,    0xBCE50864,     0x92111AEB },
        { 0,    -70,    0x971DA050,     0x74DA7BEF },
        { 0,    -74,    0xF1C90080,     0xBAF72CB1 },
        { 0,    -77,    0xC16D9A00,     0x95928A27 },
        { 0,    -80,    0x9ABE14CD,     0x44753B53 },
        { 0,    -84,    0xF79687AE,     0xD3EEC551 },
        { 0,    -87,    0xC6120625,     0x76589DDB },
        { 0,    -90,    0x9E74D1B7,     0x91E07E48 }
};

static struct EXTEND r_big_ten_powers[] = { /* representation of 10 ** -(28*i) */
        { 0,    0,      0x80000000,     0 },
        { 0,    -94,    0xFD87B5F2,     0x8300CA0E },
        { 0,    -187,   0xFB158592,     0xBE068D2F },
        { 0,    -280,   0xF8A95FCF,     0x88747D94 },
        { 0,    -373,   0xF64335BC,     0xF065D37D },
        { 0,    -466,   0xF3E2F893,     0xDEC3F126 },
        { 0,    -559,   0xF18899B1,     0xBC3F8CA2 },
        { 0,    -652,   0xEF340A98,     0x172AACE5 },
        { 0,    -745,   0xECE53CEC,     0x4A314EBE },
        { 0,    -838,   0xEA9C2277,     0x23EE8BCB },
        { 0,    -931,   0xE858AD24,     0x8F5C22CA },
        { 0,    -1024,  0xE61ACF03,     0x3D1A45DF },
        { 0,    -1117,  0xE3E27A44,     0x4D8D98B8 },
        { 0,    -1210,  0xE1AFA13A,     0xFBD14D6E },
        { 0,    -1303,  0xDF82365C,     0x497B5454 },
        { 0,    -1396,  0xDD5A2C3E,     0xAB3097CC },
        { 0,    -1489,  0xDB377599,     0xB6074245 },
        { 0,    -1582,  0xD91A0545,     0xCDB51186 },
        { 0,    -1675,  0xD701CE3B,     0xD387BF48 },
        { 0,    -1768,  0xD4EEC394,     0xD6258BF8 }
};

#define TP      (int)(sizeof(ten_powers)/sizeof(ten_powers[0]))
#define BTP     (int)(sizeof(big_ten_powers)/sizeof(big_ten_powers[0]))
#define MAX_EXP (TP * BTP - 1)

static
add_exponent(struct EXTEND *e, int exp)
{
        int neg = exp < 0;
        int divsz, modsz;
        struct EXTEND x;

        if (neg) exp = -exp;
        divsz = exp / TP;
        modsz = exp % TP;
        if (neg) {
                mul_ext(e, &r_ten_powers[modsz], &x);
                mul_ext(&x, &r_big_ten_powers[divsz], e);
        }
        else {
                mul_ext(e, &ten_powers[modsz], &x);
                mul_ext(&x, &big_ten_powers[divsz], e);
        }
}

_str_ext_cvt(const char *s, char **ss, struct EXTEND *e)
{
        /*      Like strtod, but for extended precision */
        register int    c;
        int             dotseen = 0;
        int             digitseen = 0;
        int             exp = 0;

        if (ss) *ss = (char *)s;
        while (isspace(*s)) s++;

        e->sign = 0;
        e->exp = 0;
        e->m1 = e->m2 = 0;

        c = *s;
        switch(c) {
        case '-':
                e->sign = 1;
        case '+':
                s++;
        }
        while (c = *s++, isdigit(c) || (c == '.' && ! dotseen++)) {
                if (c == '.') continue;
                digitseen = 1;
                if (e->m1 <= (unsigned long)(0xFFFFFFFF)/10) {
                        struct mantissa a1;

                        a1 = e->mantissa;
                        b64_sft(&(e->mantissa), -3);
                        b64_sft(&a1, -1);
                        b64_add(&(e->mantissa), &a1);
                        a1.h_32 = 0;
                        a1.l_32 = c - '0';
                        b64_add(&(e->mantissa), &a1);
                }
                else exp++;
                if (dotseen) exp--;
        }
        if (! digitseen) return;

        if (ss) *ss = (char *)s - 1;

        if (c == 'E' || c == 'e') {
                int     exp1 = 0;
                int     sign = 1;
                int     exp_overflow = 0;

                switch(*s) {
                case '-':
                        sign = -1;
                case '+':
                        s++;
                }
                if (c = *s, isdigit(c)) {
                        do {
                                int tmp;

                                exp1 = 10 * exp1 + (c - '0');
                                if ((tmp = sign * exp1 + exp) > MAX_EXP ||
                                     tmp < -MAX_EXP) {
                                        exp_overflow = 1;
                                }
                        } while (c = *++s, isdigit(c));
                        if (ss) *ss = (char *)s;
                }
                exp += sign * exp1;
                if (exp_overflow) {
                        exp = sign * MAX_EXP;
                        if (e->m1 != 0 || e->m2 != 0) errno = ERANGE;
                }
        }
        if (e->m1 == 0 && e->m2 == 0) return;
        e->exp = 63;
        while (! (e->m1 & 0x80000000)) {
                b64_sft(&(e->mantissa),-1);
                e->exp--;
        }
        add_exponent(e, exp);
}

#include        <math.h>

static
ten_mult(struct EXTEND *e)
{
        struct EXTEND e1 = *e;

        e1.exp++;
        e->exp += 3;
        add_ext(e, &e1, e);
}

#define NDIGITS 128
#define NSIGNIFICANT 19

char *
_ext_str_cvt(struct EXTEND *e, int ndigit, int *decpt, int *sign, int ecvtflag)
{
        /*      Like cvt(), but for extended precision */

        static char buf[NDIGITS+1];
        struct EXTEND m;
        register char *p = buf;
        register char *pe;
        int findex = 0;

        if (ndigit < 0) ndigit = 0;
        if (ndigit > NDIGITS) ndigit = NDIGITS;
        pe = &buf[ndigit];
        buf[0] = '\0';

        *sign = 0;
        if (e->sign) {
                *sign = 1;
                e->sign = 0;
        }

        *decpt = 0;
        if (e->m1 != 0) {
                register struct EXTEND *pp = &big_ten_powers[1];

                while(cmp_ext(e,pp) >= 0) {
                        pp++;
                        findex = pp - big_ten_powers;
                        if (findex >= BTP) break;
                }
                pp--;
                findex = pp - big_ten_powers;
                mul_ext(e,&r_big_ten_powers[findex],e);
                *decpt += findex * TP;
                pp = &ten_powers[1];
                while(pp < &ten_powers[TP] && cmp_ext(e, pp) >= 0) pp++;
                pp--;
                findex = pp - ten_powers;
                *decpt += findex;

                if (cmp_ext(e, &ten_powers[0]) < 0) {
                        pp = &r_big_ten_powers[1];
                        while(cmp_ext(e,pp) < 0) pp++;
                        pp--;
                        findex = pp - r_big_ten_powers;
                        mul_ext(e, &big_ten_powers[findex], e);
                        *decpt -= findex * TP;
                        /* here, value >= 10 ** -28 */
                        ten_mult(e);
                        (*decpt)--;
                        pp = &r_ten_powers[0];
                        while(cmp_ext(e, pp) < 0) pp++;
                        findex = pp - r_ten_powers;
                        mul_ext(e, &ten_powers[findex], e);
                        *decpt -= findex;
                        findex = 0;
                }
                (*decpt)++;     /* because now value in [1.0, 10.0) */
        }
        if (! ecvtflag) {
                /* for fcvt() we need ndigit digits behind the dot */
                pe += *decpt;
                if (pe > &buf[NDIGITS]) pe = &buf[NDIGITS];
        }
        m.exp = -62;
        m.sign = 0;
        m.m1 = 0xA0000000;
        m.m2 = 0;
        while (p <= pe) {
                struct EXTEND oneminm;

                if (p - pe > NSIGNIFICANT) {
                        findex = 0;
                        e->m1 = 0;
                }
                if (findex) {
                        struct EXTEND tc, oldtc;
                        int count = 0;

                        oldtc.exp = 0;
                        oldtc.sign = 0;
                        oldtc.m1 = 0;
                        oldtc.m2 = 0;
                        tc = ten_powers[findex];
                        while (cmp_ext(e, &tc) >= 0) {
                                oldtc = tc;
                                add_ext(&tc, &ten_powers[findex], &tc);
                                count++;
                        }
                        *p++ = count + '0';
                        oldtc.sign = 1;
                        add_ext(e, &oldtc, e);
                        findex--;
                        continue;
                }
                if (e->m1) {
                        m.sign = 1;
                        add_ext(&ten_powers[0], &m, &oneminm);
                        m.sign = 0;
                        if (e->exp >= 0) {
                                struct EXTEND x;

                                x.m2 = 0; x.exp = e->exp;
                                x.sign = 1;
                                x.m1 = e->m1>>(31-e->exp);
                                *p++ = (x.m1) + '0';
                                x.m1 = x.m1 << (31-e->exp);
                                add_ext(e, &x, e);
                        }
                        else *p++ = '0';
                        /* Check that remainder is still significant */
                        if (cmp_ext(&m, e) > 0 || cmp_ext(e, &oneminm) > 0) {
                                if (e->m1 && e->exp >= -1) *(p-1) += 1;
                                e->m1 = 0;
                                continue;
                        }
                        ten_mult(&m);
                        ten_mult(e);
                }
                else *p++ = '0';
        }
        if (pe >= buf) {
                p = pe;
                *p += 5;        /* round of at the end */
                while (*p > '9') {
                        *p = '0';
                        if (p > buf) ++*--p;
                        else {
                                *p = '1';
                                ++*decpt;
                                if (! ecvtflag) {
                                        /* maybe add another digit at the end,
                                           because the point was shifted right
                                        */
                                        if (pe > buf) *pe = '0';
                                        pe++;
                                }
                        }
                }
                *pe = '\0';
        }
        return buf;
}

_dbl_ext_cvt(double value, struct EXTEND *e)
{
        /*      Convert double to extended
        */
        int exponent;

        value = frexp(value, &exponent);
        e->sign = value < 0.0;
        if (e->sign) value = -value;
        e->exp = exponent - 1;
        value *= 4294967296.0;
        e->m1 = value;
        value -= e->m1;
        value *= 4294967296.0;
        e->m2 = value;
}

static struct EXTEND max_d;

double
_ext_dbl_cvt(struct EXTEND *e)
{
        /*      Convert extended to double
        */
        double f;
        int sign = e->sign;

        e->sign = 0;
        if (e->m1 == 0 && e->m2 == 0) {
                return 0.0;
        }
        if (max_d.exp == 0) {
                _dbl_ext_cvt(DBL_MAX, &max_d);
        }
        if (cmp_ext(&max_d, e) < 0) {
                f = HUGE_VAL;
                errno = ERANGE;
        }
        else    f = ldexp((double)e->m1*4294967296.0 + (double)e->m2, e->exp-63);
        if (sign) f = -f;
        if (f == 0.0 && (e->m1 != 0 || e->m2 != 0)) {
                errno = ERANGE;
        }
        return f;
}

---