Mercurial > ecos

comparison packages/language/c/libm/current/src/double/internal/k_cos.c @ 0:3111d98ba7b3 ecos-v1_1-release

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial commit of eCos version 1.1
author jlarmour
date Tue, 11 May 1999 11:16:07 +0000
parents
children 443894e2e912
comparison
equal deleted inserted replaced
-1:000000000000 0:3111d98ba7b3
1 //===========================================================================
2 //
3 // k_cos.c
4 //
5 // Part of the standard mathematical function library
6 //
7 //===========================================================================
8 //####COPYRIGHTBEGIN####
9 //
10 // -------------------------------------------
11 // The contents of this file are subject to the Cygnus eCos Public License
12 // Version 1.0 (the "License"); you may not use this file except in
13 // compliance with the License. You may obtain a copy of the License at
14 // http://sourceware.cygnus.com/ecos
15 //
16 // Software distributed under the License is distributed on an "AS IS"
17 // basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See the
18 // License for the specific language governing rights and limitations under
19 // the License.
20 //
21 // The Original Code is eCos - Embedded Cygnus Operating System, released
22 // September 30, 1998.
23 //
24 // The Initial Developer of the Original Code is Cygnus. Portions created
25 // by Cygnus are Copyright (C) 1998 Cygnus Solutions. All Rights Reserved.
26 // -------------------------------------------
27 //
28 //####COPYRIGHTEND####
29 //===========================================================================
30 //#####DESCRIPTIONBEGIN####
31 //
32 // Author(s): jlarmour@cygnus.co.uk
33 // Contributors: jlarmour@cygnus.co.uk
34 // Date: 1998-02-13
35 // Purpose:
36 // Description:
37 // Usage:
38 //
39 //####DESCRIPTIONEND####
40 //
41 //===========================================================================
42
43 // CONFIGURATION
44
45 #include <pkgconf/libm.h> // Configuration header
46
47 // Include the Math library?
48 #ifdef CYGPKG_LIBM
49
50 // Derived from code with the following copyright
51
52
53 /* @(#)k_cos.c 1.3 95/01/18 */
54 /*
55 * ====================================================
56 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
57 *
58 * Developed at SunSoft, a Sun Microsystems, Inc. business.
59 * Permission to use, copy, modify, and distribute this
60 * software is freely granted, provided that this notice
61 * is preserved.
62 * ====================================================
63 */
64
65 /*
66 * __kernel_cos( x, y )
67 * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
68 * Input x is assumed to be bounded by ~pi/4 in magnitude.
69 * Input y is the tail of x.
70 *
71 * Algorithm
72 * 1. Since cos(-x) = cos(x), we need only to consider positive x.
73 * 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
74 * 3. cos(x) is approximated by a polynomial of degree 14 on
75 * [0,pi/4]
76 * 4 14
77 * cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
78 * where the remez error is
79 *
80 * | 2 4 6 8 10 12 14 | -58
81 * |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2
82 * | |
83 *
84 * 4 6 8 10 12 14
85 * 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then
86 * cos(x) = 1 - x*x/2 + r
87 * since cos(x+y) ~ cos(x) - sin(x)*y
88 * ~ cos(x) - x*y,
89 * a correction term is necessary in cos(x) and hence
90 * cos(x+y) = 1 - (x*x/2 - (r - x*y))
91 * For better accuracy when x > 0.3, let qx = |x|/4 with
92 * the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125.
93 * Then
94 * cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)).
95 * Note that 1-qx and (x*x/2-qx) is EXACT here, and the
96 * magnitude of the latter is at least a quarter of x*x/2,
97 * thus, reducing the rounding error in the subtraction.
98 */
99
100 #include "mathincl/fdlibm.h"
101
102 static const double
103 one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
104 C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
105 C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
106 C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
107 C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
108 C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
109 C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
110
111 double __kernel_cos(double x, double y)
112 {
113 double a,hz,z,r,qx;
114 int ix;
115 ix = CYG_LIBM_HI(x)&0x7fffffff; /* ix = |x|'s high word*/
116 if(ix<0x3e400000) { /* if x < 2**27 */
117 if(((int)x)==0) return one; /* generate inexact */
118 }
119 z = x*x;
120 r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6)))));
121 if(ix < 0x3FD33333) /* if |x| < 0.3 */
122 return one - (0.5*z - (z*r - x*y));
123 else {
124 if(ix > 0x3fe90000) { /* x > 0.78125 */
125 qx = 0.28125;
126 } else {
127 CYG_LIBM_HI(qx) = ix-0x00200000; /* x/4 */
128 CYG_LIBM_LO(qx) = 0;
129 }
130 hz = 0.5*z-qx;
131 a = one-qx;
132 return a - (hz - (z*r-x*y));
133 }
134 }
135
136 #endif // ifdef CYGPKG_LIBM
137
138 // EOF k_cos.c