179 


180 
;; As above for x86. It also seems to be needed to handle 
;; As above for x86. It also seems to be needed to handle 
181 
;; constantfolding in the compiler. 
;; constantfolding in the compiler. 
182 
#+sparc 
#+(or sparc (and x86 sse2)) 
183 
(progn 
(progn 
184 
(defun %sqrt (x) 
(defun %sqrt (x) 
185 
(declare (doublefloat x) 
(declare (doublefloat x) 
293 
;;;; Power functions. 
;;;; Power functions. 
294 


295 
(defun exp (number) 
(defun exp (number) 
296 
_N"Return e raised to the power NUMBER." 
"Return e raised to the power NUMBER." 
297 
(numberdispatch ((number number)) 
(numberdispatch ((number number)) 
298 
(handlereals %exp number) 
(handlereals %exp number) 
299 
((complex) 
((complex) 
308 
((base :initarg :base :reader intexpbase) 
((base :initarg :base :reader intexpbase) 
309 
(power :initarg :power :reader intexppower)) 
(power :initarg :power :reader intexppower)) 
310 
(:report (lambda (condition stream) 
(:report (lambda (condition stream) 
311 
(format stream _"The absolute value of ~S exceeds limit ~S." 
(format stream (intl:gettext "The absolute value of ~S exceeds limit ~S.") 
312 
(intexppower condition) 
(intexppower condition) 
313 
*intexpmaximumexponent*)))) 
*intexpmaximumexponent*)))) 
314 


333 
:power power) 
:power power) 
334 
(continue () 
(continue () 
335 
:report (lambda (stream) 
:report (lambda (stream) 
336 
(writestring _"Continue with calculation" stream))) 
(writestring (intl:gettext "Continue with calculation") stream))) 
337 
(newlimit () 
(newlimit () 
338 
:report (lambda (stream) 
:report (lambda (stream) 
339 
(writestring _"Continue with calculation, update limit" stream)) 
(writestring (intl:gettext "Continue with calculation, update limit") stream)) 
340 
(setq *intexpmaximumexponent* (abs power))))) 
(setq *intexpmaximumexponent* (abs power))))) 
341 
(cond ((minusp power) 
(cond ((minusp power) 
342 
(/ (intexp base ( power)))) 
(/ (intexp base ( power)))) 
360 
;;; the complexreal and realcomplex cases from the general complex case. 
;;; the complexreal and realcomplex cases from the general complex case. 
361 
;;; 
;;; 
362 
(defun expt (base power) 
(defun expt (base power) 
363 
_N"Returns BASE raised to the POWER." 
"Returns BASE raised to the POWER." 
364 
(if (zerop power) 
(if (zerop power) 
365 
;; CLHS says that if the power is 0, the result is 1, subject to 
;; CLHS says that if the power is 0, the result is 1, subject to 
366 
;; numeric contagion. But what happens if base is infinity or 
;; numeric contagion. But what happens if base is infinity or 
676 
(+ n frac)))))) 
(+ n frac)))))) 
677 


678 
(defun log (number &optional (base nil basep)) 
(defun log (number &optional (base nil basep)) 
679 
_N"Return the logarithm of NUMBER in the base BASE, which defaults to e." 
"Return the logarithm of NUMBER in the base BASE, which defaults to e." 
680 
(if basep 
(if basep 
681 
(cond ((zerop base) 
(cond ((zerop base) 
682 
;; ANSI spec 
;; ANSI spec 
804 
(complexlog number))))) 
(complexlog number))))) 
805 


806 
(defun sqrt (number) 
(defun sqrt (number) 
807 
_N"Return the square root of NUMBER." 
"Return the square root of NUMBER." 
808 
(numberdispatch ((number number)) 
(numberdispatch ((number number)) 
809 
(((foreach fixnum bignum ratio)) 
(((foreach fixnum bignum ratio)) 
810 
(if (minusp number) 
(if (minusp number) 
829 
;;;; Trigonometic and Related Functions 
;;;; Trigonometic and Related Functions 
830 


831 
(defun abs (number) 
(defun abs (number) 
832 
_N"Returns the absolute value of the number." 
"Returns the absolute value of the number." 
833 
(numberdispatch ((number number)) 
(numberdispatch ((number number)) 
834 
(((foreach singlefloat doublefloat fixnum rational 
(((foreach singlefloat doublefloat fixnum rational 
835 
#+doubledouble doubledoublefloat)) 
#+doubledouble doubledoublefloat)) 
853 
(scalefloat (sqrt abs^2) scale)))))))) 
(scalefloat (sqrt abs^2) scale)))))))) 
854 


855 
(defun phase (number) 
(defun phase (number) 
856 
_N"Returns the angle part of the polar representation of a complex number. 
"Returns the angle part of the polar representation of a complex number. 
857 
For complex numbers, this is (atan (imagpart number) (realpart number)). 
For complex numbers, this is (atan (imagpart number) (realpart number)). 
858 
For noncomplex positive numbers, this is 0. For noncomplex negative 
For noncomplex positive numbers, this is 0. For noncomplex negative 
859 
numbers this is PI." 
numbers this is PI." 
880 


881 


882 
(defun sin (number) 
(defun sin (number) 
883 
_N"Return the sine of NUMBER." 
"Return the sine of NUMBER." 
884 
(numberdispatch ((number number)) 
(numberdispatch ((number number)) 
885 
(handlereals %sin number) 
(handlereals %sin number) 
886 
((complex) 
((complex) 
890 
(* (cos x) (sinh y))))))) 
(* (cos x) (sinh y))))))) 
891 


892 
(defun cos (number) 
(defun cos (number) 
893 
_N"Return the cosine of NUMBER." 
"Return the cosine of NUMBER." 
894 
(numberdispatch ((number number)) 
(numberdispatch ((number number)) 
895 
(handlereals %cos number) 
(handlereals %cos number) 
896 
((complex) 
((complex) 
900 
( (* (sin x) (sinh y)))))))) 
( (* (sin x) (sinh y)))))))) 
901 


902 
(defun tan (number) 
(defun tan (number) 
903 
_N"Return the tangent of NUMBER." 
"Return the tangent of NUMBER." 
904 
(numberdispatch ((number number)) 
(numberdispatch ((number number)) 
905 
(handlereals %tan number) 
(handlereals %tan number) 
906 
((complex) 
((complex) 
907 
(complextan number)))) 
(complextan number)))) 
908 


909 
(defun cis (theta) 
(defun cis (theta) 
910 
_N"Return cos(Theta) + i sin(Theta), AKA exp(i Theta)." 
"Return cos(Theta) + i sin(Theta), AKA exp(i Theta)." 
911 
(if (complexp theta) 
(if (complexp theta) 
912 
(error _"Argument to CIS is complex: ~S" theta) 
(error (intl:gettext "Argument to CIS is complex: ~S") theta) 
913 
(complex (cos theta) (sin theta)))) 
(complex (cos theta) (sin theta)))) 
914 


915 
(defun asin (number) 
(defun asin (number) 
916 
_N"Return the arc sine of NUMBER." 
"Return the arc sine of NUMBER." 
917 
(numberdispatch ((number number)) 
(numberdispatch ((number number)) 
918 
((rational) 
((rational) 
919 
(if (or (> number 1) (< number 1)) 
(if (or (> number 1) (< number 1)) 
937 
(complexasin number)))) 
(complexasin number)))) 
938 


939 
(defun acos (number) 
(defun acos (number) 
940 
_N"Return the arc cosine of NUMBER." 
"Return the arc cosine of NUMBER." 
941 
(numberdispatch ((number number)) 
(numberdispatch ((number number)) 
942 
((rational) 
((rational) 
943 
(if (or (> number 1) (< number 1)) 
(if (or (> number 1) (< number 1)) 
962 


963 


964 
(defun atan (y &optional (x nil xp)) 
(defun atan (y &optional (x nil xp)) 
965 
_N"Return the arc tangent of Y if X is omitted or Y/X if X is supplied." 
"Return the arc tangent of Y if X is omitted or Y/X if X is supplied." 
966 
(if xp 
(if xp 
967 
(flet ((atan2 (y x) 
(flet ((atan2 (y x) 
968 
(declare (type doublefloat y x) 
(declare (type doublefloat y x) 
1000 
(complexatan y))))) 
(complexatan y))))) 
1001 


1002 
(defun sinh (number) 
(defun sinh (number) 
1003 
_N"Return the hyperbolic sine of NUMBER." 
"Return the hyperbolic sine of NUMBER." 
1004 
(numberdispatch ((number number)) 
(numberdispatch ((number number)) 
1005 
(handlereals %sinh number) 
(handlereals %sinh number) 
1006 
((complex) 
((complex) 
1010 
(* (cosh x) (sin y))))))) 
(* (cosh x) (sin y))))))) 
1011 


1012 
(defun cosh (number) 
(defun cosh (number) 
1013 
_N"Return the hyperbolic cosine of NUMBER." 
"Return the hyperbolic cosine of NUMBER." 
1014 
(numberdispatch ((number number)) 
(numberdispatch ((number number)) 
1015 
(handlereals %cosh number) 
(handlereals %cosh number) 
1016 
((complex) 
((complex) 
1020 
(* (sinh x) (sin y))))))) 
(* (sinh x) (sin y))))))) 
1021 


1022 
(defun tanh (number) 
(defun tanh (number) 
1023 
_N"Return the hyperbolic tangent of NUMBER." 
"Return the hyperbolic tangent of NUMBER." 
1024 
(numberdispatch ((number number)) 
(numberdispatch ((number number)) 
1025 
(handlereals %tanh number) 
(handlereals %tanh number) 
1026 
((complex) 
((complex) 
1027 
(complextanh number)))) 
(complextanh number)))) 
1028 


1029 
(defun asinh (number) 
(defun asinh (number) 
1030 
_N"Return the hyperbolic arc sine of NUMBER." 
"Return the hyperbolic arc sine of NUMBER." 
1031 
(numberdispatch ((number number)) 
(numberdispatch ((number number)) 
1032 
(handlereals %asinh number) 
(handlereals %asinh number) 
1033 
((complex) 
((complex) 
1034 
(complexasinh number)))) 
(complexasinh number)))) 
1035 


1036 
(defun acosh (number) 
(defun acosh (number) 
1037 
_N"Return the hyperbolic arc cosine of NUMBER." 
"Return the hyperbolic arc cosine of NUMBER." 
1038 
(numberdispatch ((number number)) 
(numberdispatch ((number number)) 
1039 
((rational) 
((rational) 
1040 
;; acosh is complex if number < 1 
;; acosh is complex if number < 1 
1055 
(complexacosh number)))) 
(complexacosh number)))) 
1056 


1057 
(defun atanh (number) 
(defun atanh (number) 
1058 
_N"Return the hyperbolic arc tangent of NUMBER." 
"Return the hyperbolic arc tangent of NUMBER." 
1059 
(numberdispatch ((number number)) 
(numberdispatch ((number number)) 
1060 
((rational) 
((rational) 
1061 
;; atanh is complex if number > 1 
;; atanh is complex if number > 1 
1152 


1153 
(declaim (inline scalb)) 
(declaim (inline scalb)) 
1154 
(defun scalb (x n) 
(defun scalb (x n) 
1155 
_N"Compute 2^N * X without compute 2^N first (use properties of the 
"Compute 2^N * X without compute 2^N first (use properties of the 
1156 
underlying floatingpoint format" 
underlying floatingpoint format" 
1157 
(declare (type float x) 
(declare (type float x) 
1158 
(type doublefloatexponent n)) 
(type doublefloatexponent n)) 
1160 


1161 
(declaim (inline logbfinite)) 
(declaim (inline logbfinite)) 
1162 
(defun logbfinite (x) 
(defun logbfinite (x) 
1163 
_N"Same as logb but X is not infinity and nonzero and not a NaN, so 
"Same as logb but X is not infinity and nonzero and not a NaN, so 
1164 
that we can always return an integer" 
that we can always return an integer" 
1165 
(declare (type float x)) 
(declare (type float x)) 
1166 
(multiplevaluebind (signif expon sign) 
(multiplevaluebind (signif expon sign) 
1171 
(1 expon))) 
(1 expon))) 
1172 


1173 
(defun logb (x) 
(defun logb (x) 
1174 
_N"Compute an integer N such that 1 <= 2^(N) * x < 2. 
"Compute an integer N such that 1 <= 2^(N) * x < 2. 
1175 
For the special cases, the following values are used: 
For the special cases, the following values are used: 
1176 


1177 
x logb 
x logb 
1199 


1200 
(declaim (inline coercetocomplextype)) 
(declaim (inline coercetocomplextype)) 
1201 
(defun coercetocomplextype (x y z) 
(defun coercetocomplextype (x y z) 
1202 
_N"Create complex number with real part X and imaginary part Y such that 
"Create complex number with real part X and imaginary part Y such that 
1203 
it has the same type as Z. If Z has type (complex rational), the X 
it has the same type as Z. If Z has type (complex rational), the X 
1204 
and Y are coerced to singlefloat." 
and Y are coerced to singlefloat." 
1205 
(declare (doublefloat x y) 
(declare (doublefloat x y) 
1250 
(values rho 0))))))) 
(values rho 0))))))) 
1251 


1252 
(defun complexsqrt (z) 
(defun complexsqrt (z) 
1253 
_N"Principle square root of Z 
"Principle square root of Z 
1254 


1255 
Z may be any number, but the result is always a complex." 
Z may be any number, but the result is always a complex." 
1256 
(declare (number z)) 
(declare (number z)) 
1295 
(coercetocomplextype eta nu z))))) 
(coercetocomplextype eta nu z))))) 
1296 


1297 
(defun complexlogscaled (z j) 
(defun complexlogscaled (z j) 
1298 
_N"Compute log(2^j*z). 
"Compute log(2^j*z). 
1299 


1300 
This is for use with J /= 0 only when z is huge." 
This is for use with J /= 0 only when z is huge." 
1301 
(declare (number z) 
(declare (number z) 
1330 
z))))) 
z))))) 
1331 


1332 
(defun complexlog (z) 
(defun complexlog (z) 
1333 
_N"Log of Z = log Z + i * arg Z 
"Log of Z = log Z + i * arg Z 
1334 


1335 
Z may be any number, but the result is always a complex." 
Z may be any number, but the result is always a complex." 
1336 
(declare (number z)) 
(declare (number z)) 
1345 
;; never 0 since we have positive and negative zeroes. 
;; never 0 since we have positive and negative zeroes. 
1346 


1347 
(defun complexatanh (z) 
(defun complexatanh (z) 
1348 
_N"Compute atanh z = (log(1+z)  log(1z))/2" 
"Compute atanh z = (log(1+z)  log(1z))/2" 
1349 
(declare (number z)) 
(declare (number z)) 
1350 
#+doubledouble 
#+doubledouble 
1351 
(when (typep z '(or doubledoublefloat (complex doubledoublefloat))) 
(when (typep z '(or doubledoublefloat (complex doubledoublefloat))) 
1409 
z))))) 
z))))) 
1410 


1411 
(defun complextanh (z) 
(defun complextanh (z) 
1412 
_N"Compute tanh z = sinh z / cosh z" 
"Compute tanh z = sinh z / cosh z" 
1413 
(declare (number z)) 
(declare (number z)) 
1414 
#+doubledouble 
#+doubledouble 
1415 
(when (typep z '(or doubledoublefloat (complex doubledoublefloat))) 
(when (typep z '(or doubledoublefloat (complex doubledoublefloat))) 
1494 
(complex (+ (realpart z) 1) (imagpart z))) 
(complex (+ (realpart z) 1) (imagpart z))) 
1495 


1496 
(defun complexacos (z) 
(defun complexacos (z) 
1497 
_N"Compute acos z = pi/2  asin z 
"Compute acos z = pi/2  asin z 
1498 


1499 
Z may be any number, but the result is always a complex." 
Z may be any number, but the result is always a complex." 
1500 
(declare (number z)) 
(declare (number z)) 
1513 
sqrt1z)))))))) 
sqrt1z)))))))) 
1514 


1515 
(defun complexacosh (z) 
(defun complexacosh (z) 
1516 
_N"Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z1)/2)) 
"Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z1)/2)) 
1517 


1518 
Z may be any number, but the result is always a complex." 
Z may be any number, but the result is always a complex." 
1519 
(declare (number z)) 
(declare (number z)) 
1527 


1528 


1529 
(defun complexasin (z) 
(defun complexasin (z) 
1530 
_N"Compute asin z = asinh(i*z)/i 
"Compute asin z = asinh(i*z)/i 
1531 


1532 
Z may be any number, but the result is always a complex." 
Z may be any number, but the result is always a complex." 
1533 
(declare (number z)) 
(declare (number z)) 
1546 
sqrt1+z)))))))) 
sqrt1+z)))))))) 
1547 


1548 
(defun complexasinh (z) 
(defun complexasinh (z) 
1549 
_N"Compute asinh z = log(z + sqrt(1 + z*z)) 
"Compute asinh z = log(z + sqrt(1 + z*z)) 
1550 


1551 
Z may be any number, but the result is always a complex." 
Z may be any number, but the result is always a complex." 
1552 
(declare (number z)) 
(declare (number z)) 
1560 
( (realpart result))))) 
( (realpart result))))) 
1561 


1562 
(defun complexatan (z) 
(defun complexatan (z) 
1563 
_N"Compute atan z = atanh (i*z) / i 
"Compute atan z = atanh (i*z) / i 
1564 


1565 
Z may be any number, but the result is always a complex." 
Z may be any number, but the result is always a complex." 
1566 
(declare (number z)) 
(declare (number z)) 
1574 
( (realpart result))))) 
( (realpart result))))) 
1575 


1576 
(defun complextan (z) 
(defun complextan (z) 
1577 
_N"Compute tan z = i * tanh(i * z) 
"Compute tan z = i * tanh(i * z) 
1578 


1579 
Z may be any number, but the result is always a complex." 
Z may be any number, but the result is always a complex." 
1580 
(declare (number z)) 
(declare (number z)) 