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