code.consxy.com Git - clnl/blob - src/main/strictmath/cos.lisp

   1 (in-package #:strictmath)
   2 ; This file is taken from part of Evita Common Lisp.
   3 ;
   4 ; Copyright (C) 1996-2007 by Project Vogue.
   5 ; Written by Yoshifumi "VOGUE" INOUE. (yosi@msn.com)
   6 ;
   7 ; Before that, it was based off of fdlibm
   8 ;
   9 ;  See fdlibm (http://www.netlib.org/fdlibm/)
  10 ;  See http://sources.redhat.com/newlib/
  11 ;
  12 ; /* @(#)s_cos.c 5.1 93/09/24 */
  13 ; /*
  14 ;  * ====================================================
  15 ;  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  16 ;  *
  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
  20 ;  * is preserved.
  21 ;  * ====================================================
  22 ;  */
  23 ;
  24 ; /* cos(x)
  25 ;  * Return cosine function of x.
  26 ;  *
  27 ;  * kernel function:
  28 ;  *      __kernel_sin                ... sine function on [-pi/4,pi/4]
  29 ;  *      __kernel_cos                ... cosine function on [-pi/4,pi/4]
  30 ;  *      __ieee754_rem_pio2        ... argument reduction routine
  31 ;  *
  32 ;  * Method.
  33 ;  *      Let S,C and T denote the sin, cos and tan respectively on
  34 ;  *      [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
  35 ;  *      in [-pi/4 , +pi/4], and let n = k mod 4.
  36 ;  *      We have
  37 ;  *
  38 ;  *          n        sin(x)      cos(x)        tan(x)
  39 ;  *     ----------------------------------------------------------
  40 ;  *          0               S           C                 T
  41 ;  *          1               C          -S                -1/T
  42 ;  *          2              -S          -C                 T
  43 ;  *          3              -C           S                -1/T
  44 ;  *     ----------------------------------------------------------
  45 ;  *
  46 ;  * Special cases:
  47 ;  *      Let trig be any of sin, cos, or tan.
  48 ;  *      trig(+-INF)  is NaN, with signals;
  49 ;  *      trig(NaN)    is that NaN;
  50 ;  *
  51 ;  * Accuracy:
  52 ;  *      TRIG(x) returns trig(x) nearly rounded
  53 ;  */
  54
  55 (defun cos (x)
  56  "COS X => RESULT
  57
  58 ARGUMENTS AND VALUES:
  59
  60   X: A double representing the angle in radians bounded from [0, 2pi]
  61   RESULT: A double representing the cos of the angle
  62
  63 DESCRIPTION:
  64
  65   COS returns the cos of the angle X."
  66  (declare (values double-float))
  67  (declare (type double-float x))
  68  (let ((hx (logand (decode-float64 x) #x7fffffff)))
  69   (cond
  70    ;; |x| ~< pi/4
  71    ((<= hx #x3fe921fb) (float64-kernel-cos x 0d0))
  72
  73    ;; cos(Inf or NaN) is NaN
  74    ((>= hx #x7ff00000) (- x x))
  75
  76    ;; argument reduction needed
  77    (t
  78     (multiple-value-bind (n y0 y1) (float64-rem-pio2 x)
  79      (ecase (logand n 3)
  80       (0 (float64-kernel-cos y0 y1))
  81       (1 (- (float64-kernel-sin y0 y1 1)))
  82       (2 (- (float64-kernel-cos y0 y1)))
  83       (3 (float64-kernel-sin y0 y1 1))))))))