| |
| /* |
| * 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 GET_BITS_DP64(x, ux) \ |
| { \ |
| volatile union {double d; unsigned long long i;} _bitsy; \ |
| _bitsy.d = (x); \ |
| ux = _bitsy.i; \ |
| } |
| |
| #define PUT_BITS_DP64(ux, x) \ |
| { \ |
| volatile union {double d; unsigned long long i;} _bitsy; \ |
| _bitsy.i = (ux); \ |
| x = _bitsy.d; \ |
| } |
| |
| /* Define this to get debugging print statements activated */ |
| #define DEBUGGING_PRINT |
| #undef DEBUGGING_PRINT |
| |
| |
| #ifdef DEBUGGING_PRINT |
| #include <stdio.h> |
| char *d2b(int d, int bitsper, int point) |
| { |
| static char buff[50]; |
| int i, j; |
| j = bitsper; |
| if (point >= 0 && point <= bitsper) |
| j++; |
| buff[j] = '\0'; |
| for (i = bitsper - 1; i >= 0; i--) |
| { |
| j--; |
| if (d % 2 == 1) |
| buff[j] = '1'; |
| else |
| buff[j] = '0'; |
| if (i == point) |
| { |
| j--; |
| buff[j] = '.'; |
| } |
| d /= 2; |
| } |
| return buff; |
| } |
| #endif |
| |
| /* Given positive argument x, reduce it to the range [-pi/4,pi/4] using |
| extra precision, and return the result in r, rr. |
| 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_piby2(double x, double *r, double *rr, int *region) |
| { |
| |
| /* This method simulates multi-precision floating-point |
| arithmetic and is accurate for all 1 <= x < infinity */ |
| static const double |
| piby2_lead = 1.57079632679489655800e+00, /* 0x3ff921fb54442d18 */ |
| piby2_part1 = 1.57079631090164184570e+00, /* 0x3ff921fb50000000 */ |
| piby2_part2 = 1.58932547122958567343e-08, /* 0x3e5110b460000000 */ |
| piby2_part3 = 6.12323399573676480327e-17; /* 0x3c91a62633145c06 */ |
| const int bitsper = 10; |
| unsigned long long res[500]; |
| unsigned long long ux, u, carry, mask, mant, highbitsrr; |
| int first, last, i, rexp, xexp, resexp, ltb, determ; |
| double xx, t; |
| 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 |
| }; |
| |
| GET_BITS_DP64(x, ux); |
| |
| #ifdef DEBUGGING_PRINT |
| printf("On entry, x = %25.20e = %s\n", x, double2hex(&x)); |
| #endif |
| |
| 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 */ |
| carry = 0; |
| 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; |
| |
| #ifdef DEBUGGING_PRINT |
| printf("resexp = %d\n", resexp); |
| printf("Significant part of x * 2/pi with binary" |
| " point in correct place:\n"); |
| for (i = 0; i <= last - first; i++) |
| { |
| if (i > 0 && i % 5 == 0) |
| printf("\n "); |
| if (i == 1) |
| printf("%s ", d2b((int)res[i], bitsper, resexp)); |
| else |
| printf("%s ", d2b((int)res[i], bitsper, -1)); |
| } |
| printf("\n"); |
| #endif |
| |
| /* 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; |
| |
| #ifdef DEBUGGING_PRINT |
| printf("ltb = %d (last two bits before binary point" |
| " and first bit after)\n", ltb); |
| printf("determ = %d (1 means need to negate because the fractional\n" |
| " part of x * 2/pi is greater than 0.5)\n", determ); |
| #endif |
| |
| 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; |
| } |
| |
| #ifdef DEBUGGING_PRINT |
| printf("Normalised mantissa = 0x%016lx\n", mant); |
| printf("High bits of rest of mantissa = 0x%016lx\n", highbitsrr); |
| printf("Exponent to be inserted on mantissa = rexp = %d\n", rexp); |
| #endif |
| |
| /* 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, x); |
| |
| /* Create the bit pattern for rr */ |
| highbitsrr >>= 12; /* Note this is shifted one place too far */ |
| u = ((unsigned long long)rexp + EXPBIAS_DP64 - 53) << EXPSHIFTBITS_DP64; |
| PUT_BITS_DP64(u, t); |
| u |= highbitsrr; |
| PUT_BITS_DP64(u, xx); |
| |
| /* Subtract the implicit bit we accidentally added */ |
| xx -= t; |
| /* Set the correct sign, and double to account for the |
| "one place too far" shift */ |
| if (determ) |
| xx *= -2.0; |
| else |
| xx *= 2.0; |
| |
| #ifdef DEBUGGING_PRINT |
| printf("(lead part of x*2/pi) = %25.20e = %s\n", x, double2hex(&x)); |
| printf("(tail part of x*2/pi) = %25.20e = %s\n", xx, double2hex(&xx)); |
| #endif |
| |
| /* (x,xx) is an extra-precise version of the fractional part of |
| x * 2 / pi. Multiply (x,xx) by pi/2 in extra precision |
| to get the reduced argument (r,rr). */ |
| { |
| double hx, tx, c, cc; |
| /* Split x into hx (head) and tx (tail) */ |
| GET_BITS_DP64(x, ux); |
| ux &= 0xfffffffff8000000; |
| PUT_BITS_DP64(ux, hx); |
| tx = x - hx; |
| |
| c = piby2_lead * x; |
| cc = ((((piby2_part1 * hx - c) + piby2_part1 * tx) + |
| piby2_part2 * hx) + piby2_part2 * tx) + |
| (piby2_lead * xx + piby2_part3 * x); |
| *r = c + cc; |
| *rr = (c - *r) + cc; |
| } |
| |
| #ifdef DEBUGGING_PRINT |
| printf(" (r,rr) = lead and tail parts of frac(x*2/pi) * pi/2:\n"); |
| printf(" r = %25.20e = %s\n", *r, double2hex(r)); |
| printf("rr = %25.20e = %s\n", *rr, double2hex(rr)); |
| printf("region = (number of pi/2 subtracted from x) mod 4 = %d\n", |
| *region); |
| #endif |
| return; |
| } |
