; Copyright 2022 Frank Duncan (frank@consxy.com) under AGPL3.  See distributed LICENSE.txt.
(in-package #:strictmath)
; This file is taken from part of Evita Common Lisp.
;
; It has been updated to match the rest of the project's documentation and style
; standards.  But otherwise, the following copyright supersedes the above AGPL copyright.
;
; Copyright (C) 1996-2007 by Project Vogue.
; Written by Yoshifumi "VOGUE" INOUE. (yosi@msn.com)
;
; Before that, it was based off of fdlibm
;
;  See fdlibm (http://www.netlib.org/fdlibm/)
;  See http://sources.redhat.com/newlib/
; 
; /* @(#)s_cos.c 5.1 93/09/24 */
; /*
;  * ====================================================
;  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
;  *
;  * Developed at SunPro, a Sun Microsystems, Inc. business.
;  * Permission to use, copy, modify, and distribute this
;  * software is freely granted, provided that this notice 
;  * is preserved.
;  * ====================================================
;  */
; 
; /* cos(x)
;  * Return cosine function of x.
;  *
;  * kernel function:
;  *      __kernel_sin                ... sine function on [-pi/4,pi/4]
;  *      __kernel_cos                ... cosine function on [-pi/4,pi/4]
;  *      __ieee754_rem_pio2        ... argument reduction routine
;  *
;  * Method.
;  *      Let S,C and T denote the sin, cos and tan respectively on 
;  *      [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 
;  *      in [-pi/4 , +pi/4], and let n = k mod 4.
;  *      We have
;  *
;  *          n        sin(x)      cos(x)        tan(x)
;  *     ----------------------------------------------------------
;  *          0               S           C                 T
;  *          1               C          -S                -1/T
;  *          2              -S          -C                 T
;  *          3              -C           S                -1/T
;  *     ----------------------------------------------------------
;  *
;  * Special cases:
;  *      Let trig be any of sin, cos, or tan.
;  *      trig(+-INF)  is NaN, with signals;
;  *      trig(NaN)    is that NaN;
;  *
;  * Accuracy:
;  *      TRIG(x) returns trig(x) nearly rounded 
;  */

(defun cos (x)
 "COS X => RESULT

ARGUMENTS AND VALUES:

  X: A double representing the angle in radians bounded from [0, 2pi]
  RESULT: A double representing the cos of the angle

DESCRIPTION:

  COS returns the cos of the angle X."
 (declare (values double-float))
 (declare (type double-float x))
 (let ((hx (logand (decode-float64 x) #x7fffffff)))
  (cond
   ;; |x| ~< pi/4
   ((<= hx #x3fe921fb) (float64-kernel-cos x 0d0))

   ;; cos(Inf or NaN) is NaN
   ((>= hx #x7ff00000) (- x x))

   ;; argument reduction needed
   (t
    (multiple-value-bind (n y0 y1) (float64-rem-pio2 x)
     (ecase (logand n 3)
      (0 (float64-kernel-cos y0 y1))
      (1 (- (float64-kernel-sin y0 y1 1)))
      (2 (- (float64-kernel-cos y0 y1)))
      (3 (float64-kernel-sin y0 y1 1))))))))