1 
siliconforks 
2 
/* * Mode: C; tabwidth: 8; indenttabsmode: nil; cbasicoffset: 4 * 
2 


* 
3 


* ***** BEGIN LICENSE BLOCK ***** 
4 


* Version: MPL 1.1/GPL 2.0/LGPL 2.1 
5 


* 
6 


* The contents of this file are subject to the Mozilla Public License Version 
7 


* 1.1 (the "License"); you may not use this file except in compliance with 
8 


* the License. You may obtain a copy of the License at 
9 


* http://www.mozilla.org/MPL/ 
10 


* 
11 


* Software distributed under the License is distributed on an "AS IS" basis, 
12 


* WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License 
13 


* for the specific language governing rights and limitations under the 
14 


* License. 
15 


* 
16 


* The Original Code is Mozilla Communicator client code, released 
17 


* March 31, 1998. 
18 


* 
19 


* The Initial Developer of the Original Code is 
20 


* Sun Microsystems, Inc. 
21 


* Portions created by the Initial Developer are Copyright (C) 1998 
22 


* the Initial Developer. All Rights Reserved. 
23 


* 
24 


* Contributor(s): 
25 


* 
26 


* Alternatively, the contents of this file may be used under the terms of 
27 


* either of the GNU General Public License Version 2 or later (the "GPL"), 
28 


* or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"), 
29 


* in which case the provisions of the GPL or the LGPL are applicable instead 
30 


* of those above. If you wish to allow use of your version of this file only 
31 


* under the terms of either the GPL or the LGPL, and not to allow others to 
32 


* use your version of this file under the terms of the MPL, indicate your 
33 


* decision by deleting the provisions above and replace them with the notice 
34 


* and other provisions required by the GPL or the LGPL. If you do not delete 
35 


* the provisions above, a recipient may use your version of this file under 
36 


* the terms of any one of the MPL, the GPL or the LGPL. 
37 


* 
38 


* ***** END LICENSE BLOCK ***** */ 
39 



40 


/* @(#)k_cos.c 1.3 95/01/18 */ 
41 


/* 
42 


* ==================================================== 
43 


* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 
44 


* 
45 


* Developed at SunSoft, a Sun Microsystems, Inc. business. 
46 


* Permission to use, copy, modify, and distribute this 
47 


* software is freely granted, provided that this notice 
48 


* is preserved. 
49 


* ==================================================== 
50 


*/ 
51 



52 


/* 
53 


* __kernel_cos( x, y ) 
54 


* kernel cos function on [pi/4, pi/4], pi/4 ~ 0.785398164 
55 


* Input x is assumed to be bounded by ~pi/4 in magnitude. 
56 


* Input y is the tail of x. 
57 


* 
58 


* Algorithm 
59 


* 1. Since cos(x) = cos(x), we need only to consider positive x. 
60 


* 2. if x < 2^27 (hx<0x3e400000 0), return 1 with inexact if x!=0. 
61 


* 3. cos(x) is approximated by a polynomial of degree 14 on 
62 


* [0,pi/4] 
63 


* 4 14 
64 


* cos(x) ~ 1  x*x/2 + C1*x + ... + C6*x 
65 


* where the remez error is 
66 


* 
67 


*  2 4 6 8 10 12 14  58 
68 


* cos(x)(1.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x ) <= 2 
69 


*   
70 


* 
71 


* 4 6 8 10 12 14 
72 


* 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then 
73 


* cos(x) = 1  x*x/2 + r 
74 


* since cos(x+y) ~ cos(x)  sin(x)*y 
75 


* ~ cos(x)  x*y, 
76 


* a correction term is necessary in cos(x) and hence 
77 


* cos(x+y) = 1  (x*x/2  (r  x*y)) 
78 


* For better accuracy when x > 0.3, let qx = x/4 with 
79 


* the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125. 
80 


* Then 
81 


* cos(x+y) = (1qx)  ((x*x/2qx)  (rx*y)). 
82 


* Note that 1qx and (x*x/2qx) is EXACT here, and the 
83 


* magnitude of the latter is at least a quarter of x*x/2, 
84 


* thus, reducing the rounding error in the subtraction. 
85 


*/ 
86 



87 


#include "fdlibm.h" 
88 



89 


#ifdef __STDC__ 
90 


static const double 
91 


#else 
92 


static double 
93 


#endif 
94 


one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ 
95 


C1 = 4.16666666666666019037e02, /* 0x3FA55555, 0x5555554C */ 
96 


C2 = 1.38888888888741095749e03, /* 0xBF56C16C, 0x16C15177 */ 
97 


C3 = 2.48015872894767294178e05, /* 0x3EFA01A0, 0x19CB1590 */ 
98 


C4 = 2.75573143513906633035e07, /* 0xBE927E4F, 0x809C52AD */ 
99 


C5 = 2.08757232129817482790e09, /* 0x3E21EE9E, 0xBDB4B1C4 */ 
100 


C6 = 1.13596475577881948265e11; /* 0xBDA8FAE9, 0xBE8838D4 */ 
101 



102 


#ifdef __STDC__ 
103 


double __kernel_cos(double x, double y) 
104 


#else 
105 


double __kernel_cos(x, y) 
106 


double x,y; 
107 


#endif 
108 


{ 
109 


fd_twoints u; 
110 


double qx = 0; 
111 


double a,hz,z,r; 
112 


int ix; 
113 


u.d = x; 
114 


ix = __HI(u)&0x7fffffff; /* ix = x's high word*/ 
115 


if(ix<0x3e400000) { /* if x < 2**27 */ 
116 


if(((int)x)==0) return one; /* generate inexact */ 
117 


} 
118 


z = x*x; 
119 


r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6))))); 
120 


if(ix < 0x3FD33333) /* if x < 0.3 */ 
121 


return one  (0.5*z  (z*r  x*y)); 
122 


else { 
123 


if(ix > 0x3fe90000) { /* x > 0.78125 */ 
124 


qx = 0.28125; 
125 


} else { 
126 


u.d = qx; 
127 


__HI(u) = ix0x00200000; /* x/4 */ 
128 


__LO(u) = 0; 
129 


qx = u.d; 
130 


} 
131 


hz = 0.5*zqx; 
132 


a = oneqx; 
133 


return a  (hz  (z*rx*y)); 
134 


} 
135 


} 