Print this page
11210 libm should be cstyle(1ONBLD) clean
   1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */

  21 /*
  22  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  23  */

  24 /*
  25  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
  26  * Use is subject to license terms.
  27  */
  28 
  29 /*
  30  * __rem_pio2(x, y) passes back a better-than-double-precision
  31  * approximation to x mod pi/2 in y[0]+y[1] and returns an integer
  32  * congruent mod 8 to the integer part of x/(pi/2).
  33  *
  34  * This implementation tacitly assumes that x is finite and at
  35  * least about pi/4 in magnitude.
  36  */
  37 
  38 #include "libm.h"
  39 
  40 extern const int _TBL_ipio2_inf[];
  41 
  42 /* INDENT OFF */
  43 /*
  44  * invpio2:  53 bits of 2/pi
  45  * pio2_1:   first  33 bit of pi/2
  46  * pio2_1t:  pi/2 - pio2_1
  47  * pio2_2:   second 33 bit of pi/2
  48  * pio2_2t:  pi/2 - pio2_2
  49  * pio2_3:   third  33 bit of pi/2
  50  * pio2_3t:  pi/2 - pio2_3
  51  */
  52 static const double
  53         half    = 0.5,
  54         invpio2 = 0.636619772367581343075535,   /* 2^ -1  * 1.45F306DC9C883 */
  55         pio2_1  = 1.570796326734125614166,      /* 2^  0  * 1.921FB54400000 */
  56         pio2_1t = 6.077100506506192601475e-11,  /* 2^-34  * 1.0B4611A626331 */
  57         pio2_2  = 6.077100506303965976596e-11,  /* 2^-34  * 1.0B4611A600000 */
  58         pio2_2t = 2.022266248795950732400e-21,  /* 2^-69  * 1.3198A2E037073 */
  59         pio2_3  = 2.022266248711166455796e-21,  /* 2^-69  * 1.3198A2E000000 */
  60         pio2_3t = 8.478427660368899643959e-32;  /* 2^-104 * 1.B839A252049C1 */
  61 /* INDENT ON */
  62 
  63 int
  64 __rem_pio2(double x, double *y) {

  65         double  w, t, r, fn;
  66         double  tx[3];
  67         int     e0, i, j, nx, n, ix, hx, lx;
  68 
  69         hx = ((int *)&x)[HIWORD];
  70         ix = hx & 0x7fffffff;
  71 
  72         if (ix < 0x4002d97c) {
  73                 /* |x| < 3pi/4, special case with n=1 */
  74                 t = fabs(x) - pio2_1;

  75                 if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */
  76                         y[0] = t - pio2_1t;
  77                         y[1] = (t - y[0]) - pio2_1t;
  78                 } else {                /* near pi/2, use 33+33+53 bit pi */
  79                         t -= pio2_2;
  80                         y[0] = t - pio2_2t;
  81                         y[1] = (t - y[0]) - pio2_2t;
  82                 }

  83                 if (hx < 0) {
  84                         y[0] = -y[0];
  85                         y[1] = -y[1];
  86                         return (-1);
  87                 }

  88                 return (1);
  89         }
  90 
  91         if (ix <= 0x413921fb) {
  92                 /* |x| <= 2^19 pi */
  93                 t = fabs(x);
  94                 n = (int)(t * invpio2 + half);
  95                 fn = (double)n;
  96                 r = t - fn * pio2_1;
  97                 j = ix >> 20;
  98                 w = fn * pio2_1t;       /* 1st round good to 85 bit */
  99                 y[0] = r - w;
 100                 i = j - ((((int *)y)[HIWORD] >> 20) & 0x7ff);

 101                 if (i > 16) {        /* 2nd iteration (rare) */
 102                         /* 2nd round good to 118 bit */
 103                         if (i < 35) {
 104                                 t = r;  /* r-fn*pio2_2 may not be exact */
 105                                 w = fn * pio2_2;
 106                                 r = t - w;
 107                                 w = fn * pio2_2t - ((t - r) - w);
 108                                 y[0] = r - w;
 109                         } else {
 110                                 r -= fn * pio2_2;
 111                                 w = fn * pio2_2t;
 112                                 y[0] = r - w;
 113                                 i = j - ((((int *)y)[HIWORD] >> 20) & 0x7ff);

 114                                 if (i > 49) {
 115                                         /* 3rd iteration (extremely rare) */
 116                                         if (i < 68) {
 117                                                 t = r;
 118                                                 w = fn * pio2_3;
 119                                                 r = t - w;
 120                                                 w = fn * pio2_3t -
 121                                                     ((t - r) - w);
 122                                                 y[0] = r - w;
 123                                         } else {
 124                                                 /*
 125                                                  * 3rd round good to 151 bits;
 126                                                  * covered all possible cases
 127                                                  */
 128                                                 r -= fn * pio2_3;
 129                                                 w = fn * pio2_3t;
 130                                                 y[0] = r - w;
 131                                         }
 132                                 }
 133                         }
 134                 }

 135                 y[1] = (r - y[0]) - w;

 136                 if (hx < 0) {
 137                         y[0] = -y[0];
 138                         y[1] = -y[1];
 139                         return (-n);
 140                 }

 141                 return (n);
 142         }
 143 
 144         e0 = (ix >> 20) - 1046;   /* e0 = ilogb(x)-23; */
 145 
 146         /* break x into three 24 bit pieces */
 147         lx = ((int *)&x)[LOWORD];
 148         i = (lx & 0x1f) << 19;
 149         tx[2] = (double)i;
 150         j = (lx >> 5) & 0xffffff;
 151         tx[1] = (double)j;
 152         tx[0] = (double)((((ix & 0xfffff) | 0x100000) << 3) |
 153             ((unsigned)lx >> 29));
 154         nx = 3;

 155         if (i == 0) {
 156                 /* skip zero term */
 157                 nx--;

 158                 if (j == 0)
 159                         nx--;
 160         }

 161         n = __rem_pio2m(tx, y, e0, nx, 2, _TBL_ipio2_inf);

 162         if (hx < 0) {
 163                 y[0] = -y[0];
 164                 y[1] = -y[1];
 165                 return (-n);
 166         }

 167         return (n);
 168 }
   1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 
  26 /*
  27  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 /*
  32  * __rem_pio2(x, y) passes back a better-than-double-precision
  33  * approximation to x mod pi/2 in y[0]+y[1] and returns an integer
  34  * congruent mod 8 to the integer part of x/(pi/2).
  35  *
  36  * This implementation tacitly assumes that x is finite and at
  37  * least about pi/4 in magnitude.
  38  */
  39 
  40 #include "libm.h"
  41 
  42 extern const int _TBL_ipio2_inf[];
  43 
  44 
  45 /*
  46  * invpio2:  53 bits of 2/pi
  47  * pio2_1:   first  33 bit of pi/2
  48  * pio2_1t:  pi/2 - pio2_1
  49  * pio2_2:   second 33 bit of pi/2
  50  * pio2_2t:  pi/2 - pio2_2
  51  * pio2_3:   third  33 bit of pi/2
  52  * pio2_3t:  pi/2 - pio2_3
  53  */
  54 static const double half = 0.5,

  55         invpio2 = 0.636619772367581343075535,   /* 2^ -1  * 1.45F306DC9C883 */
  56         pio2_1 = 1.570796326734125614166,       /* 2^  0  * 1.921FB54400000 */
  57         pio2_1t = 6.077100506506192601475e-11,  /* 2^-34  * 1.0B4611A626331 */
  58         pio2_2 = 6.077100506303965976596e-11,   /* 2^-34  * 1.0B4611A600000 */
  59         pio2_2t = 2.022266248795950732400e-21,  /* 2^-69  * 1.3198A2E037073 */
  60         pio2_3 = 2.022266248711166455796e-21,   /* 2^-69  * 1.3198A2E000000 */
  61         pio2_3t = 8.478427660368899643959e-32;  /* 2^-104 * 1.B839A252049C1 */
  62 
  63 
  64 int
  65 __rem_pio2(double x, double *y)
  66 {
  67         double w, t, r, fn;
  68         double tx[3];
  69         int e0, i, j, nx, n, ix, hx, lx;
  70 
  71         hx = ((int *)&x)[HIWORD];
  72         ix = hx & 0x7fffffff;
  73 
  74         if (ix < 0x4002d97c) {
  75                 /* |x| < 3pi/4, special case with n=1 */
  76                 t = fabs(x) - pio2_1;
  77 
  78                 if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */
  79                         y[0] = t - pio2_1t;
  80                         y[1] = (t - y[0]) - pio2_1t;
  81                 } else {                /* near pi/2, use 33+33+53 bit pi */
  82                         t -= pio2_2;
  83                         y[0] = t - pio2_2t;
  84                         y[1] = (t - y[0]) - pio2_2t;
  85                 }
  86 
  87                 if (hx < 0) {
  88                         y[0] = -y[0];
  89                         y[1] = -y[1];
  90                         return (-1);
  91                 }
  92 
  93                 return (1);
  94         }
  95 
  96         if (ix <= 0x413921fb) {
  97                 /* |x| <= 2^19 pi */
  98                 t = fabs(x);
  99                 n = (int)(t * invpio2 + half);
 100                 fn = (double)n;
 101                 r = t - fn * pio2_1;
 102                 j = ix >> 20;
 103                 w = fn * pio2_1t;       /* 1st round good to 85 bit */
 104                 y[0] = r - w;
 105                 i = j - ((((int *)y)[HIWORD] >> 20) & 0x7ff);
 106 
 107                 if (i > 16) {                /* 2nd iteration (rare) */
 108                         /* 2nd round good to 118 bit */
 109                         if (i < 35) {
 110                                 t = r;  /* r-fn*pio2_2 may not be exact */
 111                                 w = fn * pio2_2;
 112                                 r = t - w;
 113                                 w = fn * pio2_2t - ((t - r) - w);
 114                                 y[0] = r - w;
 115                         } else {
 116                                 r -= fn * pio2_2;
 117                                 w = fn * pio2_2t;
 118                                 y[0] = r - w;
 119                                 i = j - ((((int *)y)[HIWORD] >> 20) & 0x7ff);
 120 
 121                                 if (i > 49) {
 122                                         /* 3rd iteration (extremely rare) */
 123                                         if (i < 68) {
 124                                                 t = r;
 125                                                 w = fn * pio2_3;
 126                                                 r = t - w;
 127                                                 w = fn * pio2_3t - ((t - r) -
 128                                                     w);
 129                                                 y[0] = r - w;
 130                                         } else {
 131                                                 /*
 132                                                  * 3rd round good to 151 bits;
 133                                                  * covered all possible cases
 134                                                  */
 135                                                 r -= fn * pio2_3;
 136                                                 w = fn * pio2_3t;
 137                                                 y[0] = r - w;
 138                                         }
 139                                 }
 140                         }
 141                 }
 142 
 143                 y[1] = (r - y[0]) - w;
 144 
 145                 if (hx < 0) {
 146                         y[0] = -y[0];
 147                         y[1] = -y[1];
 148                         return (-n);
 149                 }
 150 
 151                 return (n);
 152         }
 153 
 154         e0 = (ix >> 20) - 1046;           /* e0 = ilogb(x)-23; */
 155 
 156         /* break x into three 24 bit pieces */
 157         lx = ((int *)&x)[LOWORD];
 158         i = (lx & 0x1f) << 19;
 159         tx[2] = (double)i;
 160         j = (lx >> 5) & 0xffffff;
 161         tx[1] = (double)j;
 162         tx[0] = (double)((((ix & 0xfffff) | 0x100000) << 3) | ((unsigned)lx >>
 163             29));
 164         nx = 3;
 165 
 166         if (i == 0) {
 167                 /* skip zero term */
 168                 nx--;
 169 
 170                 if (j == 0)
 171                         nx--;
 172         }
 173 
 174         n = __rem_pio2m(tx, y, e0, nx, 2, _TBL_ipio2_inf);
 175 
 176         if (hx < 0) {
 177                 y[0] = -y[0];
 178                 y[1] = -y[1];
 179                 return (-n);
 180         }
 181 
 182         return (n);
 183 }