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

  ViewVC Help
Powered by ViewVC 1.1.5