1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2015 Free Software Foundation, Inc.
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
27 #include "libpspp/i18n.h"
34 down (struct decimal *dec)
36 if (dec->ordinate % 10 == 0 && dec->mantissa < MANT_MAX - 1)
47 up (struct decimal *dec)
49 if (llabs (dec->ordinate) < ORD_MAX / 10 && dec->mantissa > MANT_MIN)
59 /* Reduce the absolute value of the ordinate to the smallest possible,
60 without loosing precision */
62 reduce (struct decimal *dec)
64 if (dec->ordinate == 0)
70 while (dec->ordinate % 10 == 0)
77 /* Attempt to make the mantissas of BM and SM equal.
78 Prerequisite: the mantissa SM must be no greater than that of BM.
81 normalisebs (struct decimal *sm, struct decimal *bm)
83 while (sm->mantissa < bm->mantissa)
91 dec_warning = DEC_PREC;
96 while (sm->mantissa < bm->mantissa)
104 /* arrange d1 and d2 such that thier mantissas are equal */
106 normalise (struct decimal *d1, struct decimal *d2)
108 normalisebs (d1, d2);
109 normalisebs (d2, d1);
114 /* Return log base 10 of D */
116 dec_log10 (const struct decimal *d_)
118 struct decimal d = *d_;
120 while (llabs (d.ordinate) > 0)
131 /* Return the smallest integer >= d */
133 decimal_ceil_pos (const struct decimal *d)
135 mant_t m = d->mantissa;
136 ord_t o = d->ordinate;
138 assert (d->ordinate >= 0);
158 /* Return the largest integer <= d */
160 decimal_floor_pos (const struct decimal *d)
162 mant_t m = d->mantissa;
163 ord_t o = d->ordinate;
165 assert (d->ordinate >= 0);
183 /* Return the smallest integer which is no less than D.
184 (round towards minus infinity) */
186 decimal_floor (const struct decimal *d)
188 if (d->ordinate >= 0)
189 return decimal_floor_pos (d);
192 struct decimal dd = *d;
193 dd.ordinate = llabs (dd.ordinate);
194 return -decimal_ceil_pos (&dd);
198 /* Return the largest integer which is no greater than D.
199 (round towards plus infinity) */
201 decimal_ceil (const struct decimal *d)
203 if (d->ordinate >= 0)
204 return decimal_ceil_pos (d);
207 struct decimal dd = *d;
208 dd.ordinate = llabs (dd.ordinate);
209 return -decimal_floor_pos (&dd);
213 /* Add SRC onto DEST */
215 decimal_add (struct decimal *dest, const struct decimal *src_)
217 struct decimal src = *src_;
219 src.ordinate = -src.ordinate;
221 decimal_subtract (dest, &src);
224 /* Subtract SRC from DEST */
226 decimal_subtract (struct decimal *dest, const struct decimal *src_)
228 struct decimal src = *src_;
230 normalise (dest, &src);
232 bool dest_neg = dest->ordinate < 0;
233 bool src_neg = src.ordinate < 0;
235 bool expected_neg = dest_neg * src_neg;
237 if (dest->ordinate == src.ordinate)
241 else if (llabs (src.ordinate) > llabs (dest->ordinate))
243 if (dest_neg == src_neg)
244 expected_neg = !expected_neg;
247 dest->ordinate -= src.ordinate;
249 bool result_neg = dest->ordinate < 0;
251 if (expected_neg != result_neg)
253 /* The operation has resulted in an overflow.
254 To resolve this, undo the operation,
255 reduce the precision and try again */
257 dest->ordinate += src.ordinate;
259 dest->ordinate /= 10;
265 dest->ordinate -= src.ordinate;
272 /* Initialise DEC with ordinate ORD and mantissa MANT */
274 decimal_init (struct decimal *dec, ord_t ord, mant_t mant)
277 dec->mantissa = mant;
285 Returns zero if equal, +1 if D1 > D2 and -1 if D1 < D2
288 decimal_cmp (const struct decimal *d1, const struct decimal *d2)
290 struct decimal td1 = *d1;
291 struct decimal td2 = *d2;
293 normalise (&td1, &td2);
295 if (td1.ordinate < td2.ordinate)
298 return (td1.ordinate > td2.ordinate);
302 /* Multiply DEST by M */
304 decimal_int_multiply (struct decimal *dest, ord_t m)
307 while (llabs (dest->ordinate) > llabs (ORD_MAX / m))
309 dest->ordinate /= 10;
319 /* Divide DEST by M */
321 decimal_int_divide (struct decimal *dest, ord_t m)
323 while (dest->ordinate % m)
325 if (labs (dest->ordinate) > ORD_MAX / 10)
327 dec_warning = DEC_PREC;
335 /* Divide DEST by SRC */
337 decimal_divide (struct decimal *dest, const struct decimal *src)
339 while (dest->ordinate % src->ordinate)
341 if (labs (dest->ordinate) > ORD_MAX / 10)
343 dec_warning = DEC_PREC;
349 dest->ordinate /= src->ordinate;
350 dest->mantissa -= src->mantissa;
353 /* Print the value of DEC to F. Probably useful only for debugging */
355 decimal_show (const struct decimal *dec, FILE *f)
357 fprintf (f, PR_ORD " x 10^" PR_MANT "\n", dec->ordinate, dec->mantissa);
361 /* Reverse the characters in string S which has length LEN */
363 reverse (char *s, int len)
366 for (i = 0; i < len / 2; ++i)
368 char temp = s[len - i - 1];
369 s[len - i - 1] = s[i];
374 /* Return a string representation of DEC on the heap.
375 The caller is responsible for freeing the string */
377 decimal_to_string (const struct decimal *dec)
381 char *s = calloc (cap, 1);
382 ord_t ordinate = dec->ordinate;
384 while (len < dec->mantissa)
387 if (len >= cap) s = realloc (s, cap <<= 1);
392 s[len++] = labs (ordinate % 10) + '0';
393 if (len >= cap) s = realloc (s, cap <<= 1);
398 ordinate = -ordinate;
400 while (len < -dec->mantissa)
403 if (len >= cap) s = realloc (s, cap <<= 1);
406 if (dec->mantissa < 0 )
408 if (len <= -dec->mantissa)
410 s[len++] = get_system_decimal ();
411 if (len >= cap) s = realloc (s, cap <<= 1);
413 if (len >= cap) s = realloc (s, cap <<= 1);
418 if (len >= cap) s = realloc (s, cap <<= 1);
419 for (i = len - 1 ; i >= -dec->mantissa ; --i)
421 s[i + 1] = get_system_decimal ();
426 if (dec->ordinate < 0)
429 if (len >= cap) s = realloc (s, cap <<= 1);
437 if (dec->ordinate < 0)
440 while (abs_len++ <= dec->mantissa)
443 if (len >= cap) s = realloc (s, cap <<= 1);
451 /* Initialise DECIMAL from INPUT.
452 INPUT should be a convential decimal representation.
455 decimal_init_from_string (struct decimal *decimal, const char *input)
468 for (p = input; *p ; p++)
474 else if (*p == get_system_decimal ())
476 assert (point == -1);
479 else if (*p > '0' && *p <= '9')
490 fprintf (stderr, "Error: invalid character %c\n", *p);
505 for (i = lsd ; i >= fsd ; --i)
507 if (input[i] != get_system_decimal ())
509 if (ordinate > ORD_MAX - m * (input[i] - '0'))
511 fprintf (stderr, "Overflow reading value %s\n", input);
514 ordinate += m * (input[i] - '0');
520 mantissa = point - lsd;
522 mantissa = point - lsd - 1;
525 decimal->ordinate = ordinate * sign;
526 decimal->mantissa = mantissa;
531 /* Initialise DEC from the binary fp value X */
533 decimal_from_double (struct decimal *dec, double x)
537 while (trunc (x) != x)
539 if (fabs (x) > ORD_MAX / 10.0)
541 dec_warning = DEC_PREC;
551 /* Return a binary floating point value
554 decimal_to_double (const struct decimal *dec)
556 double x = dec->ordinate;
557 int mult = dec->mantissa;