/* * 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 * or http://www.opensolaris.org/os/licensing. * 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, round1 = 6755399441055744.0, /* 3 * 2^51 */ round24 = 113336795588871485128704.0, /* 3 * 2^75 */ twon24 = 5.960464477539062500E-8; int __vlibm_rem_pio2m(double *x, double *y, int e, int nx, int prec) { union { double d; int i[2]; } s; double z, t, p, q[20], r[21], *pr; int nq, ip, n, i, j, k, eq0, eqnqm1; /* determine ip and eq0; note that -48 <= eq0 <= 2 */ ip = (e - 3) / 24; if (ip < 0) ip = 0; eq0 = e - 24 * (ip + 1); /* compute q[0,...,5] = x * ipio2 and initialize nq and eqnqm1 */ if (nx == 3) { q[0] = x[0] * ipio2[ip+2] + x[1] * ipio2[ip+1] + x[2] * ipio2[ip]; q[1] = x[0] * ipio2[ip+3] + x[1] * ipio2[ip+2] + x[2] * ipio2[ip+1]; q[2] = x[0] * ipio2[ip+4] + x[1] * ipio2[ip+3] + x[2] * ipio2[ip+2]; q[3] = x[0] * ipio2[ip+5] + x[1] * ipio2[ip+4] + x[2] * ipio2[ip+3]; q[4] = x[0] * ipio2[ip+6] + x[1] * ipio2[ip+5] + x[2] * ipio2[ip+4]; q[5] = x[0] * ipio2[ip+7] + x[1] * ipio2[ip+6] + x[2] * ipio2[ip+5]; } else if (nx == 2) { q[0] = x[0] * ipio2[ip+2] + x[1] * ipio2[ip+1]; q[1] = x[0] * ipio2[ip+3] + x[1] * ipio2[ip+2]; q[2] = x[0] * ipio2[ip+4] + x[1] * ipio2[ip+3]; q[3] = x[0] * ipio2[ip+5] + x[1] * ipio2[ip+4]; q[4] = x[0] * ipio2[ip+6] + x[1] * ipio2[ip+5]; q[5] = x[0] * ipio2[ip+7] + x[1] * ipio2[ip+6]; } else { q[0] = x[0] * ipio2[ip+2]; q[1] = x[0] * ipio2[ip+3]; q[2] = x[0] * ipio2[ip+4]; q[3] = x[0] * ipio2[ip+5]; q[4] = x[0] * ipio2[ip+6]; q[5] = x[0] * ipio2[ip+7]; } nq = 5; eqnqm1 = eq0 - 96; recompute: /* propagate carries and incorporate powers of two */ s.i[HIWORD] = (0x3ff + eqnqm1) << 20; s.i[LOWORD] = 0; p = s.d; z = q[nq] * twon24; for (j = nq-1; j >= 1; j--) { z += q[j]; t = (z + round24) - round24; /* must be rounded to nearest */ r[j+1] = (z - t) * p; z = t * twon24; p *= two24; } z += q[0]; t = (z + round24) - round24; /* must be rounded to nearest */ 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]; t = (z + round1) - round1; /* must be rounded to nearest */ 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]; t = (z + round1) - round1; /* must be rounded to nearest */ r[1] -= t; z = r[0] + t; /* cut off high part of z so conversion to int doesn't overflow */ t = (z + round24) - round24; 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]; t = (z + round1) - round1; /* must be rounded to nearest */ 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 (i = nq + 1; i <= nq + k; i++) { t = x[0] * ipio2[ip+2+i]; for (j = 1; j < nx; j++) t += x[j] * ipio2[ip+2+i-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]; nq = nq - j; /* compute pio2 * pr[0,...,nq]; note that nq >= 4 here */ q[0] = pio2[0] * pr[0]; q[1] = pio2[0] * pr[1] + pio2[1] * pr[0]; q[2] = pio2[0] * pr[2] + pio2[1] * pr[1] + pio2[2] * pr[0]; q[3] = pio2[0] * pr[3] + pio2[1] * pr[2] + pio2[2] * pr[1] + pio2[3] * pr[0]; for (i = 4; i <= nq; i++) { q[i] = pio2[0] * pr[i] + pio2[1] * pr[i-1] + pio2[2] * pr[i-2] + pio2[3] * pr[i-3] + pio2[4] * pr[i-4]; } /* 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); }
