lib/trigl.c

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