lib/libm/src/k_rem_pio2.c

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