--- /dev/null
+(in-package #:strictmath)
+; This file is taken from part of Evita Common Lisp.
+;
+; 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_sin.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.
+; * ====================================================
+; */
+;
+; /*
+; FUNCTION
+; <<sin>>, <<sinf>>, <<cos>>, <<cosf>>---sine or cosine
+; INDEX
+; sin
+; INDEX
+; sinf
+; INDEX
+; cos
+; INDEX
+; cosf
+; ANSI_SYNOPSIS
+; #include <math.h>
+; double sin(double <[x]>);
+; float sinf(float <[x]>);
+; double cos(double <[x]>);
+; float cosf(float <[x]>);
+;
+; TRAD_SYNOPSIS
+; #include <math.h>
+; double sin(<[x]>)
+; double <[x]>;
+; float sinf(<[x]>)
+; float <[x]>;
+;
+; double cos(<[x]>)
+; double <[x]>;
+; float cosf(<[x]>)
+; float <[x]>;
+;
+; DESCRIPTION
+; <<sin>> and <<cos>> compute (respectively) the sine and cosine
+; of the argument <[x]>. Angles are specified in radians.
+;
+; <<sinf>> and <<cosf>> are identical, save that they take and
+; return <<float>> values.
+;
+;
+; RETURNS
+; The sine or cosine of <[x]> is returned.
+;
+; PORTABILITY
+; <<sin>> and <<cos>> are ANSI C.
+; <<sinf>> and <<cosf>> are extensions.
+;
+; QUICKREF
+; sin ansi pure
+; sinf - pure
+; */
+;
+; /* sin(x)
+; * Return sine function of x.
+; *
+; * kernel function:
+; * __kernel_sin ... sine function on [-pi/4,pi/4]
+; * __kernel_cos ... cose 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 sin (x)
+ "SIN X => RESULT
+
+ARGUMENTS AND VALUES:
+
+ X: A double representing the angle in radians bounded from [0, 2pi]
+ RESULT: A double representing the sin of the angle
+
+DESCRIPTION:
+
+ SIN returns the sin of the angle X."
+ (declare (values double-float))
+ (declare (type double-float x))
+ (let* ((hx (decode-float64 x))
+ (ix (logand hx #x7fffffff)))
+ (cond ;; |x| ~< pi/4
+ ((<= ix #x3fe921fb) (float64-kernel-sin x 0d0 0) )
+
+ ;; sin(Inf or NaN) is NaN
+ ((>= ix #x7ff00000) (- x x) )
+
+ ;; argument reduction needed
+ (t (multiple-value-bind (n y0 y1) (float64-rem-pio2 x)
+ (ecase (logand n 3) (0 (float64-kernel-sin y0 y1 1))
+ (1 (float64-kernel-cos y0 y1))
+ (2 (- (float64-kernel-sin y0 y1 1)))
+ (3 (- (float64-kernel-cos y0 y1)))))))))
projects / clnl / blobdiff