src/libm/k_cos.c
 author Ryan C. Gordon Mon, 21 May 2018 12:00:21 -0400 changeset 11994 8e094f91ddab parent 11683 48bcba563d9c permissions -rw-r--r--
 slouken@2756 ` 1` ```/* ``` slouken@2756 ` 2` ``` * ==================================================== ``` slouken@2756 ` 3` ``` * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. ``` slouken@2756 ` 4` ``` * ``` slouken@2756 ` 5` ``` * Developed at SunPro, a Sun Microsystems, Inc. business. ``` slouken@2756 ` 6` ``` * Permission to use, copy, modify, and distribute this ``` slouken@2756 ` 7` ``` * software is freely granted, provided that this notice ``` slouken@2756 ` 8` ``` * is preserved. ``` slouken@2756 ` 9` ``` * ==================================================== ``` slouken@2756 ` 10` ``` */ ``` slouken@2756 ` 11` slouken@2756 ` 12` ```/* ``` slouken@2756 ` 13` ``` * __kernel_cos( x, y ) ``` slouken@2756 ` 14` ``` * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164 ``` slouken@2756 ` 15` ``` * Input x is assumed to be bounded by ~pi/4 in magnitude. ``` slouken@2756 ` 16` ``` * Input y is the tail of x. ``` slouken@2756 ` 17` ``` * ``` slouken@2756 ` 18` ``` * Algorithm ``` slouken@2756 ` 19` ``` * 1. Since cos(-x) = cos(x), we need only to consider positive x. ``` slouken@2756 ` 20` ``` * 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0. ``` slouken@2756 ` 21` ``` * 3. cos(x) is approximated by a polynomial of degree 14 on ``` slouken@2756 ` 22` ``` * [0,pi/4] ``` slouken@2756 ` 23` ``` * 4 14 ``` slouken@2756 ` 24` ``` * cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x ``` slouken@2756 ` 25` ``` * where the remez error is ``` slouken@2756 ` 26` ``` * ``` slouken@2756 ` 27` ``` * | 2 4 6 8 10 12 14 | -58 ``` slouken@2756 ` 28` ``` * |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2 ``` slouken@2756 ` 29` ``` * | | ``` slouken@2756 ` 30` ``` * ``` slouken@2756 ` 31` ``` * 4 6 8 10 12 14 ``` slouken@2756 ` 32` ``` * 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then ``` slouken@2756 ` 33` ``` * cos(x) = 1 - x*x/2 + r ``` slouken@2756 ` 34` ``` * since cos(x+y) ~ cos(x) - sin(x)*y ``` slouken@2756 ` 35` ``` * ~ cos(x) - x*y, ``` slouken@2756 ` 36` ``` * a correction term is necessary in cos(x) and hence ``` slouken@2756 ` 37` ``` * cos(x+y) = 1 - (x*x/2 - (r - x*y)) ``` slouken@2756 ` 38` ``` * For better accuracy when x > 0.3, let qx = |x|/4 with ``` slouken@2756 ` 39` ``` * the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125. ``` slouken@2756 ` 40` ``` * Then ``` slouken@2756 ` 41` ``` * cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)). ``` slouken@2756 ` 42` ``` * Note that 1-qx and (x*x/2-qx) is EXACT here, and the ``` slouken@2756 ` 43` ``` * magnitude of the latter is at least a quarter of x*x/2, ``` slouken@2756 ` 44` ``` * thus, reducing the rounding error in the subtraction. ``` slouken@2756 ` 45` ``` */ ``` slouken@2756 ` 46` slouken@6044 ` 47` ```#include "math_libm.h" ``` slouken@2756 ` 48` ```#include "math_private.h" ``` slouken@2756 ` 49` slouken@2756 ` 50` ```static const double ``` slouken@11683 ` 51` ```one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ ``` slouken@11683 ` 52` ```C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */ ``` slouken@11683 ` 53` ```C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */ ``` slouken@11683 ` 54` ```C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */ ``` slouken@11683 ` 55` ```C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */ ``` slouken@11683 ` 56` ```C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */ ``` slouken@11683 ` 57` ```C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */ ``` slouken@2756 ` 58` slouken@11683 ` 59` ```double attribute_hidden __kernel_cos(double x, double y) ``` slouken@2756 ` 60` ```{ ``` slouken@11683 ` 61` ``` double a,hz,z,r,qx; ``` slouken@11683 ` 62` ``` int32_t ix; ``` slouken@11683 ` 63` ``` GET_HIGH_WORD(ix,x); ``` slouken@11683 ` 64` ``` ix &= 0x7fffffff; /* ix = |x|'s high word*/ ``` slouken@11683 ` 65` ``` if(ix<0x3e400000) { /* if x < 2**27 */ ``` slouken@11683 ` 66` ``` if(((int)x)==0) return one; /* generate inexact */ ``` slouken@11683 ` 67` ``` } ``` slouken@11683 ` 68` ``` z = x*x; ``` slouken@11683 ` 69` ``` r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6))))); ``` slouken@11683 ` 70` ``` if(ix < 0x3FD33333) /* if |x| < 0.3 */ ``` slouken@11683 ` 71` ``` return one - (0.5*z - (z*r - x*y)); ``` slouken@11683 ` 72` ``` else { ``` slouken@11683 ` 73` ``` if(ix > 0x3fe90000) { /* x > 0.78125 */ ``` slouken@11683 ` 74` ``` qx = 0.28125; ``` slouken@11683 ` 75` ``` } else { ``` slouken@11683 ` 76` ``` INSERT_WORDS(qx,ix-0x00200000,0); /* x/4 */ ``` slouken@11683 ` 77` ``` } ``` slouken@11683 ` 78` ``` hz = 0.5*z-qx; ``` slouken@11683 ` 79` ``` a = one-qx; ``` slouken@11683 ` 80` ``` return a - (hz - (z*r-x*y)); ``` slouken@11683 ` 81` ``` } ``` slouken@2756 ` 82` ```} ```