gcc/sreal.c

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