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