+/* PSPP - a program for statistical analysis.
+ Copyright (C) 2015 Free Software Foundation, Inc.
+
+ This program 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 3 of the License, or
+ (at your option) any later version.
+
+ This program 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 this program. If not, see <http://www.gnu.org/licenses/>. */
+
+
+#include <config.h>
+#include <stdio.h>
+#include <stdint.h>
+#include <string.h>
+#include <limits.h>
+#include <assert.h>
+#include <stdlib.h>
+#include <stdbool.h>
+#include <math.h>
+#include "libpspp/i18n.h"
+
+#include "decimal.h"
+
+int dec_warning;
+
+static bool
+down (struct decimal *dec)
+{
+ if (dec->ordinate % 10 == 0 && dec->mantissa < MANT_MAX - 1)
+ {
+ dec->ordinate /= 10;
+ dec->mantissa++;
+ return true;
+ }
+
+ return false;
+}
+
+static bool
+up (struct decimal *dec)
+{
+ if (llabs (dec->ordinate) < ORD_MAX / 10 && dec->mantissa > MANT_MIN)
+ {
+ dec->ordinate *= 10;
+ dec->mantissa--;
+ return true;
+ }
+ return false;
+}
+
+
+/* Reduce the absolute value of the ordinate to the smallest possible,
+ without loosing precision */
+static void
+reduce (struct decimal *dec)
+{
+ if (dec->ordinate == 0)
+ {
+ dec->mantissa = 0;
+ return;
+ }
+
+ while (dec->ordinate % 10 == 0)
+ {
+ if (! down (dec))
+ break;
+ }
+}
+
+/* Attempt to make the mantissas of BM and SM equal.
+ Prerequisite: the mantissa SM must be no greater than that of BM.
+ */
+static void
+normalisebs (struct decimal *sm, struct decimal *bm)
+{
+ while (sm->mantissa < bm->mantissa)
+ {
+ if (down (sm))
+ ;
+ else if (up (bm))
+ ;
+ else
+ {
+ dec_warning = DEC_PREC;
+ break;
+ }
+ }
+
+ while (sm->mantissa < bm->mantissa)
+ {
+ sm->ordinate /= 10;
+ sm->mantissa++;
+ }
+}
+
+
+/* arrange d1 and d2 such that thier mantissas are equal */
+void
+normalise (struct decimal *d1, struct decimal *d2)
+{
+ normalisebs (d1, d2);
+ normalisebs (d2, d1);
+}
+
+
+
+/* Return log base 10 of D */
+mant_t
+dec_log10 (const struct decimal *d_)
+{
+ struct decimal d = *d_;
+
+ while (llabs (d.ordinate) > 0)
+ {
+ d.ordinate /= 10;
+ d.mantissa++;
+ }
+
+ return d.mantissa;
+}
+
+
+
+/* Return the smallest integer >= d */
+static ord_t
+decimal_ceil_pos (const struct decimal *d)
+{
+ mant_t m = d->mantissa;
+ ord_t o = d->ordinate;
+
+ assert (d->ordinate >= 0);
+
+ while (m > 0)
+ {
+ o *= 10;
+ m--;
+ }
+
+ while (m < 0)
+ {
+ bool flag = o % 10;
+ o /= 10;
+ if (flag)
+ o++;
+ m++;
+ }
+
+ return o;
+}
+
+/* Return the largest integer <= d */
+static ord_t
+decimal_floor_pos (const struct decimal *d)
+{
+ mant_t m = d->mantissa;
+ ord_t o = d->ordinate;
+
+ assert (d->ordinate >= 0);
+
+ while (m > 0)
+ {
+ m--;
+ o *= 10;
+ }
+
+ while (m < 0)
+ {
+ m++;
+ o /= 10;
+ }
+
+
+ return o;
+}
+
+/* Return the smallest integer which is no less than D.
+ (round towards minus infinity) */
+ord_t
+decimal_floor (const struct decimal *d)
+{
+ if (d->ordinate >= 0)
+ return decimal_floor_pos (d);
+ else
+ {
+ struct decimal dd = *d;
+ dd.ordinate = llabs (dd.ordinate);
+ return -decimal_ceil_pos (&dd);
+ }
+}
+
+/* Return the largest integer which is no greater than D.
+ (round towards plus infinity) */
+ord_t
+decimal_ceil (const struct decimal *d)
+{
+ if (d->ordinate >= 0)
+ return decimal_ceil_pos (d);
+ else
+ {
+ struct decimal dd = *d;
+ dd.ordinate = llabs (dd.ordinate);
+ return -decimal_floor_pos (&dd);
+ }
+}
+
+/* Add SRC onto DEST */
+void
+decimal_add (struct decimal *dest, const struct decimal *src_)
+{
+ struct decimal src = *src_;
+
+ src.ordinate = -src.ordinate;
+
+ decimal_subtract (dest, &src);
+}
+
+/* Subtract SRC from DEST */
+void
+decimal_subtract (struct decimal *dest, const struct decimal *src_)
+{
+ struct decimal src = *src_;
+
+ normalise (dest, &src);
+
+ bool dest_neg = dest->ordinate < 0;
+ bool src_neg = src.ordinate < 0;
+
+ bool expected_neg = dest_neg * src_neg;
+
+ if (dest->ordinate == src.ordinate)
+ {
+ expected_neg = 0;
+ }
+ else if (llabs (src.ordinate) > llabs (dest->ordinate))
+ {
+ if (dest_neg == src_neg)
+ expected_neg = !expected_neg;
+ }
+
+ dest->ordinate -= src.ordinate;
+
+ bool result_neg = dest->ordinate < 0;
+
+ if (expected_neg != result_neg)
+ {
+ /* The operation has resulted in an overflow.
+ To resolve this, undo the operation,
+ reduce the precision and try again */
+
+ dest->ordinate += src.ordinate;
+
+ dest->ordinate /= 10;
+ src.ordinate /= 10;
+
+ dest->mantissa ++;
+ src.mantissa ++;
+
+ dest->ordinate -= src.ordinate;
+ }
+
+ reduce (dest);
+
+}
+
+/* Initialise DEC with ordinate ORD and mantissa MANT */
+void
+decimal_init (struct decimal *dec, ord_t ord, mant_t mant)
+{
+ dec->ordinate = ord;
+ dec->mantissa = mant;
+ reduce (dec);
+}
+
+
+/*
+ Compare D1 and D2.
+
+ Returns zero if equal, +1 if D1 > D2 and -1 if D1 < D2
+*/
+int
+decimal_cmp (const struct decimal *d1, const struct decimal *d2)
+{
+ struct decimal td1 = *d1;
+ struct decimal td2 = *d2;
+
+ normalise (&td1, &td2);
+
+ if (td1.ordinate < td2.ordinate)
+ return -1;
+
+ return (td1.ordinate > td2.ordinate);
+}
+
+
+/* Multiply DEST by M */
+void
+decimal_int_multiply (struct decimal *dest, ord_t m)
+{
+ if (m != 0)
+ while (llabs (dest->ordinate) > llabs (ORD_MAX / m))
+ {
+ dest->ordinate /= 10;
+ dest->mantissa++;
+ }
+
+ dest->ordinate *= m;
+
+ reduce (dest);
+}
+
+
+/* Divide DEST by M */
+void
+decimal_int_divide (struct decimal *dest, ord_t m)
+{
+ while (dest->ordinate % m)
+ {
+ if (labs (dest->ordinate) > ORD_MAX / 10)
+ {
+ dec_warning = DEC_PREC;
+ break;
+ }
+ up (dest);
+ }
+ dest->ordinate /= m;
+}
+
+/* Divide DEST by SRC */
+void
+decimal_divide (struct decimal *dest, const struct decimal *src)
+{
+ while (dest->ordinate % src->ordinate)
+ {
+ if (labs (dest->ordinate) > ORD_MAX / 10)
+ {
+ dec_warning = DEC_PREC;
+ break;
+ }
+ up (dest);
+ }
+
+ dest->ordinate /= src->ordinate;
+ dest->mantissa -= src->mantissa;
+}
+
+/* Print the value of DEC to F. Probably useful only for debugging */
+void
+decimal_show (const struct decimal *dec, FILE *f)
+{
+ fprintf (f, PR_ORD " x 10^" PR_MANT "\n", dec->ordinate, dec->mantissa);
+}
+
+
+/* Reverse the characters in string S which has length LEN */
+static void
+reverse (char *s, int len)
+{
+ int i;
+ for (i = 0; i < len / 2; ++i)
+ {
+ char temp = s[len - i - 1];
+ s[len - i - 1] = s[i];
+ s[i] = temp;
+ }
+}
+
+/* Return a string representation of DEC on the heap.
+ The caller is responsible for freeing the string */
+char *
+decimal_to_string (const struct decimal *dec)
+{
+ int cap = 16;
+ int len = 0;
+ char *s = calloc (cap, 1);
+ ord_t ordinate = dec->ordinate;
+
+ while (len < dec->mantissa)
+ {
+ s[len++] = '0';
+ if (len >= cap) s = realloc (s, cap <<= 1);
+ }
+
+ while (ordinate)
+ {
+ s[len++] = labs (ordinate % 10) + '0';
+ if (len >= cap) s = realloc (s, cap <<= 1);
+ ordinate /= 10;
+ }
+
+ if (ordinate < 0)
+ ordinate = -ordinate;
+
+ while (len < -dec->mantissa)
+ {
+ s[len++] = '0';
+ if (len >= cap) s = realloc (s, cap <<= 1);
+ }
+
+ if (dec->mantissa < 0 )
+ {
+ if (len <= -dec->mantissa)
+ {
+ s[len++] = get_system_decimal ();
+ if (len >= cap) s = realloc (s, cap <<= 1);
+ s[len++] = '0';
+ if (len >= cap) s = realloc (s, cap <<= 1);
+ }
+ else
+ {
+ int i;
+ if (len >= cap) s = realloc (s, cap <<= 1);
+ for (i = len - 1 ; i >= -dec->mantissa ; --i)
+ s[i + 1] = s[i];
+ s[i + 1] = get_system_decimal ();
+ len++;
+ }
+ }
+
+ if (dec->ordinate < 0)
+ {
+ s[len++] = '-';
+ if (len >= cap) s = realloc (s, cap <<= 1);
+ }
+
+
+ reverse (s, len);
+
+ {
+ int abs_len = len;
+ if (dec->ordinate < 0)
+ abs_len--;
+
+ while (abs_len++ <= dec->mantissa)
+ {
+ s[len++] = '0';
+ if (len >= cap) s = realloc (s, cap <<= 1);
+ }
+ }
+
+ return s;
+}
+
+
+/* Initialise DECIMAL from INPUT.
+ INPUT should be a convential decimal representation.
+ */
+void
+decimal_init_from_string (struct decimal *decimal, const char *input)
+{
+ ord_t ordinate = 0;
+
+ int point = -1;
+ int lsd = -1;
+ int fsd = -1;
+ int i = 0;
+ int len = 0;
+ int sign = 1;
+
+ const char *p;
+
+ for (p = input; *p ; p++)
+ {
+ if (*p == '-')
+ {
+ sign = -1;
+ }
+ else if (*p == get_system_decimal ())
+ {
+ assert (point == -1);
+ point = i;
+ }
+ else if (*p > '0' && *p <= '9')
+ {
+ lsd = i;
+ if (fsd == -1)
+ fsd = i;
+ }
+ else if (*p == '0')
+ /* ignore */
+ ;
+ else
+ {
+ fprintf (stderr, "Error: invalid character %c\n", *p);
+ return;
+ }
+
+ i++;
+ }
+ len = i;
+
+ if (point == -1)
+ point = len;
+
+ mant_t mantissa = 0;
+ if (fsd >= 0)
+ {
+ mant_t m = 1;
+ for (i = lsd ; i >= fsd ; --i)
+ {
+ if (input[i] != get_system_decimal ())
+ {
+ if (ordinate > ORD_MAX - m * (input[i] - '0'))
+ {
+ fprintf (stderr, "Overflow reading value %s\n", input);
+ break;
+ }
+ ordinate += m * (input[i] - '0');
+ m *= 10;
+ }
+ }
+
+ if (lsd > point)
+ mantissa = point - lsd;
+ else
+ mantissa = point - lsd - 1;
+ }
+
+ decimal->ordinate = ordinate * sign;
+ decimal->mantissa = mantissa;
+}
+
+
+
+/* Initialise DEC from the binary fp value X */
+void
+decimal_from_double (struct decimal *dec, double x)
+{
+ dec->mantissa = 0;
+
+ while (trunc (x) != x)
+ {
+ if (fabs (x) > ORD_MAX / 10.0)
+ {
+ dec_warning = DEC_PREC;
+ break;
+ }
+ x *= 10.0;
+ dec->mantissa--;
+ }
+
+ dec->ordinate = x;
+}
+
+/* Return a binary floating point value
+ approximating DEC */
+double
+decimal_to_double (const struct decimal *dec)
+{
+ double x = dec->ordinate;
+ int mult = dec->mantissa;
+
+ while (mult < 0)
+ {
+ x /= 10.0;
+ mult++;
+ }
+
+ while (mult > 0)
+ {
+ x *= 10.0;
+ mult--;
+ }
+
+ return x;
+}