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