gonzui

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