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

   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.
   4 ;
   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.
   7 ;
   8 ; Copyright (C) 1996-2007 by Project Vogue.
   9 ; Written by Yoshifumi "VOGUE" INOUE. (yosi@msn.com)
  10 ;
  11 ; Before that, it was based off of fdlibm
  12 ;
  13 ;  See fdlibm (http://www.netlib.org/fdlibm/)
  14 ;  See http://sources.redhat.com/newlib/
  15 ;
  16 ; /* @(#)s_sin.c 5.1 93/09/24 */
  17 ; /*
  18 ;  * ====================================================
  19 ;  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  20 ;  *
  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
  24 ;  * is preserved.
  25 ;  * ====================================================
  26 ;  */
  27 ;
  28 ; /*
  29 ; FUNCTION
  30 ;         <<sin>>, <<sinf>>, <<cos>>, <<cosf>>---sine or cosine
  31 ; INDEX
  32 ; sin
  33 ; INDEX
  34 ; sinf
  35 ; INDEX
  36 ; cos
  37 ; INDEX
  38 ; cosf
  39 ; ANSI_SYNOPSIS
  40 ;         #include <math.h>
  41 ;         double sin(double <[x]>);
  42 ;         float  sinf(float <[x]>);
  43 ;         double cos(double <[x]>);
  44 ;         float cosf(float <[x]>);
  45 ;
  46 ; TRAD_SYNOPSIS
  47 ;         #include <math.h>
  48 ;         double sin(<[x]>)
  49 ;         double <[x]>;
  50 ;         float  sinf(<[x]>)
  51 ;         float <[x]>;
  52 ;
  53 ;         double cos(<[x]>)
  54 ;         double <[x]>;
  55 ;         float cosf(<[x]>)
  56 ;         float <[x]>;
  57 ;
  58 ; DESCRIPTION
  59 ;         <<sin>> and <<cos>> compute (respectively) the sine and cosine
  60 ;         of the argument <[x]>.  Angles are specified in radians.
  61 ;
  62 ;         <<sinf>> and <<cosf>> are identical, save that they take and
  63 ;         return <<float>> values.
  64 ;
  65 ;
  66 ; RETURNS
  67 ;         The sine or cosine of <[x]> is returned.
  68 ;
  69 ; PORTABILITY
  70 ;         <<sin>> and <<cos>> are ANSI C.
  71 ;         <<sinf>> and <<cosf>> are extensions.
  72 ;
  73 ; QUICKREF
  74 ;         sin ansi pure
  75 ;         sinf - pure
  76 ; */
  77 ;
  78 ; /* sin(x)
  79 ;  * Return sine function of x.
  80 ;  *
  81 ;  * kernel function:
  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
  85 ;  *
  86 ;  * Method.
  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.
  90 ;  *      We have
  91 ;  *
  92 ;  *          n        sin(x)      cos(x)        tan(x)
  93 ;  *     ----------------------------------------------------------
  94 ;  *          0               S           C                 T
  95 ;  *          1               C          -S                -1/T
  96 ;  *          2              -S          -C                 T
  97 ;  *          3              -C           S                -1/T
  98 ;  *     ----------------------------------------------------------
  99 ;  *
 100 ;  * Special cases:
 101 ;  *      Let trig be any of sin, cos, or tan.
 102 ;  *      trig(+-INF)  is NaN, with signals;
 103 ;  *      trig(NaN)    is that NaN;
 104 ;  *
 105 ;  * Accuracy:
 106 ;  *      TRIG(x) returns trig(x) nearly rounded
 107 ;  */
 108
 109 (defun sin (x)
 110  "SIN X => RESULT
 111
 112 ARGUMENTS AND VALUES:
 113
 114   X: A double representing the angle in radians bounded from [0, 2pi]
 115   RESULT: A double representing the sin of the angle
 116
 117 DESCRIPTION:
 118
 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)))
 124   (cond ;; |x| ~< pi/4
 125    ((<= ix #x3fe921fb) (float64-kernel-sin x 0d0 0) )
 126
 127    ;; sin(Inf or NaN) is NaN
 128    ((>= ix #x7ff00000) (- x x) )
 129
 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)))))))))