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 void dump_sreal                 PARAMS ((FILE *, sreal *));
  60 static inline void copy         PARAMS ((sreal *, sreal *));
  61 static inline void shift_right  PARAMS ((sreal *, int));
  62 static void normalize           PARAMS ((sreal *));
  63
  64 /* Print the content of struct sreal.  */
  65
  66 void
  67 dump_sreal (file, x)
  68      FILE *file;
  69      sreal *x;
  70 {
  71 #if SREAL_PART_BITS < 32
  72   fprintf (file, "((");
  73   fprintf (file, HOST_WIDE_INT_PRINT_UNSIGNED, x->sig_hi);
  74   fprintf (file, " * 2^16 + ");
  75   fprintf (file, HOST_WIDE_INT_PRINT_UNSIGNED, x->sig_lo);
  76   fprintf (file, ") * 2^%d)", x->exp);
  77 #else
  78   fprintf (file, "(");
  79   fprintf (file, HOST_WIDE_INT_PRINT_UNSIGNED, x->sig);
  80   fprintf (file, " * 2^%d)", x->exp);
  81 #endif
  82 }
  83
  84 /* Copy the sreal number.  */
  85
  86 static inline void
  87 copy (r, a)
  88      sreal *r;
  89      sreal *a;
  90 {
  91 #if SREAL_PART_BITS < 32
  92   r->sig_lo = a->sig_lo;
  93   r->sig_hi = a->sig_hi;
  94 #else
  95   r->sig = a->sig;
  96 #endif
  97   r->exp = a->exp;
  98 }
  99
 100 /* Shift X right by S bits.  Needed: 0 < S <= SREAL_BITS.
 101    When the most significant bit shifted out is 1, add 1 to X (rounding).  */
 102
 103 static inline void
 104 shift_right (x, s)
 105      sreal *x;
 106      int s;
 107 {
 108 #ifdef ENABLE_CHECKING
 109   if (s <= 0 || s > SREAL_BITS)
 110     abort ();
 111   if (x->exp + s > SREAL_MAX_EXP)
 112     {
 113       /* Exponent should never be so large because shift_right is used only by
 114          sreal_add and sreal_sub ant thus the number cannot be shifted out from
 115          exponent range.  */
 116       abort ();
 117     }
 118 #endif
 119
 120   x->exp += s;
 121
 122 #if SREAL_PART_BITS < 32
 123   if (s > SREAL_PART_BITS)
 124     {
 125       s -= SREAL_PART_BITS;
 126       x->sig_hi += (uhwi) 1 << (s - 1);
 127       x->sig_lo = x->sig_hi >> s;
 128       x->sig_hi = 0;
 129     }
 130   else
 131     {
 132       x->sig_lo += (uhwi) 1 << (s - 1);
 133       if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
 134         {
 135           x->sig_hi++;
 136           x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
 137         }
 138       x->sig_lo >>= s;
 139       x->sig_lo |= (x->sig_hi & (((uhwi) 1 << s) - 1)) << (SREAL_PART_BITS - s);
 140       x->sig_hi >>= s;
 141     }
 142 #else
 143   x->sig += (uhwi) 1 << (s - 1);
 144   x->sig >>= s;
 145 #endif
 146 }
 147
 148 /* Normalize *X.  */
 149
 150 static void
 151 normalize (x)
 152      sreal *x;
 153 {
 154 #if SREAL_PART_BITS < 32
 155   int shift;
 156   HOST_WIDE_INT mask;
 157
 158   if (x->sig_lo == 0 && x->sig_hi == 0)
 159     {
 160       x->exp = -SREAL_MAX_EXP;
 161     }
 162   else if (x->sig_hi < SREAL_MIN_SIG)
 163     {
 164       if (x->sig_hi == 0)
 165         {
 166           /* Move lower part of significant to higher part.  */
 167           x->sig_hi = x->sig_lo;
 168           x->sig_lo = 0;
 169           x->exp -= SREAL_PART_BITS;
 170         }
 171       shift = 0;
 172       while (x->sig_hi < SREAL_MIN_SIG)
 173         {
 174           x->sig_hi <<= 1;
 175           x->exp--;
 176           shift++;
 177         }
 178       /* Check underflow.  */
 179       if (x->exp < -SREAL_MAX_EXP)
 180         {
 181           x->exp = -SREAL_MAX_EXP;
 182           x->sig_hi = 0;
 183           x->sig_lo = 0;
 184         }
 185       else if (shift)
 186         {
 187           mask = (1 << SREAL_PART_BITS) - (1 << (SREAL_PART_BITS - shift));
 188           x->sig_hi |= (x->sig_lo & mask) >> (SREAL_PART_BITS - shift);
 189           x->sig_lo = (x->sig_lo << shift) & (((uhwi) 1 << SREAL_PART_BITS) - 1);
 190         }
 191     }
 192   else if (x->sig_hi > SREAL_MAX_SIG)
 193     {
 194       unsigned HOST_WIDE_INT tmp = x->sig_hi;
 195
 196       /* Find out how many bits will be shifted.  */
 197       shift = 0;
 198       do
 199         {
 200           tmp >>= 1;
 201           shift++;
 202         }
 203       while (tmp > SREAL_MAX_SIG);
 204
 205       /* Round the number.  */
 206       x->sig_lo += (uhwi) 1 << (shift - 1);
 207
 208       x->sig_lo >>= shift;
 209       x->sig_lo += ((x->sig_hi & (((uhwi) 1 << shift) - 1))
 210                     << (SREAL_PART_BITS - shift));
 211       x->sig_hi >>= shift;
 212       x->exp += shift;
 213       if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
 214         {
 215           x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
 216           x->sig_hi++;
 217           if (x->sig_hi > SREAL_MAX_SIG)
 218             {
 219               /* x->sig_hi was SREAL_MAX_SIG before increment
 220                  so now last bit is zero.  */
 221               x->sig_hi >>= 1;
 222               x->sig_lo >>= 1;
 223               x->exp++;
 224             }
 225         }
 226
 227       /* Check overflow.  */
 228       if (x->exp > SREAL_MAX_EXP)
 229         {
 230           x->exp = SREAL_MAX_EXP;
 231           x->sig_hi = SREAL_MAX_SIG;
 232           x->sig_lo = SREAL_MAX_SIG;
 233         }
 234     }
 235 #else
 236   if (x->sig == 0)
 237     {
 238       x->exp = -SREAL_MAX_EXP;
 239     }
 240   else if (x->sig < SREAL_MIN_SIG)
 241     {
 242       do
 243         {
 244           x->sig <<= 1;
 245           x->exp--;
 246         }
 247       while (x->sig < SREAL_MIN_SIG);
 248
 249       /* Check underflow.  */
 250       if (x->exp < -SREAL_MAX_EXP)
 251         {
 252           x->exp = -SREAL_MAX_EXP;
 253           x->sig = 0;
 254         }
 255     }
 256   else if (x->sig > SREAL_MAX_SIG)
 257     {
 258       int last_bit;
 259       do
 260         {
 261           last_bit = x->sig & 1;
 262           x->sig >>= 1;
 263           x->exp++;
 264         }
 265       while (x->sig > SREAL_MAX_SIG);
 266
 267       /* Round the number.  */
 268       x->sig += last_bit;
 269       if (x->sig > SREAL_MAX_SIG)
 270         {
 271           x->sig >>= 1;
 272           x->exp++;
 273         }
 274
 275       /* Check overflow.  */
 276       if (x->exp > SREAL_MAX_EXP)
 277         {
 278           x->exp = SREAL_MAX_EXP;
 279           x->sig = SREAL_MAX_SIG;
 280         }
 281     }
 282 #endif
 283 }
 284
 285 /* Set *R to SIG * 2 ^ EXP.  Return R.  */
 286
 287 sreal *
 288 sreal_init (r, sig, exp)
 289      sreal *r;
 290      unsigned HOST_WIDE_INT sig;
 291      signed int exp;
 292 {
 293 #if SREAL_PART_BITS < 32
 294   r->sig_lo = 0;
 295   r->sig_hi = sig;
 296   r->exp = exp - 16;
 297 #else
 298   r->sig = sig;
 299   r->exp = exp;
 300 #endif
 301   normalize (r);
 302   return r;
 303 }
 304
 305 /* Return integer value of *R.  */
 306
 307 HOST_WIDE_INT
 308 sreal_to_int (r)
 309      sreal *r;
 310 {
 311 #if SREAL_PART_BITS < 32
 312   if (r->exp <= -SREAL_BITS)
 313     return 0;
 314   if (r->exp >= 0)
 315     return MAX_HOST_WIDE_INT;
 316   return ((r->sig_hi << SREAL_PART_BITS) + r->sig_lo) >> -r->exp;
 317 #else
 318   if (r->exp <= -SREAL_BITS)
 319     return 0;
 320   if (r->exp >= SREAL_PART_BITS)
 321     return MAX_HOST_WIDE_INT;
 322   if (r->exp > 0)
 323     return r->sig << r->exp;
 324   if (r->exp < 0)
 325     return r->sig >> -r->exp;
 326   return r->sig;
 327 #endif
 328 }
 329
 330 /* Compare *A and *B. Return -1 if *A < *B, 1 if *A > *B and 0 if *A == *B.  */
 331
 332 int
 333 sreal_compare (a, b)
 334      sreal *a;
 335      sreal *b;
 336 {
 337   if (a->exp > b->exp)
 338     return 1;
 339   if (a->exp < b->exp)
 340     return -1;
 341 #if SREAL_PART_BITS < 32
 342   if (a->sig_hi > b->sig_hi)
 343     return 1;
 344   if (a->sig_hi < b->sig_hi)
 345     return -1;
 346   if (a->sig_lo > b->sig_lo)
 347     return 1;
 348   if (a->sig_lo < b->sig_lo)
 349     return -1;
 350 #else
 351   if (a->sig > b->sig)
 352     return 1;
 353   if (a->sig < b->sig)
 354     return -1;
 355 #endif
 356   return 0;
 357 }
 358
 359 /* *R = *A + *B.  Return R.  */
 360
 361 sreal *
 362 sreal_add (r, a, b)
 363   sreal *r;
 364   sreal *a;
 365   sreal *b;
 366 {
 367   int dexp;
 368   sreal tmp;
 369   sreal *bb;
 370
 371   if (sreal_compare (a, b) < 0)
 372     {
 373       sreal *swap;
 374       swap = a;
 375       a = b;
 376       b = swap;
 377     }
 378
 379   dexp = a->exp - b->exp;
 380   r->exp = a->exp;
 381   if (dexp > SREAL_BITS)
 382     {
 383 #if SREAL_PART_BITS < 32
 384       r->sig_hi = a->sig_hi;
 385       r->sig_lo = a->sig_lo;
 386 #else
 387       r->sig = a->sig;
 388 #endif
 389       return r;
 390     }
 391
 392   if (dexp == 0)
 393     bb = b;
 394   else
 395     {
 396       copy (&tmp, b);
 397       shift_right (&tmp, dexp);
 398       bb = &tmp;
 399     }
 400
 401 #if SREAL_PART_BITS < 32
 402   r->sig_hi = a->sig_hi + bb->sig_hi;
 403   r->sig_lo = a->sig_lo + bb->sig_lo;
 404   if (r->sig_lo & ((uhwi) 1 << SREAL_PART_BITS))
 405     {
 406       r->sig_hi++;
 407       r->sig_lo -= (uhwi) 1 << SREAL_PART_BITS;
 408     }
 409 #else
 410   r->sig = a->sig + bb->sig;
 411 #endif
 412   normalize (r);
 413   return r;
 414 }
 415
 416 /* *R = *A - *B.  Return R.  */
 417
 418 sreal *
 419 sreal_sub (r, a, b)
 420   sreal *r;
 421   sreal *a;
 422   sreal *b;
 423 {
 424   int dexp;
 425   sreal tmp;
 426   sreal *bb;
 427
 428   if (sreal_compare (a, b) < 0)
 429     {
 430       abort ();
 431     }
 432
 433   dexp = a->exp - b->exp;
 434   r->exp = a->exp;
 435   if (dexp > SREAL_BITS)
 436     {
 437 #if SREAL_PART_BITS < 32
 438       r->sig_hi = a->sig_hi;
 439       r->sig_lo = a->sig_lo;
 440 #else
 441       r->sig = a->sig;
 442 #endif
 443       return r;
 444     }
 445   if (dexp == 0)
 446     bb = b;
 447   else
 448     {
 449       copy (&tmp, b);
 450       shift_right (&tmp, dexp);
 451       bb = &tmp;
 452     }
 453
 454 #if SREAL_PART_BITS < 32
 455   if (a->sig_lo < bb->sig_lo)
 456     {
 457       r->sig_hi = a->sig_hi - bb->sig_hi - 1;
 458       r->sig_lo = a->sig_lo + ((uhwi) 1 << SREAL_PART_BITS) - bb->sig_lo;
 459     }
 460   else
 461     {
 462       r->sig_hi = a->sig_hi - bb->sig_hi;
 463       r->sig_lo = a->sig_lo - bb->sig_lo;
 464     }
 465 #else
 466   r->sig = a->sig - bb->sig;
 467 #endif
 468   normalize (r);
 469   return r;
 470 }
 471
 472 /* *R = *A * *B.  Return R.  */
 473
 474 sreal *
 475 sreal_mul (r, a, b)
 476      sreal *r;
 477      sreal *a;
 478      sreal *b;
 479 {
 480 #if SREAL_PART_BITS < 32
 481   if (a->sig_hi < SREAL_MIN_SIG || b->sig_hi < SREAL_MIN_SIG)
 482     {
 483       r->sig_lo = 0;
 484       r->sig_hi = 0;
 485       r->exp = -SREAL_MAX_EXP;
 486     }
 487   else
 488     {
 489       unsigned HOST_WIDE_INT tmp1, tmp2, tmp3;
 490       if (sreal_compare (a, b) < 0)
 491         {
 492           sreal *swap;
 493           swap = a;
 494           a = b;
 495           b = swap;
 496         }
 497
 498       r->exp = a->exp + b->exp + SREAL_PART_BITS;
 499
 500       tmp1 = a->sig_lo * b->sig_lo;
 501       tmp2 = a->sig_lo * b->sig_hi;
 502       tmp3 = a->sig_hi * b->sig_lo + (tmp1 >> SREAL_PART_BITS);
 503
 504       r->sig_hi = a->sig_hi * b->sig_hi;
 505       r->sig_hi += (tmp2 >> SREAL_PART_BITS) + (tmp3 >> SREAL_PART_BITS);
 506       tmp2 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;
 507       tmp3 &= ((uhwi) 1 << SREAL_PART_BITS) - 1;
 508       tmp1 = tmp2 + tmp3;
 509
 510       r->sig_lo = tmp1 & (((uhwi) 1 << SREAL_PART_BITS) - 1);
 511       r->sig_hi += tmp1 >> SREAL_PART_BITS;
 512
 513       normalize (r);
 514     }
 515 #else
 516   if (a->sig < SREAL_MIN_SIG || b->sig < SREAL_MIN_SIG)
 517     {
 518       r->sig = 0;
 519       r->exp = -SREAL_MAX_EXP;
 520     }
 521   else
 522     {
 523       r->sig = a->sig * b->sig;
 524       r->exp = a->exp + b->exp;
 525       normalize (r);
 526     }
 527 #endif
 528   return r;
 529 }
 530
 531 /* *R = *A / *B.  Return R.  */
 532
 533 sreal *
 534 sreal_div (r, a, b)
 535      sreal *r;
 536      sreal *a;
 537      sreal *b;
 538 {
 539 #if SREAL_PART_BITS < 32
 540   unsigned HOST_WIDE_INT tmp, tmp1, tmp2;
 541
 542   if (b->sig_hi < SREAL_MIN_SIG)
 543     {
 544       abort ();
 545     }
 546   else if (a->sig_hi < SREAL_MIN_SIG)
 547     {
 548       r->sig_hi = 0;
 549       r->sig_lo = 0;
 550       r->exp = -SREAL_MAX_EXP;
 551     }
 552   else
 553     {
 554       /* Since division by the whole number is pretty ugly to write
 555          we are dividing by first 3/4 of bits of number.  */
 556
 557       tmp1 = (a->sig_hi << SREAL_PART_BITS) + a->sig_lo;
 558       tmp2 = ((b->sig_hi << (SREAL_PART_BITS / 2))
 559               + (b->sig_lo >> (SREAL_PART_BITS / 2)));
 560       if (b->sig_lo & ((uhwi) 1 << ((SREAL_PART_BITS / 2) - 1)))
 561         tmp2++;
 562
 563       r->sig_lo = 0;
 564       tmp = tmp1 / tmp2;
 565       tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);
 566       r->sig_hi = tmp << SREAL_PART_BITS;
 567
 568       tmp = tmp1 / tmp2;
 569       tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2);
 570       r->sig_hi += tmp << (SREAL_PART_BITS / 2);
 571
 572       tmp = tmp1 / tmp2;
 573       r->sig_hi += tmp;
 574
 575       r->exp = a->exp - b->exp - SREAL_BITS - SREAL_PART_BITS / 2;
 576       normalize (r);
 577     }
 578 #else
 579   if (b->sig == 0)
 580     {
 581       abort ();
 582     }
 583   else
 584     {
 585       r->sig = (a->sig << SREAL_PART_BITS) / b->sig;
 586       r->exp = a->exp - b->exp - SREAL_PART_BITS;
 587       normalize (r);
 588     }
 589 #endif
 590   return r;
 591 }