commit 6ad2983f94faf2374168862c2aa7f303558fe79e Author: Paolo Bonzini Date: Sat Jul 11 12:05:04 2009 +0200 fix bugs in reading floats libgst: 2009-07-11 Paolo Bonzini * 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 + + * libgst/lex.c: Use new real.c interface. + * libgst/real.c: New. + * libgst/real.h: New. + 2009-07-02 Paolo Bonzini * 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