| |
| /* |
| * Copyright (C) 2008-2009 Advanced Micro Devices, Inc. All Rights Reserved. |
| * |
| * This file is part of libacml_mv. |
| * |
| * libacml_mv is free software; you can redistribute it and/or |
| * modify it under the terms of the GNU Lesser General Public |
| * License as published by the Free Software Foundation; either |
| * version 2.1 of the License, or (at your option) any later version. |
| * |
| * libacml_mv is distributed in the hope that it will be useful, |
| * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| * Lesser General Public License for more details. |
| * |
| * You should have received a copy of the GNU Lesser General Public |
| * License along with libacml_mv. If not, see |
| * <http://www.gnu.org/licenses/>. |
| * |
| */ |
| |
| |
| |
| #define EXPBITS_DP64 0x7ff0000000000000 |
| #define EXPSHIFTBITS_DP64 52 |
| #define EXPBIAS_DP64 1023 |
| #define MANTBITS_DP64 0x000fffffffffffff |
| #define IMPBIT_DP64 0x0010000000000000 |
| #define SIGNBIT_DP64 0x8000000000000000 |
| |
| #define PUT_BITS_DP64(ux, x) \ |
| { \ |
| volatile union {double d; unsigned long long i;} _bitsy; \ |
| _bitsy.i = (ux); \ |
| x = _bitsy.d; \ |
| } |
| |
| /*Derived from static inline void __amd_remainder_piby2f_inline(unsigned long long ux, double *r, int *region) |
| in libm_inlines_amd.h. libm_inlines.h has the pure Windows one while libm_inlines_amd.h has the mixed one. |
| */ |
| /* Given positive argument x, reduce it to the range [-pi/4,pi/4] using |
| extra precision, and return the result in r. |
| Return value "region" tells how many lots of pi/2 were subtracted |
| from x to put it in the range [-pi/4,pi/4], mod 4. */ |
| void __amd_remainder_piby2d2f(unsigned long long ux, double *r, int *region) |
| { |
| /* This method simulates multi-precision floating-point |
| arithmetic and is accurate for all 1 <= x < infinity */ |
| unsigned long long u, carry, mask, mant, highbitsrr; |
| double dx; |
| unsigned long long res[500]; |
| int first, last, i, rexp, xexp, resexp, ltb, determ; |
| static const double |
| piby2 = 1.57079632679489655800e+00; /* 0x3ff921fb54442d18 */ |
| const int bitsper = 10; |
| static unsigned long long pibits[] = |
| { |
| 0, 0, 0, 0, 0, 0, |
| 162, 998, 54, 915, 580, 84, 671, 777, 855, 839, |
| 851, 311, 448, 877, 553, 358, 316, 270, 260, 127, |
| 593, 398, 701, 942, 965, 390, 882, 283, 570, 265, |
| 221, 184, 6, 292, 750, 642, 465, 584, 463, 903, |
| 491, 114, 786, 617, 830, 930, 35, 381, 302, 749, |
| 72, 314, 412, 448, 619, 279, 894, 260, 921, 117, |
| 569, 525, 307, 637, 156, 529, 504, 751, 505, 160, |
| 945, 1022, 151, 1023, 480, 358, 15, 956, 753, 98, |
| 858, 41, 721, 987, 310, 507, 242, 498, 777, 733, |
| 244, 399, 870, 633, 510, 651, 373, 158, 940, 506, |
| 997, 965, 947, 833, 825, 990, 165, 164, 746, 431, |
| 949, 1004, 287, 565, 464, 533, 515, 193, 111, 798 |
| }; |
| |
| xexp = (int)(((ux & EXPBITS_DP64) >> EXPSHIFTBITS_DP64) - EXPBIAS_DP64); |
| ux = (ux & MANTBITS_DP64) | IMPBIT_DP64; |
| |
| /* Now ux is the mantissa bit pattern of x as a long integer */ |
| mask = 1; |
| mask = (mask << bitsper) - 1; |
| |
| /* Set first and last to the positions of the first |
| and last chunks of 2/pi that we need */ |
| first = xexp / bitsper; |
| resexp = xexp - first * bitsper; |
| /* 180 is the theoretical maximum number of bits (actually |
| 175 for IEEE double precision) that we need to extract |
| from the middle of 2/pi to compute the reduced argument |
| accurately enough for our purposes */ |
| last = first + 180 / bitsper; |
| |
| /* Do a long multiplication of the bits of 2/pi by the |
| integer mantissa */ |
| /* Unroll the loop. This is only correct because we know |
| that bitsper is fixed as 10. */ |
| res[19] = 0; |
| u = pibits[last] * ux; |
| res[18] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-1] * ux + carry; |
| res[17] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-2] * ux + carry; |
| res[16] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-3] * ux + carry; |
| res[15] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-4] * ux + carry; |
| res[14] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-5] * ux + carry; |
| res[13] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-6] * ux + carry; |
| res[12] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-7] * ux + carry; |
| res[11] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-8] * ux + carry; |
| res[10] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-9] * ux + carry; |
| res[9] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-10] * ux + carry; |
| res[8] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-11] * ux + carry; |
| res[7] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-12] * ux + carry; |
| res[6] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-13] * ux + carry; |
| res[5] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-14] * ux + carry; |
| res[4] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-15] * ux + carry; |
| res[3] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-16] * ux + carry; |
| res[2] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-17] * ux + carry; |
| res[1] = u & mask; |
| carry = u >> bitsper; |
| u = pibits[last-18] * ux + carry; |
| res[0] = u & mask; |
| |
| /* Reconstruct the result */ |
| ltb = (int)((((res[0] << bitsper) | res[1]) |
| >> (bitsper - 1 - resexp)) & 7); |
| |
| /* determ says whether the fractional part is >= 0.5 */ |
| determ = ltb & 1; |
| |
| i = 1; |
| if (determ) |
| { |
| /* The mantissa is >= 0.5. We want to subtract it |
| from 1.0 by negating all the bits */ |
| *region = ((ltb >> 1) + 1) & 3; |
| mant = 1; |
| mant = ~(res[1]) & ((mant << (bitsper - resexp)) - 1); |
| while (mant < 0x0020000000000000) |
| { |
| i++; |
| mant = (mant << bitsper) | (~(res[i]) & mask); |
| } |
| highbitsrr = ~(res[i + 1]) << (64 - bitsper); |
| } |
| else |
| { |
| *region = (ltb >> 1); |
| mant = 1; |
| mant = res[1] & ((mant << (bitsper - resexp)) - 1); |
| while (mant < 0x0020000000000000) |
| { |
| i++; |
| mant = (mant << bitsper) | res[i]; |
| } |
| highbitsrr = res[i + 1] << (64 - bitsper); |
| } |
| |
| rexp = 52 + resexp - i * bitsper; |
| |
| while (mant >= 0x0020000000000000) |
| { |
| rexp++; |
| highbitsrr = (highbitsrr >> 1) | ((mant & 1) << 63); |
| mant >>= 1; |
| } |
| |
| /* Put the result exponent rexp onto the mantissa pattern */ |
| u = ((unsigned long long)rexp + EXPBIAS_DP64) << EXPSHIFTBITS_DP64; |
| ux = (mant & MANTBITS_DP64) | u; |
| if (determ) |
| /* If we negated the mantissa we negate x too */ |
| ux |= SIGNBIT_DP64; |
| PUT_BITS_DP64(ux, dx); |
| |
| /* x is a double precision version of the fractional part of |
| x * 2 / pi. Multiply x by pi/2 in double precision |
| to get the reduced argument r. */ |
| *r = dx * piby2; |
| |
| return; |
| } |
| |
| void __remainder_piby2d2f(unsigned long ux, double *r, int *region) |
| { |
| __amd_remainder_piby2d2f((unsigned long long) ux, r, region); |
| } |
| |
