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

  ViewVC Help
Powered by ViewVC 1.1.5