lib/trigl.c

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