1 /******************************** -*- C -*- ****************************
2 *
3 * Simple floating-point data type
4 *
5 *
6 ***********************************************************************/
7
8
9 /***********************************************************************
10 *
11 * Copyright 2009 Free Software Foundation, Inc.
12 * Written by Paolo Bonzini.
13 *
14 * This file is part of GNU Smalltalk.
15 *
16 * GNU Smalltalk is free software; you can redistribute it and/or modify it
17 * under the terms of the GNU General Public License as published by the Free
18 * Software Foundation; either version 2, or (at your option) any later
19 * version.
20 *
21 * Linking GNU Smalltalk statically or dynamically with other modules is
22 * making a combined work based on GNU Smalltalk. Thus, the terms and
23 * conditions of the GNU General Public License cover the whole
24 * combination.
25 *
26 * In addition, as a special exception, the Free Software Foundation
27 * give you permission to combine GNU Smalltalk with free software
28 * programs or libraries that are released under the GNU LGPL and with
29 * independent programs running under the GNU Smalltalk virtual machine.
30 *
31 * You may copy and distribute such a system following the terms of the
32 * GNU GPL for GNU Smalltalk and the licenses of the other code
33 * concerned, provided that you include the source code of that other
34 * code when and as the GNU GPL requires distribution of source code.
35 *
36 * Note that people who make modified versions of GNU Smalltalk are not
37 * obligated to grant this special exception for their modified
38 * versions; it is their choice whether to do so. The GNU General
39 * Public License gives permission to release a modified version without
40 * this exception; this exception also makes it possible to release a
41 * modified version which carries forward this exception.
42 *
43 * GNU Smalltalk is distributed in the hope that it will be useful, but WITHOUT
44 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
45 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
46 * more details.
47 *
48 * You should have received a copy of the GNU General Public License along with
49 * GNU Smalltalk; see the file COPYING. If not, write to the Free Software
50 * Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
51 *
52 ***********************************************************************/
53
54
55 #include "gstpriv.h"
56
57 #define SIG_ELEM_BITS (CHAR_BIT * sizeof (((struct real *)0)->sig[0]))
58 #define NUM_SIG_BITS (SIGSZ * SIG_ELEM_BITS)
59 #define SIG_MASK ((1 << SIG_ELEM_BITS) - 1)
60 #define SIG_MSB (1 << (SIG_ELEM_BITS - 1))
61
62 /* Shift the significant of IN by DELTA bits and store the result into
63 OUT. OUT and IN can overlap. */
64 static int rshift_sig (struct real *out, struct real *in, int delta);
65
66 /* Normalize the significant of IN so that the most significant bit is 1,
67 and store the result into OUT. OUT and IN can overlap. */
68 static void normalize (struct real *out, struct real *in);
69
70 /* Denormalize IN to increase its exponent to EXP and store the result
71 into OUT. OUT and IN can overlap. Return false if OUT would be zero,
72 in which case its contents are undefined. */
73 static int adjust_exp (struct real *out, struct real *in, int exp);
74
75 /* Sum the significands of R and S. Return the carry. */
76 static int add_significands (struct real *r, struct real *s);
77
78 /* Compare the significands of R and S and return the result like strcmp. */
79 static int cmp_significands (struct real *r, struct real *s);
80
81 /* Subtract the significands of R and S. */
82 static void sub_significands (struct real *r, struct real *s);
83
84 /* Sum S into R. */
85 static void do_add (struct real *r, struct real *s);
86
87 /* Multiply S into R. LSB is the least significant nonzero byte of the
88 significand of S, and is used to cut the iteration. */
89 static void do_mul (struct real *r, struct real *s, int lsb);
90
91 /* Divide R by S and store the result into S. OUT can overlap either R
92 or S, but R must not be the same as S. R is destroyed */
93 static void do_div (struct real *out, struct real *r, struct real *s);
94
95 /* These routines are not optimized at all. Maybe I should have bit
96 the bullet and required MPFR after all... */
97
98 static int
rshift_sig(struct real * out,struct real * in,int delta)99 rshift_sig (struct real *out, struct real *in, int delta)
100 {
101 int i, nonzero = 0;
102 int rshift = delta & (SIG_ELEM_BITS - 1);
103 int lshift = SIG_ELEM_BITS - rshift;
104
105 delta /= SIG_ELEM_BITS;
106 if (rshift)
107 {
108 for (i = 0; i + delta < SIGSZ - 1; i++)
109 {
110 out->sig[i] = ((in->sig[i + delta + 1] << lshift)
111 | (in->sig[i + delta] >> rshift));
112 nonzero |= out->sig[i];
113 }
114 out->sig[i] = in->sig[i + delta] >> rshift;
115 nonzero |= out->sig[i];
116 i++;
117 }
118 else
119 {
120 for (i = 0; i + delta < SIGSZ; i++)
121 {
122 out->sig[i] = in->sig[i + delta];
123 nonzero |= out->sig[i];
124 }
125 }
126
127 while (i < SIGSZ)
128 out->sig[i++] = 0;
129 return nonzero;
130 }
131
132 static void
normalize(struct real * out,struct real * in)133 normalize (struct real *out, struct real *in)
134 {
135 int i, msb, delta, rshift, lshift;
136 int out_exp;
137
138 out_exp = in->exp;
139 delta = 0;
140 for (i = SIGSZ; --i >= 0 && in->sig[i] == 0; )
141 {
142 delta++;
143 out_exp -= SIG_ELEM_BITS;
144 }
145
146 if (i < 0)
147 {
148 memset (out, 0, sizeof (struct real));
149 return;
150 }
151
152 /* TODO: convert this to clz. */
153 msb = in->sig[i];
154 lshift = 15;
155 if (msb & 0xFF00)
156 lshift -= 8;
157 else
158 msb <<= 8;
159 if (msb & 0xF000)
160 lshift -= 4;
161 else
162 msb <<= 4;
163 if (msb & 0xC000)
164 lshift -= 2;
165 else
166 msb <<= 2;
167 if (msb & 0x8000)
168 lshift -= 1;
169
170 rshift = 16 - lshift;
171 out->exp = out_exp - lshift;
172 out->sign = in->sign;
173 if (lshift)
174 {
175 for (i = SIGSZ; --i - delta >= 1; )
176 out->sig[i] = ((in->sig[i - delta] << lshift)
177 | (in->sig[i - delta - 1] >> rshift));
178 out->sig[i] = in->sig[0] << lshift;
179 }
180 else
181 {
182 for (i = SIGSZ; --i - delta >= 0; )
183 out->sig[i] = in->sig[i - delta];
184 }
185
186 while (--i >= 0)
187 out->sig[i] = 0;
188 }
189
190 /* Adjust IN to have exponent EXP by shifting its significand right.
191 Put the result into OUT. The shift can be done in place. */
192 static int
adjust_exp(struct real * out,struct real * in,int exp)193 adjust_exp (struct real *out, struct real *in, int exp)
194 {
195 int in_exp;
196 in_exp = in->exp;
197 assert (exp > in_exp);
198 if (exp == in_exp)
199 return true;
200 if (exp - in_exp >= NUM_SIG_BITS)
201 return false;
202
203 out->exp = exp;
204 return rshift_sig (out, in, exp - in_exp);
205 }
206
207
208 void
_gst_real_from_int(struct real * out,int s)209 _gst_real_from_int (struct real *out, int s)
210 {
211 memset (out, 0, sizeof (struct real));
212 if (s < 0)
213 {
214 out->sign = -1;
215 s = -s;
216 }
217 else
218 out->sign = 1;
219
220 /* TODO: convert this to clz. */
221 if (s & 0xFF00)
222 out->exp += 8;
223 else
224 s <<= 8;
225 if (s & 0xF000)
226 out->exp += 4;
227 else
228 s <<= 4;
229 if (s & 0xC000)
230 out->exp += 2;
231 else
232 s <<= 2;
233 if (s & 0x8000)
234 out->exp += 1;
235 else
236 s <<= 1;
237
238 out->sig[SIGSZ - 1] = s;
239 }
240
241 static int
add_significands(struct real * r,struct real * s)242 add_significands (struct real *r, struct real *s)
243 {
244 int i, carry = 0;
245 for (i = 0; i < SIGSZ; i++)
246 {
247 int result = r->sig[i] + s->sig[i] + carry;
248 carry = result >> SIG_ELEM_BITS;
249 r->sig[i] = result;
250 }
251
252 return carry;
253 }
254
255 static int
cmp_significands(struct real * r,struct real * s)256 cmp_significands (struct real *r, struct real *s)
257 {
258 int i;
259 for (i = SIGSZ; --i >= 0; )
260 if (r->sig[i] != s->sig[i])
261 return (r->sig[i] - s->sig[i]);
262
263 return 0;
264 }
265
266 static void
sub_significands(struct real * r,struct real * s)267 sub_significands (struct real *r, struct real *s)
268 {
269 int i, carry = 0;
270 for (i = 0; i < SIGSZ; i++)
271 {
272 int result = r->sig[i] - s->sig[i] + carry;
273 carry = result >> SIG_ELEM_BITS;
274 r->sig[i] = result;
275 }
276 }
277
278 static void
do_add(struct real * r,struct real * s)279 do_add (struct real *r, struct real *s)
280 {
281 struct real tmp;
282
283 if (r->exp < s->exp)
284 {
285 if (!adjust_exp (r, r, s->exp))
286 {
287 /* Underflow, R+S = S. */
288 *r = *s;
289 return;
290 }
291 }
292
293 else if (r->exp > s->exp)
294 {
295 /* We cannot modify S in place, use a temporary. */
296 if (!adjust_exp (&tmp, s, r->exp))
297 return;
298 s = &tmp;
299 }
300
301 if (add_significands (r, s))
302 {
303 /* Lose one bit of precision to fit the carry. */
304 rshift_sig (r, r, 1);
305 r->exp++;
306 r->sig[SIGSZ - 1] |= SIG_MSB;
307 }
308 }
309
310 void
_gst_real_add(struct real * r,struct real * s)311 _gst_real_add (struct real *r, struct real *s)
312 {
313 if (!s->sign)
314 return;
315 if (!r->sign)
316 memcpy (r, s, sizeof (struct real));
317 else if (s->sign == r->sign)
318 return do_add (r, s);
319 else
320 abort ();
321 }
322
323 void
_gst_real_add_int(struct real * r,int s)324 _gst_real_add_int (struct real *r, int s)
325 {
326 struct real s_real;
327 if (!s)
328 return;
329
330 _gst_real_from_int (&s_real, s);
331 if (!r->sign)
332 memcpy (r, &s_real, sizeof (struct real));
333 else if (s_real.sign == r->sign)
334 return do_add (r, &s_real);
335 else
336 abort ();
337 }
338
339 static void
do_mul(struct real * r,struct real * s,int lsb)340 do_mul (struct real *r, struct real *s, int lsb)
341 {
342 struct real rr;
343 unsigned short mask;
344 int n;
345
346 r->exp += s->exp;
347 r->sign *= s->sign;
348
349 rr = *r;
350 mask = SIG_MSB;
351 n = SIGSZ - 1;
352 assert (s->sig[n] & mask);
353 while (n > lsb || (s->sig[n] & (mask - 1)))
354 {
355 if (!(mask >>= 1))
356 {
357 mask = SIG_MSB;
358 n--;
359 }
360
361 /* Dividing rr by 2 matches the weight s->sig[n] & mask. Exit
362 early in case of underflow. */
363 if (!rshift_sig (&rr, &rr, 1))
364 break;
365
366 if (s->sig[n] & mask)
367 {
368 if (add_significands (r, &rr))
369 {
370 /* Lose one bit of precision to fit the carry. */
371 rshift_sig (r, r, 1);
372 r->sig[SIGSZ - 1] |= SIG_MSB;
373 r->exp++;
374 if (!rshift_sig (&rr, &rr, 1))
375 break;
376 rr.exp++;
377 }
378 }
379 }
380 }
381
382 void
_gst_real_mul(struct real * r,struct real * s)383 _gst_real_mul (struct real *r, struct real *s)
384 {
385 int i;
386 struct real tmp;
387 if (r->sign == 0)
388 return;
389 if (r == s)
390 {
391 tmp = *s;
392 s = &tmp;
393 }
394 if (s->sign == 0)
395 memset (r, 0, sizeof (struct real));
396
397 for (i = 0; i < SIGSZ && s->sig[i] == 0; )
398 i++;
399 do_mul (r, s, i);
400 }
401
402 void
_gst_real_mul_int(struct real * r,int s)403 _gst_real_mul_int (struct real *r, int s)
404 {
405 struct real s_real;
406
407 if (s == 0)
408 memset (r, 0, sizeof (struct real));
409 if (r->sign == 0)
410 return;
411
412 _gst_real_from_int (&s_real, s);
413 do_mul (r, &s_real, SIGSZ - 1);
414 }
415
416
417 void
_gst_real_powi(struct real * out,struct real * r,int s)418 _gst_real_powi (struct real *out, struct real *r, int s)
419 {
420 int k;
421 struct real tmp;
422 if (out == r)
423 {
424 tmp = *r;
425 r = &tmp;
426 }
427
428 assert (s > 0);
429 _gst_real_from_int (out, 1);
430 if (!s)
431 return;
432
433 for (k = 1;; k <<= 1)
434 {
435 if (s & k)
436 {
437 _gst_real_mul (out, r);
438 s ^= k;
439 if (!s)
440 break;
441 }
442 _gst_real_mul (r, r);
443 }
444 }
445
446 static void
do_div(struct real * out,struct real * r,struct real * s)447 do_div (struct real *out, struct real *r, struct real *s)
448 {
449 struct real v;
450 int msb, i;
451 int place = SIGSZ-1;
452 int bit = SIG_MSB;
453
454 memset (&v, 0, sizeof (struct real));
455 v.sign = r->sign * s->sign;
456 v.exp = r->exp - s->exp;
457 msb = 0;
458 goto start;
459 do
460 {
461 /* Get the MSB of U and shift it left by one. */
462 msb = r->sig[SIGSZ-1] & SIG_MSB;
463 for (i = SIGSZ; --i >= 1; )
464 r->sig[i] = (r->sig[i] << 1) | (r->sig[i - 1] >> 15);
465 r->sig[0] <<= 1;
466
467 start:
468 if (msb || cmp_significands (r, s) >= 0)
469 {
470 sub_significands (r, s);
471 v.sig[place] |= bit;
472 }
473 }
474 while ((bit >>= 1) || (bit = SIG_MSB, --place >= 0));
475 normalize (out, &v);
476 }
477
478 void
_gst_real_div(struct real * out,struct real * r,struct real * s)479 _gst_real_div (struct real *out, struct real *r, struct real *s)
480 {
481 assert (s->sign);
482 if (!r->sign)
483 {
484 memset (out, 0, sizeof (struct real));
485 return;
486 }
487
488 if (r == s)
489 {
490 memset (out, 0, sizeof (struct real));
491 out->sign = 1;
492 out->sig[SIGSZ-1] = SIG_MSB;
493 return;
494 }
495
496 if (out == r)
497 do_div (out, out, s);
498 else
499 {
500 /* do_div would destroy R, save it. */
501 struct real u = *r;
502 do_div (out, &u, s);
503 }
504 }
505
506
507 void
_gst_real_inv(struct real * out,struct real * s)508 _gst_real_inv (struct real *out, struct real *s)
509 {
510 struct real u;
511 assert (s->sign);
512 memset (&u, 0, sizeof (struct real));
513 u.sign = 1;
514 u.sig[SIGSZ-1] = SIG_MSB;
515 do_div (out, &u, s);
516 }
517
518
519 long double
_gst_real_get_ld(struct real * r)520 _gst_real_get_ld (struct real *r)
521 {
522 long double result = 0.0;
523 int i;
524
525 for (i = SIGSZ; --i >= 0; )
526 {
527 result *= SIG_MASK + 1;
528 result += r->sig[i];
529 }
530
531 result = ldexpl (result, r->exp - NUM_SIG_BITS + 1);
532 return r->sign == -1 ? -result : result;
533 }
534