Skip to content
Snippets Groups Projects
Commit f9c37dc7 authored by wlott's avatar wlott
Browse files

Wrote the rest of the trigonometric functions and fixed the complex case

of several others.
parent d7a9e68e
No related branches found
No related tags found
No related merge requests found
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
;;; Scott Fahlman (FAHLMAN@CMUC). ;;; Scott Fahlman (FAHLMAN@CMUC).
;;; ********************************************************************** ;;; **********************************************************************
;;; ;;;
;;; $Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/irrat.lisp,v 1.2 1990/07/21 15:33:08 wlott Exp $ ;;; $Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/irrat.lisp,v 1.3 1990/07/31 17:22:14 wlott Exp $
;;; ;;;
;;; This file contains all the irrational functions. Actually, most of the ;;; This file contains all the irrational functions. Actually, most of the
;;; work is done by calling out to C... ;;; work is done by calling out to C...
...@@ -114,7 +114,7 @@ ...@@ -114,7 +114,7 @@
;;; inverted if negative. ;;; inverted if negative.
(defun expt (base power) (defun expt (base power)
"Returns x raised to the nth power." "Returns BASE raised to the POWER."
(if (zerop power) (if (zerop power)
;; This is wrong if power isn't an integer. ;; This is wrong if power isn't an integer.
(typecase (realpart base) (typecase (realpart base)
...@@ -149,21 +149,32 @@ ...@@ -149,21 +149,32 @@
(exp (* power (log base))))))) (exp (* power (log base)))))))
(defun log (number &optional (base nil base-p)) (defun log (number &optional (base nil base-p))
"Return the logarithm of NUMBER in the base BASE, which defaults to e."
(if base-p (if base-p
(/ (log number) (log base)) (/ (log number) (log base))
(number-dispatch ((number number)) (number-dispatch ((number number))
(handle-reals %log number) (((foreach fixnum bignum ratio single-float))
((complex) (if (minusp number)
(complex (log (abs number)) (complex (log (- number)) (coerce pi 'single-float))
(phase number)))))) (coerce (%log (coerce number 'double-float)) 'single-float)))
((double-float)
(if (minusp number)
(complex (log (- number)) (coerce pi 'double-float))
(%log number)))
((complex) (complex (log (abs number)) (phase number))))))
(defun sqrt (number) (defun sqrt (number)
"Return the square root of NUMBER." "Return the square root of NUMBER."
(number-dispatch ((number number)) (number-dispatch ((number number))
(handle-reals %sqrt number) (((foreach fixnum bignum ratio single-float))
((complex) (if (minusp number)
(* (exp (/ (complex 0 (phase number)) 2)) (exp (/ (log number) 2))
(sqrt (abs number)))))) (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
((double-float)
(if (minusp number)
(exp (/ (log number) 2))
(%sqrt number)))
((complex) (exp (/ (log number) 2)))))
;;; ISQRT: Integer square root - isqrt(n)**2 <= n ;;; ISQRT: Integer square root - isqrt(n)**2 <= n
;;; Upper and lower bounds on the result are estimated using integer-length. ;;; Upper and lower bounds on the result are estimated using integer-length.
...@@ -207,27 +218,24 @@ ...@@ -207,27 +218,24 @@
(defun phase (number) (defun phase (number)
"Returns the angle part of the polar representation of a complex number. "Returns the angle part of the polar representation of a complex number.
For complex numbers, this is (atan (imagpart number) (realpart number)).
For non-complex positive numbers, this is 0. For non-complex negative For non-complex positive numbers, this is 0. For non-complex negative
numbers this is PI." numbers this is PI."
(if (zerop number) (etypecase number
(if (typep number 'double-float) 0.0d0 0.0) ((or rational single-float)
(atan (imagpart number) (realpart number)))) (if (minusp number)
(coerce pi 'single-float)
#| 0.0f0))
(if (complexp number) (double-float
(let ((ipart (imagpart number)) (if (minusp number)
(rpart (realpart number))) (coerce pi 'double-float)
(if (zerop rpart) 0.0d0))
(if (minusp ipart) (complex
(if (long-float-p ipart) (- %long-pi/2) (- %short-pi/2)) (atan (imagpart number) (realpart number)))))
(if (long-float-p ipart) %long-pi/2 %short-pi/2))
(atan (/ (imagpart number) (realpart number)))))
(if (minusp number)
(if (long-float-p number) pi %short-pi)
(if (long-float-p number) 0.0l0 0.0))))
|#
(defun sin (number) (defun sin (number)
"Return the sine of NUMBER."
(number-dispatch ((number number)) (number-dispatch ((number number))
(handle-reals %sin number) (handle-reals %sin number)
((complex) ((complex)
...@@ -236,6 +244,7 @@ ...@@ -236,6 +244,7 @@
(complex (* (sin x) (cosh y)) (* (cos x) (sinh y))))))) (complex (* (sin x) (cosh y)) (* (cos x) (sinh y)))))))
(defun cos (number) (defun cos (number)
"Return the cosine of NUMBER."
(number-dispatch ((number number)) (number-dispatch ((number number))
(handle-reals %cos number) (handle-reals %cos number)
((complex) ((complex)
...@@ -244,6 +253,7 @@ ...@@ -244,6 +253,7 @@
(complex (* (cos x) (cosh y)) (- (* (sin x) (sinh y)))))))) (complex (* (cos x) (cosh y)) (- (* (sin x) (sinh y))))))))
(defun tan (number) (defun tan (number)
"Return the tangent of NUMBER."
(number-dispatch ((number number)) (number-dispatch ((number number))
(handle-reals %tan number) (handle-reals %tan number)
((complex) ((complex)
...@@ -253,44 +263,72 @@ ...@@ -253,44 +263,72 @@
(/ num denom)))))) (/ num denom))))))
(defun cis (theta) (defun cis (theta)
"Return cos(Theta) + i sin(Theta), aka exp(i Theta)." "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
(if (complexp theta) (if (complexp theta)
(error "Argument to CIS is complex: ~S" theta) (error "Argument to CIS is complex: ~S" theta)
(complex (cos theta) (sin theta)))) (complex (cos theta) (sin theta))))
#+nil (proclaim '(inline mult-by-i))
(defun in-asin-domain (z) (defun mult-by-i (number)
(or (< (- (/ pi 2.0d0)) (realpart z) (/ pi 2.0d0)) (complex (imagpart number)
(and (= (realpart z) (- (/ pi 2.0d0))) (- (realpart number))))
(>= (imagpart z) 0))
(and (= (realpart z) (/ pi 2.0d0)) (defun complex-asin (number)
(<= (imagpart z) 0)))) (- (mult-by-i (log (+ (mult-by-i number) (sqrt (- 1 (* number number))))))))
(defun asin (number) (defun asin (number)
"Return the arc sine of NUMBER."
(number-dispatch ((number number)) (number-dispatch ((number number))
(handle-reals %asin number) ((rational)
(if (or (> number 1) (< number -1))
(complex-asin number)
(coerce (%asin (coerce number 'double-float)) 'single-float)))
(((foreach single-float double-float))
(if (or (> number (coerce 1 '(dispatch-type number)))
(< number (coerce -1 '(dispatch-type number))))
(complex-asin number)
(coerce (%asin (coerce number 'double-float))
'(dispatch-type number))))
((complex) ((complex)
(error "Can't hack complex ASIN yet: ~S" number)))) (complex-asin number))))
#+nil (defun complex-acos (number)
(defun in-acos-domain (z) (- (mult-by-i (log (+ number (mult-by-i (sqrt (- (* number number)))))))))
(or (< 0 (realpart z) pi)
(and (= 0 (realpart z))
(>= (imagpart z) 0))
(and (= (realpart z) pi)
(<= (imagpart z) 0))))
(defun acos (number) (defun acos (number)
"Return the arc cosine of NUMBER."
(number-dispatch ((number number)) (number-dispatch ((number number))
(handle-reals %acos number) ((rational)
(if (or (> number 1) (< number -1))
(complex-acos number)
(coerce (%acos (coerce number 'double-float)) 'single-float)))
(((foreach single-float double-float))
(if (or (> number (coerce 1 '(dispatch-type number)))
(< number (coerce -1 '(dispatch-type number))))
(complex-acos number)
(coerce (%acos (coerce number 'double-float))
'(dispatch-type number))))
((complex) ((complex)
(error "Can't hack complex ACOS yet: ~S" number)))) (complex-acos number))))
(defun atan (y &optional (x nil xp)) (defun atan (y &optional (x nil xp))
"Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
(if xp (if xp
(if (and (zerop x) (zerop y)) (if (and (zerop x) (zerop y))
(error "Both args to ATAN can't be zero.") (multiple-value-bind
(mag exp sign-x)
(integer-decode-float (float x))
(declare (ignore mag exp))
(if (plusp sign-x)
y
(multiple-value-bind
(mag exp sign-y)
(integer-decode-float (float y))
(declare (ignore mag exp))
(if (minusp sign-y)
(- pi)
pi))))
(number-dispatch ((y real) (x real)) (number-dispatch ((y real) (x real))
(((foreach fixnum bignum ratio single-float) (((foreach fixnum bignum ratio single-float)
(foreach fixnum bignum ratio single-float)) (foreach fixnum bignum ratio single-float))
...@@ -305,4 +343,36 @@ ...@@ -305,4 +343,36 @@
(number-dispatch ((y number)) (number-dispatch ((y number))
(handle-reals %atan y) (handle-reals %atan y)
((complex) ((complex)
(error "Can't handle complex ATAN yet: ~S" y))))) (let ((im (imagpart y))
(re (realpart y)))
(/ (- (log (complex (- 1 im) re))
(log (complex (+ 1 im) (- re))))
(complex 0 2)))))))
(defun sinh (number)
"Return the hyperbolic sine of NUMBER."
(/ (- (exp number) (exp (- number))) 2))
(defun cosh (number)
"Return the hyperbolic cosine of NUMBER."
(/ (+ (exp number) (exp (- number))) 2))
(defun tanh (number)
"Return the hyperbolic tangent of NUMBER."
(/ (- (exp number) (exp (- number)))
(+ (exp number) (exp (- number)))))
(defun asinh (number)
"Return the hyperbolic arc sine of NUMBER."
(log (+ number (sqrt (1+ (* number number))))))
(defun acosh (number)
"Return the hyperbolic arc cosine of NUMBER."
(log (+ number (* (1+ number) (sqrt (/ (1- number) (1+ number)))))))
(defun atanh (number)
"Return the hyperbolic arc tangent of NUMBER."
(log (* (1+ number) (sqrt (/ (- 1 (* number number)))))))
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment