src/cmd/gc/mparith3.c

/* [<][>][^][v][top][bottom][index][help] */
This source file includes following definitions.
sigfig
mpsetexp
mpnorm
mpaddfltflt
mpmulfltflt
mpdivfltflt
mpgetfltN
mpgetflt
mpgetflt32
mpmovecflt
mpnegflt
mptestflt
// Copyright 2009 The Go Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

#include        <u.h>
#include        <libc.h>
#include        "go.h"

/*
 * returns the leading non-zero
 * word of the number
 */
int
sigfig(Mpflt *a)
{
        int i;

        for(i=Mpprec-1; i>=0; i--)
                if(a->val.a[i] != 0)
                        break;
//print("sigfig %d %d\n", i-z+1, z);
        return i+1;
}

/*
 * sets the exponent.
 * a too large exponent is an error.
 * a too small exponent rounds the number to zero.
 */
void
mpsetexp(Mpflt *a, int exp) {
        if((short)exp != exp) {
                if(exp > 0) {
                        yyerror("float constant is too large");
                        a->exp = 0x7fff;
                }
                else {
                        mpmovecflt(a, 0);
                }
        }
        else {
                a->exp = exp;
        }
}

/*
 * shifts the leading non-zero
 * word of the number to Mpnorm
 */
void
mpnorm(Mpflt *a)
{
        int s, os;
        long x;

        os = sigfig(a);
        if(os == 0) {
                // zero
                a->exp = 0;
                a->val.neg = 0;
                return;
        }

        // this will normalize to the nearest word
        x = a->val.a[os-1];
        s = (Mpnorm-os) * Mpscale;

        // further normalize to the nearest bit
        for(;;) {
                x <<= 1;
                if(x & Mpbase)
                        break;
                s++;
                if(x == 0) {
                        // this error comes from trying to
                        // convert an Inf or something
                        // where the initial x=0x80000000
                        s = (Mpnorm-os) * Mpscale;
                        break;
                }
        }

        mpshiftfix(&a->val, s);
        mpsetexp(a, a->exp-s);
}

/// implements float arihmetic

void
mpaddfltflt(Mpflt *a, Mpflt *b)
{
        int sa, sb, s;
        Mpflt c;

        if(Mpdebug)
                print("\n%F + %F", a, b);

        sa = sigfig(a);
        if(sa == 0) {
                mpmovefltflt(a, b);
                goto out;
        }

        sb = sigfig(b);
        if(sb == 0)
                goto out;

        s = a->exp - b->exp;
        if(s > 0) {
                // a is larger, shift b right
                mpmovefltflt(&c, b);
                mpshiftfix(&c.val, -s);
                mpaddfixfix(&a->val, &c.val, 0);
                goto out;
        }
        if(s < 0) {
                // b is larger, shift a right
                mpshiftfix(&a->val, s);
                mpsetexp(a, a->exp-s);
                mpaddfixfix(&a->val, &b->val, 0);
                goto out;
        }
        mpaddfixfix(&a->val, &b->val, 0);

out:
        mpnorm(a);
        if(Mpdebug)
                print(" = %F\n\n", a);
}

void
mpmulfltflt(Mpflt *a, Mpflt *b)
{
        int sa, sb;

        if(Mpdebug)
                print("%F\n * %F\n", a, b);

        sa = sigfig(a);
        if(sa == 0) {
                // zero
                a->exp = 0;
                a->val.neg = 0;
                return;
        }

        sb = sigfig(b);
        if(sb == 0) {
                // zero
                mpmovefltflt(a, b);
                return;
        }

        mpmulfract(&a->val, &b->val);
        mpsetexp(a, (a->exp + b->exp) + Mpscale*Mpprec - Mpscale - 1);

        mpnorm(a);
        if(Mpdebug)
                print(" = %F\n\n", a);
}

void
mpdivfltflt(Mpflt *a, Mpflt *b)
{
        int sa, sb;
        Mpflt c;

        if(Mpdebug)
                print("%F\n / %F\n", a, b);

        sb = sigfig(b);
        if(sb == 0) {
                // zero and ovfl
                a->exp = 0;
                a->val.neg = 0;
                a->val.ovf = 1;
                yyerror("constant division by zero");
                return;
        }

        sa = sigfig(a);
        if(sa == 0) {
                // zero
                a->exp = 0;
                a->val.neg = 0;
                return;
        }

        // adjust b to top
        mpmovefltflt(&c, b);
        mpshiftfix(&c.val, Mpscale);

        // divide
        mpdivfract(&a->val, &c.val);
        mpsetexp(a, (a->exp-c.exp) - Mpscale*(Mpprec-1) + 1);

        mpnorm(a);
        if(Mpdebug)
                print(" = %F\n\n", a);
}

static double
mpgetfltN(Mpflt *a, int prec, int bias)
{
        int s, i, e, minexp;
        uvlong v;
        double f;

        if(a->val.ovf && nsavederrors+nerrors == 0)
                yyerror("mpgetflt ovf");

        s = sigfig(a);
        if(s == 0)
                return 0;

        if(s != Mpnorm) {
                yyerror("mpgetflt norm");
                mpnorm(a);
        }

        while((a->val.a[Mpnorm-1] & Mpsign) == 0) {
                mpshiftfix(&a->val, 1);
                mpsetexp(a, a->exp-1);  // can set 'a' to zero
                s = sigfig(a);
                if(s == 0)
                        return 0;
        }

        // pick up the mantissa, a rounding bit, and a tie-breaking bit in a uvlong
        s = prec+2;
        v = 0;
        for(i=Mpnorm-1; s>=Mpscale; i--) {
                v = (v<<Mpscale) | a->val.a[i];
                s -= Mpscale;
        }
        if(s > 0) {
                v = (v<<s) | (a->val.a[i]>>(Mpscale-s));
                if((a->val.a[i]&((1<<(Mpscale-s))-1)) != 0)
                        v |= 1;
                i--;
        }
        for(; i >= 0; i--) {
                if(a->val.a[i] != 0)
                        v |= 1;
        }

        // gradual underflow
        e = Mpnorm*Mpscale + a->exp - prec;
        minexp = bias+1-prec+1;
        if(e < minexp) {
                s = minexp - e;
                if(s > prec+1)
                        s = prec+1;
                if((v & ((1<<s)-1)) != 0)
                        v |= 1<<s;
                v >>= s;
                e = minexp;
        }
        
        // round to even
        v |= (v&4)>>2;
        v += v&1;
        v >>= 2;

        f = (double)(v);
        f = ldexp(f, e);

        if(a->val.neg)
                f = -f;

        return f;
}

double
mpgetflt(Mpflt *a)
{
        return mpgetfltN(a, 53, -1023);
}

double
mpgetflt32(Mpflt *a)
{
        return mpgetfltN(a, 24, -127);
}

void
mpmovecflt(Mpflt *a, double c)
{
        int i;
        double f;
        long l;

        if(Mpdebug)
                print("\nconst %g", c);
        mpmovecfix(&a->val, 0);
        a->exp = 0;
        if(c == 0)
                goto out;
        if(c < 0) {
                a->val.neg = 1;
                c = -c;
        }

        f = frexp(c, &i);
        a->exp = i;

        for(i=0; i<10; i++) {
                f = f*Mpbase;
                l = floor(f);
                f = f - l;
                a->exp -= Mpscale;
                a->val.a[0] = l;
                if(f == 0)
                        break;
                mpshiftfix(&a->val, Mpscale);
        }

out:
        mpnorm(a);
        if(Mpdebug)
                print(" = %F\n", a);
}

void
mpnegflt(Mpflt *a)
{
        a->val.neg ^= 1;
}

int
mptestflt(Mpflt *a)
{
        int s;

        if(Mpdebug)
                print("\n%F?", a);
        s = sigfig(a);
        if(s != 0) {
                s = +1;
                if(a->val.neg)
                        s = -1;
        }
        if(Mpdebug)
                print(" = %d\n", s);
        return s;
}
/* [<][>][^][v][top][bottom][index][help] */
root/src/cmd/gc/mparith3.c

DEFINITIONS
