__vrem_pio2m.c revision 25c28e83beb90e7c80452a7c818c5e6f73a07dc8
/*
* CDDL HEADER START
*
* The contents of this file are subject to the terms of the
* Common Development and Distribution License (the "License").
* You may not use this file except in compliance with the License.
*
* You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
* See the License for the specific language governing permissions
* and limitations under the License.
*
* When distributing Covered Code, include this CDDL HEADER in each
* file and include the License file at usr/src/OPENSOLARIS.LICENSE.
* If applicable, add the following below this CDDL HEADER, with the
* fields enclosed by brackets "[]" replaced with your own identifying
* information: Portions Copyright [yyyy] [name of copyright owner]
*
* CDDL HEADER END
*/
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
/*
* Copyright 2006 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
/*
* Given X, __vlibm_rem_pio2m finds Y and an integer n such that
* Y = X - n*pi/2 and |Y| < pi/2.
*
* On entry, X is represented by x, an array of nx 24-bit integers
* stored in double precision format, and e:
*
* X = sum (x[i] * 2^(e - 24*i))
*
* nx must be 1, 2, or 3, and e must be >= -24. For example, a
* suitable representation for the double precision number z can
* be computed as follows:
*
* e = ilogb(z)-23
* z = scalbn(z,-e)
* for i = 0,1,2
* x[i] = floor(z)
* z = (z-x[i])*2**24
*
* On exit, Y is approximated by y[0] if prec is 0 and by the un-
* evaluated sum y[0] + y[1] if prec != 0. The approximation is
* accurate to 53 bits in the former case and to at least 72 bits
* in the latter.
*
* __vlibm_rem_pio2m returns n mod 8.
*
* Notes:
*
* As n is the integer nearest X * 2/pi, we approximate the latter
* product to a precision that is determined dynamically so as to
* ensure that the final value Y is approximated accurately enough.
* We don't bother to compute terms in the product that are multiples
* of 8, so the cost of this multiplication is independent of the
* magnitude of X. The variable ip determines the offset into the
* array ipio2 of the first term we need to use. The variable eq0
* is the corresponding exponent of the first partial product.
*
* The partial products are scaled, summed, and split into an array
* of non-overlapping 24-bit terms (not necessarily having the same
* signs). Each partial product overlaps three elements of the
* resulting array:
*
* q[i] xxxxxxxxxxxxxx
* q[i+1] xxxxxxxxxxxxxx
* q[i+2] xxxxxxxxxxxxxx
* ... ...
*
*
* r[i] xxxxxx
* r[i+1] xxxxxx
* r[i+2] xxxxxx
* ... ...
*
* In order that the last element of the r array have some correct
* bits, we compute an extra term in the q array, but we don't bother
* to split this last term into 24-bit chunks; thus, the final term
* of the r array could have more than 24 bits, but this doesn't
* matter.
*
* After we subtract the nearest integer to the product, we multiply
* the remaining part of r by pi/2 to obtain Y. Before we compute
* this last product, however, we make sure that the remaining part
* of r has at least five nonzero terms, computing more if need be.
* This ensures that even if the first nonzero term is only a single
* bit and the last term is wrong in several trailing bits, we still
* have enough accuracy to obtain 72 bits of Y.
*
* IMPORTANT: This code assumes that the rounding mode is round-to-
* nearest in several key places. First, after we compute X * 2/pi,
* we round to the nearest integer by adding and subtracting a power
* of two. This step must be done in round-to-nearest mode to ensure
* that the remainder is less than 1/2 in absolute value. (Because
* we only take two adjacent terms of r into account when we perform
* this rounding, in very rare cases the remainder could be just
* barely greater than 1/2, but this shouldn't matter in practice.)
*
* Second, we also split the partial products of X * 2/pi into 24-bit
* pieces by adding and subtracting a power of two. In this step,
* round-to-nearest mode is important in order to guarantee that
* the index of the first nonzero term in the remainder gives an
* accurate indication of the number of significant terms. For
* example, suppose eq0 = -1, so that r[1] is a multiple of 1/2 and
* |r[2]| < 1/2. After we subtract the nearest integer, r[1] could
* be -1/2, and r[2] could be very nearly 1/2, so that r[1] != 0,
* yet the remainder is much smaller than the least significant bit
* corresponding to r[1]. As long as we use round-to-nearest mode,
* this can't happen; instead, the absolute value of each r[j] will
* be less than 1/2 the least significant bit corresponding to r[j-1],
* so that the entire remainder must be at least half as large as
* the first nonzero term (or perhaps just barely smaller than this).
*/
#include <sys/isa_defs.h>
#ifdef _LITTLE_ENDIAN
#define HIWORD 1
#define LOWORD 0
#else
#define HIWORD 0
#define LOWORD 1
#endif
/* 396 hex digits of 2/pi, with two leading zeroes to make life easier */
static const double ipio2[] = {
0, 0,
0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
};
/* pi/2 in 24-bit pieces */
static const double pio2[] = {
1.57079625129699707031e+00,
7.54978941586159635335e-08,
5.39030252995776476554e-15,
3.28200341580791294123e-22,
1.27065575308067607349e-29,
};
/* miscellaneous constants */
static const double
zero = 0.0,
two24 = 16777216.0,
twon24 = 5.960464477539062500E-8;
int
{
union {
double d;
int i[2];
} s;
/* determine ip and eq0; note that -48 <= eq0 <= 2 */
if (ip < 0)
ip = 0;
/* compute q[0,...,5] = x * ipio2 and initialize nq and eqnqm1 */
if (nx == 3) {
} else if (nx == 2) {
} else {
}
nq = 5;
/* propagate carries and incorporate powers of two */
s.i[LOWORD] = 0;
p = s.d;
z += q[j];
r[j+1] = (z - t) * p;
z = t * twon24;
p *= two24;
}
z += q[0];
r[1] = (z - t) * p;
r[0] = t * p;
/* form n = [r] mod 8 and leave the fractional part of r */
if (eq0 > 0) {
/* binary point lies within r[2] */
z = r[2] + r[3];
r[2] -= t;
n = (int)(r[1] + t);
r[0] = r[1] = zero;
} else if (eq0 > -24) {
/* binary point lies within or just to the right of r[1] */
z = r[1] + r[2];
r[1] -= t;
z = r[0] + t;
/* cut off high part of z so conversion to int doesn't
overflow */
n = (int)(z - t);
r[0] = zero;
} else {
/* binary point lies within or just to the right of r[0] */
z = r[0] + r[1];
r[0] -= t;
n = (int)t;
}
/* count the number of leading zeroes in r */
for (j = 0; j <= nq; j++) {
if (r[j] != zero)
break;
}
/* if fewer than 5 terms remain, add more */
if (nq - j < 4) {
k = 4 - (nq - j);
/*
* compute q[nq+1] to q[nq+k]
*
* For some reason, writing out the nx loop explicitly
* for each of the three possible values (as above) seems
* to run a little slower, so we'll leave this code as is.
*/
for (j = 1; j < nx; j++)
q[i] = t;
eqnqm1 -= 24;
}
nq += k;
goto recompute;
}
/* set pr and nq so that pr[0,...,nq] is the part of r remaining */
pr = &r[j];
/* compute pio2 * pr[0,...,nq]; note that nq >= 4 here */
for (i = 4; i <= nq; i++) {
}
/* sum q in increasing order to obtain the first term of y */
t = q[nq];
for (i = nq - 1; i >= 0; i--)
t += q[i];
y[0] = t;
if (prec) {
/* subtract and sum again in decreasing order
to obtain the second term */
t = q[0] - t;
for (i = 1; i <= nq; i++)
t += q[i];
y[1] = t;
}
return (n & 7);
}