pintos-os.org Git - pspp/blob - src/math/decimal.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2015 Free Software Foundation, Inc.
   3
   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.
   8
   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.
  13
  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/>. */
  16
  17
  18 #include <config.h>
  19 #include <stdio.h>
  20 #include <stdint.h>
  21 #include <string.h>
  22 #include <limits.h>
  23 #include <assert.h>
  24 #include <stdlib.h>
  25 #include <stdbool.h>
  26 #include <math.h>
  27 #include "libpspp/i18n.h"
  28
  29 #include "decimal.h"
  30
  31 int dec_warning;
  32
  33 static bool
  34 down (struct decimal *dec)
  35 {
  36   if (dec->ordinate % 10 == 0 &&  dec->mantissa < MANT_MAX - 1)
  37     {
  38       dec->ordinate /= 10;
  39       dec->mantissa++;
  40       return true;
  41     }
  42
  43   return false;
  44 }
  45
  46 static bool
  47 up (struct decimal *dec)
  48 {
  49   if (llabs (dec->ordinate) < ORD_MAX / 10   &&   dec->mantissa > MANT_MIN)
  50     {
  51       dec->ordinate *= 10;
  52       dec->mantissa--;
  53       return true;
  54     }
  55   return false;
  56 }
  57
  58
  59 /* Reduce the absolute value of the ordinate to the smallest possible,
  60    without loosing precision */
  61 static void
  62 reduce (struct decimal *dec)
  63 {
  64   if (dec->ordinate == 0)
  65     {
  66       dec->mantissa = 0;
  67       return;
  68     }
  69
  70   while (dec->ordinate % 10 == 0)
  71     {
  72       if (! down (dec))
  73         break;
  74     }
  75 }
  76
  77 /* Attempt to make the mantissas of BM and SM equal.
  78    Prerequisite: the mantissa SM must be no greater than that of BM.
  79  */
  80 static void
  81 normalisebs (struct decimal *sm, struct decimal *bm)
  82 {
  83   while (sm->mantissa < bm->mantissa)
  84     {
  85       if (down (sm))
  86         ;
  87       else if (up (bm))
  88         ;
  89       else
  90         {
  91           dec_warning = DEC_PREC;
  92           break;
  93         }
  94     }
  95
  96   while (sm->mantissa < bm->mantissa)
  97     {
  98       sm->ordinate /= 10;
  99       sm->mantissa++;
 100     }
 101 }
 102
 103
 104 /* arrange d1 and d2 such that thier mantissas are equal */
 105 void
 106 normalise (struct decimal *d1, struct decimal *d2)
 107 {
 108   normalisebs (d1, d2);
 109   normalisebs (d2, d1);
 110 }
 111
 112
 113
 114 /* Return log base 10 of D */
 115 mant_t
 116 dec_log10 (const struct decimal *d_)
 117 {
 118   struct decimal d = *d_;
 119
 120   while (llabs (d.ordinate) > 0)
 121     {
 122       d.ordinate /= 10;
 123       d.mantissa++;
 124     }
 125
 126   return d.mantissa;
 127 }
 128
 129
 130
 131 /* Return the smallest integer >= d */
 132 static ord_t
 133 decimal_ceil_pos (const struct decimal *d)
 134 {
 135   mant_t m = d->mantissa;
 136   ord_t o = d->ordinate;
 137
 138   assert (d->ordinate >= 0);
 139
 140   while (m > 0)
 141     {
 142       o *= 10;
 143       m--;
 144     }
 145
 146   while (m < 0)
 147     {
 148       bool flag = o % 10;
 149       o /= 10;
 150       if (flag)
 151         o++;
 152       m++;
 153     }
 154
 155   return o;
 156 }
 157
 158 /* Return the largest integer <= d */
 159 static ord_t
 160 decimal_floor_pos (const struct decimal *d)
 161 {
 162   mant_t m = d->mantissa;
 163   ord_t o = d->ordinate;
 164
 165   assert (d->ordinate >= 0);
 166
 167   while (m > 0)
 168     {
 169       m--;
 170       o *= 10;
 171     }
 172
 173   while (m < 0)
 174     {
 175       m++;
 176       o /= 10;
 177     }
 178
 179
 180   return o;
 181 }
 182
 183 /* Return the smallest integer which is no less than D.
 184    (round towards minus infinity) */
 185 ord_t
 186 decimal_floor (const struct decimal *d)
 187 {
 188   if (d->ordinate >= 0)
 189     return decimal_floor_pos (d);
 190   else
 191     {
 192       struct decimal dd = *d;
 193       dd.ordinate = llabs (dd.ordinate);
 194       return -decimal_ceil_pos (&dd);
 195     }
 196 }
 197
 198 /* Return the largest integer which is no greater than D.
 199    (round towards plus infinity) */
 200 ord_t
 201 decimal_ceil (const struct decimal *d)
 202 {
 203   if (d->ordinate >= 0)
 204     return decimal_ceil_pos (d);
 205   else
 206     {
 207       struct decimal dd = *d;
 208       dd.ordinate = llabs (dd.ordinate);
 209       return -decimal_floor_pos (&dd);
 210     }
 211 }
 212
 213 /* Add SRC onto DEST */
 214 void
 215 decimal_add (struct decimal *dest, const struct decimal *src_)
 216 {
 217   struct decimal src = *src_;
 218
 219   src.ordinate = -src.ordinate;
 220
 221   decimal_subtract (dest, &src);
 222 }
 223
 224 /* Subtract SRC from DEST */
 225 void
 226 decimal_subtract (struct decimal *dest, const struct decimal *src_)
 227 {
 228   struct decimal src = *src_;
 229
 230   normalise (dest, &src);
 231
 232   bool dest_neg = dest->ordinate < 0;
 233   bool src_neg = src.ordinate < 0;
 234
 235   bool expected_neg = dest_neg * src_neg;
 236
 237   if (dest->ordinate == src.ordinate)
 238     {
 239       expected_neg = 0;
 240     }
 241   else if (llabs (src.ordinate) > llabs (dest->ordinate))
 242     {
 243       if (dest_neg == src_neg)
 244         expected_neg = !expected_neg;
 245     }
 246
 247   dest->ordinate -= src.ordinate;
 248
 249   bool result_neg = dest->ordinate < 0;
 250
 251   if (expected_neg != result_neg)
 252     {
 253       /* The operation has resulted in an overflow.
 254          To resolve this, undo the operation,
 255          reduce the precision and try again */
 256
 257       dest->ordinate += src.ordinate;
 258
 259       dest->ordinate /= 10;
 260       src.ordinate /= 10;
 261
 262       dest->mantissa ++;
 263       src.mantissa ++;
 264
 265       dest->ordinate -= src.ordinate;
 266     }
 267
 268   reduce (dest);
 269
 270 }
 271
 272 /* Initialise DEC with ordinate ORD and mantissa MANT */
 273 void
 274 decimal_init (struct decimal *dec, ord_t ord, mant_t mant)
 275 {
 276   dec->ordinate = ord;
 277   dec->mantissa = mant;
 278   reduce (dec);
 279 }
 280
 281
 282 /*
 283   Compare D1 and D2.
 284
 285   Returns zero if equal, +1 if D1 > D2 and -1 if D1 < D2
 286 */
 287 int
 288 decimal_cmp (const struct decimal *d1, const struct decimal *d2)
 289 {
 290   struct decimal td1 = *d1;
 291   struct decimal td2 = *d2;
 292
 293   normalise (&td1, &td2);
 294
 295   if (td1.ordinate < td2.ordinate)
 296     return -1;
 297
 298   return (td1.ordinate > td2.ordinate);
 299 }
 300
 301
 302 /* Multiply DEST by M */
 303 void
 304 decimal_int_multiply (struct decimal *dest, ord_t m)
 305 {
 306   if (m != 0)
 307     while (llabs (dest->ordinate) > llabs (ORD_MAX / m))
 308       {
 309         dest->ordinate /= 10;
 310         dest->mantissa++;
 311       }
 312
 313   dest->ordinate *= m;
 314
 315   reduce (dest);
 316 }
 317
 318
 319 /* Divide DEST by M */
 320 void
 321 decimal_int_divide (struct decimal *dest, ord_t m)
 322 {
 323   while (dest->ordinate % m)
 324     {
 325       if (labs (dest->ordinate) > ORD_MAX / 10)
 326         {
 327           dec_warning = DEC_PREC;
 328           break;
 329         }
 330       up (dest);
 331     }
 332   dest->ordinate /= m;
 333 }
 334
 335 /* Divide DEST by SRC */
 336 void
 337 decimal_divide (struct decimal *dest, const struct decimal *src)
 338 {
 339   while (dest->ordinate % src->ordinate)
 340     {
 341       if (labs (dest->ordinate) > ORD_MAX / 10)
 342         {
 343           dec_warning = DEC_PREC;
 344           break;
 345         }
 346       up (dest);
 347     }
 348
 349   dest->ordinate /= src->ordinate;
 350   dest->mantissa -= src->mantissa;
 351 }
 352
 353 /* Print the value of DEC to F.  Probably useful only for debugging */
 354 void
 355 decimal_show (const struct decimal *dec, FILE *f)
 356 {
 357   fprintf (f, PR_ORD " x 10^" PR_MANT "\n", dec->ordinate, dec->mantissa);
 358 }
 359
 360
 361 /* Reverse the characters in string S which has length LEN */
 362 static void
 363 reverse (char *s, int len)
 364 {
 365   int i;
 366   for (i = 0; i < len / 2; ++i)
 367     {
 368       char temp = s[len - i - 1];
 369       s[len - i - 1] = s[i];
 370       s[i] = temp;
 371     }
 372 }
 373
 374 /* Return a string representation of DEC on the heap.
 375    The caller is responsible for freeing the string */
 376 char *
 377 decimal_to_string (const struct decimal *dec)
 378 {
 379   int cap = 16;
 380   int len = 0;
 381   char *s = calloc (cap, 1);
 382   ord_t ordinate = dec->ordinate;
 383
 384   while (len < dec->mantissa)
 385     {
 386       s[len++] = '0';
 387       if (len >= cap) s = realloc (s, cap <<= 1);
 388     }
 389
 390   while (ordinate)
 391     {
 392       s[len++] = labs (ordinate % 10) + '0';
 393       if (len >= cap) s = realloc (s, cap <<= 1);
 394       ordinate /= 10;
 395     }
 396
 397   if (ordinate < 0)
 398       ordinate = -ordinate;
 399
 400   while (len < -dec->mantissa)
 401     {
 402       s[len++] = '0';
 403       if (len >= cap) s = realloc (s, cap <<= 1);
 404     }
 405
 406   if (dec->mantissa < 0 )
 407     {
 408       if (len <= -dec->mantissa)
 409         {
 410           s[len++] = get_system_decimal ();
 411           if (len >= cap) s = realloc (s, cap <<= 1);
 412           s[len++] = '0';
 413           if (len >= cap) s = realloc (s, cap <<= 1);
 414         }
 415       else
 416         {
 417           int i;
 418           if (len >= cap) s = realloc (s, cap <<= 1);
 419           for (i = len - 1 ; i >= -dec->mantissa ; --i)
 420             s[i + 1] = s[i];
 421           s[i + 1] = get_system_decimal ();
 422           len++;
 423         }
 424     }
 425
 426   if (dec->ordinate < 0)
 427     {
 428       s[len++] = '-';
 429       if (len >= cap) s = realloc (s, cap <<= 1);
 430     }
 431
 432
 433   reverse (s, len);
 434
 435   {
 436     int abs_len = len;
 437     if (dec->ordinate < 0)
 438       abs_len--;
 439
 440     while (abs_len++ <= dec->mantissa)
 441       {
 442         s[len++] = '0';
 443         if (len >= cap) s = realloc (s, cap <<= 1);
 444       }
 445   }
 446
 447   return s;
 448 }
 449
 450
 451 /* Initialise DECIMAL from INPUT.
 452    INPUT should be a convential decimal representation.
 453  */
 454 void
 455 decimal_init_from_string (struct decimal *decimal, const char *input)
 456 {
 457   ord_t ordinate = 0;
 458
 459   int point = -1;
 460   int lsd = -1;
 461   int fsd = -1;
 462   int i = 0;
 463   int len = 0;
 464   int sign = 1;
 465
 466   const char *p;
 467
 468   for (p = input; *p ; p++)
 469     {
 470       if (*p == '-')
 471         {
 472           sign = -1;
 473         }
 474       else if (*p == get_system_decimal ())
 475         {
 476           assert (point == -1);
 477           point = i;
 478         }
 479       else if (*p > '0' && *p <= '9')
 480         {
 481           lsd = i;
 482           if (fsd == -1)
 483             fsd = i;
 484         }
 485       else if (*p == '0')
 486         /* ignore */
 487         ;
 488       else
 489         {
 490           fprintf (stderr, "Error: invalid character %c\n", *p);
 491           return;
 492         }
 493
 494       i++;
 495     }
 496   len = i;
 497
 498   if (point == -1)
 499     point = len;
 500
 501   mant_t mantissa = 0;
 502   if (fsd >= 0)
 503     {
 504       mant_t m = 1;
 505       for (i = lsd ; i >= fsd ; --i)
 506         {
 507           if (input[i] != get_system_decimal ())
 508             {
 509               if (ordinate > ORD_MAX - m * (input[i] - '0'))
 510                 {
 511                   fprintf (stderr, "Overflow reading value %s\n", input);
 512                   break;
 513                 }
 514               ordinate += m * (input[i] - '0');
 515               m *= 10;
 516             }
 517         }
 518
 519       if (lsd > point)
 520         mantissa = point - lsd;
 521       else
 522         mantissa = point - lsd - 1;
 523     }
 524
 525   decimal->ordinate = ordinate * sign;
 526   decimal->mantissa = mantissa;
 527 }
 528
 529
 530
 531 /* Initialise DEC from the binary fp value X */
 532 void
 533 decimal_from_double (struct decimal *dec, double x)
 534 {
 535   dec->mantissa = 0;
 536
 537   while (trunc (x) != x)
 538     {
 539       if (fabs (x) > ORD_MAX / 10.0)
 540         {
 541           dec_warning = DEC_PREC;
 542           break;
 543         }
 544       x *= 10.0;
 545       dec->mantissa--;
 546     }
 547
 548   dec->ordinate = x;
 549 }
 550
 551 /* Return a binary floating point value
 552    approximating DEC */
 553 double
 554 decimal_to_double (const struct decimal *dec)
 555 {
 556   double x = dec->ordinate;
 557   int mult = dec->mantissa;
 558
 559   while (mult < 0)
 560     {
 561       x /= 10.0;
 562       mult++;
 563     }
 564
 565   while (mult > 0)
 566     {
 567       x *= 10.0;
 568       mult--;
 569     }
 570
 571   return x;
 572 }