/[cmucl]/src/code/irrat-dd.lisp
ViewVC logotype

Contents of /src/code/irrat-dd.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.2.1 - (show annotations)
Thu Jun 29 01:28:02 2006 UTC (7 years, 9 months ago) by rtoy
Branch: double-double-array-branch
Changes since 1.1: +1799 -0 lines
Add implementation of special functions for double-double-float.  More
testing required, but basic functionality works.

code/irrat-dd.lisp:
o New file which implements all the required special functions for
  double-double.  Modify existing COMPLEX-<foo> functions to handle
  double-double numbers.

code/irrat.lisp:
o Update HANDLE-REALS to handle double-double float case.
o Update EXPT for double-double float.  (But negative number to
  non-integer power not working yet.)
o LOG handles double-double, but not 2-arg log yet.
o SQRT handles double-double, including complex result.
o ASIN handles double-double.
o ACOS handles double-double.
o ATAN handles double-double.
o ACOSH handles double-double.
o ATANH handles double-double.
o Adjust declaration for SQUARE, SCALB, LOGB-FINITE, and LOGB to allow
  any float type.
o COMPLEX-SQRT handles double-doubles.
o COMPLEX-LOG handles double-doubles.
o COMPLEX-ATANH handles double-doubles.
o COMPLEX-TANH handles double-doubles.
o COMPLEX-ACOS handles double-doubles.
o COMPLEX-ASIN handles double-doubles.
o COMPLEX-ASINH handles double-doubles.
o COMPLEX-ATAN handles double-doubles.
o COMPLEX-TAN handles double-doubles.

tools/worldbuild.lisp:
o Load irrat-dd.

tools/worldcom.lisp:
o Compile irrat-dd.
1 ;;; -*- Mode: Lisp; Package: KERNEL; Log: code.log -*-
2 ;;;
3 ;;; **********************************************************************
4 ;;; This code was written as part of the CMU Common Lisp project at
5 ;;; Carnegie Mellon University, and has been placed in the public domain.
6 ;;;
7 (ext:file-comment
8 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/irrat-dd.lisp,v 1.1.2.1 2006/06/29 01:28:02 rtoy Exp $")
9 ;;;
10 ;;; **********************************************************************
11 ;;;
12 ;;; This file contains all the irrational functions for double-double
13 ;;; float. The algorithms and coefficients come from the Cephes math
14 ;;; library. See http://www.moshier.net/#Cephes and
15 ;;; http://www.netlib.org/cephes/.
16 ;;;
17 ;;; Author: Raymond Toy.
18 ;;;
19
20 (in-package "KERNEL")
21
22 ;;;; Random constants, utility functions.
23
24 #+nil
25 (defconstant max-log
26 1.1356523406294143949491931077970764891253w4
27 "log(most-positive-double-double-float)")
28
29 #+nil
30 (defconstant min-log
31 -1.143276959615573793352782661133116431383730w4
32 "log(least-positive-double-double-float")
33
34 (defconstant max-log
35 7.0978271289338399678773454114191w2
36 "log(most-positive-double-double-float)")
37
38 (defconstant min-log
39 -7.4444007192138126231410729844608w2
40 "log(least-positive-double-double-float")
41
42
43 (defconstant loge2
44 0.6931471805599453094172321214581765680755001w0
45 "log(2)")
46
47 (defconstant log2e
48 1.442695040888963407359924681001892137426646w0
49 "Log base 2 of e")
50
51 (defconstant log2ea
52 4.4269504088896340735992468100189213742664595w-1
53 "log2(e)-1")
54
55 (defconstant dd-pi
56 3.141592653589793238462643383279502884197169w0
57 "Pi")
58
59 (defconstant dd-pi/2
60 1.570796326794896619231321691639751442098585w0
61 "Pi/2")
62
63 (defconstant dd-pi/4
64 0.7853981633974483096156608458198757210492923w0
65 "Pi/4")
66
67 ;; log2-c1 and log-c2 are log(2) arranged in such a way that log2-c1 +
68 ;; log2-c2 is log(2) to an accuracy greater than double-double-float.
69 (defconstant log2-c1
70 6.93145751953125w-1)
71
72 (defconstant log2-c2
73 1.428606820309417232121458176568075500134w-6)
74
75 (defconstant sqrt-1/2
76 0.7071067811865475244008443621048490392848w0
77 "Sqrt(2)")
78
79 ;; Evaluate polynomial
80 (defun poly-eval (x coef)
81 (declare (type double-double-float x)
82 (type (simple-array double-double-float (*)) coef))
83 ;; y = C0 + C1*x + C2*x^2 + ...
84 ;;
85 ;; But coefficients are stored in reverse (descending powers) order:
86 ;; coef[0] = CN, ..., coef[N] = C0.
87 (let ((y 0w0))
88 (declare (type double-double-float y))
89 (loop for c across coef do
90 (setf y (+ (* y x)
91 c)))
92 y))
93
94 (defun poly-eval-1 (x coef)
95 (declare (type double-double-float x)
96 (type (simple-array double-double-float (*)) coef))
97 ;; Like poly-eval, except it assumes coef[N] = 1 and is omitted.
98 (let ((y 1w0))
99 (declare (type double-double-float y))
100 (loop for c across coef do
101 (setf y (+ (* y x)
102 c)))
103 y))
104
105
106 ;;; dd-%expm1
107 ;;;
108 ;;; exp(x)-1
109 ;;; Range reduction is accomplished by separating the argument
110 ;;; into an integer k and fraction f such that
111 ;;;
112 ;;; x k f
113 ;;; e = 2 e.
114 ;;;
115 ;;; An expansion x + .5 x^2 + x^3 R(x) approximates exp(f) - 1
116 ;;; in the basic range [-0.5 ln 2, 0.5 ln 2].
117 ;;;
118 ;;;
119 ;;; ACCURACY:
120 ;;;
121 ;;; Relative error:
122 ;;; arithmetic domain # trials peak rms
123 ;;; IEEE -79,+MAXLOG 100,000 1.7e-34 4.5e-35
124 ;;;
125
126 (let ((p (make-array 8 :element-type 'double-double-float
127 :initial-contents '(
128 -4.888737542888633647784737721812546636240w-1
129 4.401308817383362136048032038528753151144w1
130 -1.716772506388927649032068540558788106762w3
131 4.578962475841642634225390068461943438441w4
132 -7.212432713558031519943281748462837065308w5
133 8.944630806357575461578107295909719817253w6
134 -5.722847283900608941516165725053359168840w7
135 2.943520915569954073888921213330863757240w8)))
136 (q (make-array 8 :element-type 'double-double-float
137 :initial-contents '(
138 ;; 1.000000000000000000000000000000000000000w0
139 1.766112549341972444333352727998584753865w9
140 -7.848989743695296475743081255027098295771w8
141 1.615869009634292424463780387327037251069w8
142 -2.019684072836541751428967854947019415698w7
143 1.682912729190313538934190635536631941751w6
144 -9.615511549171441430850103489315371768998w4
145 3.697714952261803935521187272204485251835w3
146 -8.802340681794263968892934703309274564037w1)))
147 ;; ln 2^-114
148 (minarg -7.9018778583833765273564461846232128760607w1))
149 (declare (type double-double-float minarg))
150 ;; exp(x) - 1 = x + 0.5*x^2+x^3*P(x)/Q(x), for -log(2)/2 < x <
151 ;; log(2)/2, where the coefficients of P and Q are given Pn and Qn
152 ;; above. Theoretical peak relative error = 8.1e-36.
153 (defun dd-%expm1 (x)
154 "exp(x) - 1"
155 (declare (type double-double-float x)
156 (optimize (speed 3)))
157 ;; Range reduction is accomplished by separating the argument
158 ;; into an integer k and fraction f such that
159 ;;
160 ;; x k f
161 ;; e = 2 e.
162 ;;
163 ;; An expansion x + .5 x^2 + x^3 R(x) approximates exp(f) - 1
164 ;; in the basic range [-0.5 ln 2, 0.5 ln 2].
165 (when (> x max-log)
166 (return-from dd-%expm1
167 (kernel:%make-double-double-float #.ext:double-float-positive-infinity
168 #.ext:double-float-positive-infinity)))
169 (when (< x minarg)
170 (return-from dd-%expm1 -1w0))
171
172 ;; Express x = ln(2)*(k+remainder), remainder not exceeding 1/2
173 (let* ((xx (+ log2-c1 log2-c2))
174 (k (floor (+ 1/2 (/ (the (double-double-float * 710w0) x) xx))))
175 (px (coerce k 'double-double-float))
176 (qx 0w0))
177 (declare (type double-double-float xx px qx))
178 ;; remainder times ln 2
179 (decf x (* px log2-c1))
180 (decf x (* px log2-c2))
181
182 ;; Approximate exp(remainder*ln(2))
183 #+nil
184 (setf px (* x
185 (+ p0
186 (* x
187 (+ p1
188 (* x
189 (+ p2
190 (* x
191 (+ p3
192 (* x
193 (+ p4
194 (* x
195 (+ p5
196 (* x
197 (+ p6
198 (* x p7))))))))))))))))
199 #+nil
200 (setf qx (+ q0
201 (* x
202 (+ q1
203 (* x
204 (+ q2
205 (* x
206 (+ q3
207 (* x
208 (+ q4
209 (* x
210 (+ q5
211 (* x
212 (+ q6
213 (* x
214 (+ x q7))))))))))))))))
215 (setf px (poly-eval x p))
216 (setf qx (poly-eval-1 x q))
217 (setf xx (* x x))
218 (setf qx (+ x (+ (* 0.5w0 xx)
219 (/ (* xx px)
220 qx))))
221 ;; exp(x) = exp(k*ln(2))*exp(remainder*ln(2)) = 2^k*exp(remainder*ln(2))
222 ;;
223 ;; We have qx = exp(remainder*ln(2))-1, so exp(x) - 1 =
224 ;; 2^k*(qx+1)-1 = 2^k*qx + 2^k - 1
225 (setf px (scale-float 1w0 k))
226 (+ (* px qx)
227 (- px 1)))))
228
229 ;;; dd-%exp
230 ;;; exp(x)
231 ;;; Range reduction is accomplished by separating the argument
232 ;;; into an integer k and fraction f such that
233 ;;;
234 ;;; x k f
235 ;;; e = 2 e.
236 ;;;
237 ;;; A Pade' form of degree 2/3 is used to approximate exp(f) - 1
238 ;;; in the basic range [-0.5 ln 2, 0.5 ln 2].
239 ;;;
240 ;;;
241 ;;; ACCURACY:
242 ;;;
243 ;;; Relative error:
244 ;;; arithmetic domain # trials peak rms
245 ;;; IEEE +-MAXLOG 100,000 2.6e-34 8.6e-35
246 ;;;
247 ;;;
248 ;;; Error amplification in the exponential function can be
249 ;;; a serious matter. The error propagation involves
250 ;;; exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
251 ;;; which shows that a 1 lsb error in representing X produces
252 ;;; a relative error of X times 1 lsb in the function.
253 ;;; While the routine gives an accurate result for arguments
254 ;;; that are exactly represented by a long double precision
255 ;;; computer number, the result contains amplified roundoff
256 ;;; error for large arguments not exactly represented.
257
258 (let ((p (make-array 5 :element-type 'double-double-float
259 :initial-contents
260 '(
261 3.279723985560247033712687707263393506266w-10
262 6.141506007208645008909088812338454698548w-7
263 2.708775201978218837374512615596512792224w-4
264 3.508710990737834361215404761139478627390w-2
265 9.999999999999999999999999999999999998502w-1
266 )))
267 (q (make-array 6 :element-type 'double-double-float
268 :initial-contents
269 '(
270 2.980756652081995192255342779918052538681w-12
271 1.771372078166251484503904874657985291164w-8
272 1.504792651814944826817779302637284053660w-5
273 3.611828913847589925056132680618007270344w-3
274 2.368408864814233538909747618894558968880w-1
275 2.000000000000000000000000000000000000150w0
276 )))
277 ;; C1 + C2 = ln 2
278 (c1 (- log2-c1))
279 (c2 (- log2-c2)))
280 (declare (type double-double-float c1 c2))
281 ;; p and q are Pade coefficients for exp(x) - 1. Theoretical peak
282 ;; relative error = 2.2e-37, relative peak error spread = 9.2e-38.
283 (defun dd-%exp (x)
284 (declare (type double-double-float x)
285 (optimize (speed 3)))
286 (when (> x max-log)
287 (return-from dd-%exp
288 (kernel:%make-double-double-float #.ext:double-float-positive-infinity
289 #.ext:double-float-positive-infinity)))
290 (when (< x min-log)
291 (return-from dd-%exp 0w0))
292 ;; Express
293 ;; exp(x) = exp(g)*2^n
294 ;; = exp(g)*exp(n*log(2))
295 ;; = exp(g + n * log2)
296 (let* ((n (floor (+ 0.5w0 (* x log2e))))
297 (px (coerce n 'double-double-float)))
298 (declare (type double-double-float px))
299 (incf x (* px c1))
300 (incf x (* px c2))
301 ;; Rational approx for exponential of fractional part:
302 ;; exp(x) = 1+2*x*p(x^2)/(q(x^2)-p(x^2))
303 ;;
304 ;; [The above comment seems not to match the code below. But
305 ;; perhaps q(x) isn't what I think it is.]
306 (let ((xx (* x x)))
307 (setf px (* x (poly-eval xx p)))
308 (setf xx (poly-eval xx q))
309 (setf x (/ px (- xx px)))
310 (setf x (+ 1 x x))
311 (scale-float x n)))))
312
313
314 ;;; dd-%log
315 ;;; log(x), natural logarithm.
316 ;;; The argument is separated into its exponent and fractional
317 ;;; parts. If the exponent is between -1 and +1, the logarithm
318 ;;; of the fraction is approximated by
319 ;;;
320 ;;; log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
321 ;;;
322 ;;; Otherwise, setting z = 2(x-1)/x+1),
323 ;;;
324 ;;; log(x) = z + z**3 P(z)/Q(z).
325 ;;;
326 ;;;
327 ;;;
328 ;;; ACCURACY:
329 ;;;
330 ;;; Relative error:
331 ;;; arithmetic domain # trials peak rms
332 ;;; IEEE exp(+-MAXLOGL) 36,000 9.5e-35 4.1e-35
333 ;;;
334
335 (let ((p (make-array 13 :element-type 'double-double-float
336 :initial-contents
337 '(1.538612243596254322971797716843006400388w-6
338 4.998469661968096229986658302195402690910w-1
339 2.321125933898420063925789532045674660756w1
340 4.114517881637811823002128927449878962058w2
341 3.824952356185897735160588078446136783779w3
342 2.128857716871515081352991964243375186031w4
343 7.594356839258970405033155585486712125861w4
344 1.797628303815655343403735250238293741397w5
345 2.854829159639697837788887080758954924001w5
346 3.007007295140399532324943111654767187848w5
347 2.014652742082537582487669938141683759923w5
348 7.771154681358524243729929227226708890930w4
349 1.313572404063446165910279910527789794488w4)))
350 (q (make-array 12 :element-type 'double-double-float
351 :initial-contents
352 '(
353 ;; 1.000000000000000000000000000000000000000w0
354 4.839208193348159620282142911143429644326w1
355 9.104928120962988414618126155557301584078w2
356 9.147150349299596453976674231612674085381w3
357 5.605842085972455027590989944010492125825w4
358 2.248234257620569139969141618556349415120w5
359 6.132189329546557743179177159925690841200w5
360 1.158019977462989115839826904108208787040w6
361 1.514882452993549494932585972882995548426w6
362 1.347518538384329112529391120390701166528w6
363 7.777690340007566932935753241556479363645w5
364 2.626900195321832660448791748036714883242w5
365 3.940717212190338497730839731583397586124w4)))
366 (r (make-array 6 :element-type 'double-double-float
367 :initial-contents
368 '(-8.828896441624934385266096344596648080902w-1
369 8.057002716646055371965756206836056074715w1
370 -2.024301798136027039250415126250455056397w3
371 2.048819892795278657810231591630928516206w4
372 -8.977257995689735303686582344659576526998w4
373 1.418134209872192732479751274970992665513w5
374 )))
375 (s (make-array 6 :element-type 'double-double-float
376 :initial-contents
377 '( ;; 1.000000000000000000000000000000000000000w0
378 -1.186359407982897997337150403816839480438w2
379 3.998526750980007367835804959888064681098w3
380 -5.748542087379434595104154610899551484314w4
381 4.001557694070773974936904547424676279307w5
382 -1.332535117259762928288745111081235577029w6
383 1.701761051846631278975701529965589676574w6))))
384 ;; p and q are coefficients in the expansion for log(1+x) = x -
385 ;; x^2/2 + x^3*p(x)/q(x), where 1/sqrt(2) <= 1+x < sqrt(2).
386 ;; Theoretical peak relative error = 5.3e-37, relative peak error
387 ;; spread = 2.3e-14.
388 ;;
389 ;; R and S are coefficients for log(x) = z + z^3*r(z^2)/s(z^2),
390 ;; where z = 2*(x-1)/(x+1), where 1/sqrt(2) <= x < sqrt(2).
391 ;; Theoretical peak relative error = 1.1e-35, relative peak error
392 ;; spread 1.1e-9.
393 ;;
394 ;;
395 (defun dd-%log (x)
396 (declare (type double-double-float x)
397 (optimize (speed 3)))
398 ;; Separate mantissa from exponent
399 (multiple-value-bind (x e)
400 (decode-float x)
401 (declare (type double-double-float x)
402 (type double-float-exponent e))
403 (let ((z 0w0)
404 (y 0w0))
405 (declare (type double-double-float z y))
406 ;; Logarithm using log(x) = z + z^3*P(z^2)/Q(z^2); where z = 2*(x-1)/x + 1;
407 (cond ((or (> e 2)
408 (< e -2))
409 (cond ((< x sqrt-1/2)
410 ;; 2*(2*x-1)/(2*x+1)
411 (decf e)
412 (setf z (- x 0.5w0))
413 (setf y (+ (* 0.5w0 z) 0.5w0)))
414 (t
415 ;; 2*(x-1)/(x+1)
416 (setf z (- x 0.5w0))
417 (decf z 0.5w0)
418 (setf y (+ (* 0.5w0 x) 0.5w0))))
419 (setf x (/ z y))
420 (setf z (* x x))
421 (setf z (* x (/ (* z (poly-eval z r))
422 (poly-eval-1 z s))))
423 (incf z (* e log2-c2))
424 (incf z x)
425 (incf z (* e log2-c1)))
426 (t
427 ;; Log using log(1+x) = x - 0.5*x^2+x^3*p(x)/q(x)
428 (cond ((< x sqrt-1/2)
429 (decf e)
430 (setf x (- (scale-float x 1) 1)))
431 (t
432 (decf x)))
433 (setf z (* x x))
434 (setf y (* x (* z (/ (poly-eval x p)
435 (poly-eval-1 x q)))))
436 (incf y (* e log2-c2))
437 ;; z = y - 0.5*z
438 (setf z (- y (scale-float z -1)))
439 (incf z x)
440 (incf z (* e log2-c1))
441 ))))))
442
443 ;;; dd-%log1p(x)
444 ;;; log(1+x)
445 ;;;
446 ;;; The argument 1+x is separated into its exponent and fractional
447 ;;; parts. If the exponent is between -1 and +1, the logarithm
448 ;;; of the fraction is approximated by
449 ;;;
450 ;;; log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
451 ;;;
452 ;;; Otherwise, setting z = 2(x-1)/x+1),
453 ;;;
454 ;;; log(x) = z + z^3 P(z)/Q(z).
455 ;;;
456 ;;;
457 ;;;
458 ;;; ACCURACY:
459 ;;;
460 ;;; Relative error:
461 ;;; arithmetic domain # trials peak rms
462 ;;; IEEE -1, 8 100000 1.9e-34 4.3e-35
463
464 (let ((p (make-array 13 :element-type 'double-double-float
465 :initial-contents '(
466 1.538612243596254322971797716843006400388w-6
467 4.998469661968096229986658302195402690910w-1
468 2.321125933898420063925789532045674660756w1
469 4.114517881637811823002128927449878962058w2
470 3.824952356185897735160588078446136783779w3
471 2.128857716871515081352991964243375186031w4
472 7.594356839258970405033155585486712125861w4
473 1.797628303815655343403735250238293741397w5
474 2.854829159639697837788887080758954924001w5
475 3.007007295140399532324943111654767187848w5
476 2.014652742082537582487669938141683759923w5
477 7.771154681358524243729929227226708890930w4
478 1.313572404063446165910279910527789794488w4
479 )))
480 (q (make-array 12 :element-type 'double-double-float
481 :initial-contents '(
482 ;; 1.000000000000000000000000000000000000000w0
483 4.839208193348159620282142911143429644326w1
484 9.104928120962988414618126155557301584078w2
485 9.147150349299596453976674231612674085381w3
486 5.605842085972455027590989944010492125825w4
487 2.248234257620569139969141618556349415120w5
488 6.132189329546557743179177159925690841200w5
489 1.158019977462989115839826904108208787040w6
490 1.514882452993549494932585972882995548426w6
491 1.347518538384329112529391120390701166528w6
492 7.777690340007566932935753241556479363645w5
493 2.626900195321832660448791748036714883242w5
494 3.940717212190338497730839731583397586124w4
495 )))
496 (r (make-array 6 :element-type 'double-double-float
497 :initial-contents '(-8.828896441624934385266096344596648080902w-1
498 8.057002716646055371965756206836056074715w1
499 -2.024301798136027039250415126250455056397w3
500 2.048819892795278657810231591630928516206w4
501 -8.977257995689735303686582344659576526998w4
502 1.418134209872192732479751274970992665513w5
503 )))
504 (s (make-array 6 :element-type 'double-double-float
505 :initial-contents '(
506 ;; 1.000000000000000000000000000000000000000w0
507 -1.186359407982897997337150403816839480438w2
508 3.998526750980007367835804959888064681098w3
509 -5.748542087379434595104154610899551484314w4
510 4.001557694070773974936904547424676279307w5
511 -1.332535117259762928288745111081235577029w6
512 1.701761051846631278975701529965589676574w6
513 )))
514 )
515 ;; Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
516 ;; 1/sqrt(2) <= 1+x < sqrt(2)
517 ;; Theoretical peak relative error = 5.3e-37,
518 ;; relative peak error spread = 2.3e-14
519 (defun dd-%log1p (xm1)
520 (declare (type double-double-float x))
521 (let ((x (+ xm1 1))
522 (z 0w0)
523 (y 0w0))
524 (declare (type double-double-float x y z))
525 (multiple-value-bind (x e)
526 (decode-float x)
527 (cond ((or (> e 2)
528 (< e -2))
529 ;; Log using log(x) = z + z^3*P(z^2)/Q(z^2)
530 (cond ((< x sqrt-1/2)
531 (decf e)
532 (setf z (- x 0.5w0))
533 (setf y (+ 0.5w0 (* 0.5w0 z))))
534 (t
535 (setf z (- x 0.5w0))
536 (decf z 0.5w0)
537 (setf y (+ 0.5w0 (* 0.5w0 x)))))
538 (setf x (/ z y))
539 (setf z (* x x))
540 (setf z (* x (/ (* z (poly-eval z r))
541 (poly-eval-1 z s))))
542 (incf z (* e log2-c2))
543 (incf z x)
544 (incf z (* e log2-c1)))
545 (t
546 ;; Log using log(1+x) = x - 0.5*x^2 + x^3*p(x)/q(x)
547 (cond ((< x sqrt-1/2)
548 (decf e)
549 (if (/= e 0)
550 (setf x (- (* 2 x) 1))
551 (setf x xm1)))
552 (t
553 (if (/= e 0)
554 (decf x 1)
555 (setf x xm1))))
556 (setf z (* x x))
557 (setf y (* x (/ (* z (poly-eval x p))
558 (poly-eval-1 x q))))
559 (incf y (* e log2-c2))
560 (setf z (- y (* 0.5w0 z)))
561 (incf z x)
562 (incf z (* e log2-c1))))))))
563
564 (let ((P (make-array 6 :element-type 'double-double-float
565 :initial-contents '(
566 1.622194395724068297909052717437740288268w3
567 1.124862584587770079742188354390171794549w6
568 3.047548980769660162696832999871894196102w8
569 3.966215348072348368191433063260384329745w10
570 2.375869584584371194838551715348965605295w12
571 6.482835792103233269752264509192030816323w13
572 )))
573 (Q (make-array 6 :element-type 'double-double-float
574 :initial-contents '(
575 ;; 1.000000000000000000000000000000000000000w0 */
576 -9.101683853129357776079049616394849086007w2
577 4.486400519836461218634448973793765123186w5
578 -1.492531313030440305095318968983514314656w8
579 3.457771488856930054902696708717192082887w10
580 -5.193289868803472640225483235513427062460w12
581 3.889701475261939961851358705515223019890w14))))
582 (defun dd-%sinh (x)
583 (declare (type double-double-float x)
584 (optimize (speed 3)))
585 (let ((a (abs x)))
586 (declare (type double-double-float a))
587 (cond ((> a 1)
588 (setf a (dd-%exp a))
589 (setf a (- (* 0.5w0 a) (/ 0.5w0 a)))
590 (if (< x 0)
591 (- a)
592 a))
593 (t
594 (setf a (* a a))
595 (let ((pp (poly-eval a p))
596 (qq (poly-eval-1 a q)))
597 (+ x (* x a
598 (/ pp qq)))))))))
599
600 #+nil
601 (defun dd-%sinh (x)
602 (declare (type double-double-float x))
603 (let ((a (abs x)))
604 (declare (type double-double-float a))
605 (cond ((> a 1)
606 (setf a (dd-%exp a))
607 (setf a (- (* 0.5w0 a) (/ 0.5w0 a)))
608 (if (< x 0)
609 (- a)
610 a))
611 (t
612 ;; Use sinh(x) = 1/2*(D(x)+D(x)/(1+D(x))), here D(x) =
613 ;; exp(x)-1.
614 (let ((d (dd-%expm1 x)))
615 (* 0.5w0 (+ d (/ d (+ 1 d)))))))))
616
617 (defun dd-%cosh (x)
618 (declare (type double-double-float x))
619 (let ((y (dd-%exp x)))
620 (scale-float (+ y (/ y)) -1)))
621
622 (let ((P (make-array 6 :element-type 'double-double-float
623 :initial-contents '(-6.505693197948351084912624750702492767503w-6
624 -9.804083860188429726356968570322356183383w-1
625 -5.055287638900473250703725789725376004355w2
626 -7.307477148073823966594990496301416814519w4
627 -3.531606586182691280701462523692471322688w6
628 -4.551377146142783468144190926206842300707w7
629 )))
630 (Q (make-array 5 :element-type 'double-double-float
631 :initial-contents '(
632 ;; 1.000000000000000000000000000000000000000w0 */
633 5.334865598460027935735737253027154828002w2
634 8.058475607422391042912151298751537172870w4
635 4.197073523796142343374222405869721575491w6
636 6.521134551226147545983467868553677881771w7
637 1.365413143842835040443257277862054198329w8))))
638 (defun dd-%tanh (x)
639 (declare (type double-double-float x)
640 (optimize (speed 3)))
641 (let ((z (abs x)))
642 (declare (type double-double-float z))
643 (cond ((> z (* 0.5w0 max-log))
644 (if (> x 0)
645 1w0
646 -1w0))
647 ((> z 0.625w0)
648 (let ((s (dd-%exp (* 2 z))))
649 (setf z (- 1 (/ 2 (+ s 1))))
650 (if (minusp x)
651 (- z)
652 z)))
653 (t
654 (let* ((s (* x x)))
655 (declare (optimize speed))
656 (setf z (/ (poly-eval s p)
657 (poly-eval-1 s q)))
658 (setf z (* x s z))
659 (+ x z)))))))
660
661 (let ((P (make-array 10 :element-type 'double-double-float
662 :initial-contents '(
663 -9.217569843805850417698565442251656375681w-1
664 5.321929116410615470118183794063211260728w1
665 -9.139522976807685333981548145417830690552w2
666 7.204314536952949779101646454146682033772w3
667 -3.097809640165146436529075324081668598891w4
668 7.865376554210973897486215630898496100534w4
669 -1.211716814094785128366087489224821937203w5
670 1.112669508789123834670923967462068457013w5
671 -5.600242872292477863751728708249167956542w4
672 1.188901082233997739779618679364295772810w4
673 )))
674 (Q (make-array 10 :element-type 'double-double-float
675 :initial-contents '(
676 ;; 1.000000000000000000000000000000000000000w0 */
677 -6.807348436010016270202879229504392062418w1
678 1.386763299649315831625106608182196351693w3
679 -1.310805752656879543134785263832907269320w4
680 6.872174720355764193772953852564737816928w4
681 -2.181008360536226513009076189881617939380w5
682 4.362736119602298592874941767284979857248w5
683 -5.535251007539393347687001489396152923502w5
684 4.321594849688346708841188057241308805551w5
685 -1.894075056489862952285849974761239845873w5
686 3.566703246701993219338856038092901974725w4
687 ))))
688 (defun dd-%atanh (x)
689 (declare (type double-double-float x)
690 (optimize (speed 3)))
691 (cond ((minusp x)
692 (- (the double-double-float (dd-%atanh (- x)))))
693 ((< x 1w-12)
694 x)
695 ((< x 0.5w0)
696 (let ((z (* x x)))
697 (+ x (* x z (/ (poly-eval z p)
698 (poly-eval-1 z q))))))
699 (t
700 (* 0.5w0 (dd-%log (/ (+ 1 x)
701 (- 1 x))))))))
702
703
704 (let ((P (make-array 9 :element-type 'double-double-float
705 :initial-contents '(
706 -8.104404283317298189545629468767571317688w-1
707 -4.954206127425209147110732546633675599008w1
708 -8.438175619831548439550086251740438689853w2
709 -6.269710069245210459536983820505214648057w3
710 -2.418935474493501382372711518024193326434w4
711 -5.208121780431312783866941311277024486498w4
712 -6.302755086521614763280617114866439227971w4
713 -4.003566436224198252093684987323233921339w4
714 -1.037690841528359305134494613113086980551w4
715 )))
716 (Q (make-array 9 :element-type 'double-double-float
717 :initial-contents '(
718 ;; 1.000000000000000000000000000000000000000w0 */
719 8.175806439951395194771977809279448392548w1
720 1.822215299975696008284027212745010251320w3
721 1.772040003462901790853111853838978236828w4
722 9.077625379864046240143413577745818879353w4
723 2.675554475070211205153169988669677418808w5
724 4.689758557916492969463473819426544383586w5
725 4.821923684550711724710891114802924039911w5
726 2.682316388947175963642524537892687560973w5
727 6.226145049170155830806967678679167550122w4))))
728 (defun dd-%asinh (x)
729 (declare (type double-double-float x)
730 (optimize (speed 3)))
731 (cond ((minusp x)
732 (- (the double-double-float (dd-%asinh (- x)))))
733 #+nil
734 ((> x 1w10)
735 (+ loge2 (dd-%log x)))
736 ((< x 0.5w0)
737 (let* ((z (* x x))
738 (a (* z (/ (poly-eval z p)
739 (poly-eval-1 z q)))))
740 (+ (* a x) x)))
741 (t
742 (dd-%log (+ x (sqrt (1+ (* x x)))))))))
743
744 (let ((P (make-array 10 :element-type 'double-double-float
745 :initial-contents '(
746 1.895467874386341763387398084072833727168w-1
747 6.443902084393244878979969557171256604767w1
748 3.914593556594721458616408528941154205393w3
749 9.164040999602964494412169748897754668733w4
750 1.065909694792026382660307834723001543839w6
751 6.899169896709615182428217047370629406305w6
752 2.599781868717579447900896150777162652518w7
753 5.663733059389964024656501196827345337766w7
754 6.606302846870644033621560858582696134512w7
755 3.190482951215438078279772140481195200593w7
756 )))
757 (Q (make-array 9 :element-type 'double-double-float
758 :initial-contents '(
759 ;; 1.000000000000000000000000000000000000000w0 */
760 1.635418024331924674147953764918262009321w2
761 7.290983678312632723073455563799692165828w3
762 1.418207894088607063257675159183397062114w5
763 1.453154285419072886840913424715826321357w6
764 8.566841438576725234955968880501739464425w6
765 3.003448667795089562511136059766833630017w7
766 6.176592872899557661256383958395266919654w7
767 6.872176426138597206811541870289420510034w7
768 3.190482951215438078279772140481195226621w7
769 ))))
770 (defun dd-%acosh (x)
771 (declare (type double-double-float x)
772 (optimize (speed 3)))
773 (cond ((> x 1w17)
774 (+ loge2 (dd-%log x)))
775 (t
776 (let ((z (- x 1)))
777 (cond ((< z 0.5w0)
778 (* (sqrt (* 2 z))
779 (/ (poly-eval z p)
780 (poly-eval-1 z q))))
781 (t
782 (dd-%log (+ x (sqrt (* z (+ 1 x))))))))))))
783
784 (let ((P (make-array 10 :element-type 'double-double-float
785 :initial-contents '(
786 -8.067112765482705313585175280952515549833w-1
787 4.845649797786849136525020822000172350977w1
788 -8.510195404865297879959793548843395926847w2
789 6.815196841370292688574521445731895826485w3
790 -2.967135182120339728996157454994675519735w4
791 7.612250656518818109652985996692466409670w4
792 -1.183360579752620455689557157684221905030w5
793 1.095432262510413338755837156377401348063w5
794 -5.554124580991113991999636773382495788705w4
795 1.187132626694762543537732514905488896985w4
796 )))
797 (Q (make-array 10 :element-type 'double-double-float
798 :initial-contents '(
799 ;; 1.000000000000000000000000000000000000000w0 */
800 -8.005471061732009595694099899234272342478w1
801 1.817324228942812880965069608562483918025w3
802 -1.867017317425756524289537002141956583706w4
803 1.048196619402464497478959760337779705622w5
804 -3.527040897897253459022458866536165564103w5
805 7.426302422018858001691440351763370029242w5
806 -9.863068411558756277454631976667880674474w5
807 8.025654653926121907774766642393757364326w5
808 -3.653000557802254281954969843055623398839w5
809 7.122795760168575261226395089432959614179w4
810 ))))
811 (defun dd-%asin (x)
812 (declare (type double-double-float x)
813 (optimize (speed 3)))
814 (cond ((minusp x)
815 (- (the double-double-float (dd-%asin (- x)))))
816 #+nil
817 ((< x 1w-8)
818 x)
819 (t
820 (let ((flag nil)
821 (z 0w0)
822 (zz 0w0))
823 (declare (type double-double-float z zz))
824 (cond ((> x 0.5w0)
825 (setf zz (- 0.5w0 x))
826 (setf zz (scale-float (+ zz 0.5w0) -1))
827 (setf z (sqrt zz))
828 (setf flag t))
829 (t
830 (setf z x)
831 (setf zz (* z z))
832 (setf flag nil)))
833 (let ((p (* zz (/ (poly-eval zz p)
834 (poly-eval-1 zz q)))))
835 (setf z (+ (* z p) z))
836 (when flag
837 (setf z (+ z z))
838 (setf z (- dd-pi/2 z)))
839 z))))))
840
841 (defun dd-%acos (x)
842 (declare (type double-double-float x)
843 (optimize (speed 3)))
844 (cond ((< x -0.5w0)
845 (- dd-pi
846 (* 2 (dd-%asin (sqrt (* 0.5w0 (+ 1 x)))))))
847 ((> x 0.5w0)
848 (* 2 (dd-%asin (sqrt (* 0.5w0 (- 1w0 x))))))
849 (t
850 (- dd-pi/2 (dd-%asin x)))))
851
852 (let ((P (make-array 9 :element-type 'double-double-float
853 :initial-contents '(
854 -6.635810778635296712545011270011752799963w-4
855 -8.768423468036849091777415076702113400070w-1
856 -2.548067867495502632615671450650071218995w1
857 -2.497759878476618348858065206895055957104w2
858 -1.148164399808514330375280133523543970854w3
859 -2.792272753241044941703278827346430350236w3
860 -3.696264445691821235400930243493001671932w3
861 -2.514829758941713674909996882101723647996w3
862 -6.880597774405940432145577545328795037141w2
863 )))
864 (Q (make-array 8 :element-type 'double-double-float
865 :initial-contents '(
866 ;; 1.000000000000000000000000000000000000000w0 */
867 3.566239794444800849656497338030115886153w1
868 4.308348370818927353321556740027020068897w2
869 2.494680540950601626662048893678584497900w3
870 7.928572347062145288093560392463784743935w3
871 1.458510242529987155225086911411015961174w4
872 1.547394317752562611786521896296215170819w4
873 8.782996876218210302516194604424986107121w3
874 2.064179332321782129643673263598686441900w3
875 )))
876
877 ;; tan( 3*pi/8 )
878 (T3P8 2.414213562373095048801688724209698078569672w0)
879 ;; tan( pi/8 )
880 (TP8 0.414213562373095048801688724209698078569672w0))
881 (declare (type double-double-float t3p8 tp8))
882 (defun dd-%atan (x)
883 (declare (type double-double-float x)
884 (optimize (speed 3)))
885 (when (minusp x)
886 (return-from dd-%atan (- (the double-double-float (dd-%atan (- x))))))
887 ;; range reduction
888 (let ((y 0w0))
889 (declare (type double-double-float y))
890 (cond ((> x t3p8)
891 (setf y dd-pi/2)
892 (setf x (/ -1 x)))
893 ((> x tp8)
894 (setf y dd-pi/4)
895 (setf x (/ (- x 1)
896 (+ x 1))))
897 (t
898 (setf y 0w0)))
899 ;; Rational form in x^2
900 (let ((z (* x x)))
901 (setf y (+ y
902 (* (/ (poly-eval z p)
903 (poly-eval-1 z q))
904 z x)
905 x))
906 y))))
907
908 (defun dd-%atan2 (y x)
909 (declare (type double-double-float x y)
910 (optimize (speed 3)))
911 (let ((code 0)
912 (w 0w0))
913 (declare (type (integer 0 3) code)
914 (type double-double-float w))
915 (when (minusp x)
916 (setf code 2))
917 (when (minusp y)
918 (setf code (logior code 1)))
919 (when (zerop x)
920 (unless (zerop (logand code 1))
921 (return-from dd-%atan2 (- dd-pi/2)))
922 (when (zerop y)
923 (return-from dd-%atan2 0w0))
924 (return-from dd-%atan2 dd-pi/2))
925 (when (zerop y)
926 (return-from dd-%atan2
927 (if (zerop (logand code 2))
928 0w0
929 dd-pi)))
930 (setf w (case code
931 (0 0w0)
932 (1 0w0)
933 (2 dd-pi)
934 (3 (- dd-pi))))
935
936 (+ w (dd-%atan (/ y x)))))
937
938
939 ;;
940 ;; Here are the fractional digits of pi, in hex.
941 ;;
942 ;; 243f6a8885 a308d31319 8a2e037073 44a4093822 299f31d008 2efa98ec4e 6c89452821 e638d01377 be5466cf34 e90c6cc0ac
943 ;; 29b7c97c50 dd3f84d5b5 b547091792 16d5d98979 fb1bd1310b a698dfb5ac 2ffd72dbd0 1adfb7b8e1 afed6a267e 96ba7c9045
944 ;;
945 ;; We want to express pi/4 in 3 parts such that the sum of the parts
946 ;; is exactly the value of pi to the higher precision. For
947 ;; double-double, we want 53 bits in each part.
948 ;;
949 ;; pi/4 in hex is
950 ;;
951 ;; C90FDAA22168C234C4C6628B80DC1CD129024E088A67CC74020BBEA63B139B22514A08798E3404DDEF9519B3CD3A431B302B0A6DF25F14374FE1356D6D51C245E485B576625E7EC6F44C42E9A637ED6B0BFF5CB6F406B7EDEE386BFB5A899FA5AE9F2411
952 ;;
953 ;;
954 #||
955 (let* ((frac #xC90FDAA22168C234C4C6628B80DC1CD129024E088A67CC74020BBEA63B139B22514A08798E3404DDEF9519B3CD3A431B302B0A6DF25F14374FE1356D6D51C245E485B576625E7EC6F44C42E9A637ED6B0BFF5CB6F406B7EDEE386BFB5A899FA5AE9F2411)
956 (len (integer-length frac)))
957 (format t "len = ~D bits~%" len)
958 (format t "pi/4 ~,,' ,106:b~%" (ldb (byte (* 4 106) (- len (* 4 106))) frac)))
959 ->
960 pi/4 11001001000011111101101010100010001000010110100011000 01000110100110001001100011001100010100010111000000011 01110000011100110100010010100100000010010011100000100 01000101001100111110011000111010000000010000010111011
961
962 ||#
963 ;; d1 = (scale-float (float #b1100100100001111110110101010001000100001011010001100001000110100110001001100011001100010100010111000000011 1w0) -106)
964 ;; d2 = (scale-float (float #b0111000001110011010001001010010000001001001110000010001000101001100111110011000111010000000010000010111011 1w0) (* 2 -106))
965 ;; d3 = (scale-float (float #b1110101001100011101100010011100110110010001001010001010010100000100001111001100011100011010000000100110111 1w0) (* 3 -106))
966 ;; d4 = (scale-float (float #b0111101111100101010001100110110011110011010011101001000011000110110011000000101011000010100110110111110010 1w0) (* 4 -106))
967
968 ;; sin(x) = x + x^3 P(x^2)
969 ;; Theoretical peak relative error = 5.6e-39
970 ;; relative peak error spread = 1.7e-9
971
972 (defconstant sincof
973 (make-array 12 :element-type 'double-double-float
974 :initial-contents
975 '(
976 6.410290407010279602425714995528976754871w-26
977 -3.868105354403065333804959405965295962871w-23
978 1.957294039628045847156851410307133941611w-20
979 -8.220635246181818130416407184286068307901w-18
980 2.811457254345322887443598804951004537784w-15
981 -7.647163731819815869711749952353081768709w-13
982 1.605904383682161459812515654720205050216w-10
983 -2.505210838544171877505034150892770940116w-8
984 2.755731922398589065255731765498970284004w-6
985 -1.984126984126984126984126984045294307281w-4
986 8.333333333333333333333333333333119885283w-3
987 -1.666666666666666666666666666666666647199w-1
988 )))
989
990 ;; cos(x) = 1 - .5 x^2 + x^2 (x^2 P(x^2))
991 ;; Theoretical peak relative error = 2.1e-37,
992 ;; relative peak error spread = 1.4e-8
993
994 (defconstant coscof
995 (make-array 11 :element-type 'double-double-float
996 :initial-contents
997 '(
998 1.601961934248327059668321782499768648351w-24
999 -8.896621117922334603659240022184527001401w-22
1000 4.110317451243694098169570731967589555498w-19
1001 -1.561920696747074515985647487260202922160w-16
1002 4.779477332386900932514186378501779328195w-14
1003 -1.147074559772972328629102981460088437917w-11
1004 2.087675698786809897637922200570559726116w-9
1005 -2.755731922398589065255365968070684102298w-7
1006 2.480158730158730158730158440896461945271w-5
1007 -1.388888888888888888888888888765724370132w-3
1008 4.166666666666666666666666666666459301466w-2
1009 )))
1010
1011 (defconstant dp1
1012 (scale-float (float #b1100100100001111110110101010001000100001011010001100001000110100110001001100011001100010100010111000000011 1w0) -106))
1013
1014 (defconstant dp2
1015 (scale-float (float #b0111000001110011010001001010010000001001001110000010001000101001100111110011000111010000000010000010111011 1w0) (* 2 -106)))
1016
1017 (defconstant dp3
1018 (scale-float (float #b1110101001100011101100010011100110110010001001010001010010100000100001111001100011100011010000000100110111 1w0) (* 3 -106)))
1019
1020 (defconstant dp4
1021 (scale-float (float #b0111101111100101010001100110110011110011010011101001000011000110110011000000101011000010100110110111110010 1w0) (* 4 -106)))
1022
1023 (defun dd-%sin (x)
1024 (declare (type double-double-float x))
1025 (when (minusp x)
1026 (return-from dd-%sin (- (dd-%sin (- x)))))
1027 ;; y = integer part of x/(pi/4).
1028 (let* ((y (float (floor (/ x dd-pi/4)) 1w0))
1029 (z (scale-float y -4)))
1030 (declare (type double-double-float y z))
1031 (setf z (float (floor z) 1w0)) ; integer part of y/8
1032 (setf z (- y (scale-float z 4))) ; y - 16*(y/16)
1033
1034 (let ((j (truncate z))
1035 (sign 1))
1036 (unless (zerop (logand j 1))
1037 (incf j)
1038 (incf y))
1039 (setf j (logand j 7))
1040
1041 (when (> j 3)
1042 (setf sign (- sign))
1043 (decf j 4))
1044
1045 ;; Extended precision modular arithmetic
1046 (setf z (- (- (- x (* y dp1))
1047 (* y dp2))
1048 (* y dp3)))
1049 (let ((zz (* z z)))
1050 (if (or (= j 1)
1051 (= j 2))
1052 (setf y (+ (- 1 (scale-float zz -1))
1053 (* zz zz (poly-eval zz coscof))))
1054 (setf y (+ z (* z (* zz (poly-eval zz sincof))))))
1055 (if (< sign 0)
1056 (- y)
1057 y)))))
1058
1059 (defun dd-%cos (x)
1060 (declare (type double-double-float x)
1061 (optimize (speed 3)))
1062 (when (minusp x)
1063 (return-from dd-%cos (dd-%cos (- x))))
1064 ;; y = integer part of x/(pi/4).
1065 (let* ((y (float (floor (/ x dd-pi/4)) 1w0))
1066 (z (scale-float y -4)))
1067 (declare (type double-double-float y z))
1068 (setf z (float (floor z) 1w0)) ; integer part of y/8
1069 (setf z (- y (scale-float z 4))) ; y - 16*(y/16)
1070
1071 (let ((i (truncate z))
1072 (j 0)
1073 (sign 1))
1074 (unless (zerop (logand i 1))
1075 (incf i)
1076 (incf y))
1077 (setf j (logand i 7))
1078
1079 (when (> j 3)
1080 (setf sign (- sign))
1081 (decf j 4))
1082 (when (> j 1)
1083 (setf sign (- sign)))
1084
1085 ;; Extended precision modular arithmetic
1086 (setf z (- (- (- x (* y dp1))
1087 (* y dp2))
1088 (* y dp3)))
1089 (let ((zz (* z z)))
1090 (if (or (= j 1)
1091 (= j 2))
1092 (setf y (+ z (* z (* zz (poly-eval zz sincof)))))
1093 (setf y (+ (- 1 (scale-float zz -1))
1094 (* zz zz (poly-eval zz coscof)))))
1095 (if (< sign 0)
1096 (- y)
1097 y)))))
1098
1099 (let ((P (make-array 6 :element-type 'double-double-float
1100 :initial-contents
1101 '(
1102 -9.889929415807650724957118893791829849557w-1
1103 1.272297782199996882828849455156962260810w3
1104 -4.249691853501233575668486667664718192660w5
1105 5.160188250214037865511600561074819366815w7
1106 -2.307030822693734879744223131873392503321w9
1107 2.883414728874239697964612246732416606301w10
1108 )))
1109 (Q (make-array 6 :element-type 'double-double-float
1110 :initial-contents
1111 '(
1112 ;; 1.000000000000000000000000000000000000000w0 */
1113 -1.317243702830553658702531997959756728291w3
1114 4.529422062441341616231663543669583527923w5
1115 -5.733709132766856723608447733926138506824w7
1116 2.758476078803232151774723646710890525496w9
1117 -4.152206921457208101480801635640958361612w10
1118 8.650244186622719093893836740197250197602w10
1119 ))))
1120 (defun dd-tancot (xx cotflag)
1121 (declare (type double-double-float xx)
1122 (optimize (speed 3)))
1123 (let ((x 0w0)
1124 (sign 1))
1125 (declare (type double-double-float x)
1126 (type (integer -1 1) sign))
1127 (cond ((minusp xx)
1128 (setf x (- xx))
1129 (setf sign -1))
1130 (t
1131 (setf x xx)))
1132 (let* ((y (float (floor (/ x dd-pi/4)) 1w0))
1133 (z (scale-float y -4))
1134 (j 0))
1135 (declare (type double-double-float y z)
1136 (type fixnum j))
1137 (setf z (float (floor z) 1w0))
1138 (setf z (- y (scale-float z 4)))
1139
1140 (setf j (truncate z))
1141
1142 (unless (zerop (logand j 1))
1143 (incf j)
1144 (incf y))
1145
1146 (setf z (- (- (- x (* y dp1))
1147 (* y dp2))
1148 (* y dp3)))
1149 (let ((zz (* z z)))
1150 (if (> zz 1w-40)
1151 (setf y (+ z
1152 (* z (* zz (/ (poly-eval zz p)
1153 (poly-eval-1 zz q))))))
1154 (setf y z))
1155 (if (not (zerop (logand j 2)))
1156 (if cotflag
1157 (setf y (- y))
1158 (setf y (/ -1 y)))
1159 (if cotflag
1160 (setf y (/ y))))
1161 (if (< sign 0)
1162 (- y)
1163 y))))))
1164
1165 (defun dd-%tan (x)
1166 (declare (type double-double-float x))
1167 (dd-tancot x nil))
1168
1169 ;;; dd-%log2
1170 ;;; Base 2 logarithm.
1171
1172 (let ((P (make-array 13 :element-type 'double-double-float
1173 :initial-contents
1174 '(
1175 1.538612243596254322971797716843006400388w-6
1176 4.998469661968096229986658302195402690910w-1
1177 2.321125933898420063925789532045674660756w1
1178 4.114517881637811823002128927449878962058w2
1179 3.824952356185897735160588078446136783779w3
1180 2.128857716871515081352991964243375186031w4
1181 7.594356839258970405033155585486712125861w4
1182 1.797628303815655343403735250238293741397w5
1183 2.854829159639697837788887080758954924001w5
1184 3.007007295140399532324943111654767187848w5
1185 2.014652742082537582487669938141683759923w5
1186 7.771154681358524243729929227226708890930w4
1187 1.313572404063446165910279910527789794488w4
1188 )))
1189 (Q (make-array 12 :element-type 'double-double-float
1190 :initial-contents
1191 '(
1192 ;; 1.000000000000000000000000000000000000000w0
1193 4.839208193348159620282142911143429644326w1
1194 9.104928120962988414618126155557301584078w2
1195 9.147150349299596453976674231612674085381w3
1196 5.605842085972455027590989944010492125825w4
1197 2.248234257620569139969141618556349415120w5
1198 6.132189329546557743179177159925690841200w5
1199 1.158019977462989115839826904108208787040w6
1200 1.514882452993549494932585972882995548426w6
1201 1.347518538384329112529391120390701166528w6
1202 7.777690340007566932935753241556479363645w5
1203 2.626900195321832660448791748036714883242w5
1204 3.940717212190338497730839731583397586124w4
1205 )))
1206 ;; Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
1207 ;; where z = 2(x-1)/(x+1)
1208 ;; 1/sqrt(2) <= x < sqrt(2)
1209 ;; Theoretical peak relative error = 1.1e-35,
1210 ;; relative peak error spread 1.1e-9
1211 (R (make-array 6 :element-type 'double-double-float
1212 :initial-contents
1213 '(
1214 -8.828896441624934385266096344596648080902w-1
1215 8.057002716646055371965756206836056074715w1
1216 -2.024301798136027039250415126250455056397w3
1217 2.048819892795278657810231591630928516206w4
1218 -8.977257995689735303686582344659576526998w4
1219 1.418134209872192732479751274970992665513w5
1220 )))
1221 (S (make-array 6 :element-type 'double-double-float
1222 :initial-contents
1223 '(
1224 ;; 1.000000000000000000000000000000000000000w0 */
1225 -1.186359407982897997337150403816839480438w2
1226 3.998526750980007367835804959888064681098w3
1227 -5.748542087379434595104154610899551484314w4
1228 4.001557694070773974936904547424676279307w5
1229 -1.332535117259762928288745111081235577029w6
1230 1.701761051846631278975701529965589676574w6
1231 ))))
1232 (defun dd-%log2 (x)
1233 (declare (type double-double-float x)
1234 (optimize (speed 3)))
1235 (multiple-value-bind (x e)
1236 (decode-float x)
1237 (declare (type double-double-float x)
1238 (type double-float-exponent e))
1239 (let ((z 0w0)
1240 (y 0w0))
1241 (declare (type double-double-float z y))
1242 (cond ((or (> e 2)
1243 (< e -2))
1244 (cond ((< x sqrt-1/2)
1245 ;; 2*(2*x-1)/(2*x+1)
1246 (decf e)
1247 (setf z (- x 0.5w0))
1248 (setf y (+ (* 0.5w0 z) 0.5w0)))
1249 (t
1250 ;; 2*(x-1)/(x+1)
1251 (setf z (- x 0.5w0))
1252 (decf z 0.5w0)
1253 (setf y (+ (* 0.5w0 z) 0.5w0))))
1254 (setf x (/ z y))
1255 (setf z (* x x))
1256 (setf y (* x (/ (* z (poly-eval z r))
1257 (poly-eval-1 z s)))))
1258 (t
1259 (cond ((< x sqrt-1/2)
1260 (decf e)
1261 (setf x (- (scale-float x 1) 1)))
1262 (t
1263 (decf x)))
1264 (setf z (* x x))
1265 (setf y (* x (/ (* z (poly-eval x p))
1266 (poly-eval-1 x q))))
1267 (decf y (scale-float z -1))))
1268 ;; Multiply log of fraction by log2(e) and base 2 exponent by 1
1269 ;;
1270 ;; This sequence of operations is critical
1271 (setf z (* y log2ea))
1272 (setf z (+ z (* x log2ea)))
1273 (setf z (+ z y))
1274 (setf z (+ z x))
1275 (setf z (+ z e))
1276 z))))
1277
1278 ;;; dd-%exp2
1279 ;;; 2^x
1280
1281 (let ((P (make-array 5 :element-type 'double-double-float
1282 :initial-contents
1283 '(
1284 1.587171580015525194694938306936721666031w2
1285 6.185032670011643762127954396427045467506w5
1286 5.677513871931844661829755443994214173883w8
1287 1.530625323728429161131811299626419117557w11
1288 9.079594442980146270952372234833529694788w12
1289 )))
1290 (Q (make-array 5 :element-type 'double-double-float
1291 :initial-contents
1292 '(
1293 ;; 1.000000000000000000000000000000000000000w0 */
1294 1.236602014442099053716561665053645270207w4
1295 2.186249607051644894762167991800811827835w7
1296 1.092141473886177435056423606755843616331w10
1297 1.490560994263653042761789432690793026977w12
1298 2.619817175234089411411070339065679229869w13
1299 ))))
1300 (defun dd-%exp2 (x)
1301 (declare (type double-double-float x)
1302 (optimize (speed 3)))
1303 (when (>= x 1024w0)
1304 (return-from dd-%exp2
1305 (%make-double-double-float ext:double-float-positive-infinity
1306 ext:double-float-positive-infinity)))
1307 (when (<= x -1024w0)
1308 (return-from dd-%exp2 0w0))
1309 (multiple-value-bind (n x)
1310 (floor (the (double-double-float -1024w0 1024w0) x))
1311 (declare (type double-double-float x))
1312 (let* ((xx (* x x))
1313 (px (* x (poly-eval xx p))))
1314 (setf x (/ px (- (poly-eval-1 xx q) px)))
1315 (setf x (+ 1 (scale-float x 1)))
1316 (scale-float x n)))))
1317
1318 ;;; dd-%powil
1319 ;;; x^n where n is an integer
1320
1321 (defun dd-%powil (x nn)
1322 (declare (type double-double-float x)
1323 (fixnum nn)
1324 (optimize (speed 3)))
1325 (when (zerop x)
1326 (return-from dd-%powil
1327 (cond ((zerop nn)
1328 1w0)
1329 ((minusp nn)
1330 (%make-double-double-float ext:double-float-positive-infinity
1331 ext:double-float-positive-infinity))
1332 (t
1333 0w0))))
1334 (when (zerop nn)
1335 (return-from dd-%powil 1w0))
1336
1337 (let ((asign 0)
1338 (sign 0)
1339 (n nn))
1340 (declare (type (integer -1 0) asign sign)
1341 (fixnum n))
1342 (when (minusp x)
1343 (setf asign -1)
1344 (setf x (- x)))
1345 (cond ((minusp nn)
1346 (setf sign -1)
1347 (setf n (- nn)))
1348 (t
1349 (setf sign 0)
1350 (setf n nn)))
1351 ;; Overflow detection
1352
1353 ;; Calculate approximate log of answer
1354 (multiple-value-bind (s lx)
1355 (decode-float x)
1356 (declare (type double-double-float s)
1357 (type double-float-exponent lx))
1358 (let ((e (* n (1- lx))))
1359 (cond ((or (zerop e)
1360 (> e 64)
1361 (< e -64))
1362 (setf s (/ (- s 7.0710678118654752w-1)
1363 (+ s 7.0710678118654752w-1)))
1364 (setf s (* (+ (- (* 2.9142135623730950w0 s)
1365 0.5w0)
1366 lx)
1367 nn log2e)))
1368 (t
1369 (setf s (* loge2 e))))
1370 (when (> s max-log)
1371 ;; Overflow. What to do?
1372 (error "Overflow"))
1373 (when (< s min-log)
1374 (return-from dd-%powil 0w0))
1375
1376 ;; Handle denormal answer
1377 (when (< s (- 2 max-log))
1378 (setf x (/ x))
1379 (setf sign 0))
1380
1381 ;; First bit of the power
1382 (let ((y 0w0))
1383 (declare (type double-double-float y))
1384 (cond ((zerop (logand n 1))
1385 (setf y 1w0)
1386 (setf asign 0))
1387 (t
1388 (setf y x)))
1389 (let ((w x))
1390 (setf n (ash n -1))
1391 (loop while (not (zerop n))
1392 do
1393 (setf w (* w w))
1394 (unless (zerop (logand n 1))
1395 (setf y (* y w)))
1396 (setf n (ash n -1))))
1397 ;; Odd power of a negative number
1398 (unless (zerop asign)
1399 (setf y (- y)))
1400
1401 (unless (zerop sign)
1402 (setf y (/ y)))
1403 y)))))
1404
1405
1406 (defun dd-complex-pow (x y)
1407 (error "Not implemented yet"))
1408
1409 ;;; dd-real-pow
1410 ;;; x^y, for x and y real, and real result.
1411
1412 (defun dd-real-pow (x y)
1413 (declare (type double-double-float x y)
1414 (optimize (speed 3)))
1415 (let ((nflg 0)
1416 (w (floor y)))
1417 ;; nflg = 1 if x < 0 raised to integer power
1418 (when (and (= w y)
1419 (< (abs w) 32768))
1420 (return-from dd-real-pow (dd-%powil x w)))
1421
1422 (when (<= x 0)
1423 (cond ((zerop x)
1424 (if (zerop y)
1425 ;; 0^0
1426 (return-from dd-real-pow 1w0)
1427 ;; 0^y
1428 (return-from dd-real-pow 0w0)))
1429 (t
1430 (when (/= w y)
1431 ;; noninteger power of negative number
1432 (return-from dd-real-pow (dd-complex-pow x y)))
1433
1434 ;; For negative x, find out if the integer exponent is odd or even.
1435 (let ((w (scale-float y -1)))
1436 (declare (type double-double-float w))
1437 (setf w (ffloor w))
1438 (setf w (scale-float w 1))
1439 (when (/= w y)
1440 (setf nflg 1))
1441 (setf x (abs x))))))
1442 (let ((z (dd-%log2 x)))
1443 (declare (type double-double-float z))
1444 (setf z (dd-%exp2 (* y z)))
1445 (unless (zerop nflg)
1446 (setf z (- z)))
1447 z)))
1448
1449 (defun dd-%pow (x y)
1450 (declare (type double-double-float x y))
1451 (dd-real-pow x y))
1452
1453
1454 (defun dd-cssqs (z)
1455 ;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
1456 ;; result is r + i*k, where k is an integer.
1457
1458 ;; Save all FP flags
1459 (let ((x (float (realpart z) 1w0))
1460 (y (float (imagpart z) 1w0)))
1461 ;; Would this be better handled using an exception handler to
1462 ;; catch the overflow or underflow signal? For now, we turn all
1463 ;; traps off and look at the accrued exceptions to see if any
1464 ;; signal would have been raised.
1465 (with-float-traps-masked (:underflow :overflow)
1466 (let ((rho (+ (square x) (square y))))
1467 (declare (optimize (speed 3) (space 0)))
1468 (cond ((and (or (float-nan-p rho)
1469 (float-infinity-p rho))
1470 (or (float-infinity-p (abs x))
1471 (float-infinity-p (abs y))))
1472 (values (%make-double-double-float ext:double-float-positive-infinity 0w0)
1473 0))
1474 ((let ((threshold #.(/ least-positive-double-float
1475 double-float-epsilon))
1476 (traps (ldb vm::float-sticky-bits
1477 (vm:floating-point-modes))))
1478 ;; Overflow raised or (underflow raised and rho <
1479 ;; lambda/eps)
1480 (or (not (zerop (logand vm:float-overflow-trap-bit traps)))
1481 (and (not (zerop (logand vm:float-underflow-trap-bit traps)))
1482 (< rho threshold))))
1483 ;; If we're here, neither x nor y are infinity and at
1484 ;; least one is non-zero.. Thus logb returns a nice
1485 ;; integer.
1486 (let ((k (- (logb-finite (max (abs x) (abs y))))))
1487 (values (+ (square (scalb x k))
1488 (square (scalb y k)))
1489 (- k))))
1490 (t
1491 (values rho 0)))))))
1492
1493 (defun dd-complex-sqrt (z)
1494 "Principle square root of Z
1495
1496 Z may be any number, but the result is always a complex."
1497 (declare (number z))
1498 (multiple-value-bind (rho k)
1499 (dd-cssqs z)
1500 (declare (type (or (member 0w0) (double-double-float 0w0)) rho)
1501 (type fixnum k))
1502 (let ((x (float (realpart z) 1.0w0))
1503 (y (float (imagpart z) 1.0w0))
1504 (eta 0w0)
1505 (nu 0w0))
1506 (declare (type double-double-float x y eta nu))
1507
1508 (locally
1509 ;; space 0 to get maybe-inline functions inlined.
1510 (declare (optimize (speed 3) (space 0)))
1511
1512 (if (not (float-nan-p x))
1513 (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
1514
1515 (cond ((oddp k)
1516 (setf k (ash k -1)))
1517 (t
1518 (setf k (1- (ash k -1)))
1519 (setf rho (+ rho rho))))
1520
1521 (setf rho (scalb (sqrt rho) k))
1522
1523 (setf eta rho)
1524 (setf nu y)
1525
1526 (when (/= rho 0d0)
1527 (when (not (float-infinity-p (abs nu)))
1528 (setf nu (/ (/ nu rho) 2d0)))
1529 (when (< x 0d0)
1530 (setf eta (abs nu))
1531 (setf nu (float-sign y rho))))
1532 (complex eta nu)))))
1533
1534 (defun dd-complex-log-scaled (z j)
1535 "Compute log(2^j*z).
1536
1537 This is for use with J /= 0 only when |z| is huge."
1538 (declare (number z)
1539 (fixnum j))
1540 ;; The constants t0, t1, t2 should be evaluated to machine
1541 ;; precision. In addition, Kahan says the accuracy of log1p
1542 ;; influences the choices of these constants but doesn't say how to
1543 ;; choose them. We'll just assume his choices matches our
1544 ;; implementation of log1p.
1545 (let ((t0 #.(/ 1 (sqrt 2.0w0)))
1546 (t1 1.2w0)
1547 (t2 3w0)
1548 (ln2 #.(log 2w0))
1549 (x (float (realpart z) 1.0w0))
1550 (y (float (imagpart z) 1.0w0)))
1551 (multiple-value-bind (rho k)
1552 (dd-cssqs z)
1553 (declare (optimize (speed 3)))
1554 (let ((beta (max (abs x) (abs y)))
1555 (theta (min (abs x) (abs y))))
1556 (complex (if (and (zerop k)
1557 (< t0 beta)
1558 (or (<= beta t1)
1559 (< rho t2)))
1560 (/ (dd-%log1p (+ (* (- beta 1.0d0)
1561 (+ beta 1.0d0))
1562 (* theta theta)))
1563 2d0)
1564 (+ (/ (log rho) 2d0)
1565 (* (+ k j) ln2)))
1566 (atan y x))))))
1567
1568 (defun dd-complex-log (z)
1569 "Log of Z = log |Z| + i * arg Z
1570
1571 Z may be any number, but the result is always a complex."
1572 (declare (number z))
1573 (dd-complex-log-scaled z 0))
1574
1575 ;; Let us note the following "strange" behavior. atanh 1.0d0 is
1576 ;; +infinity, but the following code returns approx 176 + i*pi/4. The
1577 ;; reason for the imaginary part is caused by the fact that arg i*y is
1578 ;; never 0 since we have positive and negative zeroes.
1579
1580 (defun dd-complex-atanh (z)
1581 "Compute atanh z = (log(1+z) - log(1-z))/2"
1582 (declare (number z))
1583 (if (and (realp z) (< z -1))
1584 ;; atanh is continuous in quadrant III in this case.
1585 (dd-complex-atanh (complex z -0f0))
1586 (let* ( ;; Constants
1587 (theta (/ (sqrt most-positive-double-float) 4.0w0))
1588 (rho (/ 4.0w0 (sqrt most-positive-double-float)))
1589 (half-pi dd-pi/2)
1590 (rp (float (realpart z) 1.0w0))
1591 (beta (float-sign rp 1.0w0))
1592 (x (* beta rp))
1593 (y (* beta (- (float (imagpart z) 1.0w0))))
1594 (eta 0.0w0)
1595 (nu 0.0w0))
1596 ;; Shouldn't need this declare.
1597 (declare (double-double-float x y))
1598 (locally
1599 (declare (optimize (speed 3)))
1600 (cond ((or (> x theta)
1601 (> (abs y) theta))
1602 ;; To avoid overflow...
1603 (setf nu (float-sign y half-pi))
1604 ;; eta is real part of 1/(x + iy). This is x/(x^2+y^2),
1605 ;; which can cause overflow. Arrange this computation so
1606 ;; that it won't overflow.
1607 (setf eta (let* ((x-bigger (> x (abs y)))
1608 (r (if x-bigger (/ y x) (/ x y)))
1609 (d (+ 1.0d0 (* r r))))
1610 (if x-bigger
1611 (/ (/ x) d)
1612 (/ (/ r y) d)))))
1613 ((= x 1.0w0)
1614 ;; Should this be changed so that if y is zero, eta is set
1615 ;; to +infinity instead of approx 176? In any case
1616 ;; tanh(176) is 1.0d0 within working precision.
1617 (let ((t1 (+ 4w0 (square y)))
1618 (t2 (+ (abs y) rho)))
1619 (setf eta (dd-%log (/ (sqrt (sqrt t1))
1620 (sqrt t2))))
1621 (setf nu (* 0.5d0
1622 (float-sign y
1623 (+ half-pi (dd-%atan (* 0.5d0 t2))))))))
1624 (t
1625 (let ((t1 (+ (abs y) rho)))
1626 ;; Normal case using log1p(x) = log(1 + x)
1627 (setf eta (* 0.25d0
1628 (dd-%log1p (/ (* 4.0d0 x)
1629 (+ (square (- 1.0d0 x))
1630 (square t1))))))
1631 (setf nu (* 0.5d0
1632 (dd-%atan2 (* 2.0d0 y)
1633 (- (* (- 1.0d0 x)
1634 (+ 1.0d0 x))
1635 (square t1))))))))
1636 (format t "beta = ~A~%" beta)
1637 (format t "eta = ~A~%" eta)
1638 (format t "nu = ~A~%" nu)
1639 (complex (* beta eta)
1640 (- (* beta nu)))))))
1641
1642 (defun dd-complex-tanh (z)
1643 "Compute tanh z = sinh z / cosh z"
1644 (declare (number z))
1645 (let ((x (float (realpart z) 1.0w0))
1646 (y (float (imagpart z) 1.0w0)))
1647 (locally
1648 ;; space 0 to get maybe-inline functions inlined
1649 (declare (optimize (speed 3) (space 0)))
1650 (cond ((> (abs x)
1651 #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)
1652 ;; This is more accurate under linux.
1653 #+(or linux hpux) #.(/ (+ (%log 2.0d0)
1654 (%log most-positive-double-float)) 4d0))
1655 (complex (float-sign x)
1656 (float-sign y)))
1657 (t
1658 (let* ((tv (dd-%tan y))
1659 (beta (+ 1.0d0 (* tv tv)))
1660 (s (sinh x))
1661 (rho (sqrt (+ 1.0d0 (* s s)))))
1662 (if (float-infinity-p (abs tv))
1663 (complex (/ rho s)
1664 (/ tv))
1665 (let ((den (+ 1.0d0 (* beta s s))))
1666 (complex (/ (* beta rho s)
1667 den)
1668 (/ tv den))))))))))
1669
1670 ;; Kahan says we should only compute the parts needed. Thus, the
1671 ;; realpart's below should only compute the real part, not the whole
1672 ;; complex expression. Doing this can be important because we may get
1673 ;; spurious signals that occur in the part that we are not using.
1674 ;;
1675 ;; However, we take a pragmatic approach and just use the whole
1676 ;; expression.
1677
1678 ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
1679 ;; it's the conjugate of the square root or the square root of the
1680 ;; conjugate. This needs to be checked.
1681
1682 ;; I checked. It doesn't matter because (conjugate (sqrt z)) is the
1683 ;; same as (sqrt (conjugate z)) for all z. This follows because
1684 ;;
1685 ;; (conjugate (sqrt z)) = exp(0.5*log |z|)*exp(-0.5*j*arg z).
1686 ;;
1687 ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
1688 ;;
1689 ;; and these two expressions are equal if and only if arg conj z =
1690 ;; -arg z, which is clearly true for all z.
1691
1692 ;; NOTE: The rules of Common Lisp says that if you mix a real with a
1693 ;; complex, the real is converted to a complex before performing the
1694 ;; operation. However, Kahan says in this paper (pg 176):
1695 ;;
1696 ;; (iii) Careless handling can turn infinity or the sign of zero into
1697 ;; misinformation that subsequently disappears leaving behind
1698 ;; only a plausible but incorrect result. That is why compilers
1699 ;; must not transform z-1 into z-(1+i*0), as we have seen above,
1700 ;; nor -(-x-x^2) into (x+x^2), as we shall see below, lest a
1701 ;; subsequent logarithm or square root produce a non-zero
1702 ;; imaginary part whose sign is opposite to what was intended.
1703 ;;
1704 ;; The interesting examples are too long and complicated to reproduce
1705 ;; here. We refer the reader to his paper.
1706 ;;
1707 ;; The functions below are intended to handle the cases where a real
1708 ;; is mixed with a complex and we don't want CL complex contagion to
1709 ;; occur..
1710
1711 (declaim (inline 1+z 1-z z-1 z+1))
1712 (defun dd-1+z (z)
1713 (complex (+ 1 (realpart z)) (imagpart z)))
1714 (defun dd-1-z (z)
1715 (complex (- 1 (realpart z)) (- (imagpart z))))
1716 (defun dd-z-1 (z)
1717 (complex (- (realpart z) 1) (imagpart z)))
1718 (defun dd-z+1 (z)
1719 (complex (+ (realpart z) 1) (imagpart z)))
1720
1721 (defun dd-complex-acos (z)
1722 "Compute acos z = pi/2 - asin z
1723
1724 Z may be any number, but the result is always a complex."
1725 (declare (number z))
1726 (if (and (realp z) (> z 1))
1727 ;; acos is continuous in quadrant IV in this case.
1728 (complex-acos (complex z -0f0))
1729 (let ((sqrt-1+z (complex-sqrt (1+z z)))
1730 (sqrt-1-z (complex-sqrt (1-z z))))
1731 (with-float-traps-masked (:divide-by-zero)
1732 (complex (* 2 (atan (/ (realpart sqrt-1-z)
1733 (realpart sqrt-1+z))))
1734 (asinh (imagpart (* (conjugate sqrt-1+z)
1735 sqrt-1-z))))))))
1736
1737 (defun dd-complex-acosh (z)
1738 "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
1739
1740 Z may be any number, but the result is always a complex."
1741 (declare (number z))
1742 (let ((sqrt-z-1 (complex-sqrt (z-1 z)))
1743 (sqrt-z+1 (complex-sqrt (z+1 z))))
1744 (with-float-traps-masked (:divide-by-zero)
1745 (complex (asinh (realpart (* (conjugate sqrt-z-1)
1746 sqrt-z+1)))
1747 (* 2 (atan (/ (imagpart sqrt-z-1)
1748 (realpart sqrt-z+1))))))))
1749
1750
1751 (defun dd-complex-asin (z)
1752 "Compute asin z = asinh(i*z)/i
1753
1754 Z may be any number, but the result is always a complex."
1755 (declare (number z))
1756 (if (and (realp z) (> z 1))
1757 ;; asin is continuous in quadrant IV in this case.
1758 (dd-complex-asin (complex z -0f0))
1759 (let ((sqrt-1-z (complex-sqrt (1-z z)))
1760 (sqrt-1+z (complex-sqrt (1+z z))))
1761 (with-float-traps-masked (:divide-by-zero)
1762 ;; We get a invalid operation here when z is real and |z| > 1.
1763 (complex (atan (/ (realpart z)
1764 (realpart (* sqrt-1-z sqrt-1+z))))
1765 (asinh (imagpart (* (conjugate sqrt-1-z)
1766 sqrt-1+z))))))))
1767
1768 (defun dd-complex-asinh (z)
1769 "Compute asinh z = log(z + sqrt(1 + z*z))
1770
1771 Z may be any number, but the result is always a complex."
1772 (declare (number z))
1773 ;; asinh z = -i * asin (i*z)
1774 (let* ((iz (complex (- (imagpart z)) (realpart z)))
1775 (result (complex-asin iz)))
1776 (complex (imagpart result)
1777 (- (realpart result)))))
1778
1779 (defun dd-complex-atan (z)
1780 "Compute atan z = atanh (i*z) / i
1781
1782 Z may be any number, but the result is always a complex."
1783 (declare (number z))
1784 ;; atan z = -i * atanh (i*z)
1785 (let* ((iz (complex (- (imagpart z)) (realpart z)))
1786 (result (complex-atanh iz)))
1787 (complex (imagpart result)
1788 (- (realpart result)))))
1789
1790 (defun dd-complex-tan (z)
1791 "Compute tan z = -i * tanh(i * z)
1792
1793 Z may be any number, but the result is always a complex."
1794 (declare (number z))
1795 ;; tan z = -i * tanh(i*z)
1796 (let* ((iz (complex (- (imagpart z)) (realpart z)))
1797 (result (complex-tanh iz)))
1798 (complex (imagpart result)
1799 (- (realpart result)))))

  ViewVC Help
Powered by ViewVC 1.1.5