1*a2142627SBen Jolitz /*
2*a2142627SBen Jolitz * Copyright (c) 1987 Regents of the University of California.
3*a2142627SBen Jolitz * All rights reserved.
4*a2142627SBen Jolitz *
5*a2142627SBen Jolitz * Redistribution and use in source and binary forms, with or without
6*a2142627SBen Jolitz * modification, are permitted provided that the following conditions
7*a2142627SBen Jolitz * are met:
8*a2142627SBen Jolitz * 1. Redistributions of source code must retain the above copyright
9*a2142627SBen Jolitz * notice, this list of conditions and the following disclaimer.
10*a2142627SBen Jolitz * 2. Redistributions in binary form must reproduce the above copyright
11*a2142627SBen Jolitz * notice, this list of conditions and the following disclaimer in the
12*a2142627SBen Jolitz * documentation and/or other materials provided with the distribution.
13*a2142627SBen Jolitz * 3. All advertising materials mentioning features or use of this software
14*a2142627SBen Jolitz * must display the following acknowledgement:
15*a2142627SBen Jolitz * This product includes software developed by the University of
16*a2142627SBen Jolitz * California, Berkeley and its contributors.
17*a2142627SBen Jolitz * 4. Neither the name of the University nor the names of its contributors
18*a2142627SBen Jolitz * may be used to endorse or promote products derived from this software
19*a2142627SBen Jolitz * without specific prior written permission.
20*a2142627SBen Jolitz *
21*a2142627SBen Jolitz * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
22*a2142627SBen Jolitz * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23*a2142627SBen Jolitz * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24*a2142627SBen Jolitz * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
25*a2142627SBen Jolitz * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26*a2142627SBen Jolitz * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
27*a2142627SBen Jolitz * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
28*a2142627SBen Jolitz * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29*a2142627SBen Jolitz * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
30*a2142627SBen Jolitz * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
31*a2142627SBen Jolitz * SUCH DAMAGE.
32*a2142627SBen Jolitz */
33*a2142627SBen Jolitz
34*a2142627SBen Jolitz #ifndef lint
35*a2142627SBen Jolitz static char sccsid[] = "@(#)sincos.c 5.5 (Berkeley) 10/9/90";
36*a2142627SBen Jolitz #endif /* not lint */
37*a2142627SBen Jolitz
38*a2142627SBen Jolitz #include "trig.h"
39*a2142627SBen Jolitz double
sin(x)40*a2142627SBen Jolitz sin(x)
41*a2142627SBen Jolitz double x;
42*a2142627SBen Jolitz {
43*a2142627SBen Jolitz double a,c,z;
44*a2142627SBen Jolitz
45*a2142627SBen Jolitz if(!finite(x)) /* sin(NaN) and sin(INF) must be NaN */
46*a2142627SBen Jolitz return x-x;
47*a2142627SBen Jolitz x=drem(x,PI2); /* reduce x into [-PI,PI] */
48*a2142627SBen Jolitz a=copysign(x,one);
49*a2142627SBen Jolitz if (a >= PIo4) {
50*a2142627SBen Jolitz if(a >= PI3o4) /* ... in [3PI/4,PI] */
51*a2142627SBen Jolitz x = copysign((a = PI-a),x);
52*a2142627SBen Jolitz else { /* ... in [PI/4,3PI/4] */
53*a2142627SBen Jolitz a = PIo2-a; /* rtn. sign(x)*C(PI/2-|x|) */
54*a2142627SBen Jolitz z = a*a;
55*a2142627SBen Jolitz c = cos__C(z);
56*a2142627SBen Jolitz z *= half;
57*a2142627SBen Jolitz a = (z >= thresh ? half-((z-half)-c) : one-(z-c));
58*a2142627SBen Jolitz return copysign(a,x);
59*a2142627SBen Jolitz }
60*a2142627SBen Jolitz }
61*a2142627SBen Jolitz
62*a2142627SBen Jolitz if (a < small) { /* rtn. S(x) */
63*a2142627SBen Jolitz big+a;
64*a2142627SBen Jolitz return x;
65*a2142627SBen Jolitz }
66*a2142627SBen Jolitz return x+x*sin__S(x*x);
67*a2142627SBen Jolitz }
68*a2142627SBen Jolitz
69*a2142627SBen Jolitz double
cos(x)70*a2142627SBen Jolitz cos(x)
71*a2142627SBen Jolitz double x;
72*a2142627SBen Jolitz {
73*a2142627SBen Jolitz double a,c,z,s = 1.0;
74*a2142627SBen Jolitz
75*a2142627SBen Jolitz if(!finite(x)) /* cos(NaN) and cos(INF) must be NaN */
76*a2142627SBen Jolitz return x-x;
77*a2142627SBen Jolitz x=drem(x,PI2); /* reduce x into [-PI,PI] */
78*a2142627SBen Jolitz a=copysign(x,one);
79*a2142627SBen Jolitz if (a >= PIo4) {
80*a2142627SBen Jolitz if (a >= PI3o4) { /* ... in [3PI/4,PI] */
81*a2142627SBen Jolitz a = PI-a;
82*a2142627SBen Jolitz s = negone;
83*a2142627SBen Jolitz }
84*a2142627SBen Jolitz else { /* ... in [PI/4,3PI/4] */
85*a2142627SBen Jolitz a = PIo2-a;
86*a2142627SBen Jolitz return a+a*sin__S(a*a); /* rtn. S(PI/2-|x|) */
87*a2142627SBen Jolitz }
88*a2142627SBen Jolitz }
89*a2142627SBen Jolitz if (a < small) {
90*a2142627SBen Jolitz big+a;
91*a2142627SBen Jolitz return s; /* rtn. s*C(a) */
92*a2142627SBen Jolitz }
93*a2142627SBen Jolitz z = a*a;
94*a2142627SBen Jolitz c = cos__C(z);
95*a2142627SBen Jolitz z *= half;
96*a2142627SBen Jolitz a = (z >= thresh ? half-((z-half)-c) : one-(z-c));
97*a2142627SBen Jolitz return copysign(a,s);
98*a2142627SBen Jolitz }
99