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_cos.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 ; * ====================================================
29 ; * Return cosine function of x.
32 ; * __kernel_sin ... sine function on [-pi/4,pi/4]
33 ; * __kernel_cos ... cosine function on [-pi/4,pi/4]
34 ; * __ieee754_rem_pio2 ... argument reduction routine
37 ; * Let S,C and T denote the sin, cos and tan respectively on
38 ; * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
39 ; * in [-pi/4 , +pi/4], and let n = k mod 4.
42 ; * n sin(x) cos(x) tan(x)
43 ; * ----------------------------------------------------------
48 ; * ----------------------------------------------------------
51 ; * Let trig be any of sin, cos, or tan.
52 ; * trig(+-INF) is NaN, with signals;
53 ; * trig(NaN) is that NaN;
56 ; * TRIG(x) returns trig(x) nearly rounded
64 X: A double representing the angle in radians bounded from [0, 2pi]
65 RESULT: A double representing the cos of the angle
69 COS returns the cos of the angle X."
70 (declare (values double-float))
71 (declare (type double-float x))
72 (let ((hx (logand (decode-float64 x) #x7fffffff)))
75 ((<= hx #x3fe921fb) (float64-kernel-cos x 0d0))
77 ;; cos(Inf or NaN) is NaN
78 ((>= hx #x7ff00000) (- x x))
80 ;; argument reduction needed
82 (multiple-value-bind (n y0 y1) (float64-rem-pio2 x)
84 (0 (float64-kernel-cos y0 y1))
85 (1 (- (float64-kernel-sin y0 y1 1)))
86 (2 (- (float64-kernel-cos y0 y1)))
87 (3 (float64-kernel-sin y0 y1 1))))))))