1132718Skan/* Simple data type for positive real numbers for the GNU compiler. 2169689Skan Copyright (C) 2002, 2003, 2004 Free Software Foundation, Inc. 3132718Skan 4132718SkanThis file is part of GCC. 5132718Skan 6132718SkanGCC is free software; you can redistribute it and/or modify it under 7132718Skanthe terms of the GNU General Public License as published by the Free 8132718SkanSoftware Foundation; either version 2, or (at your option) any later 9132718Skanversion. 10132718Skan 11132718SkanGCC is distributed in the hope that it will be useful, but WITHOUT ANY 12132718SkanWARRANTY; without even the implied warranty of MERCHANTABILITY or 13132718SkanFITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 14132718Skanfor more details. 15132718Skan 16132718SkanYou should have received a copy of the GNU General Public License 17132718Skanalong with GCC; see the file COPYING. If not, write to the Free 18169689SkanSoftware Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 19169689Skan02110-1301, USA. */ 20132718Skan 21132718Skan/* This library supports positive real numbers and 0; 22132718Skan inf and nan are NOT supported. 23132718Skan It is written to be simple and fast. 24132718Skan 25132718Skan Value of sreal is 26132718Skan x = sig * 2 ^ exp 27132718Skan where 28132718Skan sig = significant 29132718Skan (for < 64-bit machines sig = sig_lo + sig_hi * 2 ^ SREAL_PART_BITS) 30132718Skan exp = exponent 31132718Skan 32132718Skan One HOST_WIDE_INT is used for the significant on 64-bit (and more than 33132718Skan 64-bit) machines, 34132718Skan otherwise two HOST_WIDE_INTs are used for the significant. 35132718Skan Only a half of significant bits is used (in normalized sreals) so that we do 36132718Skan not have problems with overflow, for example when c->sig = a->sig * b->sig. 37132718Skan So the precision for 64-bit and 32-bit machines is 32-bit. 38132718Skan 39132718Skan Invariant: The numbers are normalized before and after each call of sreal_*. 40132718Skan 41132718Skan Normalized sreals: 42132718Skan All numbers (except zero) meet following conditions: 43132718Skan SREAL_MIN_SIG <= sig && sig <= SREAL_MAX_SIG 44132718Skan -SREAL_MAX_EXP <= exp && exp <= SREAL_MAX_EXP 45132718Skan 46132718Skan If the number would be too large, it is set to upper bounds of these 47132718Skan conditions. 48132718Skan 49132718Skan If the number is zero or would be too small it meets following conditions: 50132718Skan sig == 0 && exp == -SREAL_MAX_EXP 51132718Skan*/ 52132718Skan 53132718Skan#include "config.h" 54132718Skan#include "system.h" 55132718Skan#include "coretypes.h" 56132718Skan#include "tm.h" 57132718Skan#include "sreal.h" 58132718Skan 59132718Skanstatic inline void copy (sreal *, sreal *); 60132718Skanstatic inline void shift_right (sreal *, int); 61132718Skanstatic void normalize (sreal *); 62132718Skan 63132718Skan/* Print the content of struct sreal. */ 64132718Skan 65132718Skanvoid 66132718Skandump_sreal (FILE *file, sreal *x) 67132718Skan{ 68132718Skan#if SREAL_PART_BITS < 32 69132718Skan fprintf (file, "((" HOST_WIDE_INT_PRINT_UNSIGNED " * 2^16 + " 70132718Skan HOST_WIDE_INT_PRINT_UNSIGNED ") * 2^%d)", 71132718Skan x->sig_hi, x->sig_lo, x->exp); 72132718Skan#else 73132718Skan fprintf (file, "(" HOST_WIDE_INT_PRINT_UNSIGNED " * 2^%d)", x->sig, x->exp); 74132718Skan#endif 75132718Skan} 76132718Skan 77132718Skan/* Copy the sreal number. */ 78132718Skan 79132718Skanstatic inline void 80132718Skancopy (sreal *r, sreal *a) 81132718Skan{ 82132718Skan#if SREAL_PART_BITS < 32 83132718Skan r->sig_lo = a->sig_lo; 84132718Skan r->sig_hi = a->sig_hi; 85132718Skan#else 86132718Skan r->sig = a->sig; 87132718Skan#endif 88132718Skan r->exp = a->exp; 89132718Skan} 90132718Skan 91132718Skan/* Shift X right by S bits. Needed: 0 < S <= SREAL_BITS. 92132718Skan When the most significant bit shifted out is 1, add 1 to X (rounding). */ 93132718Skan 94132718Skanstatic inline void 95132718Skanshift_right (sreal *x, int s) 96132718Skan{ 97169689Skan gcc_assert (s > 0); 98169689Skan gcc_assert (s <= SREAL_BITS); 99169689Skan /* Exponent should never be so large because shift_right is used only by 100169689Skan sreal_add and sreal_sub ant thus the number cannot be shifted out from 101169689Skan exponent range. */ 102169689Skan gcc_assert (x->exp + s <= SREAL_MAX_EXP); 103132718Skan 104132718Skan x->exp += s; 105132718Skan 106132718Skan#if SREAL_PART_BITS < 32 107132718Skan if (s > SREAL_PART_BITS) 108132718Skan { 109132718Skan s -= SREAL_PART_BITS; 110132718Skan x->sig_hi += (uhwi) 1 << (s - 1); 111132718Skan x->sig_lo = x->sig_hi >> s; 112132718Skan x->sig_hi = 0; 113132718Skan } 114132718Skan else 115132718Skan { 116132718Skan x->sig_lo += (uhwi) 1 << (s - 1); 117132718Skan if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS)) 118132718Skan { 119132718Skan x->sig_hi++; 120132718Skan x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS; 121132718Skan } 122132718Skan x->sig_lo >>= s; 123132718Skan x->sig_lo |= (x->sig_hi & (((uhwi) 1 << s) - 1)) << (SREAL_PART_BITS - s); 124132718Skan x->sig_hi >>= s; 125132718Skan } 126132718Skan#else 127132718Skan x->sig += (uhwi) 1 << (s - 1); 128132718Skan x->sig >>= s; 129132718Skan#endif 130132718Skan} 131132718Skan 132132718Skan/* Normalize *X. */ 133132718Skan 134132718Skanstatic void 135132718Skannormalize (sreal *x) 136132718Skan{ 137132718Skan#if SREAL_PART_BITS < 32 138132718Skan int shift; 139132718Skan HOST_WIDE_INT mask; 140132718Skan 141132718Skan if (x->sig_lo == 0 && x->sig_hi == 0) 142132718Skan { 143132718Skan x->exp = -SREAL_MAX_EXP; 144132718Skan } 145132718Skan else if (x->sig_hi < SREAL_MIN_SIG) 146132718Skan { 147132718Skan if (x->sig_hi == 0) 148132718Skan { 149132718Skan /* Move lower part of significant to higher part. */ 150132718Skan x->sig_hi = x->sig_lo; 151132718Skan x->sig_lo = 0; 152132718Skan x->exp -= SREAL_PART_BITS; 153132718Skan } 154132718Skan shift = 0; 155132718Skan while (x->sig_hi < SREAL_MIN_SIG) 156132718Skan { 157132718Skan x->sig_hi <<= 1; 158132718Skan x->exp--; 159132718Skan shift++; 160132718Skan } 161132718Skan /* Check underflow. */ 162132718Skan if (x->exp < -SREAL_MAX_EXP) 163132718Skan { 164132718Skan x->exp = -SREAL_MAX_EXP; 165132718Skan x->sig_hi = 0; 166132718Skan x->sig_lo = 0; 167132718Skan } 168132718Skan else if (shift) 169132718Skan { 170132718Skan mask = (1 << SREAL_PART_BITS) - (1 << (SREAL_PART_BITS - shift)); 171132718Skan x->sig_hi |= (x->sig_lo & mask) >> (SREAL_PART_BITS - shift); 172132718Skan x->sig_lo = (x->sig_lo << shift) & (((uhwi) 1 << SREAL_PART_BITS) - 1); 173132718Skan } 174132718Skan } 175132718Skan else if (x->sig_hi > SREAL_MAX_SIG) 176132718Skan { 177132718Skan unsigned HOST_WIDE_INT tmp = x->sig_hi; 178132718Skan 179132718Skan /* Find out how many bits will be shifted. */ 180132718Skan shift = 0; 181132718Skan do 182132718Skan { 183132718Skan tmp >>= 1; 184132718Skan shift++; 185132718Skan } 186132718Skan while (tmp > SREAL_MAX_SIG); 187132718Skan 188132718Skan /* Round the number. */ 189132718Skan x->sig_lo += (uhwi) 1 << (shift - 1); 190132718Skan 191132718Skan x->sig_lo >>= shift; 192132718Skan x->sig_lo += ((x->sig_hi & (((uhwi) 1 << shift) - 1)) 193132718Skan << (SREAL_PART_BITS - shift)); 194132718Skan x->sig_hi >>= shift; 195132718Skan x->exp += shift; 196132718Skan if (x->sig_lo & ((uhwi) 1 << SREAL_PART_BITS)) 197132718Skan { 198132718Skan x->sig_lo -= (uhwi) 1 << SREAL_PART_BITS; 199132718Skan x->sig_hi++; 200132718Skan if (x->sig_hi > SREAL_MAX_SIG) 201132718Skan { 202132718Skan /* x->sig_hi was SREAL_MAX_SIG before increment 203132718Skan so now last bit is zero. */ 204132718Skan x->sig_hi >>= 1; 205132718Skan x->sig_lo >>= 1; 206132718Skan x->exp++; 207132718Skan } 208132718Skan } 209132718Skan 210132718Skan /* Check overflow. */ 211132718Skan if (x->exp > SREAL_MAX_EXP) 212132718Skan { 213132718Skan x->exp = SREAL_MAX_EXP; 214132718Skan x->sig_hi = SREAL_MAX_SIG; 215132718Skan x->sig_lo = SREAL_MAX_SIG; 216132718Skan } 217132718Skan } 218132718Skan#else 219132718Skan if (x->sig == 0) 220132718Skan { 221132718Skan x->exp = -SREAL_MAX_EXP; 222132718Skan } 223132718Skan else if (x->sig < SREAL_MIN_SIG) 224132718Skan { 225132718Skan do 226132718Skan { 227132718Skan x->sig <<= 1; 228132718Skan x->exp--; 229132718Skan } 230132718Skan while (x->sig < SREAL_MIN_SIG); 231132718Skan 232132718Skan /* Check underflow. */ 233132718Skan if (x->exp < -SREAL_MAX_EXP) 234132718Skan { 235132718Skan x->exp = -SREAL_MAX_EXP; 236132718Skan x->sig = 0; 237132718Skan } 238132718Skan } 239132718Skan else if (x->sig > SREAL_MAX_SIG) 240132718Skan { 241132718Skan int last_bit; 242132718Skan do 243132718Skan { 244132718Skan last_bit = x->sig & 1; 245132718Skan x->sig >>= 1; 246132718Skan x->exp++; 247132718Skan } 248132718Skan while (x->sig > SREAL_MAX_SIG); 249132718Skan 250132718Skan /* Round the number. */ 251132718Skan x->sig += last_bit; 252132718Skan if (x->sig > SREAL_MAX_SIG) 253132718Skan { 254132718Skan x->sig >>= 1; 255132718Skan x->exp++; 256132718Skan } 257132718Skan 258132718Skan /* Check overflow. */ 259132718Skan if (x->exp > SREAL_MAX_EXP) 260132718Skan { 261132718Skan x->exp = SREAL_MAX_EXP; 262132718Skan x->sig = SREAL_MAX_SIG; 263132718Skan } 264132718Skan } 265132718Skan#endif 266132718Skan} 267132718Skan 268132718Skan/* Set *R to SIG * 2 ^ EXP. Return R. */ 269132718Skan 270132718Skansreal * 271132718Skansreal_init (sreal *r, unsigned HOST_WIDE_INT sig, signed int exp) 272132718Skan{ 273132718Skan#if SREAL_PART_BITS < 32 274132718Skan r->sig_lo = 0; 275132718Skan r->sig_hi = sig; 276132718Skan r->exp = exp - 16; 277132718Skan#else 278132718Skan r->sig = sig; 279132718Skan r->exp = exp; 280132718Skan#endif 281132718Skan normalize (r); 282132718Skan return r; 283132718Skan} 284132718Skan 285132718Skan/* Return integer value of *R. */ 286132718Skan 287132718SkanHOST_WIDE_INT 288132718Skansreal_to_int (sreal *r) 289132718Skan{ 290132718Skan#if SREAL_PART_BITS < 32 291132718Skan if (r->exp <= -SREAL_BITS) 292132718Skan return 0; 293132718Skan if (r->exp >= 0) 294132718Skan return MAX_HOST_WIDE_INT; 295132718Skan return ((r->sig_hi << SREAL_PART_BITS) + r->sig_lo) >> -r->exp; 296132718Skan#else 297132718Skan if (r->exp <= -SREAL_BITS) 298132718Skan return 0; 299132718Skan if (r->exp >= SREAL_PART_BITS) 300132718Skan return MAX_HOST_WIDE_INT; 301132718Skan if (r->exp > 0) 302132718Skan return r->sig << r->exp; 303132718Skan if (r->exp < 0) 304132718Skan return r->sig >> -r->exp; 305132718Skan return r->sig; 306132718Skan#endif 307132718Skan} 308132718Skan 309132718Skan/* Compare *A and *B. Return -1 if *A < *B, 1 if *A > *B and 0 if *A == *B. */ 310132718Skan 311132718Skanint 312132718Skansreal_compare (sreal *a, sreal *b) 313132718Skan{ 314132718Skan if (a->exp > b->exp) 315132718Skan return 1; 316132718Skan if (a->exp < b->exp) 317132718Skan return -1; 318132718Skan#if SREAL_PART_BITS < 32 319132718Skan if (a->sig_hi > b->sig_hi) 320132718Skan return 1; 321132718Skan if (a->sig_hi < b->sig_hi) 322132718Skan return -1; 323132718Skan if (a->sig_lo > b->sig_lo) 324132718Skan return 1; 325132718Skan if (a->sig_lo < b->sig_lo) 326132718Skan return -1; 327132718Skan#else 328132718Skan if (a->sig > b->sig) 329132718Skan return 1; 330132718Skan if (a->sig < b->sig) 331132718Skan return -1; 332132718Skan#endif 333132718Skan return 0; 334132718Skan} 335132718Skan 336132718Skan/* *R = *A + *B. Return R. */ 337132718Skan 338132718Skansreal * 339132718Skansreal_add (sreal *r, sreal *a, sreal *b) 340132718Skan{ 341132718Skan int dexp; 342132718Skan sreal tmp; 343132718Skan sreal *bb; 344132718Skan 345132718Skan if (sreal_compare (a, b) < 0) 346132718Skan { 347132718Skan sreal *swap; 348132718Skan swap = a; 349132718Skan a = b; 350132718Skan b = swap; 351132718Skan } 352132718Skan 353132718Skan dexp = a->exp - b->exp; 354132718Skan r->exp = a->exp; 355132718Skan if (dexp > SREAL_BITS) 356132718Skan { 357132718Skan#if SREAL_PART_BITS < 32 358132718Skan r->sig_hi = a->sig_hi; 359132718Skan r->sig_lo = a->sig_lo; 360132718Skan#else 361132718Skan r->sig = a->sig; 362132718Skan#endif 363132718Skan return r; 364132718Skan } 365132718Skan 366132718Skan if (dexp == 0) 367132718Skan bb = b; 368132718Skan else 369132718Skan { 370132718Skan copy (&tmp, b); 371132718Skan shift_right (&tmp, dexp); 372132718Skan bb = &tmp; 373132718Skan } 374132718Skan 375132718Skan#if SREAL_PART_BITS < 32 376132718Skan r->sig_hi = a->sig_hi + bb->sig_hi; 377132718Skan r->sig_lo = a->sig_lo + bb->sig_lo; 378132718Skan if (r->sig_lo & ((uhwi) 1 << SREAL_PART_BITS)) 379132718Skan { 380132718Skan r->sig_hi++; 381132718Skan r->sig_lo -= (uhwi) 1 << SREAL_PART_BITS; 382132718Skan } 383132718Skan#else 384132718Skan r->sig = a->sig + bb->sig; 385132718Skan#endif 386132718Skan normalize (r); 387132718Skan return r; 388132718Skan} 389132718Skan 390132718Skan/* *R = *A - *B. Return R. */ 391132718Skan 392132718Skansreal * 393132718Skansreal_sub (sreal *r, sreal *a, sreal *b) 394132718Skan{ 395132718Skan int dexp; 396132718Skan sreal tmp; 397132718Skan sreal *bb; 398132718Skan 399169689Skan gcc_assert (sreal_compare (a, b) >= 0); 400132718Skan 401132718Skan dexp = a->exp - b->exp; 402132718Skan r->exp = a->exp; 403132718Skan if (dexp > SREAL_BITS) 404132718Skan { 405132718Skan#if SREAL_PART_BITS < 32 406132718Skan r->sig_hi = a->sig_hi; 407132718Skan r->sig_lo = a->sig_lo; 408132718Skan#else 409132718Skan r->sig = a->sig; 410132718Skan#endif 411132718Skan return r; 412132718Skan } 413132718Skan if (dexp == 0) 414132718Skan bb = b; 415132718Skan else 416132718Skan { 417132718Skan copy (&tmp, b); 418132718Skan shift_right (&tmp, dexp); 419132718Skan bb = &tmp; 420132718Skan } 421132718Skan 422132718Skan#if SREAL_PART_BITS < 32 423132718Skan if (a->sig_lo < bb->sig_lo) 424132718Skan { 425132718Skan r->sig_hi = a->sig_hi - bb->sig_hi - 1; 426132718Skan r->sig_lo = a->sig_lo + ((uhwi) 1 << SREAL_PART_BITS) - bb->sig_lo; 427132718Skan } 428132718Skan else 429132718Skan { 430132718Skan r->sig_hi = a->sig_hi - bb->sig_hi; 431132718Skan r->sig_lo = a->sig_lo - bb->sig_lo; 432132718Skan } 433132718Skan#else 434132718Skan r->sig = a->sig - bb->sig; 435132718Skan#endif 436132718Skan normalize (r); 437132718Skan return r; 438132718Skan} 439132718Skan 440132718Skan/* *R = *A * *B. Return R. */ 441132718Skan 442132718Skansreal * 443132718Skansreal_mul (sreal *r, sreal *a, sreal *b) 444132718Skan{ 445132718Skan#if SREAL_PART_BITS < 32 446132718Skan if (a->sig_hi < SREAL_MIN_SIG || b->sig_hi < SREAL_MIN_SIG) 447132718Skan { 448132718Skan r->sig_lo = 0; 449132718Skan r->sig_hi = 0; 450132718Skan r->exp = -SREAL_MAX_EXP; 451132718Skan } 452132718Skan else 453132718Skan { 454132718Skan unsigned HOST_WIDE_INT tmp1, tmp2, tmp3; 455132718Skan if (sreal_compare (a, b) < 0) 456132718Skan { 457132718Skan sreal *swap; 458132718Skan swap = a; 459132718Skan a = b; 460132718Skan b = swap; 461132718Skan } 462132718Skan 463132718Skan r->exp = a->exp + b->exp + SREAL_PART_BITS; 464132718Skan 465132718Skan tmp1 = a->sig_lo * b->sig_lo; 466132718Skan tmp2 = a->sig_lo * b->sig_hi; 467132718Skan tmp3 = a->sig_hi * b->sig_lo + (tmp1 >> SREAL_PART_BITS); 468132718Skan 469132718Skan r->sig_hi = a->sig_hi * b->sig_hi; 470132718Skan r->sig_hi += (tmp2 >> SREAL_PART_BITS) + (tmp3 >> SREAL_PART_BITS); 471132718Skan tmp2 &= ((uhwi) 1 << SREAL_PART_BITS) - 1; 472132718Skan tmp3 &= ((uhwi) 1 << SREAL_PART_BITS) - 1; 473132718Skan tmp1 = tmp2 + tmp3; 474132718Skan 475132718Skan r->sig_lo = tmp1 & (((uhwi) 1 << SREAL_PART_BITS) - 1); 476132718Skan r->sig_hi += tmp1 >> SREAL_PART_BITS; 477132718Skan 478132718Skan normalize (r); 479132718Skan } 480132718Skan#else 481132718Skan if (a->sig < SREAL_MIN_SIG || b->sig < SREAL_MIN_SIG) 482132718Skan { 483132718Skan r->sig = 0; 484132718Skan r->exp = -SREAL_MAX_EXP; 485132718Skan } 486132718Skan else 487132718Skan { 488132718Skan r->sig = a->sig * b->sig; 489132718Skan r->exp = a->exp + b->exp; 490132718Skan normalize (r); 491132718Skan } 492132718Skan#endif 493132718Skan return r; 494132718Skan} 495132718Skan 496132718Skan/* *R = *A / *B. Return R. */ 497132718Skan 498132718Skansreal * 499132718Skansreal_div (sreal *r, sreal *a, sreal *b) 500132718Skan{ 501132718Skan#if SREAL_PART_BITS < 32 502132718Skan unsigned HOST_WIDE_INT tmp, tmp1, tmp2; 503132718Skan 504169689Skan gcc_assert (b->sig_hi >= SREAL_MIN_SIG); 505169689Skan if (a->sig_hi < SREAL_MIN_SIG) 506132718Skan { 507132718Skan r->sig_hi = 0; 508132718Skan r->sig_lo = 0; 509132718Skan r->exp = -SREAL_MAX_EXP; 510132718Skan } 511132718Skan else 512132718Skan { 513132718Skan /* Since division by the whole number is pretty ugly to write 514132718Skan we are dividing by first 3/4 of bits of number. */ 515132718Skan 516132718Skan tmp1 = (a->sig_hi << SREAL_PART_BITS) + a->sig_lo; 517132718Skan tmp2 = ((b->sig_hi << (SREAL_PART_BITS / 2)) 518132718Skan + (b->sig_lo >> (SREAL_PART_BITS / 2))); 519132718Skan if (b->sig_lo & ((uhwi) 1 << ((SREAL_PART_BITS / 2) - 1))) 520132718Skan tmp2++; 521132718Skan 522132718Skan r->sig_lo = 0; 523132718Skan tmp = tmp1 / tmp2; 524132718Skan tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2); 525132718Skan r->sig_hi = tmp << SREAL_PART_BITS; 526132718Skan 527132718Skan tmp = tmp1 / tmp2; 528132718Skan tmp1 = (tmp1 % tmp2) << (SREAL_PART_BITS / 2); 529132718Skan r->sig_hi += tmp << (SREAL_PART_BITS / 2); 530132718Skan 531132718Skan tmp = tmp1 / tmp2; 532132718Skan r->sig_hi += tmp; 533132718Skan 534132718Skan r->exp = a->exp - b->exp - SREAL_BITS - SREAL_PART_BITS / 2; 535132718Skan normalize (r); 536132718Skan } 537132718Skan#else 538169689Skan gcc_assert (b->sig != 0); 539169689Skan r->sig = (a->sig << SREAL_PART_BITS) / b->sig; 540169689Skan r->exp = a->exp - b->exp - SREAL_PART_BITS; 541169689Skan normalize (r); 542132718Skan#endif 543132718Skan return r; 544132718Skan} 545