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