lib/msun/src/k_rem_pio2.c

   1
   2 /* @(#)k_rem_pio2.c 1.3 95/01/18 */
   3 /*
   4  * ====================================================
   5  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   6  *
   7  * Developed at SunSoft, a Sun Microsystems, Inc. business.
   8  * Permission to use, copy, modify, and distribute this
   9  * software is freely granted, provided that this notice
  10  * is preserved.
  11  * ====================================================
  12  */
  13
  14 #include <sys/cdefs.h>
  15 __FBSDID("$FreeBSD$");
  16
  17 /*
  18  * __kernel_rem_pio2(x,y,e0,nx,prec)
  19  * double x[],y[]; int e0,nx,prec;
  20  *
  21  * __kernel_rem_pio2 return the last three digits of N with
  22  *              y = x - N*pi/2
  23  * so that |y| < pi/2.
  24  *
  25  * The method is to compute the integer (mod 8) and fraction parts of
  26  * (2/pi)*x without doing the full multiplication. In general we
  27  * skip the part of the product that are known to be a huge integer (
  28  * more accurately, = 0 mod 8 ). Thus the number of operations are
  29  * independent of the exponent of the input.
  30  *
  31  * (2/pi) is represented by an array of 24-bit integers in ipio2[].
  32  *
  33  * Input parameters:
  34  *      x[]     The input value (must be positive) is broken into nx
  35  *              pieces of 24-bit integers in double precision format.
  36  *              x[i] will be the i-th 24 bit of x. The scaled exponent
  37  *              of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
  38  *              match x's up to 24 bits.
  39  *
  40  *              Example of breaking a double positive z into x[0]+x[1]+x[2]:
  41  *                      e0 = ilogb(z)-23
  42  *                      z  = scalbn(z,-e0)
  43  *              for i = 0,1,2
  44  *                      x[i] = floor(z)
  45  *                      z    = (z-x[i])*2**24
  46  *
  47  *
  48  *      y[]     output result in an array of double precision numbers.
  49  *              The dimension of y[] is:
  50  *                      24-bit  precision       1
  51  *                      53-bit  precision       2
  52  *                      64-bit  precision       2
  53  *                      113-bit precision       3
  54  *              The actual value is the sum of them. Thus for 113-bit
  55  *              precison, one may have to do something like:
  56  *
  57  *              long double t,w,r_head, r_tail;
  58  *              t = (long double)y[2] + (long double)y[1];
  59  *              w = (long double)y[0];
  60  *              r_head = t+w;
  61  *              r_tail = w - (r_head - t);
  62  *
  63  *      e0      The exponent of x[0]. Must be <= 16360 or you need to
  64  *              expand the ipio2 table.
  65  *
  66  *      nx      dimension of x[]
  67  *
  68  *      prec    an integer indicating the precision:
  69  *                      0       24  bits (single)
  70  *                      1       53  bits (double)
  71  *                      2       64  bits (extended)
  72  *                      3       113 bits (quad)
  73  *
  74  * External function:
  75  *      double scalbn(), floor();
  76  *
  77  *
  78  * Here is the description of some local variables:
  79  *
  80  *      jk      jk+1 is the initial number of terms of ipio2[] needed
  81  *              in the computation. The minimum and recommended value
  82  *              for jk is 3,4,4,6 for single, double, extended, and quad.
  83  *              jk+1 must be 2 larger than you might expect so that our
  84  *              recomputation test works. (Up to 24 bits in the integer
  85  *              part (the 24 bits of it that we compute) and 23 bits in
  86  *              the fraction part may be lost to cancelation before we
  87  *              recompute.)
  88  *
  89  *      jz      local integer variable indicating the number of
  90  *              terms of ipio2[] used.
  91  *
  92  *      jx      nx - 1
  93  *
  94  *      jv      index for pointing to the suitable ipio2[] for the
  95  *              computation. In general, we want
  96  *                      ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
  97  *              is an integer. Thus
  98  *                      e0-3-24*jv >= 0 or (e0-3)/24 >= jv
  99  *              Hence jv = max(0,(e0-3)/24).
 100  *
 101  *      jp      jp+1 is the number of terms in PIo2[] needed, jp = jk.
 102  *
 103  *      q[]     double array with integral value, representing the
 104  *              24-bits chunk of the product of x and 2/pi.
 105  *
 106  *      q0      the corresponding exponent of q[0]. Note that the
 107  *              exponent for q[i] would be q0-24*i.
 108  *
 109  *      PIo2[]  double precision array, obtained by cutting pi/2
 110  *              into 24 bits chunks.
 111  *
 112  *      f[]     ipio2[] in floating point
 113  *
 114  *      iq[]    integer array by breaking up q[] in 24-bits chunk.
 115  *
 116  *      fq[]    final product of x*(2/pi) in fq[0],..,fq[jk]
 117  *
 118  *      ih      integer. If >0 it indicates q[] is >= 0.5, hence
 119  *              it also indicates the *sign* of the result.
 120  *
 121  */
 122
 123
 124 /*
 125  * Constants:
 126  * The hexadecimal values are the intended ones for the following
 127  * constants. The decimal values may be used, provided that the
 128  * compiler will convert from decimal to binary accurately enough
 129  * to produce the hexadecimal values shown.
 130  */
 131
 132 #include <float.h>
 133
 134 #include "math.h"
 135 #include "math_private.h"
 136
 137 static const int init_jk[] = {3,4,4,6}; /* initial value for jk */
 138
 139 /*
 140  * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
 141  *
 142  *              integer array, contains the (24*i)-th to (24*i+23)-th
 143  *              bit of 2/pi after binary point. The corresponding
 144  *              floating value is
 145  *
 146  *                      ipio2[i] * 2^(-24(i+1)).
 147  *
 148  * NB: This table must have at least (e0-3)/24 + jk terms.
 149  *     For quad precision (e0 <= 16360, jk = 6), this is 686.
 150  */
 151 static const int32_t ipio2[] = {
 152 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
 153 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
 154 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
 155 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
 156 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
 157 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
 158 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
 159 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
 160 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
 161 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
 162 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
 163
 164 #if LDBL_MAX_EXP > 1024
 165 #if LDBL_MAX_EXP > 16384
 166 #error "ipio2 table needs to be expanded"
 167 #endif
 168 0x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6,
 169 0xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2,
 170 0xDE4F98, 0x327DBB, 0xC33D26, 0xEF6B1E, 0x5EF89F, 0x3A1F35,
 171 0xCAF27F, 0x1D87F1, 0x21907C, 0x7C246A, 0xFA6ED5, 0x772D30,
 172 0x433B15, 0xC614B5, 0x9D19C3, 0xC2C4AD, 0x414D2C, 0x5D000C,
 173 0x467D86, 0x2D71E3, 0x9AC69B, 0x006233, 0x7CD2B4, 0x97A7B4,
 174 0xD55537, 0xF63ED7, 0x1810A3, 0xFC764D, 0x2A9D64, 0xABD770,
 175 0xF87C63, 0x57B07A, 0xE71517, 0x5649C0, 0xD9D63B, 0x3884A7,
 176 0xCB2324, 0x778AD6, 0x23545A, 0xB91F00, 0x1B0AF1, 0xDFCE19,
 177 0xFF319F, 0x6A1E66, 0x615799, 0x47FBAC, 0xD87F7E, 0xB76522,
 178 0x89E832, 0x60BFE6, 0xCDC4EF, 0x09366C, 0xD43F5D, 0xD7DE16,
 179 0xDE3B58, 0x929BDE, 0x2822D2, 0xE88628, 0x4D58E2, 0x32CAC6,
 180 0x16E308, 0xCB7DE0, 0x50C017, 0xA71DF3, 0x5BE018, 0x34132E,
 181 0x621283, 0x014883, 0x5B8EF5, 0x7FB0AD, 0xF2E91E, 0x434A48,
 182 0xD36710, 0xD8DDAA, 0x425FAE, 0xCE616A, 0xA4280A, 0xB499D3,
 183 0xF2A606, 0x7F775C, 0x83C2A3, 0x883C61, 0x78738A, 0x5A8CAF,
 184 0xBDD76F, 0x63A62D, 0xCBBFF4, 0xEF818D, 0x67C126, 0x45CA55,
 185 0x36D9CA, 0xD2A828, 0x8D61C2, 0x77C912, 0x142604, 0x9B4612,
 186 0xC459C4, 0x44C5C8, 0x91B24D, 0xF31700, 0xAD43D4, 0xE54929,
 187 0x10D5FD, 0xFCBE00, 0xCC941E, 0xEECE70, 0xF53E13, 0x80F1EC,
 188 0xC3E7B3, 0x28F8C7, 0x940593, 0x3E71C1, 0xB3092E, 0xF3450B,
 189 0x9C1288, 0x7B20AB, 0x9FB52E, 0xC29247, 0x2F327B, 0x6D550C,
 190 0x90A772, 0x1FE76B, 0x96CB31, 0x4A1679, 0xE27941, 0x89DFF4,
 191 0x9794E8, 0x84E6E2, 0x973199, 0x6BED88, 0x365F5F, 0x0EFDBB,
 192 0xB49A48, 0x6CA467, 0x427271, 0x325D8D, 0xB8159F, 0x09E5BC,
 193 0x25318D, 0x3974F7, 0x1C0530, 0x010C0D, 0x68084B, 0x58EE2C,
 194 0x90AA47, 0x02E774, 0x24D6BD, 0xA67DF7, 0x72486E, 0xEF169F,
 195 0xA6948E, 0xF691B4, 0x5153D1, 0xF20ACF, 0x339820, 0x7E4BF5,
 196 0x6863B2, 0x5F3EDD, 0x035D40, 0x7F8985, 0x295255, 0xC06437,
 197 0x10D86D, 0x324832, 0x754C5B, 0xD4714E, 0x6E5445, 0xC1090B,
 198 0x69F52A, 0xD56614, 0x9D0727, 0x50045D, 0xDB3BB4, 0xC576EA,
 199 0x17F987, 0x7D6B49, 0xBA271D, 0x296996, 0xACCCC6, 0x5414AD,
 200 0x6AE290, 0x89D988, 0x50722C, 0xBEA404, 0x940777, 0x7030F3,
 201 0x27FC00, 0xA871EA, 0x49C266, 0x3DE064, 0x83DD97, 0x973FA3,
 202 0xFD9443, 0x8C860D, 0xDE4131, 0x9D3992, 0x8C70DD, 0xE7B717,
 203 0x3BDF08, 0x2B3715, 0xA0805C, 0x93805A, 0x921110, 0xD8E80F,
 204 0xAF806C, 0x4BFFDB, 0x0F9038, 0x761859, 0x15A562, 0xBBCB61,
 205 0xB989C7, 0xBD4010, 0x04F2D2, 0x277549, 0xF6B6EB, 0xBB22DB,
 206 0xAA140A, 0x2F2689, 0x768364, 0x333B09, 0x1A940E, 0xAA3A51,
 207 0xC2A31D, 0xAEEDAF, 0x12265C, 0x4DC26D, 0x9C7A2D, 0x9756C0,
 208 0x833F03, 0xF6F009, 0x8C402B, 0x99316D, 0x07B439, 0x15200C,
 209 0x5BC3D8, 0xC492F5, 0x4BADC6, 0xA5CA4E, 0xCD37A7, 0x36A9E6,
 210 0x9492AB, 0x6842DD, 0xDE6319, 0xEF8C76, 0x528B68, 0x37DBFC,
 211 0xABA1AE, 0x3115DF, 0xA1AE00, 0xDAFB0C, 0x664D64, 0xB705ED,
 212 0x306529, 0xBF5657, 0x3AFF47, 0xB9F96A, 0xF3BE75, 0xDF9328,
 213 0x3080AB, 0xF68C66, 0x15CB04, 0x0622FA, 0x1DE4D9, 0xA4B33D,
 214 0x8F1B57, 0x09CD36, 0xE9424E, 0xA4BE13, 0xB52333, 0x1AAAF0,
 215 0xA8654F, 0xA5C1D2, 0x0F3F0B, 0xCD785B, 0x76F923, 0x048B7B,
 216 0x721789, 0x53A6C6, 0xE26E6F, 0x00EBEF, 0x584A9B, 0xB7DAC4,
 217 0xBA66AA, 0xCFCF76, 0x1D02D1, 0x2DF1B1, 0xC1998C, 0x77ADC3,
 218 0xDA4886, 0xA05DF7, 0xF480C6, 0x2FF0AC, 0x9AECDD, 0xBC5C3F,
 219 0x6DDED0, 0x1FC790, 0xB6DB2A, 0x3A25A3, 0x9AAF00, 0x9353AD,
 220 0x0457B6, 0xB42D29, 0x7E804B, 0xA707DA, 0x0EAA76, 0xA1597B,
 221 0x2A1216, 0x2DB7DC, 0xFDE5FA, 0xFEDB89, 0xFDBE89, 0x6C76E4,
 222 0xFCA906, 0x70803E, 0x156E85, 0xFF87FD, 0x073E28, 0x336761,
 223 0x86182A, 0xEABD4D, 0xAFE7B3, 0x6E6D8F, 0x396795, 0x5BBF31,
 224 0x48D784, 0x16DF30, 0x432DC7, 0x356125, 0xCE70C9, 0xB8CB30,
 225 0xFD6CBF, 0xA200A4, 0xE46C05, 0xA0DD5A, 0x476F21, 0xD21262,
 226 0x845CB9, 0x496170, 0xE0566B, 0x015299, 0x375550, 0xB7D51E,
 227 0xC4F133, 0x5F6E13, 0xE4305D, 0xA92E85, 0xC3B21D, 0x3632A1,
 228 0xA4B708, 0xD4B1EA, 0x21F716, 0xE4698F, 0x77FF27, 0x80030C,
 229 0x2D408D, 0xA0CD4F, 0x99A520, 0xD3A2B3, 0x0A5D2F, 0x42F9B4,
 230 0xCBDA11, 0xD0BE7D, 0xC1DB9B, 0xBD17AB, 0x81A2CA, 0x5C6A08,
 231 0x17552E, 0x550027, 0xF0147F, 0x8607E1, 0x640B14, 0x8D4196,
 232 0xDEBE87, 0x2AFDDA, 0xB6256B, 0x34897B, 0xFEF305, 0x9EBFB9,
 233 0x4F6A68, 0xA82A4A, 0x5AC44F, 0xBCF82D, 0x985AD7, 0x95C7F4,
 234 0x8D4D0D, 0xA63A20, 0x5F57A4, 0xB13F14, 0x953880, 0x0120CC,
 235 0x86DD71, 0xB6DEC9, 0xF560BF, 0x11654D, 0x6B0701, 0xACB08C,
 236 0xD0C0B2, 0x485551, 0x0EFB1E, 0xC37295, 0x3B06A3, 0x3540C0,
 237 0x7BDC06, 0xCC45E0, 0xFA294E, 0xC8CAD6, 0x41F3E8, 0xDE647C,
 238 0xD8649B, 0x31BED9, 0xC397A4, 0xD45877, 0xC5E369, 0x13DAF0,
 239 0x3C3ABA, 0x461846, 0x5F7555, 0xF5BDD2, 0xC6926E, 0x5D2EAC,
 240 0xED440E, 0x423E1C, 0x87C461, 0xE9FD29, 0xF3D6E7, 0xCA7C22,
 241 0x35916F, 0xC5E008, 0x8DD7FF, 0xE26A6E, 0xC6FDB0, 0xC10893,
 242 0x745D7C, 0xB2AD6B, 0x9D6ECD, 0x7B723E, 0x6A11C6, 0xA9CFF7,
 243 0xDF7329, 0xBAC9B5, 0x5100B7, 0x0DB2E2, 0x24BA74, 0x607DE5,
 244 0x8AD874, 0x2C150D, 0x0C1881, 0x94667E, 0x162901, 0x767A9F,
 245 0xBEFDFD, 0xEF4556, 0x367ED9, 0x13D9EC, 0xB9BA8B, 0xFC97C4,
 246 0x27A831, 0xC36EF1, 0x36C594, 0x56A8D8, 0xB5A8B4, 0x0ECCCF,
 247 0x2D8912, 0x34576F, 0x89562C, 0xE3CE99, 0xB920D6, 0xAA5E6B,
 248 0x9C2A3E, 0xCC5F11, 0x4A0BFD, 0xFBF4E1, 0x6D3B8E, 0x2C86E2,
 249 0x84D4E9, 0xA9B4FC, 0xD1EEEF, 0xC9352E, 0x61392F, 0x442138,
 250 0xC8D91B, 0x0AFC81, 0x6A4AFB, 0xD81C2F, 0x84B453, 0x8C994E,
 251 0xCC2254, 0xDC552A, 0xD6C6C0, 0x96190B, 0xB8701A, 0x649569,
 252 0x605A26, 0xEE523F, 0x0F117F, 0x11B5F4, 0xF5CBFC, 0x2DBC34,
 253 0xEEBC34, 0xCC5DE8, 0x605EDD, 0x9B8E67, 0xEF3392, 0xB817C9,
 254 0x9B5861, 0xBC57E1, 0xC68351, 0x103ED8, 0x4871DD, 0xDD1C2D,
 255 0xA118AF, 0x462C21, 0xD7F359, 0x987AD9, 0xC0549E, 0xFA864F,
 256 0xFC0656, 0xAE79E5, 0x362289, 0x22AD38, 0xDC9367, 0xAAE855,
 257 0x382682, 0x9BE7CA, 0xA40D51, 0xB13399, 0x0ED7A9, 0x480569,
 258 0xF0B265, 0xA7887F, 0x974C88, 0x36D1F9, 0xB39221, 0x4A827B,
 259 0x21CF98, 0xDC9F40, 0x5547DC, 0x3A74E1, 0x42EB67, 0xDF9DFE,
 260 0x5FD45E, 0xA4677B, 0x7AACBA, 0xA2F655, 0x23882B, 0x55BA41,
 261 0x086E59, 0x862A21, 0x834739, 0xE6E389, 0xD49EE5, 0x40FB49,
 262 0xE956FF, 0xCA0F1C, 0x8A59C5, 0x2BFA94, 0xC5C1D3, 0xCFC50F,
 263 0xAE5ADB, 0x86C547, 0x624385, 0x3B8621, 0x94792C, 0x876110,
 264 0x7B4C2A, 0x1A2C80, 0x12BF43, 0x902688, 0x893C78, 0xE4C4A8,
 265 0x7BDBE5, 0xC23AC4, 0xEAF426, 0x8A67F7, 0xBF920D, 0x2BA365,
 266 0xB1933D, 0x0B7CBD, 0xDC51A4, 0x63DD27, 0xDDE169, 0x19949A,
 267 0x9529A8, 0x28CE68, 0xB4ED09, 0x209F44, 0xCA984E, 0x638270,
 268 0x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5,
 269 0x4D7E6F, 0x5119A5, 0xABF9B5, 0xD6DF82, 0x61DD96, 0x023616,
 270 0x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B,
 271 0x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901, 0x8071E0,
 272 #endif
 273
 274 };
 275
 276 static const double PIo2[] = {
 277   1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
 278   7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
 279   5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
 280   3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
 281   1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
 282   1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
 283   2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
 284   2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
 285 };
 286
 287 static const double
 288 zero   = 0.0,
 289 one    = 1.0,
 290 two24   =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
 291 twon24  =  5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
 292
 293 int
 294 __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec)
 295 {
 296         int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
 297         double z,fw,f[20],fq[20],q[20];
 298
 299     /* initialize jk*/
 300         jk = init_jk[prec];
 301         jp = jk;
 302
 303     /* determine jx,jv,q0, note that 3>q0 */
 304         jx =  nx-1;
 305         jv = (e0-3)/24; if(jv<0) jv=0;
 306         q0 =  e0-24*(jv+1);
 307
 308     /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
 309         j = jv-jx; m = jx+jk;
 310         for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
 311
 312     /* compute q[0],q[1],...q[jk] */
 313         for (i=0;i<=jk;i++) {
 314             for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;
 315         }
 316
 317         jz = jk;
 318 recompute:
 319     /* distill q[] into iq[] reversingly */
 320         for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
 321             fw    =  (double)((int32_t)(twon24* z));
 322             iq[i] =  (int32_t)(z-two24*fw);
 323             z     =  q[j-1]+fw;
 324         }
 325
 326     /* compute n */
 327         z  = scalbn(z,q0);              /* actual value of z */
 328         z -= 8.0*floor(z*0.125);                /* trim off integer >= 8 */
 329         n  = (int32_t) z;
 330         z -= (double)n;
 331         ih = 0;
 332         if(q0>0) {      /* need iq[jz-1] to determine n */
 333             i  = (iq[jz-1]>>(24-q0)); n += i;
 334             iq[jz-1] -= i<<(24-q0);
 335             ih = iq[jz-1]>>(23-q0);
 336         }
 337         else if(q0==0) ih = iq[jz-1]>>23;
 338         else if(z>=0.5) ih=2;
 339
 340         if(ih>0) {      /* q > 0.5 */
 341             n += 1; carry = 0;
 342             for(i=0;i<jz ;i++) {        /* compute 1-q */
 343                 j = iq[i];
 344                 if(carry==0) {
 345                     if(j!=0) {
 346                         carry = 1; iq[i] = 0x1000000- j;
 347                     }
 348                 } else  iq[i] = 0xffffff - j;
 349             }
 350             if(q0>0) {          /* rare case: chance is 1 in 12 */
 351                 switch(q0) {
 352                 case 1:
 353                    iq[jz-1] &= 0x7fffff; break;
 354                 case 2:
 355                    iq[jz-1] &= 0x3fffff; break;
 356                 }
 357             }
 358             if(ih==2) {
 359                 z = one - z;
 360                 if(carry!=0) z -= scalbn(one,q0);
 361             }
 362         }
 363
 364     /* check if recomputation is needed */
 365         if(z==zero) {
 366             j = 0;
 367             for (i=jz-1;i>=jk;i--) j |= iq[i];
 368             if(j==0) { /* need recomputation */
 369                 for(k=1;iq[jk-k]==0;k++);   /* k = no. of terms needed */
 370
 371                 for(i=jz+1;i<=jz+k;i++) {   /* add q[jz+1] to q[jz+k] */
 372                     f[jx+i] = (double) ipio2[jv+i];
 373                     for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
 374                     q[i] = fw;
 375                 }
 376                 jz += k;
 377                 goto recompute;
 378             }
 379         }
 380
 381     /* chop off zero terms */
 382         if(z==0.0) {
 383             jz -= 1; q0 -= 24;
 384             while(iq[jz]==0) { jz--; q0-=24;}
 385         } else { /* break z into 24-bit if necessary */
 386             z = scalbn(z,-q0);
 387             if(z>=two24) {
 388                 fw = (double)((int32_t)(twon24*z));
 389                 iq[jz] = (int32_t)(z-two24*fw);
 390                 jz += 1; q0 += 24;
 391                 iq[jz] = (int32_t) fw;
 392             } else iq[jz] = (int32_t) z ;
 393         }
 394
 395     /* convert integer "bit" chunk to floating-point value */
 396         fw = scalbn(one,q0);
 397         for(i=jz;i>=0;i--) {
 398             q[i] = fw*(double)iq[i]; fw*=twon24;
 399         }
 400
 401     /* compute PIo2[0,...,jp]*q[jz,...,0] */
 402         for(i=jz;i>=0;i--) {
 403             for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
 404             fq[jz-i] = fw;
 405         }
 406
 407     /* compress fq[] into y[] */
 408         switch(prec) {
 409             case 0:
 410                 fw = 0.0;
 411                 for (i=jz;i>=0;i--) fw += fq[i];
 412                 y[0] = (ih==0)? fw: -fw;
 413                 break;
 414             case 1:
 415             case 2:
 416                 fw = 0.0;
 417                 for (i=jz;i>=0;i--) fw += fq[i];
 418                 STRICT_ASSIGN(double,fw,fw);
 419                 y[0] = (ih==0)? fw: -fw;
 420                 fw = fq[0]-fw;
 421                 for (i=1;i<=jz;i++) fw += fq[i];
 422                 y[1] = (ih==0)? fw: -fw;
 423                 break;
 424             case 3:     /* painful */
 425                 for (i=jz;i>0;i--) {
 426                     fw      = fq[i-1]+fq[i];
 427                     fq[i]  += fq[i-1]-fw;
 428                     fq[i-1] = fw;
 429                 }
 430                 for (i=jz;i>1;i--) {
 431                     fw      = fq[i-1]+fq[i];
 432                     fq[i]  += fq[i-1]-fw;
 433                     fq[i-1] = fw;
 434                 }
 435                 for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
 436                 if(ih==0) {
 437                     y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
 438                 } else {
 439                     y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
 440                 }
 441         }
 442         return n&7;
 443 }