/* $Header$ */
/*
* Copyright © 1988-2004 Keith Packard and Bart Massey.
* All Rights Reserved. See the file COPYING in this directory
* for licensing information.
*/
#include <math.h>
#include "nickle.h"
Fpart *zero_fpart, *one_fpart;
#if 0
#define DebugV(s,v) FilePrintf (FileStdout, "%s %v\n", s, v)
#define DebugN(s,n) FilePrintf (FileStdout, "%s %n\n", s, n)
#define DebugFp(s,f) FilePrintf (FileStdout, "%s %s%n\n", s, \
(f)->sign == Negative ? "-" : "", \
(f)->mag)
#define DebugF(s,f) { \
DebugFp(s,(f)->mant); \
DebugFp(" e ", (f)->exp); \
}
#else
#define DebugV(s,v)
#define DebugN(s,n)
#define DebugFp(s,f)
#define DebugF(s,f)
#endif
static void
FpartMark (void *object)
{
Fpart *f = object;
MemReference (f->mag);
}
DataType FpartType = { FpartMark, 0, "FpartType" };
static Fpart *
NewFpart (Sign sign, Natural *mag)
{
ENTER ();
Fpart *ret;
if (NaturalZero (mag))
sign = Positive;
ret = ALLOCATE (&FpartType, sizeof (Fpart));
ret->sign = sign;
ret->mag = mag;
RETURN (ret);
}
static Fpart *
NewIntFpart (int i)
{
ENTER ();
Sign sign;
unsigned long mag;
if (i < 0)
{
sign = Negative;
mag = -i;
}
else
{
sign = Positive;
mag = i;
}
RETURN (NewFpart (sign, NewNatural (mag)));
}
static Fpart *
NewValueFpart (Value v)
{
ENTER ();
Fpart *ret;
switch (ValueTag(v)) {
case rep_int:
ret = NewIntFpart (ValueInt(v));
break;
case rep_integer:
ret = NewFpart (IntegerSign(v), IntegerMag(v));
break;
default:
ret = zero_fpart;
break;
}
RETURN (ret);
}
static Fpart *
FpartAdd (Fpart *a, Fpart *b, Bool negate)
{
ENTER ();
Fpart *ret;
switch (catagorize_signs(a->sign, negate ? SignNegate (b->sign):b->sign)) {
default:
case BothPositive:
ret = NewFpart (Positive, NaturalPlus (a->mag, b->mag));
break;
case FirstPositive:
if (NaturalLess (a->mag, b->mag))
ret = NewFpart (Negative, NaturalMinus (b->mag, a->mag));
else
ret = NewFpart (Positive, NaturalMinus (a->mag, b->mag));
break;
case SecondPositive:
if (NaturalLess (a->mag, b->mag))
ret = NewFpart (Positive, NaturalMinus (b->mag, a->mag));
else
ret = NewFpart (Negative, NaturalMinus (a->mag, b->mag));
break;
case BothNegative:
ret = NewFpart (Negative, NaturalPlus (a->mag, b->mag));
break;
}
RETURN (ret);
}
static Fpart *
FpartMult (Fpart *a, Fpart *b)
{
ENTER ();
Sign sign;
sign = Positive;
if (a->sign != b->sign)
sign = Negative;
RETURN (NewFpart (sign, NaturalTimes (a->mag, b->mag)));
}
static Fpart *
FpartDivide (Fpart *a, Fpart *b)
{
ENTER ();
Natural *rem;
Natural *quo;
Sign sign;
sign = Positive;
if (a->sign != b->sign)
sign = Negative;
quo = NaturalDivide (a->mag, b->mag, &rem);
RETURN (NewFpart (sign, quo));
}
static Fpart *
FpartRsl (Fpart *a, int shift)
{
ENTER ();
RETURN (NewFpart (a->sign, NaturalRsl (a->mag, shift)));
}
static Fpart *
FpartLsl (Fpart *a, int shift)
{
ENTER ();
RETURN (NewFpart (a->sign, NaturalLsl (a->mag, shift)));
}
static Bool
FpartLess (Fpart *a, Fpart *b)
{
switch (catagorize_signs(a->sign, b->sign)) {
default:
case BothPositive:
return NaturalLess (a->mag, b->mag);
case FirstPositive:
return False;
case SecondPositive:
return True;
case BothNegative:
return NaturalLess (b->mag, a->mag);
}
}
static Bool
FpartEqual (Fpart *a, Fpart *b)
{
switch (catagorize_signs(a->sign, b->sign)) {
default:
case BothPositive:
case BothNegative:
return NaturalEqual (a->mag, b->mag);
case FirstPositive:
return False;
case SecondPositive:
return False;
}
}
static Bool
FpartZero (Fpart *a)
{
return NaturalZero (a->mag);
}
unsigned
FpartLength (Fpart *a)
{
unsigned bits;
digit top;
if (a->mag->length == 0)
return 0;
bits = (a->mag->length - 1) * LBASE2;
top = NaturalDigits(a->mag)[a->mag->length - 1];
while (top)
{
bits++;
top >>= 1;
}
return bits;
}
static unsigned
FpartZeros (Fpart *a)
{
int i;
int zeros = 0;
digit top;
if (a->mag->length == 0)
return 0;
for (i = 0; i < a->mag->length - 1; i++)
{
if (NaturalDigits(a->mag)[i] != 0)
break;
zeros += LBASE2;
}
top = NaturalDigits(a->mag)[i];
while ((top & 1) == 0)
{
zeros++;
top >>= 1;
}
return zeros;
}
static Fpart *
FpartNegate (Fpart *a)
{
ENTER ();
RETURN (NewFpart (SignNegate (a->sign), a->mag));
}
int
FpartInit (void)
{
ENTER ();
zero_fpart = NewFpart (Positive, zero_natural);
MemAddRoot (zero_fpart);
one_fpart = NewFpart (Positive, one_natural);
MemAddRoot (one_fpart);
EXIT ();
return 1;
}
static Value
FloatAdd (Value av, Value bv, int expandOk, Bool negate)
{
ENTER ();
Value ret;
Float *a = &av->floats, *b = &bv->floats;
Fpart *dist;
Fpart *amant, *bmant, *mant;
Fpart *exp;
int d;
unsigned prec;
int alen, blen;
dist = FpartAdd (a->exp, b->exp, True);
ret = 0;
if (NaturalLess (dist->mag, max_int_natural))
{
d = NaturalToInt (dist->mag);
if (dist->sign == Negative)
d = -d;
amant = a->mant;
bmant = b->mant;
alen = FpartLength (amant);
blen = FpartLength (bmant);
prec = 0;
exp = 0;
if (d >= 0)
{
if (alen + d <= blen + a->prec)
{
amant = FpartLsl (amant, d);
exp = b->exp;
prec = b->prec;
if (a->prec + d < prec)
prec = a->prec + d;
}
}
else
{
d = -d;
if (blen + d <= alen + b->prec)
{
bmant = FpartLsl (bmant, d);
exp = a->exp;
prec = a->prec;
if (b->prec + d < prec)
prec = b->prec + d;
}
}
if (prec)
{
mant = FpartAdd (amant, bmant, negate);
ret = NewFloat (mant, exp, prec);
}
}
if (!ret)
{
if (dist->sign == Negative)
{
ret = bv;
if (negate)
ret = NewFloat (NewFpart (SignNegate (bv->floats.mant->sign),
bv->floats.mant->mag),
bv->floats.exp,
bv->floats.prec);
}
else
ret = av;
}
RETURN (ret);
}
static Value
FloatPlus (Value av, Value bv, int expandOk)
{
return FloatAdd (av, bv, expandOk, False);
}
static Value
FloatMinus (Value av, Value bv, int expandOk)
{
return FloatAdd (av, bv, expandOk, True);
}
static Value
FloatTimes (Value av, Value bv, int expandOk)
{
ENTER ();
Float *a = &av->floats, *b = &bv->floats;
Fpart *mant;
Fpart *exp;
unsigned prec;
mant = FpartMult (a->mant, b->mant);
exp = FpartAdd (a->exp, b->exp, False);
if (a->prec < b->prec)
prec = a->prec;
else
prec = b->prec;
RETURN (NewFloat (mant, exp, prec));
}
static Value
FloatDivide (Value av, Value bv, int expandOk)
{
ENTER ();
Float *a = &av->floats, *b = &bv->floats;
Fpart *mant;
Fpart *amant = a->mant, *bmant = b->mant;
Fpart *exp;
unsigned prec;
unsigned iprec, alen;
if (FpartZero (b->mant))
{
RaiseStandardException (exception_divide_by_zero,
"real divide by zero",
2, av, bv);
RETURN (Void);
}
DebugF ("Dividend ", a);
DebugF ("Divisor ", b);
if (a->prec < b->prec)
prec = a->prec;
else
prec = b->prec;
iprec = prec + FpartLength (bmant) + 1;
alen = FpartLength (amant);
exp = b->exp;
if (alen < iprec)
{
amant = FpartLsl (amant, iprec - alen);
exp = FpartAdd (NewIntFpart (iprec-alen), exp, False);
}
exp = FpartAdd (a->exp, exp, True);
DebugFp ("amant ", amant);
DebugFp ("bmant ", bmant);
mant = FpartDivide (amant, bmant);
DebugFp ("mant ", mant);
DebugFp ("exp ", exp);
RETURN (NewFloat (mant, exp, prec));
}
static Value
FloatMod (Value av, Value bv, int expandOk)
{
ENTER ();
Value q;
q = Floor (Divide (av, bv));
av = Minus (av, Times (q, bv));
RETURN (av);
}
static Value
FloatLess (Value av, Value bv, int expandOk)
{
ENTER ();
Value ret;
Float *a = &av->floats, *b = &bv->floats;
if (FpartEqual (a->mant, zero_fpart))
{
if (b->mant->sign == Positive &&
!FpartEqual (b->mant, zero_fpart))
ret = TrueVal;
else
ret = FalseVal;
}
else if (FpartEqual (b->mant, zero_fpart))
{
if (a->mant->sign == Negative)
ret = TrueVal;
else
ret = FalseVal;
}
else if (FpartEqual (a->exp, b->exp))
{
ret = FalseVal;
if (FpartLess (a->mant, b->mant))
ret = TrueVal;
}
else
{
av = FloatMinus (av, bv, expandOk);
ret = FalseVal;
if (av->floats.mant->sign == Negative)
ret = TrueVal;
}
RETURN (ret);
}
static Value
FloatEqual (Value av, Value bv, int expandOk)
{
ENTER ();
Value ret;
Float *a = &av->floats, *b = &bv->floats;
if (FpartEqual (a->exp, b->exp))
{
ret = FalseVal;
if (FpartEqual (a->mant, b->mant))
ret = TrueVal;
}
else
{
av = FloatMinus (av, bv, expandOk);
ret = FalseVal;
if (NaturalZero (av->floats.mant->mag))
ret = TrueVal;
}
RETURN (ret);
}
static Value
FloatNegate (Value av, int expandOk)
{
ENTER ();
Float *a = &av->floats;
RETURN (NewFloat (FpartNegate (a->mant), a->exp, a->prec));
}
static Value
FloatInteger (Value av)
{
ENTER ();
Float *a = &av->floats;
Natural *mag;
int dist;
/*
* Can only reduce floats that are integral
*
* Ensure that the precision of the number holds every bit
* This requires that the precision of the representation
* be no smaller than length of the numbers plus the
* number of implied zeros.
*
* precision >= length (mant) + exponent
* precision - length (mant) >= exponent
* !(precision - length (mant) < exponent)
*
* The canonical representation ensures that length <= prec
*/
if (a->exp->sign == Positive &&
!NaturalLess (NewNatural (a->prec - FpartLength (a->mant)),
a->exp->mag))
{
mag = a->mant->mag;
dist = NaturalToInt (a->exp->mag);
if (dist)
mag = NaturalLsl (mag, dist);
av = Reduce (NewInteger (a->mant->sign, mag));
}
else
{
RaiseStandardException (exception_invalid_unop_value,
"ambiguous conversion to int",
1,
av);
}
RETURN (av);
}
static Value
FloatFloor (Value av, int expandOk)
{
ENTER ();
Float *a = &av->floats;
Fpart *mant;
Fpart *exp;
int d;
if (a->exp->sign == Positive)
RETURN (FloatInteger (av));
if (NaturalLess (NewNatural (a->prec), a->exp->mag))
RETURN (Zero);
d = NaturalToInt (a->exp->mag);
mant = FpartRsl (a->mant, d);
if (d && a->mant->sign == Negative)
{
mant = FpartAdd (mant, one_fpart, True);
d--;
}
exp = zero_fpart;
RETURN (FloatInteger (NewFloat (mant, exp, a->prec - d)));
}
static Value
FloatCeil (Value av, int expandOk)
{
ENTER ();
Float *a = &av->floats;
Fpart *mant;
Fpart *exp;
int d;
if (a->exp->sign == Positive)
RETURN (FloatInteger (av));
if (NaturalLess (NewNatural (a->prec), a->exp->mag))
RETURN (Zero);
d = NaturalToInt (a->exp->mag);
mant = FpartRsl (a->mant, d);
if (d && a->mant->sign == Positive)
{
mant = FpartAdd (mant, one_fpart, False);
d--;
}
exp = zero_fpart;
RETURN (FloatInteger (NewFloat (mant, exp, a->prec - d)));
}
static Value
FloatPromote (Value av, Value bv)
{
ENTER ();
int prec;
if (!ValueIsFloat(av))
{
if (bv && ValueIsFloat(bv))
prec = bv->floats.prec;
else
prec = DEFAULT_FLOAT_PREC;
av = NewValueFloat (av, prec);
}
RETURN (av);
}
static Value
FloatReduce (Value av)
{
return av;
}
/*
* 1/2 <= value / 2^exp2 < 1
* 1/base <= value / base^expbase < 1
*
* 2^(exp2-1) <= value < 2^exp2
*
* assign value = 2^(exp2-1)
*
* then
*
* 1/base <= 2^(exp2-1) / base^expbase < 1
*
* 1 <= 2^(exp2-1) / (base^(expbase-1)) < base
*
* -log(base) <= (exp2-1) * log(2) - expbase * log(base) < 1
*
* ignoring the right inequality
*
* 0 <= (exp2 - 1) * log(2) - (expbase-1) * log(base)
* (expbase - 1) * log(base) <= (exp2 - 1) * log(2)
* (expbase - 1) <= (exp2 - 1) * log(2) / log(base);
* expbase <= (exp2 - 1) * log(2) / log(base) + 1;
* expbase = floor ((exp2 - 1) * log(2) / log(base) + 1);
*
* Depending on value, expbase may need an additional digit
*/
#if 0
static Bool
NaturalBitSet (Natural *n, int i)
{
int d = i / LBASE2;
int b = i & LBASE2;
return d < NaturalLength (n) && (NaturalDigits(n)[d] & 1 << b);
}
#endif
static Value
FloatExp (Value exp2, Value *ratio, int ibase, unsigned prec)
{
ENTER ();
double dscale;
Value scale;
Value r;
Value min, max, mean, nmean;
Value pow2;
Value base;
Value base_f;
Value two;
Value two_f;
Bool done;
DebugV ("exp2", exp2);
two = NewInt (2);
two_f = NewIntFloat (2, prec + 32);
base = NewInt (ibase);
base_f = NewIntFloat (ibase, prec + 32);
/*
* Compute expbase, this is a bit tricky as log is only
* available in floats
*/
dscale = log(2) / log(ibase) * MAX_NICKLE_INT;
scale = Divide (NewInt ((int) dscale),
NewInt (MAX_NICKLE_INT));
/*
* min = floor (((exp2 - 1) * scale) + 1);
*/
min = Floor (Plus (Times (Minus (exp2, One), scale), One));
if (Negativep (min))
max = Div (min, two);
else
max = Times (min, two);
/*
* pow2 = 2 ** (exp2-1)
*/
pow2 = Pow (two_f, Minus (exp2, One));
mean = 0;
done = False;
do
{
if (aborting)
{
EXIT ();
*ratio = Void;
return Void;
}
nmean = Div (Plus (min, max), two);
if (mean && True(Equal (nmean, mean)))
{
nmean = Plus (nmean, One);
done = True;
}
mean = nmean;
DebugV ("min ", min);
DebugV ("mean", mean);
DebugV ("max ", max);
/*
* r = 2 ** (exp2-1) / (base ** (mean - 1))
*/
r = Divide (pow2, Pow (base_f, Minus (mean, One)));
if (done)
break;
if (True (Less (One, r)))
min = mean;
else
max = mean;
} while (False (Equal (max, min)));
mean = Minus (mean, One);
/* r = Divide (pow2, Pow (base, Minus (mean, One)));*/
r = Divide (Pow (two_f, exp2), Pow (base_f, mean));
/* r = Divide (pow2, Pow (base, mean)); */
EXIT ();
REFERENCE (mean);
REFERENCE (r);
*ratio = r;
return mean;
}
static Bool
FloatPrint (Value f, Value fv, char format, int base, int width, int prec, int fill)
{
ENTER ();
Float *a = &fv->floats;
Value expbase;
Fpart *exp;
Natural *int_n;
Natural *frac_n;
Value ratio;
Value down;
Value fratio;
Value m;
Value int_part;
Value frac_part;
unsigned length;
int orig_prec = prec;
int mant_prec;
int frac_prec;
int dig_max;
int exp_width;
int int_width;
int frac_width;
int print_width;
Bool negative;
char *int_buffer;
char *int_string;
char *frac_buffer;
char *frac_string;
char *exp_string = 0;
if (base <= 0)
base = 10;
if (prec == DEFAULT_OUTPUT_PRECISION)
prec = 15;
mant_prec = a->prec * log(2) / log(base);
DebugFp ("mant", a->mant);
DebugFp ("exp ", a->exp);
length = FpartLength (a->mant);
expbase = FloatExp (Plus (NewInt (length),
NewInteger (a->exp->sign,
a->exp->mag)),
&ratio,
base,
a->prec);
if (aborting)
{
EXIT ();
return False;
}
DebugV ("expbase", expbase);
DebugF ("ratio ", &ratio->floats);
down = Pow (NewInt (2),
NewInt ((int) length));
DebugV ("down ", down);
fratio = Divide (ratio, down);
DebugF ("fratio ", &fratio->floats);
negative = a->mant->sign == Negative;
m = NewInteger (Positive, a->mant->mag);
/*
* Round the mantissa up by adding a bit at the extreme of the precision
*/
m = Plus (m, NewFloat (one_fpart,
NewIntFpart (length - a->prec), a->prec + 2));
m = Times (m, fratio);
if (True (Less (m, One)))
{
m = Times (m, NewInt (base));
expbase = Minus (expbase, One);
}
else if (False (Less (m, NewInt (base))))
{
m = Divide (m, NewInt (base));
expbase = Plus (expbase, One);
}
exp = NewValueFpart (expbase);
switch (format) {
case 'e':
case 'E':
case 'f':
break;
default:
dig_max = prec;
if ((exp->sign == Positive &&
!NaturalLess (exp->mag, NewNatural (dig_max))) ||
(exp->sign == Negative &&
NaturalLess (NewNatural (4), exp->mag)))
{
format = 'e';
}
else
{
format = 'f';
}
}
if (format == 'f')
{
m = Times (m, Pow (NewInt (base), expbase));
exp_width = 0;
if (prec == INFINITE_OUTPUT_PRECISION)
{
prec = mant_prec;
if (ValueIsInt(expbase))
{
if (ValueInt(expbase) < 0)
prec -= ValueInt(expbase);
else if (ValueInt(expbase) > prec)
prec = ValueInt(expbase);
}
}
}
else
{
exp_string = NaturalSprint (0, exp->mag, base, &exp_width);
if (aborting)
{
EXIT ();
return True;
}
exp_width++;
if (exp->sign == Negative)
exp_width++;
if (prec == INFINITE_OUTPUT_PRECISION)
prec = mant_prec;
}
int_part = Floor (m);
frac_part = Minus (m, int_part);
if (ValueIsInteger(int_part))
int_n = IntegerMag(int_part);
else
int_n = NewNatural (ValueInt(int_part));
int_width = NaturalEstimateLength (int_n, base);
if (negative)
int_width++;
int_buffer = malloc (int_width + 1);
int_string = NaturalSprint (int_buffer + int_width + 1,
int_n, base, &int_width);
if (aborting)
{
EXIT ();
return True;
}
frac_prec = mant_prec - int_width;
if (*int_string == '0')
frac_prec++;
if (negative)
{
*--int_string = '-';
int_width++;
}
if (width)
{
if (width > 0)
frac_width = width - int_width - exp_width;
else
frac_width = -width - int_width - exp_width;
if (prec > 0)
if (frac_width > prec + 1)
frac_width = prec + 1;
}
else
{
if (prec == INFINITE_OUTPUT_PRECISION)
frac_width = frac_prec + 1;
else
frac_width = prec + 1;
}
/*
* Limit fraction to available precision
*/
if (frac_width > frac_prec + 1)
frac_width = frac_prec + 1;
if (frac_width < 2)
frac_width = 0;
frac_buffer = 0;
frac_string = 0;
if (frac_width)
frac_part = Floor (Times (frac_part, Pow (NewInt (base),
NewInt (frac_width - 1))));
if (frac_width && (!Zerop (frac_part) || orig_prec > 0))
{
int frac_wrote;
if (ValueIsInteger(frac_part))
frac_n = IntegerMag(frac_part);
else
frac_n = NewNatural (ValueInt(frac_part));
frac_buffer = malloc (frac_width + 1);
frac_string = NaturalSprint (frac_buffer + frac_width + 1,
frac_n, base, &frac_wrote);
if (aborting)
{
EXIT ();
return True;
}
while (frac_wrote < frac_width - 1)
{
*--frac_string = '0';
frac_wrote++;
}
*--frac_string = '.';
if (orig_prec < 0)
while (frac_buffer[frac_width-1] == '0')
frac_buffer[--frac_width] = '\0';
}
else
frac_width = 0;
print_width = int_width + frac_width + exp_width;
while (width > print_width)
{
FileOutchar (f, fill);
width--;
}
FilePuts (f, int_string);
if (frac_string)
FilePuts (f, frac_string);
if (exp_width)
{
FilePuts (f, "e");
if (exp->sign == Negative)
FilePuts (f, "-");
FilePuts (f, exp_string);
}
while (-width > print_width)
{
FileOutchar (f, fill);
width++;
}
free (int_buffer);
if (frac_buffer)
free (frac_buffer);
EXIT ();
return True;
}
static HashValue
FloatHash (Value av)
{
Float *a = &av->floats;
return (NaturalHash(a->mant->mag) ^ a->mant->sign ^
NaturalHash(a->exp->mag) ^ a->exp->sign);
}
static void
FloatMark (void *object)
{
Float *f = object;
MemReference (f->mant);
MemReference (f->exp);
}
ValueRep FloatRep = {
{ FloatMark, 0, "FloatRep" }, /* base */
rep_float, /* tag */
{ /* binary */
FloatPlus,
FloatMinus,
FloatTimes,
FloatDivide,
NumericDiv,
FloatMod,
FloatLess,
FloatEqual,
0,
0,
},
{ /* unary */
FloatNegate,
FloatFloor,
FloatCeil,
},
FloatPromote,
FloatReduce,
FloatPrint,
0,
FloatHash,
};
Value
NewFloat (Fpart *mant, Fpart *exp, unsigned prec)
{
ENTER ();
unsigned bits, dist;
Value ret;
DebugFp ("New mant", mant);
DebugFp ("New exp ", exp);
/*
* Trim to specified precision
*/
bits = FpartLength (mant);
if (bits > prec)
{
dist = bits - prec;
exp = FpartAdd (exp, NewIntFpart (dist), False);
mant = FpartRsl (mant, dist);
}
/*
* Canonicalize by shifting to a 1 in the LSB
*/
dist = FpartZeros (mant);
if (dist)
{
exp = FpartAdd (exp, NewIntFpart (dist), False);
mant = FpartRsl (mant, dist);
}
bits = FpartLength (mant);
if (bits == 0)
exp = mant = zero_fpart;
DebugFp ("Can mant", mant);
DebugFp ("Can exp ", exp);
ret = ALLOCATE (&FloatRep.data, sizeof (Float));
ret->floats.mant = mant;
ret->floats.exp = exp;
ret->floats.prec = prec;
RETURN (ret);
}
Value
NewIntFloat (int i, unsigned prec)
{
ENTER ();
RETURN (NewFloat (NewIntFpart (i), zero_fpart, prec));
}
Value
NewIntegerFloat (Integer *i, unsigned prec)
{
ENTER ();
Fpart *mant;
mant = NewFpart (IntegerSign((Value) i), IntegerMag((Value) i));
RETURN (NewFloat (mant, zero_fpart, prec));
}
Value
NewNaturalFloat (Sign sign, Natural *n, unsigned prec)
{
ENTER ();
Fpart *mant;
mant = NewFpart (sign, n);
RETURN (NewFloat (mant, zero_fpart, prec));
}
Value
NewRationalFloat (Rational *r, unsigned prec)
{
ENTER ();
Value num, den;
num = NewNaturalFloat (r->sign, r->num, prec);
den = NewNaturalFloat (Positive, r->den, prec);
RETURN (FloatDivide (num, den, 1));
}
#define SCALE_BITS 52
#define SCALE 4503599627370496.0 /* 2 ** 52 */
Value
NewDoubleFloat (double d)
{
ENTER ();
int e;
double m;
Sign ms;
double_digit dd;
if (d == 0.0) RETURN (Zero);
e = ilogb (d);
m = significand (d);
ms = Positive;
if (m < 0)
{
ms = Negative;
m = -m;
}
e = e - SCALE_BITS;
dd = (double_digit) (m * SCALE + 0.5);
RETURN (NewFloat (NewFpart (ms, NewDoubleDigitNatural (dd)),
NewIntFpart (e), SCALE_BITS));
}
Value
NewValueFloat (Value av, unsigned prec)
{
ENTER ();
switch (ValueTag(av)) {
case rep_int:
av = NewIntFloat (ValueInt(av), prec);
break;
case rep_integer:
av = NewIntegerFloat (&av->integer, prec);
break;
case rep_rational:
av = NewRationalFloat (&av->rational, prec);
break;
case rep_float:
av = NewFloat (av->floats.mant, av->floats.exp, prec);
break;
default:
break;
}
RETURN (av);
}
double
DoublePart (Value av, char *error)
{
double mantissa;
int i;
int e;
digit *mt;
double div;
av = NewValueFloat (av, 64);
if (!ValueIsFloat (av))
{
RaiseStandardException (exception_invalid_argument, error,
2, NewInt (0), av);
return 0.0;
}
if (NaturalLess (av->floats.exp->mag, max_int_natural))
e = NaturalToInt (av->floats.exp->mag);
else
e = MAX_NICKLE_INT;
if (e > 1023)
{
RaiseStandardException (exception_invalid_argument, error,
2, NewInt (0), av);
return 0.0;
}
if (av->floats.exp->sign == Negative)
e = -e;
mantissa = 0.0;
i = av->floats.mant->mag->length;
e += DIGITBITS * i;
mt = NaturalDigits (av->floats.mant->mag) + i;
div = 1.0 / (double) BASE;
while (i--)
{
mantissa = mantissa + (double) *--mt * div;
div *= 1.0 / (double) BASE;
}
if (av->floats.mant->sign == Negative)
mantissa = -mantissa;
return mantissa * pow (2.0, (double) e);
}
syntax highlighted by Code2HTML, v. 0.9.1