## slow-portable-floating-point-constants-g.pkg
# Compiled by:
# src/lib/compiler/core.sublib# slow_portable_floating_point_constants_g: generate Mythryl floating-point constants.
# slow_portable_floating_point_constants_g uses long multiplication to find the correct bit pattern for
# the real. This method is slow, but accurate, and works to any precision,
# which means that floats can be cross-compiled correctly.
#
# The function emitreal should take (Int * Bool Rw_Vector * Int) which represents
# a real value as (sign * fraction * exponent).
# The sign is 0 if the real is positive, 1 if negative.
# The fraction is a boolean rw_vector representing the bits; note that the most
# significant bit is in position 0.
# The exponent is the binary exponent of the normalized fraction.
# "Normalized" here means a number between 0 and 1.
#
# The algorithm is inefficient on forms like 10000000.0; forms like 1E7 (with
# no bogus zeros) are better. Also inefficient on forms like 0E23 or 1E-212.
api Primitive_Floating_Point { # Unreferenced outside this file.
#
significant: Int;
minexp: Int;
maxexp: Int;
transreal: ((Int, ((Int,Int) -> Int), Int)) -> String;
};
api Floating_Point_Constant { # Unreferenced outside this file.
#
exception BAD_FLOAT String;
realconst: String -> String;
};
generic package slow_portable_floating_point_constants_g (
p: Primitive_Floating_Point # Primitive_Floating_Point is from src/lib/compiler/src/fconst/slow-portable-floating-point-constants-g.pkg)
: (weak) Floating_Point_Constant # Floating_Point_Constant is from src/lib/compiler/src/fconst/slow-portable-floating-point-constants-g.pkg{
include package p;
exception BAD_FLOAT String;
# Use more than the required precision, then round at the end.
# This criterion works well enough for the 53 bits required by
# Vax G format and IEEE double format, but has not been tested
# with other values of significant.
fun log2 0 => 1;
log2 i => 1 + log2 (i / 2);
end;
precision = significant + log2 (maxexp-minexp) + 3;
# A float is a WHOLE "fraction" and an exponent base TWO.
Float = { frac: multiword_int::Int, exp: Int };
bigint = multiword_int::from_int;
plus = multiword_int::(+);
times = multiword_int::(*);
# Size of a bigint in bits:
# bigsize 0 = 1
# bigsize i = floor (log2 (i))+1
# This means that:
# bigsize i = floor (log2 (2*i+1))
bigone = multiword_int::from_int 1;
fun bigsize x
=
multiword_int::log2 (plus (plus (x, x), bigone));
fun getbit (b, w)
=
multiword_int::compare (multiword_int::bitwise_and (multiword_int::(>>>) (b, w), bigone), bigone) == EQUAL;
infix my plus times ;
# Take a bigint and return a function that will represent the
# fraction. The function is called with two integers (start, width), and returns
# an integer represented by the bits from start to start+width-1.
# The high (1/2) bit is in position 0. Assumes that
# the bigint is positive. This will work if the bigint is smaller than
# the rw_vector or vice versa; however, the number will be truncated, not
# rounded.
exception BITS;
fun makebits frac (start, width)
=
{ s = bigsize frac;
fun onebit b
=
getbit (frac, unt::from_int (s - 1-b));
fun b TRUE => 1;
b FALSE => 0; end;
fun f 0 => b (onebit start);
f i => b (onebit (start+i)) + 2 * f (i - 1); end;
if ( start < 0 or
width < 0
)
raise exception BITS;
else
f (width - 1);
fi;
};
# Round a float to n significant digits
stipulate
one = bigint 1;
herein
fun round (float as { frac=>f, exp=>e }, n)
=
{ shift = bigsize f + 1 - n;
if (shift <= 0)
float;
else
{ exp => e + shift,
frac => (getbit (f, unt::from_int (shift - 1)))
?? multiword_int::(>>>) (f, unt::from_int shift) plus one
:: multiword_int::(>>>) (f, unt::from_int shift)
};
fi;
};
end;
# maketenth: create the float of one tenth,
# to any number of significant
# digits, with no rounding on the last digit.
#
stipulate
zero = bigint 0;
one = bigint 1;
two = bigint 2;
herein
fun maketenth 1
=>
{ frac=>one, exp=> -4 };
maketenth n
=>
{ my { frac, exp } = maketenth (n - 1);
recursive my tenthbit
=
\\ 0 => zero;
1 => one;
2 => one;
3 => zero;
n => tenthbit (n % 4);
end;
f = (frac times two) plus tenthbit n;
e = exp - 1;
{ frac => f,
exp => e
};
};
end;
end;
# Float values ten and one tenth, to the correct precision.
ten = { frac => bigint 5, exp => 1 };
tenth = round (maketenth (precision+1), precision);
# Multiply two floats together to the correct precision:
fun mult { frac => f1, exp => e1 }
{ frac => f2, exp => e2 }
=
{ f = f1 times f2;
my e: Int = e1 + e2; # Shouldn't need the type constraint, our comp bug XXX BUGGO FIXME
round ( { frac => f, exp => e }, precision );
};
# Create a dynamic rw_vector of powers of ten
package dfa
=
expanding_rw_vector_g (
package {
include package rw_vector;
Float = { frac: multiword_int::Int, exp: Int };
Element = Void -> Float;
Rw_Vector = Rw_Vector( Element );
Vector = Vector( Element );
}
);
stipulate
include package rw_vector;
include package list;
include package dfa;
infix my 9 sub ;
exception UNKNOWN;
fun makelem e
=
(\\ () = e);
one = { frac=>bigint 1, exp=>0 };
herein
pos10 = dfa::rw_vector (0, \\ () = raise exception UNKNOWN); /* 10^2^n */ my _ =
dfa::set (pos10, 0, makelem ten);
neg10 = dfa::rw_vector (0, \\ () = raise exception UNKNOWN); /* 10^-2^n */ my _ =
dfa::set (neg10, 0, makelem tenth);
fun access (arr, n)
=
(dfa::get (arr, n)) ()
except
UNKNOWN
=
{ last = access (arr, n - 1);
new = mult last last;
dfa::set (arr, n, makelem new);
new;
};
fun pow10_2 0 => one;
pow10_2 n => if (n > 0 ) access (pos10, n - 1);
else access (neg10, -n - 1); fi;
end;
fun raisepower (f, 0)
=>
f;
raisepower (f, e)
=>
{ sign = if (e < 0) -1;
else 1;
fi;
fun power (f, p)
=
mult f (pow10_2 (sign*p));
fun raisep (f, 0u0, _) => f;
raisep (f, e, p)
=>
if (unt::bitwise_and (e, 0u1) == 0u1)
raisep ( power (f, p), unt::(>>) (e, 0u1), p + 1 );
else raisep ( f, unt::(>>) (e, 0u1), p + 1 );
fi;
end;
raisep (f, unt::from_int (abs e), 1);
};
end;
end;
# Take a string list of the form { digit*.[digit*] },
# and return a bigint and the exponent base 10.
# Requires that the list contain a decimal point and
# no trailing zeros (useless zeros after the decimal point).
#
stipulate
ten = bigint 10;
zero = bigint 0;
herein
fun reducefrac f
=
{ fun getexp NIL => 0;
getexp ('.' ! _) => 0;
getexp (_ ! tl) => getexp tl - 1;
end;
fun getwhole NIL => zero;
getwhole ('.' ! tl) => getwhole tl;
getwhole ('0' ! tl) => ten times getwhole tl;
getwhole (n ! tl)
=>
bigint (char::to_int n - char::to_int '0') plus (ten times getwhole tl);
end;
backwards = reverse f;
whole = getwhole backwards;
exp = getexp backwards;
(whole, exp);
};
end;
# Take a legal ML float string and return an (Int * bigint * Int)
# which is the sign, whole "fraction", and power of ten exponent
#
fun getparts s
=
{ Trailing = SIGNIFICANT | TRAILING;
# Separate the fraction from the exponent,
# adding a decimal point if there is none
# and eliminating trailing zeros:
#
fun separate (NIL, s) => (NIL, NIL, s);
separate (('E' | 'e') ! tl, SIGNIFICANT) => (['.'], tl, SIGNIFICANT);
separate (('E' | 'e') ! tl, TRAILING ) => (NIL, tl, TRAILING );
separate ('0' ! tl, s)
=>
{ my (r, e, s) = separate (tl, s);
case s
TRAILING => ( r, e, TRAILING );
SIGNIFICANT => ('0' ! r, e, SIGNIFICANT);
esac;
};
separate ('.' ! tl, _)
=>
{ my (r, e, _) = separate (tl, TRAILING);
('.' ! r, e, SIGNIFICANT);
};
separate (hd ! tl, s)
=>
{ my (r, e, _) = separate (tl, s);
(hd ! r, e, SIGNIFICANT);
};
end;
my (unsigned, sign)
=
case (explode s)
('-' ! more) => (more, 1);
other => (other, 0);
esac;
my (frac_s, exp_s, _)
=
separate (unsigned, SIGNIFICANT);
fun atoi strlist
=
{ numlist
=
map (\\ n = char::to_int n - char::to_int '0')
strlist;
list::fold_forward
(\\ (a: Int, b) = b*10 + a)
0
numlist;
};
exp10
=
case exp_s
NIL => 0;
'-' ! more => -(atoi more);
other => atoi other;
esac
except
OVERFLOW = raise exception BAD_FLOAT s;
my (frac, exp)
=
reducefrac frac_s;
(sign, frac, exp10 + exp);
};
fun realconst f
=
{ my (sign, frac10, exp10) = getparts f;
float
=
raisepower (
round (
{ frac => frac10,
exp => 0
},
precision
),
exp10
);
my (newf as { frac, exp } )
=
round (float, significant+1);
size = bigsize frac;
bits = makebits frac;
exp = exp + size;
transreal (
case size
0 => (0, bits, 0);
_ => if (exp < minexp or
exp > maxexp
)
raise exception BAD_FLOAT f;
else
(sign, bits, exp);
fi;
esac
);
};
}; # generic package slow_portable_floating_point_constants_g
## Copyright 1989 by AT&T Bell Laboratories
## Subsequent changes by Jeff Prothero Copyright (c) 2010-2015,
## released per terms of SMLNJ-COPYRIGHT.