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

  ViewVC Help
Powered by ViewVC 1.1.5