gcc/sreal.c

   1 /* Simple data type for positive real numbers for the GNU compiler.
   2    Copyright (C) 2002, 2003 Free Software Foundation, Inc.
   3
   4 This file is part of GCC.
   5
   6 GCC is free software; you can redistribute it and/or modify it under
   7 the terms of the GNU General Public License as published by the Free
   8 Software Foundation; either version 2, or (at your option) any later
   9 version.
  10
  11 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  12 WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  14 for more details.
  15
  16 You should have received a copy of the GNU General Public License
  17 along with GCC; see the file COPYING.  If not, write to the Free
  18 Software Foundation, 59 Temple Place - Suite 330, Boston, MA
  19 02111-1307, USA.  */
  20
  21 /* This library supports positive real numbers and 0;
  22    inf and nan are NOT supported.
  23    It is written to be simple and fast.
  24
  25    Value of sreal is
  26         x = sig * 2 ^ exp
  27    where
  28         sig = significant
  29           (for < 64-bit machines sig = sig_lo + sig_hi * 2 ^ SREAL_PART_BITS)
  30         exp = exponent
  31
  32    One HOST_WIDE_INT is used for the significant on 64-bit (and more than
  33    64-bit) machines,
  34    otherwise two HOST_WIDE_INTs are used for the significant.
  35    Only a half of significant bits is used (in normalized sreals) so that we do
  36    not have problems with overflow, for example when c->sig = a->sig * b->sig.
  37    So the precission for 64-bit and 32-bit machines is 32-bit.
  38
  39    Invariant: The numbers are normalized before and after each call of sreal_*.
  40
  41    Normalized sreals:
  42    All numbers (except zero) meet following conditions:
  43          SREAL_MIN_SIG <= sig && sig <= SREAL_MAX_SIG
  44         -SREAL_MAX_EXP <= exp && exp <= SREAL_MAX_EXP
  45
  46    If the number would be too large, it is set to upper bounds of these
  47    conditions.
  48
  49    If the number is zero or would be too small it meets following conditions:
  50         sig == 0 && exp == -SREAL_MAX_EXP
  51 */
  52
  53 #include "config.h"
  54 #include "system.h"
  55 #include "coretypes.h"
  56 #include "tm.h"
  57 #include "sreal.h"
  58
  59 static inline void copy         PARAMS ((sreal *, sreal *));
  60 static inline void shift_right  PARAMS ((sreal *, int));
  61 static void normalize           PARAMS ((sreal *));
  62
  63 /* Print the content of struct sreal.  */
  64
  65 void
  66 dump_sreal (file, x)
  67      FILE *file;
  68      sreal *x;
  69 {
  70 #if SREAL_PART_BITS < 32
  71   fprintf (file, "((" HOST_WIDE_INT_PRINT_UNSIGNED " * 2^16 + "
  72            HOST_WIDE_INT_PRINT_UNSIGNED ") * 2^%d)",
  73            x->sig_hi, x->sig_lo, x->exp);
  74 #else
  75   fprintf (file, "(" HOST_WIDE_INT_PRINT_UNSIGNED " * 2^%d)", x->sig, x->exp);
  76 #endif
  77 }
  78
  79 /* Copy the sreal number.  */
  80
  81 static inline void
  82 copy (r, a)
  83      sreal *r;
  84      sreal *a;
  85 {
  86 #if SREAL_PART_BITS < 32
  87   r->sig_lo = a->sig_lo;
  88   r->sig_hi = a->sig_hi;
  89 #else
  90   r->sig = a->sig;
  91 #endif
  92   r->exp = a->exp;
  93 }
  94
  95 /* Shift X right by S bits.  Needed: 0 < S <= SREAL_BITS.
  96    When the most significant bit shifted out is 1, add 1 to X (rounding).  */
  97
  98 static inline void
  99 shift_right (x, s)
 100      sreal *x;
 101      int s;
 102 {
 103 #ifdef ENABLE_CHECKING
 104   if (s <= 0 || s > SREAL_BITS)
 105     abort ();
 106   if (x->exp + s > SREAL_MAX_EXP)
 107     {
 108       /* Exponent should never be so large because shift_right is used only by
 109          sreal_add and sreal_sub ant thus the number cannot be shifted out from
 110          exponent range.  */
 111       abort ();
 112     }
 113 #endif
 114
 115   x->exp += s;
 116
 117 #if SREAL_PART_BITS < 32
 118   if (s > SREAL_PART_BITS)
 119     {
 120       s -= SREAL_PART_BITS;
 121       x->sig_hi += (uhwi) 1 << (s - 1);
 122       x->sig_lo = x->sig_hi >> s;
 123       x->sig_hi = 0;
 124     }
 125   else
 126     {
 127       x->sig_lo += (uhwi) 1 << (s - 1);
 128       if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
 129         {
 130           x->sig_hi++;
 131           x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
 132         }
 133       x->sig_lo >>= s;
 134       x->sig_lo |= (x->sig_hi & (((uhwi) 1 << s) - 1)) << (SREAL_PART_BITS - s);
 135       x->sig_hi >>= s;
 136     }
 137 #else
 138   x->sig += (uhwi) 1 << (s - 1);
 139   x->sig >>= s;
 140 #endif
 141 }
 142
 143 /* Normalize *X.  */
 144
 145 static void
 146 normalize (x)
 147      sreal *x;
 148 {
 149 #if SREAL_PART_BITS < 32
 150   int shift;
 151   HOST_WIDE_INT mask;
 152
 153   if (x->sig_lo == 0 && x->sig_hi == 0)
 154     {
 155       x->exp = -SREAL_MAX_EXP;
 156     }
 157   else if (x->sig_hi < SREAL_MIN_SIG)
 158     {
 159       if (x->sig_hi == 0)
 160         {
 161           /* Move lower part of significant to higher part.  */
 162           x->sig_hi = x->sig_lo;
 163           x->sig_lo = 0;
 164           x->exp -= SREAL_PART_BITS;
 165         }
 166       shift = 0;
 167       while (x->sig_hi < SREAL_MIN_SIG)
 168         {
 169           x->sig_hi <<= 1;
 170           x->exp--;
 171           shift++;
 172         }
 173       /* Check underflow.  */
 174       if (x->exp < -SREAL_MAX_EXP)
 175         {
 176           x->exp = -SREAL_MAX_EXP;
 177           x->sig_hi = 0;
 178           x->sig_lo = 0;
 179         }
 180       else if (shift)
 181         {
 182           mask = (1 << SREAL_PART_BITS) - (1 << (SREAL_PART_BITS - shift));
 183           x->sig_hi |= (x->sig_lo & mask) >> (SREAL_PART_BITS - shift);
 184           x->sig_lo = (x->sig_lo << shift) & (((uhwi) 1 << SREAL_PART_BITS) - 1);
 185         }
 186     }
 187   else if (x->sig_hi > SREAL_MAX_SIG)
 188     {
 189       unsigned HOST_WIDE_INT tmp = x->sig_hi;
 190
 191       /* Find out how many bits will be shifted.  */
 192       shift = 0;
 193       do
 194         {
 195           tmp >>= 1;
 196           shift++;
 197         }
 198       while (tmp > SREAL_MAX_SIG);
 199
 200       /* Round the number.  */
 201       x->sig_lo += (uhwi) 1 << (shift - 1);
 202
 203       x->sig_lo >>= shift;
 204       x->sig_lo += ((x->sig_hi & (((uhwi) 1 << shift) - 1))
 205                     << (SREAL_PART_BITS - shift));
 206       x->sig_hi >>= shift;
 207       x->exp += shift;
 208       if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
 209         {
 210           x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
 211           x->sig_hi++;
 212           if (x->sig_hi > SREAL_MAX_SIG)
 213             {
 214               /* x->sig_hi was SREAL_MAX_SIG before increment
 215                  so now last bit is zero.  */
 216               x->sig_hi >>= 1;
 217               x->sig_lo >>= 1;
 218               x->exp++;
 219             }
 220         }
 221
 222       /* Check overflow.  */
 223       if (x->exp > SREAL_MAX_EXP)
 224         {
 225           x->exp = SREAL_MAX_EXP;
 226           x->sig_hi = SREAL_MAX_SIG;
 227           x->sig_lo = SREAL_MAX_SIG;
 228         }
 229     }
 230 #else
 231   if (x->sig == 0)
 232     {
 233       x->exp = -SREAL_MAX_EXP;
 234     }
 235   else if (x->sig < SREAL_MIN_SIG)
 236     {
 237       do
 238         {
 239           x->sig <<= 1;
 240           x->exp--;
 241         }
 242       while (x->sig < SREAL_MIN_SIG);
 243
 244       /* Check underflow.  */
 245       if (x->exp < -SREAL_MAX_EXP)
 246         {
 247           x->exp = -SREAL_MAX_EXP;
 248           x->sig = 0;
 249         }
 250     }
 251   else if (x->sig > SREAL_MAX_SIG)
 252     {
 253       int last_bit;
 254       do
 255         {
 256           last_bit = x->sig & 1;
 257           x->sig >>= 1;
 258           x->exp++;
 259         }
 260       while (x->sig > SREAL_MAX_SIG);
 261
 262       /* Round the number.  */
 263       x->sig += last_bit;
 264       if (x->sig > SREAL_MAX_SIG)
 265         {
 266           x->sig >>= 1;
 267           x->exp++;
 268         }
 269
 270       /* Check overflow.  */
 271       if (x->exp > SREAL_MAX_EXP)
 272         {
 273           x->exp = SREAL_MAX_EXP;
 274           x->sig = SREAL_MAX_SIG;
 275         }
 276     }
 277 #endif
 278 }
 279
 280 /* Set *R to SIG * 2 ^ EXP.  Return R.  */
 281
 282 sreal *
 283 sreal_init (r, sig, exp)
 284      sreal *r;
 285      unsigned HOST_WIDE_INT sig;
 286      signed int exp;
 287 {
 288 #if SREAL_PART_BITS < 32
 289   r->sig_lo = 0;
 290   r->sig_hi = sig;
 291   r->exp = exp - 16;
 292 #else
 293   r->sig = sig;
 294   r->exp = exp;
 295 #endif
 296   normalize (r);
 297   return r;
 298 }
 299
 300 /* Return integer value of *R.  */
 301
 302 HOST_WIDE_INT
 303 sreal_to_int (r)
 304      sreal *r;
 305 {
 306 #if SREAL_PART_BITS < 32
 307   if (r->exp <= -SREAL_BITS)
 308     return 0;
 309   if (r->exp >= 0)
 310     return MAX_HOST_WIDE_INT;
 311   return ((r->sig_hi << SREAL_PART_BITS) + r->sig_lo) >> -r->exp;
 312 #else
 313   if (r->exp <= -SREAL_BITS)
 314     return 0;
 315   if (r->exp >= SREAL_PART_BITS)
 316     return MAX_HOST_WIDE_INT;
 317   if (r->exp > 0)
 318     return r->sig << r->exp;
 319   if (r->exp < 0)
 320     return r->sig >> -r->exp;
 321   return r->sig;
 322 #endif
 323 }
 324
 325 /* Compare *A and *B. Return -1 if *A < *B, 1 if *A > *B and 0 if *A == *B.  */
 326
 327 int
 328 sreal_compare (a, b)
 329      sreal *a;
 330      sreal *b;
 331 {
 332   if (a->exp > b->exp)
 333     return 1;
 334   if (a->exp < b->exp)
 335     return -1;
 336 #if SREAL_PART_BITS < 32
 337   if (a->sig_hi > b->sig_hi)
 338     return 1;
 339   if (a->sig_hi < b->sig_hi)
 340     return -1;
 341   if (a->sig_lo > b->sig_lo)
 342     return 1;
 343   if (a->sig_lo < b->sig_lo)
 344     return -1;
 345 #else
 346   if (a->sig > b->sig)
 347     return 1;
 348   if (a->sig < b->sig)
 349     return -1;
 350 #endif
 351   return 0;
 352 }
 353
 354 /* *R = *A + *B.  Return R.  */
 355
 356 sreal *
 357 sreal_add (r, a, b)
 358   sreal *r;
 359   sreal *a;
 360   sreal *b;
 361 {
 362   int dexp;
 363   sreal tmp;
 364   sreal *bb;
 365
 366   if (sreal_compare (a, b) < 0)
 367     {
 368       sreal *swap;
 369       swap = a;
 370       a = b;
 371       b = swap;
 372     }
 373
 374   dexp = a->exp - b->exp;
 375   r->exp = a->exp;
 376   if (dexp > SREAL_BITS)
 377     {
 378 #if SREAL_PART_BITS < 32
 379       r->sig_hi = a->sig_hi;
 380       r->sig_lo = a->sig_lo;
 381 #else
 382       r->sig = a->sig;
 383 #endif
 384       return r;
 385     }
 386
 387   if (dexp == 0)
 388     bb = b;
 389   else
 390     {
 391       copy (&tmp, b);
 392       shift_right (&tmp, dexp);
 393       bb = &tmp;
 394     }
 395
 396 #if SREAL_PART_BITS < 32
 397   r->sig_hi = a->sig_hi + bb->sig_hi;
 398   r->sig_lo = a->sig_lo + bb->sig_lo;
 399   if (r->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
 400     {
 401       r->sig_hi++;
 402       r->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
 403     }
 404 #else
 405   r->sig = a->sig + bb->sig;
 406 #endif
 407   normalize (r);
 408   return r;
 409 }
 410
 411 /* *R = *A - *B.  Return R.  */
 412
 413 sreal *
 414 sreal_sub (r, a, b)
 415   sreal *r;
 416   sreal *a;
 417   sreal *b;
 418 {
 419   int dexp;
 420   sreal tmp;
 421   sreal *bb;
 422
 423   if (sreal_compare (a, b) < 0)
 424     {
 425       abort ();
 426     }
 427
 428   dexp = a->exp - b->exp;
 429   r->exp = a->exp;
 430   if (dexp > SREAL_BITS)
 431     {
 432 #if SREAL_PART_BITS < 32
 433       r->sig_hi = a->sig_hi;
 434       r->sig_lo = a->sig_lo;
 435 #else
 436       r->sig = a->sig;
 437 #endif
 438       return r;
 439     }
 440   if (dexp == 0)
 441     bb = b;
 442   else
 443     {
 444       copy (&tmp, b);
 445       shift_right (&tmp, dexp);
 446       bb = &tmp;
 447     }
 448
 449 #if SREAL_PART_BITS < 32
 450   if (a->sig_lo < bb->sig_lo)
 451     {
 452       r->sig_hi = a->sig_hi - bb->sig_hi - 1;
 453       r->sig_lo = a->sig_lo + ((uhwi) 1 << SREAL_PART_BITS) - bb->sig_lo;
 454     }
 455   else
 456     {
 457       r->sig_hi = a->sig_hi - bb->sig_hi;
 458       r->sig_lo = a->sig_lo - bb->sig_lo;
 459     }
 460 #else
 461   r->sig = a->sig - bb->sig;
 462 #endif
 463   normalize (r);
 464   return r;
 465 }
 466
 467 /* *R = *A * *B.  Return R.  */
 468
 469 sreal *
 470 sreal_mul (r, a, b)
 471      sreal *r;
 472      sreal *a;
 473      sreal *b;
 474 {
 475 #if SREAL_PART_BITS < 32
 476   if (a->sig_hi < SREAL_MIN_SIG || b->sig_hi < SREAL_MIN_SIG)
 477     {
 478       r->sig_lo = 0;
 479       r->sig_hi = 0;
 480       r->exp = -SREAL_MAX_EXP;
 481     }
 482   else
 483     {
 484       unsigned HOST_WIDE_INT tmp1, tmp2, tmp3;
 485       if (sreal_compare (a, b) < 0)
 486         {
 487           sreal *swap;
 488           swap = a;
 489           a = b;
 490           b = swap;
 491         }
 492
 493       r->exp = a->exp + b->exp + SREAL_PART_BITS;
 494
 495       tmp1 = a->sig_lo * b->sig_lo;
 496       tmp2 = a->sig_lo * b->sig_hi;
 497       tmp3 = a->sig_hi * b->sig_lo + (tmp1 >> SREAL_PART_BITS);
 498
 499       r->sig_hi = a->sig_hi * b->sig_hi;
 500       r->sig_hi += (tmp2 >> SREAL_PART_BITS) + (tmp3 >> SREAL_PART_BITS);
 501       tmp2 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;
 502       tmp3 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;
 503       tmp1 = tmp2 + tmp3;
 504
 505       r->sig_lo = tmp1 & (((uhwi) 1 << SREAL_PART_BITS) - 1);
 506       r->sig_hi += tmp1 >> SREAL_PART_BITS;
 507
 508       normalize (r);
 509     }
 510 #else
 511   if (a->sig < SREAL_MIN_SIG || b->sig < SREAL_MIN_SIG)
 512     {
 513       r->sig = 0;
 514       r->exp = -SREAL_MAX_EXP;
 515     }
 516   else
 517     {
 518       r->sig = a->sig * b->sig;
 519       r->exp = a->exp + b->exp;
 520       normalize (r);
 521     }
 522 #endif
 523   return r;
 524 }
 525
 526 /* *R = *A / *B.  Return R.  */
 527
 528 sreal *
 529 sreal_div (r, a, b)
 530      sreal *r;
 531      sreal *a;
 532      sreal *b;
 533 {
 534 #if SREAL_PART_BITS < 32
 535   unsigned HOST_WIDE_INT tmp, tmp1, tmp2;
 536
 537   if (b->sig_hi < SREAL_MIN_SIG)
 538     {
 539       abort ();
 540     }
 541   else if (a->sig_hi < SREAL_MIN_SIG)
 542     {
 543       r->sig_hi = 0;
 544       r->sig_lo = 0;
 545       r->exp = -SREAL_MAX_EXP;
 546     }
 547   else
 548     {
 549       /* Since division by the whole number is pretty ugly to write
 550          we are dividing by first 3/4 of bits of number.  */
 551
 552       tmp1 = (a->sig_hi << SREAL_PART_BITS) + a->sig_lo;
 553       tmp2 = ((b->sig_hi << (SREAL_PART_BITS / 2))
 554               + (b->sig_lo >> (SREAL_PART_BITS / 2)));
 555       if (b->sig_lo & ((uhwi) 1 << ((SREAL_PART_BITS / 2) - 1)))
 556         tmp2++;
 557
 558       r->sig_lo = 0;
 559       tmp = tmp1 / tmp2;
 560       tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);
 561       r->sig_hi = tmp << SREAL_PART_BITS;
 562
 563       tmp = tmp1 / tmp2;
 564       tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);
 565       r->sig_hi += tmp << (SREAL_PART_BITS / 2);
 566
 567       tmp = tmp1 / tmp2;
 568       r->sig_hi += tmp;
 569
 570       r->exp = a->exp - b->exp - SREAL_BITS - SREAL_PART_BITS / 2;
 571       normalize (r);
 572     }
 573 #else
 574   if (b->sig == 0)
 575     {
 576       abort ();
 577     }
 578   else
 579     {
 580       r->sig = (a->sig << SREAL_PART_BITS) / b->sig;
 581       r->exp = a->exp - b->exp - SREAL_PART_BITS;
 582       normalize (r);
 583     }
 584 #endif
 585   return r;
 586 }