sincospil.c revision 25c28e83beb90e7c80452a7c818c5e6f73a07dc8
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER START
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * The contents of this file are subject to the terms of the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Common Development and Distribution License (the "License").
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You may not use this file except in compliance with the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * See the License for the specific language governing permissions
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * and limitations under the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER END
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Use is subject to license terms.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * void sincospil(long double x, long double *s, long double *c)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * *s = sinl(pi*x); *c = cosl(pi*x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Algorithm, 10/17/2002, K.C. Ng
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * ------------------------------
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Let y = |4x|, z = floor(y), and n = (int)(z mod 8.0) (displayed in binary).
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1. If y == z, then x is a multiple of pi/4. Return the following values:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * ---------------------------------------------------
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * n x mod 2 sin(x*pi) cos(x*pi) tan(x*pi)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * ---------------------------------------------------
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 000 0.00 +0 ___ +1 ___ +0
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 001 0.25 +\/0.5 +\/0.5 +1
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 010 0.50 +1 ___ +0 ___ +inf
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 011 0.75 +\/0.5 -\/0.5 -1
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 100 1.00 -0 ___ -1 ___ +0
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 101 1.25 -\/0.5 -\/0.5 +1
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 110 1.50 -1 ___ -0 ___ +inf
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 111 1.75 -\/0.5 +\/0.5 -1
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * ---------------------------------------------------
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 2. Otherwise,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * ---------------------------------------------------
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * n t sin(x*pi) cos(x*pi) tan(x*pi)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * ---------------------------------------------------
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 000 (y-z)/4 sinpi(t) cospi(t) tanpi(t)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 001 (z+1-y)/4 cospi(t) sinpi(t) 1/tanpi(t)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 010 (y-z)/4 cospi(t) -sinpi(t) -1/tanpi(t)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 011 (z+1-y)/4 sinpi(t) -cospi(t) -tanpi(t)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 100 (y-z)/4 -sinpi(t) -cospi(t) tanpi(t)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 101 (z+1-y)/4 -cospi(t) -sinpi(t) 1/tanpi(t)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 110 (y-z)/4 -cospi(t) sinpi(t) -1/tanpi(t)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 111 (z+1-y)/4 -sinpi(t) cospi(t) -tanpi(t)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * ---------------------------------------------------
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * NOTE. This program compute sinpi/cospi(t<0.25) by __k_sin/cos(pi*t, 0.0).
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * This will return a result with error slightly more than one ulp (but less
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * than 2 ulp). If one wants accurate result, one may break up pi*t in
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * high (tpi_h) and low (tpi_l) parts and call __k_sin/cos(tip_h, tip_lo)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define I(q, m) ((int *) &(q))[m]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define U(q, m) ((unsigned *) &(q))[m]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define LDBL_MOST_SIGNIF_I(ld) ((I(ld, 2) << 16) | (0xffff & (I(ld, 1) >> 15)))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const long double twoPRECM2 = 9.223372036854775808000000000000000e+18L;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define LDBL_LEAST_SIGNIF_U(ld) U(ld, sizeof(long double) / sizeof(int) - 1)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const long double twoPRECM2 = 5.192296858534827628530496329220096e+33L;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const long double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtispi = 3.141592653589793238462643383279502884197e+0000L,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtissqrth = 0.707106781186547524400844362104849039284835937688474,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtissincospil(long double x, long double *s, long double *c) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis long double y, z, t;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (k >= PRECM2) { /* |x| >= 2**(Prec-2) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (k >= 16384) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *s = *c = x - x;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else if (k == PRECM1) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else { /* k = Prec - 2 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* y = |4x|, z = floor(y), and n = (int)(z mod 8.0) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis n = LDBL_LEAST_SIGNIF_U(z) & 7; /* 3 LSb of z */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else if (t > y) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = (y - t) * quater;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else { /* k = Prec-3 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis n = LDBL_LEAST_SIGNIF_U(y) & 7; /* 3 LSb of z */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (k) { /* x = N/4 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((n & 1) != 0)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((n & 2) == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((n & 4) != 0)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((n & 1) != 0)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((n & 4) != 0)