gcc/sreal.c

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