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 } 169