Hi all and especially Nicolas,
here is my shotgun approach to fixing the reading of FloatQ. The
problem lies in Smalltalk's support for arbitrary bases, which does not
allow me to use the optimized strtod function from the C library.
I decided to write my own tiny floating-point math library (I couldn't
decide whether using mpf/mpfr was *more* or *less* shotgun) so that I
can read the FP numbers in 160 bit precision. This is somewhat dual to
what the Smalltalk code does, which uses fractions and multi-precision
integers, and it is also how GCC works by the way (though they have more
stuff to do because of handling denormals, NaNs and infinities).
While I tried it a bit, there might be bugs, so I ask you all to figure
out some stress tests. These are needed because at some point all
floats were returned multiplied by 2 :-) and "make" still built... I
tried on the ones from http://smalltalk.gnu.org/project/issue/313, as
well as a few others like 1q3000/1q2999 or 3.6451995318824746025q-4951
(FloatQ fmin), and they seem to work. Ideas on how to incorporate these
into floatmath.st as portably as possible are also welcome.
The code is already in git.
Paolo
commit 6ad2983f94faf2374168862c2aa7f303558fe79e
Author: Paolo Bonzini <[email protected]>
Date: Sat Jul 11 12:05:04 2009 +0200
fix bugs in reading floats
libgst:
2009-07-11 Paolo Bonzini <[email protected]>
* libgst/lex.c: Use new real.c interface.
* libgst/real.c: New.
* libgst/real.h: New.
diff --git a/Makefile.am b/Makefile.am
index 374d31c..1a5dc43 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -102,7 +102,7 @@ bin_PROGRAMS = gst
gst_SOURCES = main.c
gst_LDADD = libgst/libgst.la lib-src/library.la @ICON@
gst_DEPENDENCIES = libgst/libgst.la lib-src/library.la @ICON@
-gst_LDFLAGS = -export-dynamic $(RELOC_LDFLAGS)
+gst_LDFLAGS = -export-dynamic $(RELOC_LDFLAGS) -static
if ENABLE_DISASSEMBLER
gst_LDADD += opcode/libdisass.la
diff --git a/libgst/ChangeLog b/libgst/ChangeLog
index e637d0d..68d63fe 100644
--- a/libgst/ChangeLog
+++ b/libgst/ChangeLog
@@ -1,3 +1,9 @@
+2009-07-11 Paolo Bonzini <[email protected]>
+
+ * libgst/lex.c: Use new real.c interface.
+ * libgst/real.c: New.
+ * libgst/real.h: New.
+
2009-07-02 Paolo Bonzini <[email protected]>
* libgst/input.c: Include all symbols after popular demand for
diff --git a/libgst/Makefile.am b/libgst/Makefile.am
index 404cfff..a7c3ad8 100644
--- a/libgst/Makefile.am
+++ b/libgst/Makefile.am
@@ -36,7 +36,7 @@ libgst_la_SOURCES = \
save.c cint.c heap.c input.c \
sysdep.c callin.c xlat.c events.c \
mpz.c print.c alloc.c security.c \
- re.c interp.c
+ re.c interp.c real.c
# definitions for genprims
diff --git a/libgst/gstpriv.h b/libgst/gstpriv.h
index 36d71a5..f86d37c 100644
--- a/libgst/gstpriv.h
+++ b/libgst/gstpriv.h
@@ -612,6 +612,7 @@ extern OOP _gst_nil_oop
#include "mpz.h"
#include "print.h"
#include "security.h"
+#include "real.h"
/* Include this last, it has the bad habit of #defining printf
and this fools gcc's __attribute__ (format) */
diff --git a/libgst/lex.c b/libgst/lex.c
index 2835ac7..f644d80 100644
--- a/libgst/lex.c
+++ b/libgst/lex.c
@@ -119,9 +119,11 @@ static mst_Boolean is_base_digit (int c,
the digits are stored in an obstack, and LARGEINTEGER is set to true
if numPtr does not have sufficient precision. */
static int scan_fraction (int c,
- int base,
- long double *numPtr,
- mst_Boolean *largeInteger);
+ mst_Boolean negative,
+ unsigned base,
+ uintptr_t *intNumPtr,
+ struct real *numPtr,
+ mst_Boolean *largeInteger);
/* Parse a numeric constant and return it. Read numbers in
base-BASE, the first one being C. If a - was parsed, NEGATIVE
@@ -129,10 +131,11 @@ static int scan_fraction (int c,
If LARGEINTEGER is not NULL, the digits are stored in an obstack,
and LARGEINTEGER is set to true if the return value does not have
sufficient precision. */
-static long double scan_digits (int c,
- mst_Boolean negative,
- int base,
- mst_Boolean * largeInteger);
+static uintptr_t scan_digits (int c,
+ mst_Boolean negative,
+ unsigned base,
+ struct real * n,
+ mst_Boolean * largeInteger);
/* Parse the large integer constant stored as base-BASE
digits in the buffer maintained by str.c, adjusting
@@ -210,11 +213,6 @@ static int scan_symbol (int c,
static int digit_to_int (int c,
int base);
-/* Raise BASE to the N-th power and multiply MANT by the result. */
-static inline long double mul_powl (long double mant,
- long double base,
- int n);
-
#ifdef LEXDEBUG
static void print_token (int token,
YYSTYPE *yylval);
@@ -745,45 +743,19 @@ scan_ident (int c,
}
-long double
-mul_powl (long double mant,
- long double base,
- int n)
-{
- long double result = 1.0;
- int n_abs = n < 0 ? -n : n;
- int exp;
- int k = 1;
-
- /* Restrict base to the range from 1.0 to 2.0. */
- base = frexpl (base, &exp);
- base *= 2;
- exp--;
-
- while (n_abs)
- {
- if (n_abs & k)
- {
- result *= base;
- n_abs ^= k;
- }
- base *= base;
- k <<= 1;
- }
-
- if (n < 0)
- return ldexpl (mant, exp * n) / result;
- else
- return ldexpl (mant * result, exp * n);
-}
+/* TODO: We track the number in *three* formats: struct real, uintptr_t,
+ * and just save the bytes for large integers. We should just save
+ * the bytes and work on those. */
int
scan_number (int c,
YYSTYPE * lvalp)
{
- OOP numOOP;
+ OOP intNumOOP;
int base, exponent, ic;
- long double num, floatExponent;
+ uintptr_t intNum;
+ struct real num, dummy;
+ int floatExponent;
mst_Boolean isNegative = false, largeInteger = false;
int float_type = 0;
@@ -792,19 +764,20 @@ scan_number (int c,
ic = c;
assert (ic != '-');
- num = scan_digits (ic, false, 10, &largeInteger);
+ intNum = scan_digits (ic, false, 10, &num, &largeInteger);
ic = _gst_next_char ();
if (ic == 'r')
{
char *p = obstack_finish (_gst_compilation_obstack);
obstack_free (_gst_compilation_obstack, p);
- base = (int) num;
- if (base > 36 || largeInteger)
+ if (intNum > 36 || largeInteger)
{
_gst_errorf ("Numeric base too large %d", base);
_gst_had_error = true;
}
+ else
+ base = intNum;
ic = _gst_next_char ();
/* Having to support things like 16r-123 is a pity :-) because we
@@ -815,7 +788,7 @@ scan_number (int c,
ic = _gst_next_char ();
}
- num = scan_digits (ic, isNegative, base, &largeInteger);
+ intNum = scan_digits (ic, isNegative, base, &num, &largeInteger);
ic = _gst_next_char ();
}
@@ -834,7 +807,7 @@ scan_number (int c,
else
{
float_type = FLOATD_LITERAL;
- exponent = scan_fraction (ic, base, &num, &largeInteger);
+ exponent = scan_fraction (ic, isNegative, base, &intNum, &num,
&largeInteger);
ic = _gst_next_char ();
}
}
@@ -852,7 +825,7 @@ scan_number (int c,
else if (CHAR_TAB (ic)->char_class & DIGIT)
{
/* 123s4 format -- parse the exponent */
- floatExponent = scan_digits (ic, false, 10, NULL);
+ floatExponent = scan_digits (ic, false, 10, &dummy, NULL);
}
else if (CHAR_TAB (ic)->char_class & ID_CHAR)
{
@@ -864,29 +837,26 @@ scan_number (int c,
else
_gst_unread_char (ic);
- if (isNegative)
- num = -num;
-
- if (largeInteger || num < MIN_ST_INT || num > MAX_ST_INT)
+ if (largeInteger)
{
/* Make a LargeInteger constant and create an object out of
it. */
byte_object bo = scan_large_integer (isNegative, base);
- gst_object result = instantiate_with (bo->class, bo->size, &numOOP);
+ gst_object result = instantiate_with (bo->class, bo->size,
&intNumOOP);
memcpy (result->data, bo->body, bo->size);
}
else
- numOOP = FROM_INT(num);
+ intNumOOP = FROM_INT((intptr_t) (isNegative ? -intNum : intNum));
/* too much of a chore to create a Fraction, so we call-in. We
lose the ability to create ScaledDecimals during the very
first phases of bootstrapping, but who cares?...
This is equivalent to
- (num * (10 raisedToInteger: exponent)
+ (intNumOOP * (10 raisedToInteger: exponent)
asScaledDecimal: floatExponent) */
lvalp->oval =
- _gst_msg_send (numOOP, _gst_as_scaled_decimal_radix_scale_symbol,
+ _gst_msg_send (intNumOOP, _gst_as_scaled_decimal_radix_scale_symbol,
FROM_INT (exponent),
FROM_INT (base),
FROM_INT ((int) floatExponent),
@@ -914,12 +884,12 @@ scan_number (int c,
;
else if (ic == '-') {
floatExponent =
- scan_digits (_gst_next_char (), true, 10, NULL);
+ scan_digits (_gst_next_char (), true, 10, &dummy, NULL);
exponent -= (int) floatExponent;
}
else if (CHAR_TAB (ic)->char_class & DIGIT)
{
- floatExponent = scan_digits (ic, false, 10, NULL);
+ floatExponent = scan_digits (ic, false, 10, &dummy, NULL);
exponent += (int) floatExponent;
}
else if (CHAR_TAB (ic)->char_class & ID_CHAR)
@@ -935,24 +905,25 @@ scan_number (int c,
else
_gst_unread_char (ic);
- if (exponent != 0)
- num = mul_powl (num, (long double) base, exponent);
-
- if (isNegative)
- num = -num;
-
-#if defined(__FreeBSD__)
- fpsetmask (0);
-#endif
-
if (float_type)
{
char *p = obstack_finish (_gst_compilation_obstack);
obstack_free (_gst_compilation_obstack, p);
- lvalp->fval = num;
+
+ if (exponent)
+ {
+ struct real r;
+ _gst_real_from_int (&r, base);
+ _gst_real_powi (&r, &r, exponent < 0 ? -exponent : exponent);
+ if (exponent < 0)
+ _gst_real_div (&num, &num, &r);
+ else
+ _gst_real_mul (&num, &r);
+ }
+ lvalp->fval = _gst_real_get_ld (&num);
return (float_type);
}
- else if (largeInteger || num < MIN_ST_INT || num > MAX_ST_INT)
+ else if (largeInteger)
{
lvalp->boval = scan_large_integer (isNegative, base);
return (LARGE_INTEGER_LITERAL);
@@ -961,36 +932,44 @@ scan_number (int c,
{
char *p = obstack_finish (_gst_compilation_obstack);
obstack_free (_gst_compilation_obstack, p);
- lvalp->ival = (intptr_t) num;
+ lvalp->ival = (intptr_t) (isNegative ? -intNum : intNum);
return (INTEGER_LITERAL);
}
}
-long double
+uintptr_t
scan_digits (int c,
- mst_Boolean negative,
- int base,
- mst_Boolean * largeInteger)
+ mst_Boolean negative,
+ unsigned base,
+ struct real * n,
+ mst_Boolean * largeInteger)
{
- long double result;
+ uintptr_t result;
mst_Boolean oneDigit = false;
while (c == '_')
c = _gst_next_char ();
+ memset (n, 0, sizeof (*n));
for (result = 0.0; is_base_digit (c, base); )
{
- result *= base;
- oneDigit = true;
- result += digit_to_int (c, base);
+ unsigned value = digit_to_int (c, base);
if (largeInteger)
{
obstack_1grow (_gst_compilation_obstack, digit_to_int (c, base));
- if (result > -MIN_ST_INT
- || (!negative && result > MAX_ST_INT))
+ if (result >
+ (negative
+ /* We want (uintptr_t) -MIN_ST_INT, but it's the same. */
+ ? (uintptr_t) MIN_ST_INT - value
+ : (uintptr_t) MAX_ST_INT - value) / base)
*largeInteger = true;
}
+ _gst_real_mul_int (n, base);
+ _gst_real_add_int (n, value);
+ oneDigit = true;
+ result *= base;
+ result += value;
do
c = _gst_next_char ();
while (c == '_');
@@ -1009,30 +988,39 @@ scan_digits (int c,
int
scan_fraction (int c,
- int base,
- long double *numPtr,
- mst_Boolean *largeInteger)
+ mst_Boolean negative,
+ unsigned base,
+ uintptr_t *intNumPtr,
+ struct real *numPtr,
+ mst_Boolean *largeInteger)
{
+ uintptr_t intNum;
int scale;
- long double num;
scale = 0;
while (c == '_')
c = _gst_next_char ();
- for (num = *numPtr; is_base_digit (c, base); )
+ for (intNum = *intNumPtr; is_base_digit (c, base); )
{
- num *= base;
- num += digit_to_int (c, base);
- scale--;
-
+ unsigned value = digit_to_int (c, base);
if (largeInteger)
- {
- obstack_1grow (_gst_compilation_obstack, digit_to_int (c, base));
- if (num > MAX_ST_INT)
- *largeInteger = true;
- }
+ {
+ obstack_1grow (_gst_compilation_obstack, digit_to_int (c, base));
+ if (intNum >
+ (negative
+ /* We want (uintptr_t) -MIN_ST_INT, but it's the same. */
+ ? (uintptr_t) MIN_ST_INT - value
+ : (uintptr_t) MAX_ST_INT - value) / base)
+ *largeInteger = true;
+ }
+
+ _gst_real_mul_int (numPtr, base);
+ _gst_real_add_int (numPtr, value);
+ intNum *= base;
+ intNum += value;
+ scale--;
do
c = _gst_next_char ();
@@ -1041,7 +1029,7 @@ scan_fraction (int c,
_gst_unread_char (c);
- *numPtr = num;
+ *intNumPtr = intNum;
return scale;
}
diff --git a/libgst/real.c b/libgst/real.c
new file mode 100644
index 0000000..a7aaab2
--- /dev/null
+++ b/libgst/real.c
@@ -0,0 +1,533 @@
+/******************************** -*- C -*- ****************************
+ *
+ * Simple floating-point data type
+ *
+ *
+ ***********************************************************************/
+
+
+/***********************************************************************
+ *
+ * Copyright 2009 Free Software Foundation, Inc.
+ * Written by Paolo Bonzini.
+ *
+ * This file is part of GNU Smalltalk.
+ *
+ * GNU Smalltalk is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by the Free
+ * Software Foundation; either version 2, or (at your option) any later
+ * version.
+ *
+ * Linking GNU Smalltalk statically or dynamically with other modules is
+ * making a combined work based on GNU Smalltalk. Thus, the terms and
+ * conditions of the GNU General Public License cover the whole
+ * combination.
+ *
+ * In addition, as a special exception, the Free Software Foundation
+ * give you permission to combine GNU Smalltalk with free software
+ * programs or libraries that are released under the GNU LGPL and with
+ * independent programs running under the GNU Smalltalk virtual machine.
+ *
+ * You may copy and distribute such a system following the terms of the
+ * GNU GPL for GNU Smalltalk and the licenses of the other code
+ * concerned, provided that you include the source code of that other
+ * code when and as the GNU GPL requires distribution of source code.
+ *
+ * Note that people who make modified versions of GNU Smalltalk are not
+ * obligated to grant this special exception for their modified
+ * versions; it is their choice whether to do so. The GNU General
+ * Public License gives permission to release a modified version without
+ * this exception; this exception also makes it possible to release a
+ * modified version which carries forward this exception.
+ *
+ * GNU Smalltalk is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
+ * more details.
+ *
+ * You should have received a copy of the GNU General Public License along with
+ * GNU Smalltalk; see the file COPYING. If not, write to the Free Software
+ * Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ ***********************************************************************/
+
+
+#include "gstpriv.h"
+
+#define SIG_ELEM_BITS (CHAR_BIT * sizeof (((struct real *)0)->sig[0]))
+#define NUM_SIG_BITS (SIGSZ * SIG_ELEM_BITS)
+#define SIG_MASK ((1 << SIG_ELEM_BITS) - 1)
+#define SIG_MSB (1 << (SIG_ELEM_BITS - 1))
+
+/* Shift the significant of IN by DELTA bits and store the result into
+ OUT. OUT and IN can overlap. */
+static int rshift_sig (struct real *out, struct real *in, int delta);
+
+/* Normalize the significant of IN so that the most significant bit is 1,
+ and store the result into OUT. OUT and IN can overlap. */
+static void normalize (struct real *out, struct real *in);
+
+/* Denormalize IN to increase its exponent to EXP and store the result
+ into OUT. OUT and IN can overlap. Return false if OUT would be zero,
+ in which case its contents are undefined. */
+static int adjust_exp (struct real *out, struct real *in, int exp);
+
+/* Sum the significands of R and S. Return the carry. */
+static int add_significands (struct real *r, struct real *s);
+
+/* Compare the significands of R and S and return the result like strcmp. */
+static int cmp_significands (struct real *r, struct real *s);
+
+/* Subtract the significands of R and S. */
+static void sub_significands (struct real *r, struct real *s);
+
+/* Sum S into R. */
+static void do_add (struct real *r, struct real *s);
+
+/* Multiply S into R. LSB is the least significant nonzero byte of the
+ significand of S, and is used to cut the iteration. */
+static void do_mul (struct real *r, struct real *s, int lsb);
+
+/* Divide R by S and store the result into S. OUT can overlap either R
+ or S, but R must not be the same as S. R is destroyed */
+static void do_div (struct real *out, struct real *r, struct real *s);
+
+/* These routines are not optimized at all. Maybe I should have bit
+ the bullet and required MPFR after all... */
+
+static int
+rshift_sig (struct real *out, struct real *in, int delta)
+{
+ int i, nonzero = 0;
+ int rshift = delta & (SIG_ELEM_BITS - 1);
+ int lshift = SIG_ELEM_BITS - rshift;
+
+ delta /= SIG_ELEM_BITS;
+ if (rshift)
+ {
+ for (i = 0; i + delta < SIGSZ - 1; i++)
+ {
+ out->sig[i] = ((in->sig[i + delta + 1] << lshift)
+ | (in->sig[i + delta] >> rshift));
+ nonzero |= out->sig[i];
+ }
+ out->sig[i] = in->sig[i + delta] >> rshift;
+ nonzero |= out->sig[i];
+ i++;
+ }
+ else
+ {
+ for (i = 0; i + delta < SIGSZ; i++)
+ {
+ out->sig[i] = in->sig[i + delta];
+ nonzero |= out->sig[i];
+ }
+ }
+
+ while (i < SIGSZ)
+ out->sig[i++] = 0;
+ return nonzero;
+}
+
+static void
+normalize (struct real *out, struct real *in)
+{
+ int i, msb, delta, rshift, lshift;
+ int out_exp;
+
+ out_exp = in->exp;
+ delta = 0;
+ for (i = SIGSZ; --i >= 0 && in->sig[i] == 0; )
+ {
+ delta++;
+ out_exp -= SIG_ELEM_BITS;
+ }
+
+ if (i < 0)
+ {
+ memset (out, 0, sizeof (struct real));
+ return;
+ }
+
+ /* TODO: convert this to clz. */
+ msb = in->sig[i];
+ lshift = 15;
+ if (msb & 0xFF00)
+ lshift -= 8;
+ else
+ msb <<= 8;
+ if (msb & 0xF000)
+ lshift -= 4;
+ else
+ msb <<= 4;
+ if (msb & 0xC000)
+ lshift -= 2;
+ else
+ msb <<= 2;
+ if (msb & 0x8000)
+ lshift -= 1;
+
+ rshift = 16 - lshift;
+ out->exp = out_exp - lshift;
+ out->sign = in->sign;
+ if (lshift)
+ {
+ for (i = SIGSZ; --i - delta >= 1; )
+ out->sig[i] = ((in->sig[i - delta] << lshift)
+ | (in->sig[i - delta - 1] >> rshift));
+ out->sig[i] = in->sig[0] << lshift;
+ }
+ else
+ {
+ for (i = SIGSZ; --i - delta >= 0; )
+ out->sig[i] = in->sig[i - delta];
+ }
+
+ while (--i >= 0)
+ out->sig[i] = 0;
+}
+
+/* Adjust IN to have exponent EXP by shifting its significand right.
+ Put the result into OUT. The shift can be done in place. */
+static int
+adjust_exp (struct real *out, struct real *in, int exp)
+{
+ int in_exp;
+ in_exp = in->exp;
+ assert (exp > in_exp);
+ if (exp == in_exp)
+ return true;
+ if (exp - in_exp >= NUM_SIG_BITS)
+ return false;
+
+ out->exp = exp;
+ return rshift_sig (out, in, exp - in_exp);
+}
+
+
+void
+_gst_real_from_int (struct real *out, int s)
+{
+ memset (out, 0, sizeof (struct real));
+ if (s < 0)
+ {
+ out->sign = -1;
+ s = -s;
+ }
+ else
+ out->sign = 1;
+
+ /* TODO: convert this to clz. */
+ if (s & 0xFF00)
+ out->exp += 8;
+ else
+ s <<= 8;
+ if (s & 0xF000)
+ out->exp += 4;
+ else
+ s <<= 4;
+ if (s & 0xC000)
+ out->exp += 2;
+ else
+ s <<= 2;
+ if (s & 0x8000)
+ out->exp += 1;
+ else
+ s <<= 1;
+
+ out->sig[SIGSZ - 1] = s;
+}
+
+static int
+add_significands (struct real *r, struct real *s)
+{
+ int i, carry = 0;
+ for (i = 0; i < SIGSZ; i++)
+ {
+ int result = r->sig[i] + s->sig[i] + carry;
+ carry = result >> SIG_ELEM_BITS;
+ r->sig[i] = result;
+ }
+
+ return carry;
+}
+
+static int
+cmp_significands (struct real *r, struct real *s)
+{
+ int i;
+ for (i = SIGSZ; --i >= 0; )
+ if (r->sig[i] != s->sig[i])
+ return (r->sig[i] - s->sig[i]);
+
+ return 0;
+}
+
+static void
+sub_significands (struct real *r, struct real *s)
+{
+ int i, carry = 0;
+ for (i = 0; i < SIGSZ; i++)
+ {
+ int result = r->sig[i] - s->sig[i] + carry;
+ carry = result >> SIG_ELEM_BITS;
+ r->sig[i] = result;
+ }
+}
+
+static void
+do_add (struct real *r, struct real *s)
+{
+ struct real tmp;
+
+ if (r->exp < s->exp)
+ {
+ if (!adjust_exp (r, r, s->exp))
+ {
+ /* Underflow, R+S = S. */
+ *r = *s;
+ return;
+ }
+ }
+
+ else if (r->exp > s->exp)
+ {
+ /* We cannot modify S in place, use a temporary. */
+ if (!adjust_exp (&tmp, s, r->exp))
+ return;
+ s = &tmp;
+ }
+
+ if (add_significands (r, s))
+ {
+ /* Lose one bit of precision to fit the carry. */
+ rshift_sig (r, r, 1);
+ r->exp++;
+ r->sig[SIGSZ - 1] |= SIG_MSB;
+ }
+}
+
+void
+_gst_real_add (struct real *r, struct real *s)
+{
+ if (!s->sign)
+ return;
+ if (!r->sign)
+ memcpy (r, s, sizeof (struct real));
+ else if (s->sign == r->sign)
+ return do_add (r, s);
+ else
+ abort ();
+}
+
+void
+_gst_real_add_int (struct real *r, int s)
+{
+ struct real s_real;
+ if (!s)
+ return;
+
+ _gst_real_from_int (&s_real, s);
+ if (!r->sign)
+ memcpy (r, &s_real, sizeof (struct real));
+ else if (s_real.sign == r->sign)
+ return do_add (r, &s_real);
+ else
+ abort ();
+}
+
+static void
+do_mul (struct real *r, struct real *s, int lsb)
+{
+ struct real rr;
+ unsigned short mask;
+ int n;
+
+ r->exp += s->exp;
+ r->sign *= s->sign;
+
+ rr = *r;
+ mask = SIG_MSB;
+ n = SIGSZ - 1;
+ assert (s->sig[n] & mask);
+ while (n > lsb || (s->sig[n] & (mask - 1)))
+ {
+ if (!(mask >>= 1))
+ {
+ mask = SIG_MSB;
+ n--;
+ }
+
+ /* Dividing rr by 2 matches the weight s->sig[n] & mask. Exit
+ early in case of underflow. */
+ if (!rshift_sig (&rr, &rr, 1))
+ break;
+
+ if (s->sig[n] & mask)
+ {
+ if (add_significands (r, &rr))
+ {
+ /* Lose one bit of precision to fit the carry. */
+ rshift_sig (r, r, 1);
+ r->sig[SIGSZ - 1] |= SIG_MSB;
+ r->exp++;
+ if (!rshift_sig (&rr, &rr, 1))
+ break;
+ rr.exp++;
+ }
+ }
+ }
+}
+
+void
+_gst_real_mul (struct real *r, struct real *s)
+{
+ int i;
+ struct real tmp;
+ if (r->sign == 0)
+ return;
+ if (r == s)
+ {
+ tmp = *s;
+ s = &tmp;
+ }
+ if (s->sign == 0)
+ memset (r, 0, sizeof (struct real));
+
+ for (i = 0; i <= SIGSZ && s->sig[i] == 0; )
+ i++;
+ do_mul (r, s, i);
+}
+
+void
+_gst_real_mul_int (struct real *r, int s)
+{
+ struct real s_real;
+
+ if (s == 0)
+ memset (r, 0, sizeof (struct real));
+ if (r->sign == 0)
+ return;
+
+ _gst_real_from_int (&s_real, s);
+ do_mul (r, &s_real, SIGSZ - 1);
+}
+
+
+void
+_gst_real_powi (struct real *out, struct real *r, int s)
+{
+ int k;
+ struct real tmp;
+ if (out == r)
+ {
+ tmp = *r;
+ r = &tmp;
+ }
+
+ assert (s > 0);
+ _gst_real_from_int (out, 1);
+ if (!s)
+ return;
+
+ for (k = 1;; k <<= 1)
+ {
+ if (s & k)
+ {
+ _gst_real_mul (out, r);
+ s ^= k;
+ if (!s)
+ break;
+ }
+ _gst_real_mul (r, r);
+ }
+}
+
+static void
+do_div (struct real *out, struct real *r, struct real *s)
+{
+ struct real v;
+ int msb, i;
+ int place = SIGSZ-1;
+ int bit = SIG_MSB;
+
+ memset (&v, 0, sizeof (struct real));
+ v.sign = r->sign * s->sign;
+ v.exp = r->exp - s->exp;
+ msb = 0;
+ goto start;
+ do
+ {
+ /* Get the MSB of U and shift it left by one. */
+ msb = r->sig[SIGSZ-1] & SIG_MSB;
+ for (i = SIGSZ; --i >= 1; )
+ r->sig[i] = (r->sig[i] << 1) | (r->sig[i - 1] >> 15);
+ r->sig[0] <<= 1;
+
+ start:
+ if (msb || cmp_significands (r, s) >= 0)
+ {
+ sub_significands (r, s);
+ v.sig[place] |= bit;
+ }
+ }
+ while ((bit >>= 1) || (bit = SIG_MSB, --place >= 0));
+ normalize (out, &v);
+}
+
+void
+_gst_real_div (struct real *out, struct real *r, struct real *s)
+{
+ assert (s->sign);
+ if (!r->sign)
+ {
+ memset (out, 0, sizeof (struct real));
+ return;
+ }
+
+ if (r == s)
+ {
+ memset (out, 0, sizeof (struct real));
+ out->sign = 1;
+ out->sig[SIGSZ-1] = SIG_MSB;
+ return;
+ }
+
+ if (out == r)
+ do_div (out, out, s);
+ else
+ {
+ /* do_div would destroy R, save it. */
+ struct real u = *r;
+ do_div (out, &u, s);
+ }
+}
+
+
+void
+_gst_real_inv (struct real *out, struct real *s)
+{
+ struct real u;
+ assert (s->sign);
+ memset (&u, 0, sizeof (struct real));
+ u.sign = 1;
+ u.sig[SIGSZ-1] = SIG_MSB;
+ do_div (out, &u, s);
+}
+
+
+long double
+_gst_real_get_ld (struct real *r)
+{
+ long double result = 0.0;
+ int i;
+
+ for (i = SIGSZ; --i >= 0; )
+ {
+ result *= SIG_MASK + 1;
+ result += r->sig[i];
+ }
+
+ result = ldexpl (result, r->exp - NUM_SIG_BITS + 1);
+ return r->sign == -1 ? -result : result;
+}
diff --git a/libgst/real.h b/libgst/real.h
new file mode 100644
index 0000000..588b658
--- /dev/null
+++ b/libgst/real.h
@@ -0,0 +1,105 @@
+/******************************** -*- C -*- ****************************
+ *
+ * Simple floating-point data type
+ *
+ *
+ ***********************************************************************/
+
+
+/***********************************************************************
+ *
+ * Copyright 2009 Free Software Foundation, Inc.
+ * Written by Paolo Bonzini.
+ *
+ * This file is part of GNU Smalltalk.
+ *
+ * GNU Smalltalk is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by the Free
+ * Software Foundation; either version 2, or (at your option) any later
+ * version.
+ *
+ * Linking GNU Smalltalk statically or dynamically with other modules is
+ * making a combined work based on GNU Smalltalk. Thus, the terms and
+ * conditions of the GNU General Public License cover the whole
+ * combination.
+ *
+ * In addition, as a special exception, the Free Software Foundation
+ * give you permission to combine GNU Smalltalk with free software
+ * programs or libraries that are released under the GNU LGPL and with
+ * independent programs running under the GNU Smalltalk virtual machine.
+ *
+ * You may copy and distribute such a system following the terms of the
+ * GNU GPL for GNU Smalltalk and the licenses of the other code
+ * concerned, provided that you include the source code of that other
+ * code when and as the GNU GPL requires distribution of source code.
+ *
+ * Note that people who make modified versions of GNU Smalltalk are not
+ * obligated to grant this special exception for their modified
+ * versions; it is their choice whether to do so. The GNU General
+ * Public License gives permission to release a modified version without
+ * this exception; this exception also makes it possible to release a
+ * modified version which carries forward this exception.
+ *
+ * GNU Smalltalk is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
+ * more details.
+ *
+ * You should have received a copy of the GNU General Public License along with
+ * GNU Smalltalk; see the file COPYING. If not, write to the Free Software
+ * Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ ***********************************************************************/
+
+
+#ifndef GST_REAL_H
+#define GST_REAL_H
+
+#define SIGSZ 10
+
+struct real {
+ /* Little-endian significant. sig[0] is the least significant part.
+ The 1 is not implicit, so in a normalized real sig[9] == 0 means
+ the value is zero. sig[9]'s MSB has weight 2^exp. */
+ unsigned short sig[SIGSZ];
+
+ int exp;
+ int sign;
+};
+
+/* Convert an integer number S to floating point and move it to OUT. */
+extern void _gst_real_from_int (struct real *out, int s);
+
+/* Sum S to R and store the result into R. */
+extern void _gst_real_add_int (struct real *r, int s)
+ ATTRIBUTE_HIDDEN;
+
+/* Multiply R by S and store the result into R. */
+extern void _gst_real_mul_int (struct real *r, int s)
+ ATTRIBUTE_HIDDEN;
+
+/* Compute R^S and store the result into R. */
+extern void _gst_real_powi (struct real *out, struct real *r, int s)
+ ATTRIBUTE_HIDDEN;
+
+/* Sum S to R and store the result into R. */
+extern void _gst_real_add (struct real *r, struct real *s)
+ ATTRIBUTE_HIDDEN;
+
+/* Multiply R by S and store the result into R. */
+extern void _gst_real_mul (struct real *r, struct real *s)
+ ATTRIBUTE_HIDDEN;
+
+/* Divide R by S and store the result into OUT. */
+extern void _gst_real_div (struct real *out, struct real *r, struct real *s)
+ ATTRIBUTE_HIDDEN;
+
+/* Compute the inverse of R and store it into OUT (OUT and R can overlap). */
+extern void _gst_real_inv (struct real *out, struct real *r)
+ ATTRIBUTE_HIDDEN;
+
+/* Convert R to a long double and return it. */
+extern long double _gst_real_get_ld (struct real *r)
+ ATTRIBUTE_HIDDEN;
+
+#endif
_______________________________________________
help-smalltalk mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/help-smalltalk
