gcc/sreal.c

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