Mercurial > hg > CbC > CbC_gcc
annotate gcc/sreal.c @ 67:f6334be47118
update gcc from gcc-4.6-20100522 to gcc-4.6-20110318
| author | nobuyasu <dimolto@cr.ie.u-ryukyu.ac.jp> |
|---|---|
| date | Tue, 22 Mar 2011 17:18:12 +0900 |
| parents | a06113de4d67 |
| children | 04ced10e8804 |
| rev | line source |
|---|---|
| 0 | 1 /* Simple data type for positive real numbers for the GNU compiler. |
|
67
f6334be47118
update gcc from gcc-4.6-20100522 to gcc-4.6-20110318
nobuyasu <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
0
diff
changeset
|
2 Copyright (C) 2002, 2003, 2004, 2007, 2010 Free Software Foundation, Inc. |
| 0 | 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 } |