lib/trigl.c

   1 /* Quad-precision floating point argument reduction.
   2    Copyright (C) 1999, 2007 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4    Contributed by Jakub Jelinek <jj@ultra.linux.cz>
   5
   6    This program is free software; you can redistribute it and/or modify
   7    it under the terms of the GNU General Public License as published by
   8    the Free Software Foundation; either version 2, or (at your option)
   9    any later version.
  10
  11    This program is distributed in the hope that it will be useful,
  12    but WITHOUT ANY WARRANTY; without even the implied warranty of
  13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14    GNU General Public License for more details.
  15
  16    You should have received a copy of the GNU General Public License along
  17    with this program; if not, write to the Free Software Foundation,
  18    Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.  */
  19
  20 #include <config.h>
  21
  22 /* Specification.  */
  23 #include <math.h>
  24
  25 #include <float.h>
  26
  27 /* Table of constants for 2/pi, 5628 hexadecimal digits of 2/pi */
  28 static const int two_over_pi[] = {
  29   0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62,
  30   0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a,
  31   0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129,
  32   0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41,
  33   0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8,
  34   0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf,
  35   0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5,
  36   0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08,
  37   0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3,
  38   0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880,
  39   0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b,
  40   0x47c419, 0xc367cd, 0xdce809, 0x2a8359, 0xc4768b, 0x961ca6,
  41   0xddaf44, 0xd15719, 0x053ea5, 0xff0705, 0x3f7e33, 0xe832c2,
  42   0xde4f98, 0x327dbb, 0xc33d26, 0xef6b1e, 0x5ef89f, 0x3a1f35,
  43   0xcaf27f, 0x1d87f1, 0x21907c, 0x7c246a, 0xfa6ed5, 0x772d30,
  44   0x433b15, 0xc614b5, 0x9d19c3, 0xc2c4ad, 0x414d2c, 0x5d000c,
  45   0x467d86, 0x2d71e3, 0x9ac69b, 0x006233, 0x7cd2b4, 0x97a7b4,
  46   0xd55537, 0xf63ed7, 0x1810a3, 0xfc764d, 0x2a9d64, 0xabd770,
  47   0xf87c63, 0x57b07a, 0xe71517, 0x5649c0, 0xd9d63b, 0x3884a7,
  48   0xcb2324, 0x778ad6, 0x23545a, 0xb91f00, 0x1b0af1, 0xdfce19,
  49   0xff319f, 0x6a1e66, 0x615799, 0x47fbac, 0xd87f7e, 0xb76522,
  50   0x89e832, 0x60bfe6, 0xcdc4ef, 0x09366c, 0xd43f5d, 0xd7de16,
  51   0xde3b58, 0x929bde, 0x2822d2, 0xe88628, 0x4d58e2, 0x32cac6,
  52   0x16e308, 0xcb7de0, 0x50c017, 0xa71df3, 0x5be018, 0x34132e,
  53   0x621283, 0x014883, 0x5b8ef5, 0x7fb0ad, 0xf2e91e, 0x434a48,
  54   0xd36710, 0xd8ddaa, 0x425fae, 0xce616a, 0xa4280a, 0xb499d3,
  55   0xf2a606, 0x7f775c, 0x83c2a3, 0x883c61, 0x78738a, 0x5a8caf,
  56   0xbdd76f, 0x63a62d, 0xcbbff4, 0xef818d, 0x67c126, 0x45ca55,
  57   0x36d9ca, 0xd2a828, 0x8d61c2, 0x77c912, 0x142604, 0x9b4612,
  58   0xc459c4, 0x44c5c8, 0x91b24d, 0xf31700, 0xad43d4, 0xe54929,
  59   0x10d5fd, 0xfcbe00, 0xcc941e, 0xeece70, 0xf53e13, 0x80f1ec,
  60   0xc3e7b3, 0x28f8c7, 0x940593, 0x3e71c1, 0xb3092e, 0xf3450b,
  61   0x9c1288, 0x7b20ab, 0x9fb52e, 0xc29247, 0x2f327b, 0x6d550c,
  62   0x90a772, 0x1fe76b, 0x96cb31, 0x4a1679, 0xe27941, 0x89dff4,
  63   0x9794e8, 0x84e6e2, 0x973199, 0x6bed88, 0x365f5f, 0x0efdbb,
  64   0xb49a48, 0x6ca467, 0x427271, 0x325d8d, 0xb8159f, 0x09e5bc,
  65   0x25318d, 0x3974f7, 0x1c0530, 0x010c0d, 0x68084b, 0x58ee2c,
  66   0x90aa47, 0x02e774, 0x24d6bd, 0xa67df7, 0x72486e, 0xef169f,
  67   0xa6948e, 0xf691b4, 0x5153d1, 0xf20acf, 0x339820, 0x7e4bf5,
  68   0x6863b2, 0x5f3edd, 0x035d40, 0x7f8985, 0x295255, 0xc06437,
  69   0x10d86d, 0x324832, 0x754c5b, 0xd4714e, 0x6e5445, 0xc1090b,
  70   0x69f52a, 0xd56614, 0x9d0727, 0x50045d, 0xdb3bb4, 0xc576ea,
  71   0x17f987, 0x7d6b49, 0xba271d, 0x296996, 0xacccc6, 0x5414ad,
  72   0x6ae290, 0x89d988, 0x50722c, 0xbea404, 0x940777, 0x7030f3,
  73   0x27fc00, 0xa871ea, 0x49c266, 0x3de064, 0x83dd97, 0x973fa3,
  74   0xfd9443, 0x8c860d, 0xde4131, 0x9d3992, 0x8c70dd, 0xe7b717,
  75   0x3bdf08, 0x2b3715, 0xa0805c, 0x93805a, 0x921110, 0xd8e80f,
  76   0xaf806c, 0x4bffdb, 0x0f9038, 0x761859, 0x15a562, 0xbbcb61,
  77   0xb989c7, 0xbd4010, 0x04f2d2, 0x277549, 0xf6b6eb, 0xbb22db,
  78   0xaa140a, 0x2f2689, 0x768364, 0x333b09, 0x1a940e, 0xaa3a51,
  79   0xc2a31d, 0xaeedaf, 0x12265c, 0x4dc26d, 0x9c7a2d, 0x9756c0,
  80   0x833f03, 0xf6f009, 0x8c402b, 0x99316d, 0x07b439, 0x15200c,
  81   0x5bc3d8, 0xc492f5, 0x4badc6, 0xa5ca4e, 0xcd37a7, 0x36a9e6,
  82   0x9492ab, 0x6842dd, 0xde6319, 0xef8c76, 0x528b68, 0x37dbfc,
  83   0xaba1ae, 0x3115df, 0xa1ae00, 0xdafb0c, 0x664d64, 0xb705ed,
  84   0x306529, 0xbf5657, 0x3aff47, 0xb9f96a, 0xf3be75, 0xdf9328,
  85   0x3080ab, 0xf68c66, 0x15cb04, 0x0622fa, 0x1de4d9, 0xa4b33d,
  86   0x8f1b57, 0x09cd36, 0xe9424e, 0xa4be13, 0xb52333, 0x1aaaf0,
  87   0xa8654f, 0xa5c1d2, 0x0f3f0b, 0xcd785b, 0x76f923, 0x048b7b,
  88   0x721789, 0x53a6c6, 0xe26e6f, 0x00ebef, 0x584a9b, 0xb7dac4,
  89   0xba66aa, 0xcfcf76, 0x1d02d1, 0x2df1b1, 0xc1998c, 0x77adc3,
  90   0xda4886, 0xa05df7, 0xf480c6, 0x2ff0ac, 0x9aecdd, 0xbc5c3f,
  91   0x6dded0, 0x1fc790, 0xb6db2a, 0x3a25a3, 0x9aaf00, 0x9353ad,
  92   0x0457b6, 0xb42d29, 0x7e804b, 0xa707da, 0x0eaa76, 0xa1597b,
  93   0x2a1216, 0x2db7dc, 0xfde5fa, 0xfedb89, 0xfdbe89, 0x6c76e4,
  94   0xfca906, 0x70803e, 0x156e85, 0xff87fd, 0x073e28, 0x336761,
  95   0x86182a, 0xeabd4d, 0xafe7b3, 0x6e6d8f, 0x396795, 0x5bbf31,
  96   0x48d784, 0x16df30, 0x432dc7, 0x356125, 0xce70c9, 0xb8cb30,
  97   0xfd6cbf, 0xa200a4, 0xe46c05, 0xa0dd5a, 0x476f21, 0xd21262,
  98   0x845cb9, 0x496170, 0xe0566b, 0x015299, 0x375550, 0xb7d51e,
  99   0xc4f133, 0x5f6e13, 0xe4305d, 0xa92e85, 0xc3b21d, 0x3632a1,
 100   0xa4b708, 0xd4b1ea, 0x21f716, 0xe4698f, 0x77ff27, 0x80030c,
 101   0x2d408d, 0xa0cd4f, 0x99a520, 0xd3a2b3, 0x0a5d2f, 0x42f9b4,
 102   0xcbda11, 0xd0be7d, 0xc1db9b, 0xbd17ab, 0x81a2ca, 0x5c6a08,
 103   0x17552e, 0x550027, 0xf0147f, 0x8607e1, 0x640b14, 0x8d4196,
 104   0xdebe87, 0x2afdda, 0xb6256b, 0x34897b, 0xfef305, 0x9ebfb9,
 105   0x4f6a68, 0xa82a4a, 0x5ac44f, 0xbcf82d, 0x985ad7, 0x95c7f4,
 106   0x8d4d0d, 0xa63a20, 0x5f57a4, 0xb13f14, 0x953880, 0x0120cc,
 107   0x86dd71, 0xb6dec9, 0xf560bf, 0x11654d, 0x6b0701, 0xacb08c,
 108   0xd0c0b2, 0x485551, 0x0efb1e, 0xc37295, 0x3b06a3, 0x3540c0,
 109   0x7bdc06, 0xcc45e0, 0xfa294e, 0xc8cad6, 0x41f3e8, 0xde647c,
 110   0xd8649b, 0x31bed9, 0xc397a4, 0xd45877, 0xc5e369, 0x13daf0,
 111   0x3c3aba, 0x461846, 0x5f7555, 0xf5bdd2, 0xc6926e, 0x5d2eac,
 112   0xed440e, 0x423e1c, 0x87c461, 0xe9fd29, 0xf3d6e7, 0xca7c22,
 113   0x35916f, 0xc5e008, 0x8dd7ff, 0xe26a6e, 0xc6fdb0, 0xc10893,
 114   0x745d7c, 0xb2ad6b, 0x9d6ecd, 0x7b723e, 0x6a11c6, 0xa9cff7,
 115   0xdf7329, 0xbac9b5, 0x5100b7, 0x0db2e2, 0x24ba74, 0x607de5,
 116   0x8ad874, 0x2c150d, 0x0c1881, 0x94667e, 0x162901, 0x767a9f,
 117   0xbefdfd, 0xef4556, 0x367ed9, 0x13d9ec, 0xb9ba8b, 0xfc97c4,
 118   0x27a831, 0xc36ef1, 0x36c594, 0x56a8d8, 0xb5a8b4, 0x0ecccf,
 119   0x2d8912, 0x34576f, 0x89562c, 0xe3ce99, 0xb920d6, 0xaa5e6b,
 120   0x9c2a3e, 0xcc5f11, 0x4a0bfd, 0xfbf4e1, 0x6d3b8e, 0x2c86e2,
 121   0x84d4e9, 0xa9b4fc, 0xd1eeef, 0xc9352e, 0x61392f, 0x442138,
 122   0xc8d91b, 0x0afc81, 0x6a4afb, 0xd81c2f, 0x84b453, 0x8c994e,
 123   0xcc2254, 0xdc552a, 0xd6c6c0, 0x96190b, 0xb8701a, 0x649569,
 124   0x605a26, 0xee523f, 0x0f117f, 0x11b5f4, 0xf5cbfc, 0x2dbc34,
 125   0xeebc34, 0xcc5de8, 0x605edd, 0x9b8e67, 0xef3392, 0xb817c9,
 126   0x9b5861, 0xbc57e1, 0xc68351, 0x103ed8, 0x4871dd, 0xdd1c2d,
 127   0xa118af, 0x462c21, 0xd7f359, 0x987ad9, 0xc0549e, 0xfa864f,
 128   0xfc0656, 0xae79e5, 0x362289, 0x22ad38, 0xdc9367, 0xaae855,
 129   0x382682, 0x9be7ca, 0xa40d51, 0xb13399, 0x0ed7a9, 0x480569,
 130   0xf0b265, 0xa7887f, 0x974c88, 0x36d1f9, 0xb39221, 0x4a827b,
 131   0x21cf98, 0xdc9f40, 0x5547dc, 0x3a74e1, 0x42eb67, 0xdf9dfe,
 132   0x5fd45e, 0xa4677b, 0x7aacba, 0xa2f655, 0x23882b, 0x55ba41,
 133   0x086e59, 0x862a21, 0x834739, 0xe6e389, 0xd49ee5, 0x40fb49,
 134   0xe956ff, 0xca0f1c, 0x8a59c5, 0x2bfa94, 0xc5c1d3, 0xcfc50f,
 135   0xae5adb, 0x86c547, 0x624385, 0x3b8621, 0x94792c, 0x876110,
 136   0x7b4c2a, 0x1a2c80, 0x12bf43, 0x902688, 0x893c78, 0xe4c4a8,
 137   0x7bdbe5, 0xc23ac4, 0xeaf426, 0x8a67f7, 0xbf920d, 0x2ba365,
 138   0xb1933d, 0x0b7cbd, 0xdc51a4, 0x63dd27, 0xdde169, 0x19949a,
 139   0x9529a8, 0x28ce68, 0xb4ed09, 0x209f44, 0xca984e, 0x638270,
 140   0x237c7e, 0x32b90f, 0x8ef5a7, 0xe75614, 0x08f121, 0x2a9db5,
 141   0x4d7e6f, 0x5119a5, 0xabf9b5, 0xd6df82, 0x61dd96, 0x023616,
 142   0x9f3ac4, 0xa1a283, 0x6ded72, 0x7a8d39, 0xa9b882, 0x5c326b,
 143   0x5b2746, 0xed3400, 0x7700d2, 0x55f4fc, 0x4d5901, 0x8071e0,
 144   0xe13f89, 0xb295f3, 0x64a8f1, 0xaea74b, 0x38fc4c, 0xeab2bb,
 145   0x47270b, 0xabc3a7, 0x34ba60, 0x52dd34, 0xf8563a, 0xeb7e8a,
 146   0x31bb36, 0x5895b7, 0x47f7a9, 0x94c3aa, 0xd39225, 0x1e7f3e,
 147   0xd8974e, 0xbba94f, 0xd8ae01, 0xe661b4, 0x393d8e, 0xa523aa,
 148   0x33068e, 0x1633b5, 0x3bb188, 0x1d3a9d, 0x4013d0, 0xcc1be5,
 149   0xf862e7, 0x3bf28f, 0x39b5bf, 0x0bc235, 0x22747e, 0xa247c0,
 150   0xd52d1f, 0x19add3, 0x9094df, 0x9311d0, 0xb42b25, 0x496db2,
 151   0xe264b2, 0x5ef135, 0x3bc6a4, 0x1a4ad0, 0xaac92e, 0x64e886,
 152   0x573091, 0x982cfb, 0x311b1a, 0x08728b, 0xbdcee1, 0x60e142,
 153   0xeb641d, 0xd0bba3, 0xe559d4, 0x597b8c, 0x2a4483, 0xf332ba,
 154   0xf84867, 0x2c8d1b, 0x2fa9b0, 0x50f3dd, 0xf9f573, 0xdb61b4,
 155   0xfe233e, 0x6c41a6, 0xeea318, 0x775a26, 0xbc5e5c, 0xcea708,
 156   0x94dc57, 0xe20196, 0xf1e839, 0xbe4851, 0x5d2d2f, 0x4e9555,
 157   0xd96ec2, 0xe7d755, 0x6304e0, 0xc02e0e, 0xfc40a0, 0xbbf9b3,
 158   0x7125a7, 0x222dfb, 0xf619d8, 0x838c1c, 0x6619e6, 0xb20d55,
 159   0xbb5137, 0x79e809, 0xaf9149, 0x0d73de, 0x0b0da5, 0xce7f58,
 160   0xac1934, 0x724667, 0x7a1a13, 0x9e26bc, 0x4555e7, 0x585cb5,
 161   0x711d14, 0x486991, 0x480d60, 0x56adab, 0xd62f64, 0x96ee0c,
 162   0x212ff3, 0x5d6d88, 0xa67684, 0x95651e, 0xab9e0a, 0x4ddefe,
 163   0x571010, 0x836a39, 0xf8ea31, 0x9e381d, 0xeac8b1, 0xcac96b,
 164   0x37f21e, 0xd505e9, 0x984743, 0x9fc56c, 0x0331b7, 0x3b8bf8,
 165   0x86e56a, 0x8dc343, 0x6230e7, 0x93cfd5, 0x6a8f2d, 0x733005,
 166   0x1af021, 0xa09fcb, 0x7415a1, 0xd56b23, 0x6ff725, 0x2f4bc7,
 167   0xb8a591, 0x7fac59, 0x5c55de, 0x212c38, 0xb13296, 0x5cff50,
 168   0x366262, 0xfa7b16, 0xf4d9a6, 0x2acfe7, 0xf07403, 0xd4d604,
 169   0x6fd916, 0x31b1bf, 0xcbb450, 0x5bd7c8, 0x0ce194, 0x6bd643,
 170   0x4fd91c, 0xdf4543, 0x5f3453, 0xe2b5aa, 0xc9aec8, 0x131485,
 171   0xf9d2bf, 0xbadb9e, 0x76f5b9, 0xaf15cf, 0xca3182, 0x14b56d,
 172   0xe9fe4d, 0x50fc35, 0xf5aed5, 0xa2d0c1, 0xc96057, 0x192eb6,
 173   0xe91d92, 0x07d144, 0xaea3c6, 0x343566, 0x26d5b4, 0x3161e2,
 174   0x37f1a2, 0x209eff, 0x958e23, 0x493798, 0x35f4a6, 0x4bdc02,
 175   0xc2be13, 0xbe80a0, 0x0b72a3, 0x115c5f, 0x1e1bd1, 0x0db4d3,
 176   0x869e85, 0x96976b, 0x2ac91f, 0x8a26c2, 0x3070f0, 0x041412,
 177   0xfc9fa5, 0xf72a38, 0x9c6878, 0xe2aa76, 0x50cfe1, 0x559274,
 178   0x934e38, 0x0a92f7, 0x5533f0, 0xa63db4, 0x399971, 0xe2b755,
 179   0xa98a7c, 0x008f19, 0xac54d2, 0x2ea0b4, 0xf5f3e0, 0x60c849,
 180   0xffd269, 0xae52ce, 0x7a5fdd, 0xe9ce06, 0xfb0ae8, 0xa50cce,
 181   0xea9d3e, 0x3766dd, 0xb834f5, 0x0da090, 0x846f88, 0x4ae3d5,
 182   0x099a03, 0x2eae2d, 0xfcb40a, 0xfb9b33, 0xe281dd, 0x1b16ba,
 183   0xd8c0af, 0xd96b97, 0xb52dc9, 0x9c277f, 0x5951d5, 0x21ccd6,
 184   0xb6496b, 0x584562, 0xb3baf2, 0xa1a5c4, 0x7ca2cf, 0xa9b93d,
 185   0x7b7b89, 0x483d38,
 186 };
 187
 188 static const long double c[] = {
 189 /* 93 bits of pi/2 */
 190 #define PI_2_1 c[0]
 191   1.57079632679489661923132169155131424e+00L,   /* 3fff921fb54442d18469898cc5100000 */
 192
 193 /* pi/2 - PI_2_1 */
 194 #define PI_2_1t c[1]
 195   8.84372056613570112025531863263659260e-29L,   /* 3fa1c06e0e68948127044533e63a0106 */
 196 };
 197
 198 static int kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec,
 199                             const int *ipio2);
 200
 201 int
 202 ieee754_rem_pio2l (long double x, long double *y)
 203 {
 204   long double z, w, t;
 205   double tx[8];
 206   int exp, n;
 207
 208   if (x >= -0.78539816339744830961566084581987572104929234984377
 209       && x < 0.78539816339744830961566084581987572104929234984377)
 210     /* x in <-pi/4, pi/4> */
 211     {
 212       y[0] = x;
 213       y[1] = 0;
 214       return 0;
 215     }
 216
 217   if (x >= 2.35619449019234492884698253745962716314787704953131
 218       && x < 2.35619449019234492884698253745962716314787704953131)
 219     if (x > 0)
 220       {
 221         /* 113 + 93 bit PI is ok */
 222         z = x - PI_2_1;
 223         y[0] = z - PI_2_1t;
 224         y[1] = (z - y[0]) - PI_2_1t;
 225         return 1;
 226       }
 227     else
 228       {
 229         /* 113 + 93 bit PI is ok */
 230         z = x + PI_2_1;
 231         y[0] = z + PI_2_1t;
 232         y[1] = (z - y[0]) + PI_2_1t;
 233         return -1;
 234       }
 235
 236   if (x + x == x || x != x)     /* x is +=oo or NaN */
 237     {
 238       y[0] = x - x;
 239       y[1] = y[0];
 240       return 0;
 241     }
 242
 243   /* Handle large arguments.
 244      We split the 113 bits of the mantissa into 5 24bit integers
 245      stored in a double array.  */
 246   z = frexp (x, &exp);
 247   tx[0] = floorl (z *= 16777216.0);
 248   z -= tx[0];
 249   tx[1] = floorl (z *= 16777216.0);
 250   z -= tx[1];
 251   tx[2] = floorl (z *= 16777216.0);
 252   z -= tx[2];
 253   tx[3] = floorl (z *= 16777216.0);
 254   z -= tx[3];
 255   tx[4] = floorl (z *= 16777216.0);
 256
 257   n = kernel_rem_pio2 (tx, tx + 5, exp - 24, tx[4] ? 5 : 4, 3, two_over_pi);
 258
 259   /* The result is now stored in 3 double values, we need to convert it into
 260      two long double values.  */
 261   t = (long double) tx[6] + (long double) tx[7];
 262   w = (long double) tx[5];
 263
 264   if (x > 0)
 265     {
 266       y[0] = w + t;
 267       y[1] = t - (y[0] - w);
 268       return n;
 269     }
 270   else
 271     {
 272       y[0] = -(w + t);
 273       y[1] = -t - (y[0] + w);
 274       return -n;
 275     }
 276 }
 277
 278 /* @(#)k_rem_pio2.c 5.1 93/09/24 */
 279 /*
 280  * ====================================================
 281  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 282  *
 283  * Developed at SunPro, a Sun Microsystems, Inc. business.
 284  * Permission to use, copy, modify, and distribute this
 285  * software is freely granted, provided that this notice
 286  * is preserved.
 287  * ====================================================
 288  */
 289
 290 #if defined(LIBM_SCCS) && !defined(lint)
 291 static char rcsid[] =
 292   "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $";
 293 #endif
 294
 295 /*
 296  * kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
 297  * double x[],y[]; int e0,nx,prec; int ipio2[];
 298  *
 299  * kernel_rem_pio2 return the last three digits of N with
 300  *              y = x - N*pi/2
 301  * so that |y| < pi/2.
 302  *
 303  * The method is to compute the integer (mod 8) and fraction parts of
 304  * (2/pi)*x without doing the full multiplication. In general we
 305  * skip the part of the product that are known to be a huge integer (
 306  * more accurately, = 0 mod 8 ). Thus the number of operations are
 307  * independent of the exponent of the input.
 308  *
 309  * (2/pi) is represented by an array of 24-bit integers in ipio2[].
 310  *
 311  * Input parameters:
 312  *      x[]     The input value (must be positive) is broken into nx
 313  *              pieces of 24-bit integers in double precision format.
 314  *              x[i] will be the i-th 24 bit of x. The scaled exponent
 315  *              of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
 316  *              match x's up to 24 bits.
 317  *
 318  *              Example of breaking a double positive z into x[0]+x[1]+x[2]:
 319  *                      e0 = ilogb(z)-23
 320  *                      z  = scalbn(z,-e0)
 321  *              for i = 0,1,2
 322  *                      x[i] = floor(z)
 323  *                      z    = (z-x[i])*2**24
 324  *
 325  *
 326  *      y[]     ouput result in an array of double precision numbers.
 327  *              The dimension of y[] is:
 328  *                      24-bit  precision       1
 329  *                      53-bit  precision       2
 330  *                      64-bit  precision       2
 331  *                      113-bit precision       3
 332  *              The actual value is the sum of them. Thus for 113-bit
 333  *              precision, one may have to do something like:
 334  *
 335  *              long double t,w,r_head, r_tail;
 336  *              t = (long double)y[2] + (long double)y[1];
 337  *              w = (long double)y[0];
 338  *              r_head = t+w;
 339  *              r_tail = w - (r_head - t);
 340  *
 341  *      e0      The exponent of x[0]
 342  *
 343  *      nx      dimension of x[]
 344  *
 345  *      prec    an integer indicating the precision:
 346  *                      0       24  bits (single)
 347  *                      1       53  bits (double)
 348  *                      2       64  bits (extended)
 349  *                      3       113 bits (quad)
 350  *
 351  *      ipio2[]
 352  *              integer array, contains the (24*i)-th to (24*i+23)-th
 353  *              bit of 2/pi after binary point. The corresponding
 354  *              floating value is
 355  *
 356  *                      ipio2[i] * 2^(-24(i+1)).
 357  *
 358  * External function:
 359  *      double scalbn(), floor();
 360  *
 361  *
 362  * Here is the description of some local variables:
 363  *
 364  *      jk      jk+1 is the initial number of terms of ipio2[] needed
 365  *              in the computation. The recommended value is 2,3,4,
 366  *              6 for single, double, extended,and quad.
 367  *
 368  *      jz      local integer variable indicating the number of
 369  *              terms of ipio2[] used.
 370  *
 371  *      jx      nx - 1
 372  *
 373  *      jv      index for pointing to the suitable ipio2[] for the
 374  *              computation. In general, we want
 375  *                      ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
 376  *              is an integer. Thus
 377  *                      e0-3-24*jv >= 0 or (e0-3)/24 >= jv
 378  *              Hence jv = max(0,(e0-3)/24).
 379  *
 380  *      jp      jp+1 is the number of terms in PIo2[] needed, jp = jk.
 381  *
 382  *      q[]     double array with integral value, representing the
 383  *              24-bits chunk of the product of x and 2/pi.
 384  *
 385  *      q0      the corresponding exponent of q[0]. Note that the
 386  *              exponent for q[i] would be q0-24*i.
 387  *
 388  *      PIo2[]  double precision array, obtained by cutting pi/2
 389  *              into 24 bits chunks.
 390  *
 391  *      f[]     ipio2[] in floating point
 392  *
 393  *      iq[]    integer array by breaking up q[] in 24-bits chunk.
 394  *
 395  *      fq[]    final product of x*(2/pi) in fq[0],..,fq[jk]
 396  *
 397  *      ih      integer. If >0 it indicates q[] is >= 0.5, hence
 398  *              it also indicates the *sign* of the result.
 399  *
 400  */
 401
 402
 403 /*
 404  * Constants:
 405  * The hexadecimal values are the intended ones for the following
 406  * constants. The decimal values may be used, provided that the
 407  * compiler will convert from decimal to binary accurately enough
 408  * to produce the hexadecimal values shown.
 409  */
 410
 411 static const int init_jk[] = { 2, 3, 4, 6 };    /* initial value for jk */
 412 static const double PIo2[] = {
 413   1.57079625129699707031e+00,   /* 0x3FF921FB, 0x40000000 */
 414   7.54978941586159635335e-08,   /* 0x3E74442D, 0x00000000 */
 415   5.39030252995776476554e-15,   /* 0x3CF84698, 0x80000000 */
 416   3.28200341580791294123e-22,   /* 0x3B78CC51, 0x60000000 */
 417   1.27065575308067607349e-29,   /* 0x39F01B83, 0x80000000 */
 418   1.22933308981111328932e-36,   /* 0x387A2520, 0x40000000 */
 419   2.73370053816464559624e-44,   /* 0x36E38222, 0x80000000 */
 420   2.16741683877804819444e-51,   /* 0x3569F31D, 0x00000000 */
 421 };
 422
 423 static const double zero = 0.0, one = 1.0, two24 = 1.67772160000000000000e+07,  /* 0x41700000, 0x00000000 */
 424   twon24 = 5.96046447753906250000e-08;  /* 0x3E700000, 0x00000000 */
 425
 426 static int
 427 kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec,
 428                  const int *ipio2)
 429 {
 430   int jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih;
 431   double z, fw, f[20], fq[20], q[20];
 432
 433   /* initialize jk */
 434   jk = init_jk[prec];
 435   jp = jk;
 436
 437   /* determine jx,jv,q0, note that 3>q0 */
 438   jx = nx - 1;
 439   jv = (e0 - 3) / 24;
 440   if (jv < 0)
 441     jv = 0;
 442   q0 = e0 - 24 * (jv + 1);
 443
 444   /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
 445   j = jv - jx;
 446   m = jx + jk;
 447   for (i = 0; i <= m; i++, j++)
 448     f[i] = (j < 0) ? zero : (double) ipio2[j];
 449
 450   /* compute q[0],q[1],...q[jk] */
 451   for (i = 0; i <= jk; i++)
 452     {
 453       for (j = 0, fw = 0.0; j <= jx; j++)
 454         fw += x[j] * f[jx + i - j];
 455       q[i] = fw;
 456     }
 457
 458   jz = jk;
 459 recompute:
 460   /* distill q[] into iq[] reversingly */
 461   for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--)
 462     {
 463       fw = (double) ((int) (twon24 * z));
 464       iq[i] = (int) (z - two24 * fw);
 465       z = q[j - 1] + fw;
 466     }
 467
 468   /* compute n */
 469   z = ldexp (z, q0);            /* actual value of z */
 470   z -= 8.0 * floor (z * 0.125); /* trim off integer >= 8 */
 471   n = (int) z;
 472   z -= (double) n;
 473   ih = 0;
 474   if (q0 > 0)
 475     {                           /* need iq[jz-1] to determine n */
 476       i = (iq[jz - 1] >> (24 - q0));
 477       n += i;
 478       iq[jz - 1] -= i << (24 - q0);
 479       ih = iq[jz - 1] >> (23 - q0);
 480     }
 481   else if (q0 == 0)
 482     ih = iq[jz - 1] >> 23;
 483   else if (z >= 0.5)
 484     ih = 2;
 485
 486   if (ih > 0)
 487     {                           /* q > 0.5 */
 488       n += 1;
 489       carry = 0;
 490       for (i = 0; i < jz; i++)
 491         {                       /* compute 1-q */
 492           j = iq[i];
 493           if (carry == 0)
 494             {
 495               if (j != 0)
 496                 {
 497                   carry = 1;
 498                   iq[i] = 0x1000000 - j;
 499                 }
 500             }
 501           else
 502             iq[i] = 0xffffff - j;
 503         }
 504       if (q0 > 0)
 505         {                       /* rare case: chance is 1 in 12 */
 506           switch (q0)
 507             {
 508             case 1:
 509               iq[jz - 1] &= 0x7fffff;
 510               break;
 511             case 2:
 512               iq[jz - 1] &= 0x3fffff;
 513               break;
 514             }
 515         }
 516       if (ih == 2)
 517         {
 518           z = one - z;
 519           if (carry != 0)
 520             z -= ldexp (one, q0);
 521         }
 522     }
 523
 524   /* check if recomputation is needed */
 525   if (z == zero)
 526     {
 527       j = 0;
 528       for (i = jz - 1; i >= jk; i--)
 529         j |= iq[i];
 530       if (j == 0)
 531         {                       /* need recomputation */
 532           for (k = 1; iq[jk - k] == 0; k++);    /* k = no. of terms needed */
 533
 534           for (i = jz + 1; i <= jz + k; i++)
 535             {                   /* add q[jz+1] to q[jz+k] */
 536               f[jx + i] = (double) ipio2[jv + i];
 537               for (j = 0, fw = 0.0; j <= jx; j++)
 538                 fw += x[j] * f[jx + i - j];
 539               q[i] = fw;
 540             }
 541           jz += k;
 542           goto recompute;
 543         }
 544     }
 545
 546   /* chop off zero terms */
 547   if (z == 0.0)
 548     {
 549       jz -= 1;
 550       q0 -= 24;
 551       while (iq[jz] == 0)
 552         {
 553           jz--;
 554           q0 -= 24;
 555         }
 556     }
 557   else
 558     {                           /* break z into 24-bit if necessary */
 559       z = ldexp (z, -q0);
 560       if (z >= two24)
 561         {
 562           fw = (double) ((int) (twon24 * z));
 563           iq[jz] = (int) (z - two24 * fw);
 564           jz += 1;
 565           q0 += 24;
 566           iq[jz] = (int) fw;
 567         }
 568       else
 569         iq[jz] = (int) z;
 570     }
 571
 572   /* convert integer "bit" chunk to floating-point value */
 573   fw = ldexp (one, q0);
 574   for (i = jz; i >= 0; i--)
 575     {
 576       q[i] = fw * (double) iq[i];
 577       fw *= twon24;
 578     }
 579
 580   /* compute PIo2[0,...,jp]*q[jz,...,0] */
 581   for (i = jz; i >= 0; i--)
 582     {
 583       for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++)
 584         fw += PIo2[k] * q[i + k];
 585       fq[jz - i] = fw;
 586     }
 587
 588   /* compress fq[] into y[] */
 589   switch (prec)
 590     {
 591     case 0:
 592       fw = 0.0;
 593       for (i = jz; i >= 0; i--)
 594         fw += fq[i];
 595       y[0] = (ih == 0) ? fw : -fw;
 596       break;
 597     case 1:
 598     case 2:
 599       fw = 0.0;
 600       for (i = jz; i >= 0; i--)
 601         fw += fq[i];
 602       y[0] = (ih == 0) ? fw : -fw;
 603       fw = fq[0] - fw;
 604       for (i = 1; i <= jz; i++)
 605         fw += fq[i];
 606       y[1] = (ih == 0) ? fw : -fw;
 607       break;
 608     case 3:                     /* painful */
 609       for (i = jz; i > 0; i--)
 610         {
 611           fw = fq[i - 1] + fq[i];
 612           fq[i] += fq[i - 1] - fw;
 613           fq[i - 1] = fw;
 614         }
 615       for (i = jz; i > 1; i--)
 616         {
 617           fw = fq[i - 1] + fq[i];
 618           fq[i] += fq[i - 1] - fw;
 619           fq[i - 1] = fw;
 620         }
 621       for (fw = 0.0, i = jz; i >= 2; i--)
 622         fw += fq[i];
 623       if (ih == 0)
 624         {
 625           y[0] = fq[0];
 626           y[1] = fq[1];
 627           y[2] = fw;
 628         }
 629       else
 630         {
 631           y[0] = -fq[0];
 632           y[1] = -fq[1];
 633           y[2] = -fw;
 634         }
 635     }
 636   return n & 7;
 637 }