1 ; Copyright 2022 Frank Duncan (frank@consxy.com) under AGPL3. See distributed LICENSE.txt.
2 (in-package #:strictmath)
3 ; This file is taken from part of Evita Common Lisp.
5 ; It has been updated to match the rest of the project's documentation and style
6 ; standards. But otherwise, the following copyright supersedes the above AGPL copyright.
8 ; Copyright (C) 1996-2007 by Project Vogue.
9 ; Written by Yoshifumi "VOGUE" INOUE. (yosi@msn.com)
11 ; Before that, it was based off of fdlibm
13 ; See fdlibm (http://www.netlib.org/fdlibm/)
14 ; See http://sources.redhat.com/newlib/
16 ; /* @(#)s_sin.c 5.1 93/09/24 */
18 ; * ====================================================
19 ; * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
21 ; * Developed at SunPro, a Sun Microsystems, Inc. business.
22 ; * Permission to use, copy, modify, and distribute this
23 ; * software is freely granted, provided that this notice
25 ; * ====================================================
30 ; <<sin>>, <<sinf>>, <<cos>>, <<cosf>>---sine or cosine
41 ; double sin(double <[x]>);
42 ; float sinf(float <[x]>);
43 ; double cos(double <[x]>);
44 ; float cosf(float <[x]>);
59 ; <<sin>> and <<cos>> compute (respectively) the sine and cosine
60 ; of the argument <[x]>. Angles are specified in radians.
62 ; <<sinf>> and <<cosf>> are identical, save that they take and
63 ; return <<float>> values.
67 ; The sine or cosine of <[x]> is returned.
70 ; <<sin>> and <<cos>> are ANSI C.
71 ; <<sinf>> and <<cosf>> are extensions.
79 ; * Return sine function of x.
82 ; * __kernel_sin ... sine function on [-pi/4,pi/4]
83 ; * __kernel_cos ... cose function on [-pi/4,pi/4]
84 ; * __ieee754_rem_pio2 ... argument reduction routine
87 ; * Let S,C and T denote the sin, cos and tan respectively on
88 ; * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
89 ; * in [-pi/4 , +pi/4], and let n = k mod 4.
92 ; * n sin(x) cos(x) tan(x)
93 ; * ----------------------------------------------------------
98 ; * ----------------------------------------------------------
101 ; * Let trig be any of sin, cos, or tan.
102 ; * trig(+-INF) is NaN, with signals;
103 ; * trig(NaN) is that NaN;
106 ; * TRIG(x) returns trig(x) nearly rounded
112 ARGUMENTS AND VALUES:
114 X: A double representing the angle in radians bounded from [0, 2pi]
115 RESULT: A double representing the sin of the angle
119 SIN returns the sin of the angle X."
120 (declare (values double-float))
121 (declare (type double-float x))
122 (let* ((hx (decode-float64 x))
123 (ix (logand hx #x7fffffff)))
125 ((<= ix #x3fe921fb) (float64-kernel-sin x 0d0 0) )
127 ;; sin(Inf or NaN) is NaN
128 ((>= ix #x7ff00000) (- x x) )
130 ;; argument reduction needed
131 (t (multiple-value-bind (n y0 y1) (float64-rem-pio2 x)
132 (ecase (logand n 3) (0 (float64-kernel-sin y0 y1 1))
133 (1 (float64-kernel-cos y0 y1))
134 (2 (- (float64-kernel-sin y0 y1 1)))
135 (3 (- (float64-kernel-cos y0 y1)))))))))