## eight-byte-float-guts.pkg
# Compiled by:
# src/lib/std/src/standard-core.sublib### "Science is what you know,
### philosophy is what you don't know."
###
### -- Bertrand Russell
stipulate
package lms = list_mergesort; # list_mergesort is from src/lib/src/list-mergesort.pkgherein
package eight_byte_float_guts
: (weak) Float # Float is from src/lib/std/src/float.api {
package i= inline_t::default_int; # inline_t is from src/lib/core/init/built-in.pkg package math= math64; # math64 is from src/lib/std/src/math64-intel32.pkg infix my 50 ==== !=;
Float = Float;
fun *+(a: Float, b, c) = a*b+c;
fun *-(a: Float, b, c) = a*b-c;
my (====) = inline_t::f64::(====);
my (!=) = inline_t::f64::(!=);
fun unordered (x: Float, y)
=
bool::not (x>y or x <= y);
fun ?=== (x, y)
=
(x ==== y) or unordered (x, y);
fun is_normal x
=
case (runtime::asm::logb x)
#
-1023 => FALSE; # 0.0 or subnormal
1024 => FALSE; # inf or nan
_ => TRUE;
esac;
w31_r = inline_t::f64::from_int1 o inline_t::i1::copy_tagged_unt;
rbase = w31_r core_multiword_int::base;
base_bits = inline_t::tu::copyt_tagged_int core_multiword_int::base_bits;
intbound = w31_r 0ux40000000; # not necessarily the same as rbase
negintbound = -intbound;
# The next three values are computed laboriously, partly to
# avoid problems with inaccurate String->float conversions
# in the compiler itself. XXX BUGGO FIXME
max_finite
=
g (0.0, y, 53)
where
fun f (x, i)
=
if (i==1023 ) x; else f (x*2.0, i + 1);fi;
y = f (1.0, 0);
fun g (z, y, 0) => z;
g (z, y, i) => g (z+y, y*0.5, i - 1);
end;
end;
min_normal_pos
=
f 1.0
where
fun f (x)
=
{ y = x * 0.5;
#
if (is_normal y) f y;
else x;
fi;
};
end;
stipulate
# The intel32 uses extended precision (80 bits) internally, therefore
# it is necessary to write out the result of r * 0.5 to get
# 64 bit precision.
mem = inline_t::poly_rw_vector::make_nonempty_rw_vector (1, min_normal_pos);
set = inline_t::poly_rw_vector::set;
get = inline_t::poly_rw_vector::get_with_boundscheck;
fun f ()
=
{ r = get (mem, 0);
#
y = r * 0.5;
#
set (mem, 0, y);
#
if (get (mem, 0) ==== 0.0) r;
else f ();
fi;
};
herein
min_pos = f();
end;
pos_inf = max_finite * max_finite;
neg_inf = -pos_inf;
fun is_finite x
=
neg_inf < x and x < pos_inf;
fun is_nan x
=
bool::not (x====x);
fun floor x
=
if (x < intbound and
x >= negintbound
)
runtime::asm::floor x;
else
if (is_nan x) raise exception exceptions_guts::DOMAIN; # exceptions_guts is from src/lib/std/src/exceptions-guts.pkg else raise exception exceptions_guts::OVERFLOW; # exceptions_guts is from src/lib/std/src/exceptions-guts.pkg fi;
fi;
fun ceil n = -1 - floor (-1.0 - n);
fun truncate n = if (n < 0.0 ) ceil n; else floor n;fi;
fun round x
=
# Ties round to the nearest even number:
#
{ fl = floor (x + 0.5);
cl = ceil (x - 0.5);
if (fl == cl)
fl;
else
if (tagged_unt_guts::bitwise_and (tagged_unt_guts::from_int fl, 0u1) == 0u1)
cl;
else
fl;
fi;
fi;
};
# This is the IEEE double-
# precision maxint:
#
max_int = 4503599627370496.0;
stipulate
# realround mode x returns x rounded
# to the nearest integer using the
# given rounding mode.
# May be applied to inf's and nan's.
#
fun floatround mode x
=
{ save_mode = ieee_float::get_rounding_mode ();
ieee_float::set_rounding_mode mode;
(if (x>=0.0 ) x+max_int-max_int; else x-max_int+max_int;fi)
then
ieee_float::set_rounding_mode save_mode;
};
herein
float_floor = floatround ieee_float::TO_NEGINF;
float_ceil = floatround ieee_float::TO_POSINF;
float_truncate = floatround ieee_float::TO_ZERO;
float_round = floatround ieee_float::TO_NEAREST;
end;
my abs: Float -> Float
= inline_t::f64::abs;
my from_int: Int -> Float
= inline_t::f64::from_tagged_int;
fun to_int ieee_float::TO_NEGINF => floor;
to_int ieee_float::TO_POSINF => ceil;
to_int ieee_float::TO_ZERO => truncate;
to_int ieee_float::TO_NEAREST => round;
end;
fun to_eight_byte_float x = x;
fun from_eight_byte_float _ x = x;
fun sign x
=
if (x < 0.0) -1;
elif (x > 0.0) 1;
elif (is_nan x) raise exception DOMAIN;
else 0;
fi;
fun sign_bit x # Bug: negative zero not handled properly XXX BUGGO FIXME
=
runtime::asm::scalb (x, -(runtime::asm::logb x)) < 0.0;
fun same_sign (x, y)
=
sign_bit x == sign_bit y;
fun copy_sign (x, y)
= # May not work if x is Nan.
if (same_sign (x, y)) x;
else -x;
fi;
fun compare (x, y)
=
if (x < y ) exceptions_guts::LESS; # exceptions_guts is from src/lib/std/src/exceptions-guts.pkg elif (x > y ) exceptions_guts::GREATER;
elif (x ==== y) exceptions_guts::EQUAL;
else raise exception ieee_float::UNORDERED_EXCEPTION;
fi;
fun compare_real (x, y)
=
if (x<y ) ieee_float::LESS;
elif (x>y ) ieee_float::GREATER;
elif (x ==== y) ieee_float::EQUAL;
else ieee_float::UNORDERED;
fi;
# * This probably needs to be reorganized * XXX BUGGO FIXME
#
fun ilk x
= # Does not distinguish between quiet and signalling NaN
if (sign_bit x)
if (x>neg_inf) if (x ==== 0.0) ieee_float::ZERO;
elif (runtime::asm::logb x == -1023) ieee_float::SUBNORMAL;
else ieee_float::NORMAL;
fi;
elif (x====x) ieee_float::INF;
else ieee_float::NAN ieee_float::QUIET;
fi;
elif (x<pos_inf) if (x ==== 0.0) ieee_float::ZERO;
elif (runtime::asm::logb x == -1023) ieee_float::SUBNORMAL;
else ieee_float::NORMAL;
fi;
elif (x====x ) ieee_float::INF;
else ieee_float::NAN ieee_float::QUIET;
fi;
radix = 2;
precision = 53; # hidden bit gets counted, too
two_to_the_neg_1000
=
f (1000, 1.0)
where
fun f (i, x)
=
if (i == 0) x;
else f (i - 1, x*0.5); fi;
end;
# AARGH! Our version of logb gives a value that's one less than the
# rest of the world's logb functions.
# We should fix this systematically some time. XXX BUGGO FIXME
fun to_mantissa_exponent x
=
case (runtime::asm::logb x + 1)
-1023 => if (x====0.0 ) { mantissa => x, exponent => 0 };
else my { mantissa => m, exponent => e } = to_mantissa_exponent (x*1048576.0);
{ mantissa => m, exponent => e - 20 };
fi;
1024 => { mantissa => x, exponent => 0 };
i => { mantissa => runtime::asm::scalb (x, -i), exponent => i };
esac;
fun from_mantissa_exponent { mantissa => m, exponent => e: Int }
=
if (m >= 0.5 and m <= 1.0 or m <= -0.5 and m >= -1.0)
if (e > 1020)
if (e > 1050)
if (m>0.0) pos_inf;
else neg_inf;
fi;
else
f (e - 1020, runtime::asm::scalb (m, 1020))
where
fun f (i, x)
=
if (i==0) x;
else f (i - 1, x+x);
fi;
end;
fi;
elif (e < -1020)
if (e < -1200 ) 0.0;
else
fun f (i, x)
=
if (i==0) x;
else f (i - 1, x*0.5);
fi;
f (1020-e, runtime::asm::scalb (m, -1020));
fi;
else
runtime::asm::scalb (m, e); # This is the common case!
fi;
else
(to_mantissa_exponent m) -> { mantissa => m', exponent => e' };
from_mantissa_exponent { mantissa => m', exponent => e' + e };
fi;
# Some protection against insanity...
my _ =
if (base_bits < 18 ) # i.e., 3 * base_bits < 53
raise exception DIE
"big digits in integer implementation do not have enough bits";
fi;
fun from_multiword_int (x: multiword_int::Int)
=
{
my core_multiword_int::BI { negative, digits } = core_multiword_int::concrete x;
w2r = from_int o inline_t::tu::copyt_tagged_int;
# Looking at at most 3 "big digits" is always
# enough to get 53 bits of precision...
# (See insanity insurance above.)
fun dosign (x: Float)
=
if negative -x;
else x;
fi;
fun calc (k, d1, d2, d3, [])
=>
dosign (runtime::asm::scalb (w2r d1 +
rbase * (w2r d2 + rbase * w2r d3),
k));
calc (k, _, d1, d2, d3 ! r)
=>
calc (k + base_bits, d1, d2, d3, r);
end;
case digits
[] => 0.0;
[d] => dosign (w2r d);
[d1, d2] => dosign (rbase * w2r d2 + w2r d1);
d1 ! d2 ! d3 ! r => calc (0, d1, d2, d3, r);
esac;
};
# whole and split could be implemented more efficiently if we had
# control over the rounding mode; but for now we don't. XXX BUGGO FIXME
fun whole x
=
if (x > 0.0)
if (x > 0.5)
x - 0.5+max_int-max_int;
else
whole (x+1.0) - 1.0;
fi;
else
if (x < 0.0)
if (x < -0.5)
x + 0.5 - max_int + max_int;
else
whole (x - 1.0)+1.0;
fi;
else
x;
fi;
fi;
fun split x
=
{ w = whole x;
f = x-w;
if (abs f ==== 1.0) { whole=>w+f, frac=>0.0 };
else { whole=>w, frac=>f };
fi;
};
fun float_mod x
=
{ f = x - whole x;
if (abs f ==== 1.0) 0.0;
else f;
fi;
};
fun rem (x, y)
=
y * .frac (split (x/y));
fun check_float x
=
if (x>neg_inf and x<pos_inf) x;
elif (is_nan x) raise exception exceptions_guts::DIVIDE_BY_ZERO; # exceptions_guts is from src/lib/std/src/exceptions-guts.pkg else raise exception exceptions_guts::OVERFLOW;
fi;
fun to_multiword_int mode x
=
if (is_nan x ) raise exception DOMAIN;
elif (x ==== pos_inf or x ==== neg_inf ) raise exception OVERFLOW;
else
my (negative, x) =
if (x < 0.0 ) (TRUE, -x); else (FALSE, x);fi;
fun feven x = .frac (split (x / 2.0)) ==== 0.0;
# If the magnitute is less than 1.0, then
# we just have to figure out whether to return -1, 0, or 1
#
if (x < 1.0 )
case mode
ieee_float::TO_ZERO => 0;
ieee_float::TO_POSINF =>
if negative 0; else 1;fi;
ieee_float::TO_NEGINF =>
if negative -1; else 0;fi;
ieee_float::TO_NEAREST =>
if (x < 0.5 ) 0;
elif (x > 0.5 ) if negative -1; else 1;fi;
else 0;
fi;
esac; # 0 is even
else
# Otherwise we start with an integral value,
# suitably adjusted according to fractional part
# and rounding mode:
my { whole, frac } = split x;
start =
case mode
ieee_float::TO_NEGINF =>
if (frac > 0.0 and negative)
whole + 1.0;
else whole;fi;
ieee_float::TO_POSINF =>
if (frac > 0.0 and not negative )
whole + 1.0;
else whole;fi;
ieee_float::TO_ZERO => whole;
ieee_float::TO_NEAREST =>
if (frac > 0.5 ) whole + 1.0;
elif (frac < 0.5 ) whole;
elif (feven whole) whole;
else whole + 1.0;
fi;
esac;
# Now, for efficiency, we construct a
# fairly "small" whole number with
# all the significant bits. First
# we get mantissa and exponent:
my { mantissa, exponent } = to_mantissa_exponent start;
# Then we adjust both to make sure the mantissa
# is whole:
# We know that man is between .5 and 1, so
# multiplying man by 2^53 will guarantee wholeness.
# However, exp might be < 53 -- which would be
# bad. The correct solution is to multiply
# by 2^min (exp, 53) and adjust exp by subtracting
# min (exp, 53):
adj = int_guts::min (precision, exponent);
man = from_mantissa_exponent { mantissa, exponent => adj };
exponent = exponent - adj;
# Now we can construct our bignum digits by
# repeated div/mod using the bignum base.
# This loop will terminate after two rounds at
# the most because we chop off 30 bits each
# time:
fun loop x
=
if (x ==== 0.0)
#
[];
else
my { whole, frac } = split (x / rbase);
dig = inline_t::tu::copyf_tagged_int
(runtime::asm::floor
(frac * rbase));
dig ! loop whole;
fi;
# Now we make a bignum out of those digits:
iman =
core_multiword_int::abstract
(core_multiword_int::BI { negative,
digits => loop man } );
# Finally, we have to put the exponent back
# into the picture:
multiword_int_guts::(<<) (iman, inline_t::tu::copyf_tagged_int exponent);
fi;
fi;
fun next_after _
=
raise exception DIE "float::nextAfter unimplemented";
my min: (Float, Float) -> Float = inline_t::f64::min;
my max: (Float, Float) -> Float = inline_t::f64::max;
fun to_decimal _ = raise exception DIE "float::toDecimal unimplemented";
fun from_decimal _ = raise exception DIE "float::fromDecimal unimplemented";
format = float_format::format_float;
to_string = format (number_string::GEN NULL);
scan = number_scan::scan_real;
from_string = number_string::scan_string scan;
my (-_) = inline_t::f64::neg;
my neg = inline_t::f64::neg;
my (+) = inline_t::f64::(+);
my (-) = inline_t::f64::(-);
my (*) = inline_t::f64::(*);
my (/) = inline_t::f64::(/);
my (>) = inline_t::f64::(>);
my (<) = inline_t::f64::(<);
my (>=) = inline_t::f64::(>=);
my (<=) = inline_t::f64::(<=);
fun sum floats # Should maybe recode this at some point to first sum pairs, then sum pairs in the resultlist, etc: Less prone to floating point underflow.
=
sum' (floats, 0.0)
where
fun sum' ( [], result) => result;
sum' (i ! rest, result) => sum' (rest, i + result);
end;
end;
fun product floats
=
product' (floats, 1.0)
where
fun product' ( [], result) => result;
product' (i ! rest, result) => product' (rest, i * result);
end;
end;
fun list_min [] => raise exception DIE "Cannot do float::list_min on empty list";
#
list_min (f ! floats)
=>
min' (floats, f: Float)
where
fun min' ( [], result) => result;
min' (f ! rest, result) => min' (rest, f < result ?? f :: result);
end;
end;
end;
fun list_max [] => raise exception DIE "Cannot do float::list_max on empty list";
#
list_max (f ! floats)
=>
min' (floats, f: Float)
where
fun min' ( [], result) => result;
min' (f ! rest, result) => min' (rest, f > result ?? f :: result);
end;
end;
end;
fun sort floats
=
lms::sort_list (>) floats;
fun sort_and_drop_duplicates floats
=
lms::sort_list_and_drop_duplicates compare floats;
fun mean [] => 0.0; # Would throwing an exception be better? In graphics, at least, often it is better to just gloss over the occasional special case...
mean floats => sum floats / from_int (length floats);
end;
fun median []
=>
0.0; # As above, arbitrary, possibly should throw exception.
median floats
=>
{ len = length floats;
floats = lms::sort_list (>) floats;
#
i1 = tagged_int_guts::(/) (len,2); # Regular / and + defs are not yet in effect at this level of the library.
i2 = tagged_int_guts::(-) (i1, 1);
if (is_odd(len))
#
# Return middle element:
#
list::nth (floats, i1);
else
# Return average of the two middle elements:
#
f1 = list::nth (floats, i1);
f2 = list::nth (floats, i2);
(f1 + f2) * 0.5;
fi;
}
where
fun is_odd(i) = (i & 1 == 1);
end;
end;
};
end;