go/src/cmd/gc/mparith3.c

// 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;
}

/*
 * 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);
	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);
		goto out;
	}
	if(s < 0) {
		// b is larger, shift a right
		mpshiftfix(&a->val, s);
		a->exp -= s;
		mpaddfixfix(&a->val, &b->val);
		goto out;
	}
	mpaddfixfix(&a->val, &b->val);

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);
	a->exp = (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("mpdivfltflt divide 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);
	a->exp = (a->exp-c.exp) - Mpscale*(Mpprec-1) + 1;

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

double
mpgetflt(Mpflt *a)
{
	int s, i, e;
	uvlong v, vm;
	double f;

	if(a->val.ovf)
		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);
		a->exp -= 1;
	}

	// the magic numbers (64, 63, 53, 10, -1074) are
	// IEEE specific. this should be done machine
	// independently or in the 6g half of the compiler

	// pick up the mantissa and a rounding bit in a uvlong
	s = 53+1;
	v = 0;
	for(i=Mpnorm-1; s>=Mpscale; i--) {
		v = (v<<Mpscale) | a->val.a[i];
		s -= Mpscale;
	}
	vm = v;
	if(s > 0)
		vm = (vm<<s) | (a->val.a[i]>>(Mpscale-s));

	// continue with 64 more bits
	s += 64;
	for(; s>=Mpscale; i--) {
		v = (v<<Mpscale) | a->val.a[i];
		s -= Mpscale;
	}
	if(s > 0)
		v = (v<<s) | (a->val.a[i]>>(Mpscale-s));

	// gradual underflow
	e = Mpnorm*Mpscale + a->exp - 53;
	if(e < -1074) {
		s = -e - 1074;
		if(s > 54)
			s = 54;
		v |= vm & ((1ULL<<s) - 1);
		vm >>= s;
		e = -1074;
	}

//print("vm=%.16llux v=%.16llux\n", vm, v);
	// round toward even
	if(v != 0 || (vm&2ULL) != 0)
		vm = (vm>>1) + (vm&1ULL);
	else
		vm >>= 1;

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

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

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;
}