code.consxy.com Git - clnl/blob - src/main/strictmath/sin.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_sin.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 ; /*
  25 ; FUNCTION
  26 ;         <<sin>>, <<sinf>>, <<cos>>, <<cosf>>---sine or cosine
  27 ; INDEX
  28 ; sin
  29 ; INDEX
  30 ; sinf
  31 ; INDEX
  32 ; cos
  33 ; INDEX
  34 ; cosf
  35 ; ANSI_SYNOPSIS
  36 ;         #include <math.h>
  37 ;         double sin(double <[x]>);
  38 ;         float  sinf(float <[x]>);
  39 ;         double cos(double <[x]>);
  40 ;         float cosf(float <[x]>);
  41 ;
  42 ; TRAD_SYNOPSIS
  43 ;         #include <math.h>
  44 ;         double sin(<[x]>)
  45 ;         double <[x]>;
  46 ;         float  sinf(<[x]>)
  47 ;         float <[x]>;
  48 ;
  49 ;         double cos(<[x]>)
  50 ;         double <[x]>;
  51 ;         float cosf(<[x]>)
  52 ;         float <[x]>;
  53 ;
  54 ; DESCRIPTION
  55 ;         <<sin>> and <<cos>> compute (respectively) the sine and cosine
  56 ;         of the argument <[x]>.  Angles are specified in radians.
  57 ;
  58 ;         <<sinf>> and <<cosf>> are identical, save that they take and
  59 ;         return <<float>> values.
  60 ;
  61 ;
  62 ; RETURNS
  63 ;         The sine or cosine of <[x]> is returned.
  64 ;
  65 ; PORTABILITY
  66 ;         <<sin>> and <<cos>> are ANSI C.
  67 ;         <<sinf>> and <<cosf>> are extensions.
  68 ;
  69 ; QUICKREF
  70 ;         sin ansi pure
  71 ;         sinf - pure
  72 ; */
  73 ;
  74 ; /* sin(x)
  75 ;  * Return sine function of x.
  76 ;  *
  77 ;  * kernel function:
  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
  81 ;  *
  82 ;  * Method.
  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.
  86 ;  *      We have
  87 ;  *
  88 ;  *          n        sin(x)      cos(x)        tan(x)
  89 ;  *     ----------------------------------------------------------
  90 ;  *          0               S           C                 T
  91 ;  *          1               C          -S                -1/T
  92 ;  *          2              -S          -C                 T
  93 ;  *          3              -C           S                -1/T
  94 ;  *     ----------------------------------------------------------
  95 ;  *
  96 ;  * Special cases:
  97 ;  *      Let trig be any of sin, cos, or tan.
  98 ;  *      trig(+-INF)  is NaN, with signals;
  99 ;  *      trig(NaN)    is that NaN;
 100 ;  *
 101 ;  * Accuracy:
 102 ;  *      TRIG(x) returns trig(x) nearly rounded
 103 ;  */
 104
 105 (defun sin (x)
 106  "SIN X => RESULT
 107
 108 ARGUMENTS AND VALUES:
 109
 110   X: A double representing the angle in radians bounded from [0, 2pi]
 111   RESULT: A double representing the sin of the angle
 112
 113 DESCRIPTION:
 114
 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)))
 120   (cond ;; |x| ~< pi/4
 121    ((<= ix #x3fe921fb) (float64-kernel-sin x 0d0 0) )
 122
 123    ;; sin(Inf or NaN) is NaN
 124    ((>= ix #x7ff00000) (- x x) )
 125
 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)))))))))