1 (in-package #:strictmath)
2 ; This file is taken from part of Evita Common Lisp.
4 ; Copyright (C) 1996-2007 by Project Vogue.
5 ; Written by Yoshifumi "VOGUE" INOUE. (yosi@msn.com)
7 ; Before that, it was based off of fdlibm
9 ; See fdlibm (http://www.netlib.org/fdlibm/)
10 ; See http://sources.redhat.com/newlib/
12 ; /* @(#)s_sin.c 5.1 93/09/24 */
14 ; * ====================================================
15 ; * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
17 ; * Developed at SunPro, a Sun Microsystems, Inc. business.
18 ; * Permission to use, copy, modify, and distribute this
19 ; * software is freely granted, provided that this notice
21 ; * ====================================================
26 ; <<sin>>, <<sinf>>, <<cos>>, <<cosf>>---sine or cosine
37 ; double sin(double <[x]>);
38 ; float sinf(float <[x]>);
39 ; double cos(double <[x]>);
40 ; float cosf(float <[x]>);
55 ; <<sin>> and <<cos>> compute (respectively) the sine and cosine
56 ; of the argument <[x]>. Angles are specified in radians.
58 ; <<sinf>> and <<cosf>> are identical, save that they take and
59 ; return <<float>> values.
63 ; The sine or cosine of <[x]> is returned.
66 ; <<sin>> and <<cos>> are ANSI C.
67 ; <<sinf>> and <<cosf>> are extensions.
75 ; * Return sine function of x.
78 ; * __kernel_sin ... sine function on [-pi/4,pi/4]
79 ; * __kernel_cos ... cose function on [-pi/4,pi/4]
80 ; * __ieee754_rem_pio2 ... argument reduction routine
83 ; * Let S,C and T denote the sin, cos and tan respectively on
84 ; * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
85 ; * in [-pi/4 , +pi/4], and let n = k mod 4.
88 ; * n sin(x) cos(x) tan(x)
89 ; * ----------------------------------------------------------
94 ; * ----------------------------------------------------------
97 ; * Let trig be any of sin, cos, or tan.
98 ; * trig(+-INF) is NaN, with signals;
99 ; * trig(NaN) is that NaN;
102 ; * TRIG(x) returns trig(x) nearly rounded
108 ARGUMENTS AND VALUES:
110 X: A double representing the angle in radians bounded from [0, 2pi]
111 RESULT: A double representing the sin of the angle
115 SIN returns the sin of the angle X."
116 (declare (values double-float))
117 (declare (type double-float x))
118 (let* ((hx (decode-float64 x))
119 (ix (logand hx #x7fffffff)))
121 ((<= ix #x3fe921fb) (float64-kernel-sin x 0d0 0) )
123 ;; sin(Inf or NaN) is NaN
124 ((>= ix #x7ff00000) (- x x) )
126 ;; argument reduction needed
127 (t (multiple-value-bind (n y0 y1) (float64-rem-pio2 x)
128 (ecase (logand n 3) (0 (float64-kernel-sin y0 y1 1))
129 (1 (float64-kernel-cos y0 y1))
130 (2 (- (float64-kernel-sin y0 y1 1)))
131 (3 (- (float64-kernel-cos y0 y1)))))))))