/[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.10 - (show annotations)
Thu May 24 19:13:17 2007 UTC (6 years, 10 months ago) by rtoy
Branch: MAIN
Changes since 1.9: +64 -26 lines
o Make DD-%ATAN2 handle signed zeroes.  This is needed to get branch
  cuts right for some special functions.
o DD-COMPLEX-ATANH needs to be careful about producing signed zeroes
  in the right spots.  Add a CAREFUL-MUL to make sure a zero of the
  correct sign is produced when needed.  (What other places need such
  care?)
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.10 2007/05/24 19:13:17 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 (neg-x (minusp (float-sign x)))
888 (neg-y (minusp (float-sign y))))
889 (declare (type (integer 0 3) code))
890 (when neg-x
891 (setf code 2))
892 (when neg-y
893 (setf code (logior code 1)))
894 ;; Handle the case where x or y is zero, taking into account the
895 ;; sign of the zero.
896 (cond ((zerop x)
897 ;; x = 0.
898 (cond ((zerop y)
899 ;; x = 0, y = 0
900 ;;
901 ;; CLHS Figure 12-15 shows what the results should
902 ;; be.
903 (if neg-x
904 (if neg-y
905 (- dd-pi)
906 dd-pi)
907 (if neg-y
908 -0w0
909 0w0)))
910 (t
911 ;; x = 0, y /= 0
912 (if neg-y
913 (- dd-pi/2)
914 dd-pi/2))))
915 ((zerop y)
916 ;; y = 0, x /= 0
917 (cond (neg-x
918 (if neg-y
919 (- dd-pi)
920 dd-pi))
921 (t
922 (if neg-y
923 -0w0
924 0w0))))
925 (t
926 (let ((w (ecase code
927 (0
928 ;; x > 0, y > 0
929 0w0)
930 (1
931 ;; x > 0, y < 0
932 -0w0)
933 (2
934 ;; x < 0, y > 0
935 dd-pi)
936 (3
937 ;; x < 0, y < 0
938 (- dd-pi)))))
939
940 (+ w (dd-%atan (/ y x))))))))
941
942
943 ;;
944 ;; Here are the fractional digits of pi, in hex.
945 ;;
946 ;; 243f6a8885 a308d31319 8a2e037073 44a4093822 299f31d008 2efa98ec4e 6c89452821 e638d01377 be5466cf34 e90c6cc0ac
947 ;; 29b7c97c50 dd3f84d5b5 b547091792 16d5d98979 fb1bd1310b a698dfb5ac 2ffd72dbd0 1adfb7b8e1 afed6a267e 96ba7c9045
948 ;;
949 ;; We want to express pi/4 in 3 parts such that the sum of the parts
950 ;; is exactly the value of pi to the higher precision. For
951 ;; double-double, we want 53 bits in each part.
952 ;;
953 ;; pi/4 in hex is
954 ;;
955 ;; C90FDAA22168C234C4C6628B80DC1CD129024E088A67CC74020BBEA63B139B22514A08798E3404DDEF9519B3CD3A431B302B0A6DF25F14374FE1356D6D51C245E485B576625E7EC6F44C42E9A637ED6B0BFF5CB6F406B7EDEE386BFB5A899FA5AE9F2411
956 ;;
957 ;;
958 #||
959 (let* ((frac #xC90FDAA22168C234C4C6628B80DC1CD129024E088A67CC74020BBEA63B139B22514A08798E3404DDEF9519B3CD3A431B302B0A6DF25F14374FE1356D6D51C245E485B576625E7EC6F44C42E9A637ED6B0BFF5CB6F406B7EDEE386BFB5A899FA5AE9F2411)
960 (len (integer-length frac)))
961 (format t "len = ~D bits~%" len)
962 (format t "pi/4 ~,,' ,106:b~%" (ldb (byte (* 4 106) (- len (* 4 106))) frac)))
963 ->
964 pi/4 11001001000011111101101010100010001000010110100011000 01000110100110001001100011001100010100010111000000011 01110000011100110100010010100100000010010011100000100 01000101001100111110011000111010000000010000010111011
965
966 ||#
967 ;; d1 = (scale-float (float #b1100100100001111110110101010001000100001011010001100001000110100110001001100011001100010100010111000000011 1w0) -106)
968 ;; d2 = (scale-float (float #b0111000001110011010001001010010000001001001110000010001000101001100111110011000111010000000010000010111011 1w0) (* 2 -106))
969 ;; d3 = (scale-float (float #b1110101001100011101100010011100110110010001001010001010010100000100001111001100011100011010000000100110111 1w0) (* 3 -106))
970 ;; d4 = (scale-float (float #b0111101111100101010001100110110011110011010011101001000011000110110011000000101011000010100110110111110010 1w0) (* 4 -106))
971
972 ;; sin(x) = x + x^3 P(x^2)
973 ;; Theoretical peak relative error = 5.6e-39
974 ;; relative peak error spread = 1.7e-9
975
976 (defconstant sincof
977 (make-array 12 :element-type 'double-double-float
978 :initial-contents
979 '(
980 6.410290407010279602425714995528976754871w-26
981 -3.868105354403065333804959405965295962871w-23
982 1.957294039628045847156851410307133941611w-20
983 -8.220635246181818130416407184286068307901w-18
984 2.811457254345322887443598804951004537784w-15
985 -7.647163731819815869711749952353081768709w-13
986 1.605904383682161459812515654720205050216w-10
987 -2.505210838544171877505034150892770940116w-8
988 2.755731922398589065255731765498970284004w-6
989 -1.984126984126984126984126984045294307281w-4
990 8.333333333333333333333333333333119885283w-3
991 -1.666666666666666666666666666666666647199w-1
992 )))
993
994 ;; cos(x) = 1 - .5 x^2 + x^2 (x^2 P(x^2))
995 ;; Theoretical peak relative error = 2.1e-37,
996 ;; relative peak error spread = 1.4e-8
997
998 (defconstant coscof
999 (make-array 11 :element-type 'double-double-float
1000 :initial-contents
1001 '(
1002 1.601961934248327059668321782499768648351w-24
1003 -8.896621117922334603659240022184527001401w-22
1004 4.110317451243694098169570731967589555498w-19
1005 -1.561920696747074515985647487260202922160w-16
1006 4.779477332386900932514186378501779328195w-14
1007 -1.147074559772972328629102981460088437917w-11
1008 2.087675698786809897637922200570559726116w-9
1009 -2.755731922398589065255365968070684102298w-7
1010 2.480158730158730158730158440896461945271w-5
1011 -1.388888888888888888888888888765724370132w-3
1012 4.166666666666666666666666666666459301466w-2
1013 )))
1014
1015 (defconstant dp1
1016 (scale-float (float #b1100100100001111110110101010001000100001011010001100001000110100110001001100011001100010100010111000000011 1w0) -106))
1017
1018 (defconstant dp2
1019 (scale-float (float #b0111000001110011010001001010010000001001001110000010001000101001100111110011000111010000000010000010111011 1w0) (* 2 -106)))
1020
1021 (defconstant dp3
1022 (scale-float (float #b1110101001100011101100010011100110110010001001010001010010100000100001111001100011100011010000000100110111 1w0) (* 3 -106)))
1023
1024 (defconstant dp4
1025 (scale-float (float #b0111101111100101010001100110110011110011010011101001000011000110110011000000101011000010100110110111110010 1w0) (* 4 -106)))
1026
1027 (defun dd-%%sin (x)
1028 (declare (type double-double-float x)
1029 (optimize (speed 3) (space 0)))
1030 (when (minusp x)
1031 (return-from dd-%%sin (- (the double-double-float (dd-%%sin (- x))))))
1032 ;; y = integer part of x/(pi/4).
1033 (let* ((y (float (floor (/ x dd-pi/4)) 1w0))
1034 (z (scale-float y -4)))
1035 (declare (type double-double-float y z))
1036 (setf z (float (floor z) 1w0)) ; integer part of y/8
1037 (setf z (- y (scale-float z 4))) ; y - 16*(y/16)
1038
1039 (let ((j (truncate z))
1040 (sign 1))
1041 (unless (zerop (logand j 1))
1042 (incf j)
1043 (incf y))
1044 (setf j (logand j 7))
1045
1046 (when (> j 3)
1047 (setf sign (- sign))
1048 (decf j 4))
1049
1050 ;; Extended precision modular arithmetic
1051 (setf z (- (- (- x (* y dp1))
1052 (* y dp2))
1053 (* y dp3)))
1054 (let ((zz (* z z)))
1055 (if (or (= j 1)
1056 (= j 2))
1057 (setf y (+ (- 1 (scale-float zz -1))
1058 (* zz zz (poly-eval zz coscof))))
1059 (setf y (+ z (* z (* zz (poly-eval zz sincof))))))
1060 (if (< sign 0)
1061 (- y)
1062 y)))))
1063
1064 (defun dd-%%cos (x)
1065 (declare (type double-double-float x)
1066 (optimize (speed 3) (space 0)))
1067 (when (minusp x)
1068 (return-from dd-%%cos (dd-%%cos (- x))))
1069 ;; y = integer part of x/(pi/4).
1070 (let* ((y (float (floor (/ x dd-pi/4)) 1w0))
1071 (z (scale-float y -4)))
1072 (declare (type double-double-float y z))
1073 (setf z (float (floor z) 1w0)) ; integer part of y/8
1074 (setf z (- y (scale-float z 4))) ; y - 16*(y/16)
1075
1076 (let ((i (truncate z))
1077 (j 0)
1078 (sign 1))
1079 (declare (type (integer 0 7) j)
1080 (type (integer -1 1) sign))
1081 (unless (zerop (logand i 1))
1082 (incf i)
1083 (incf y))
1084 (setf j (logand i 7))
1085
1086 (when (> j 3)
1087 (setf sign (- sign))
1088 (decf j 4))
1089 (when (> j 1)
1090 (setf sign (- sign)))
1091
1092 ;; Extended precision modular arithmetic. This is basically
1093 ;; computing x - y*(pi/4) accurately so that |z| < pi/4.
1094 (setf z (- (- (- x (* y dp1))
1095 (* y dp2))
1096 (* y dp3)))
1097 (let ((zz (* z z)))
1098 (if (or (= j 1)
1099 (= j 2))
1100 (setf y (+ z (* z (* zz (poly-eval zz sincof)))))
1101 (setf y (+ (- 1 (scale-float zz -1))
1102 (* zz (poly-eval zz coscof) zz))))
1103 (if (< sign 0)
1104 (- y)
1105 y)))))
1106
1107 (let ((P (make-array 6 :element-type 'double-double-float
1108 :initial-contents
1109 '(
1110 -9.889929415807650724957118893791829849557w-1
1111 1.272297782199996882828849455156962260810w3
1112 -4.249691853501233575668486667664718192660w5
1113 5.160188250214037865511600561074819366815w7
1114 -2.307030822693734879744223131873392503321w9
1115 2.883414728874239697964612246732416606301w10
1116 )))
1117 (Q (make-array 6 :element-type 'double-double-float
1118 :initial-contents
1119 '(
1120 ;; 1.000000000000000000000000000000000000000w0 */
1121 -1.317243702830553658702531997959756728291w3
1122 4.529422062441341616231663543669583527923w5
1123 -5.733709132766856723608447733926138506824w7
1124 2.758476078803232151774723646710890525496w9
1125 -4.152206921457208101480801635640958361612w10
1126 8.650244186622719093893836740197250197602w10
1127 ))))
1128 (defun dd-tancot (xx cotflag)
1129 (declare (type double-double-float xx)
1130 (optimize (speed 3) (space 0)))
1131 (let ((x 0w0)
1132 (sign 1))
1133 (declare (type double-double-float x)
1134 (type (integer -1 1) sign))
1135 (cond ((minusp xx)
1136 (setf x (- xx))
1137 (setf sign -1))
1138 (t
1139 (setf x xx)))
1140 (let* ((y (float (floor (/ x dd-pi/4)) 1w0))
1141 (z (scale-float y -4))
1142 (j 0))
1143 (declare (type double-double-float y z)
1144 (type fixnum j))
1145 (setf z (float (floor z) 1w0))
1146 (setf z (- y (scale-float z 4)))
1147
1148 (setf j (truncate z))
1149
1150 (unless (zerop (logand j 1))
1151 (incf j)
1152 (incf y))
1153
1154 (setf z (- (- (- x (* y dp1))
1155 (* y dp2))
1156 (* y dp3)))
1157 (let ((zz (* z z)))
1158 (if (> zz 1w-40)
1159 (setf y (+ z
1160 (* z (* zz (/ (poly-eval zz p)
1161 (poly-eval-1 zz q))))))
1162 (setf y z))
1163 (if (not (zerop (logand j 2)))
1164 (if cotflag
1165 (setf y (- y))
1166 (setf y (/ -1 y)))
1167 (if cotflag
1168 (setf y (/ y))))
1169 (if (< sign 0)
1170 (- y)
1171 y))))))
1172
1173 (defun dd-%%tan (x)
1174 (declare (type double-double-float x))
1175 (dd-tancot x nil))
1176
1177 (declaim (inline %kernel-rem-pi/2))
1178 (alien:def-alien-routine ("__kernel_rem_pio2" %kernel-rem-pi/2) c-call:int
1179 (x (* double-float))
1180 (y (* double-float))
1181 (e0 c-call:int)
1182 (nx c-call:int)
1183 (prec c-call:int)
1184 (ipio2 (* c-call:int)))
1185
1186 ;; This is taken from two_over_pi in fdlibm's e_rem_pio2.c. We do
1187 ;; this here so that the Sparc version doesn't have to compile in
1188 ;; e_rem_pio2, which we don't need. (But x86 and ppc do.)
1189 (defconstant two-over-pi
1190 (make-array 66 :element-type '(unsigned-byte 32)
1191 :initial-contents
1192 '(#xA2F983 #x6E4E44 #x1529FC #x2757D1 #xF534DD #xC0DB62
1193 #x95993C #x439041 #xFE5163 #xABDEBB #xC561B7 #x246E3A
1194 #x424DD2 #xE00649 #x2EEA09 #xD1921C #xFE1DEB #x1CB129
1195 #xA73EE8 #x8235F5 #x2EBB44 #x84E99C #x7026B4 #x5F7E41
1196 #x3991D6 #x398353 #x39F49C #x845F8B #xBDF928 #x3B1FF8
1197 #x97FFDE #x05980F #xEF2F11 #x8B5A0A #x6D1F6D #x367ECF
1198 #x27CB09 #xB74F46 #x3F669E #x5FEA2D #x7527BA #xC7EBE5
1199 #xF17B3D #x0739F7 #x8A5292 #xEA6BFB #x5FB11F #x8D5D08
1200 #x560330 #x46FC7B #x6BABF0 #xCFBC20 #x9AF436 #x1DA9E3
1201 #x91615E #xE61B08 #x659985 #x5F14A0 #x68408D #xFFD880
1202 #x4D7327 #x310606 #x1556CA #x73A8C9 #x60E27B #xC08C6B
1203 ))
1204 "396 (hex) digits of 2/pi")
1205
1206
1207 (let ((y (make-array 3 :element-type 'double-float))
1208 (parts (make-array 5 :element-type 'double-float)))
1209 (declare (type (simple-array double-float (3)) y)
1210 (type (simple-array double-float (5)) parts))
1211 ;; Take the double-double-float number and break it into 24-bit
1212 ;; chunks. Each chunk is an integer, which is coerced to a
1213 ;; double-float and stored in PARTS.
1214 (defun dd-expand (x)
1215 (declare (double-double-float x)
1216 (optimize (speed 3) (space 0)))
1217 (multiple-value-bind (frac exp)
1218 (decode-float x)
1219 (declare (double-double-float frac)
1220 (type (signed-byte 16) exp))
1221 (setf frac (scale-float frac 24))
1222 (decf exp 24)
1223 (dotimes (k 5)
1224 (setf (aref parts k) (coerce (ffloor frac) 'double-float))
1225 (setf frac (scale-float (- frac (aref parts k)) 24)))
1226 exp))
1227 (defun reduce-arg (x)
1228 (declare (double-double-float x)
1229 (optimize (speed 3)))
1230 (let* ((e0 (dd-expand x))
1231 (n (sys:without-gcing
1232 (%kernel-rem-pi/2 (vector-sap parts)
1233 (vector-sap y)
1234 e0
1235 (length parts)
1236 3
1237 (vector-sap two-over-pi))))
1238 (sum (+ (coerce (aref y 2) 'double-double-float)
1239 (coerce (aref y 1) 'double-double-float)
1240 (coerce (aref y 0) 'double-double-float))))
1241 (values n sum))))
1242
1243
1244 (defun dd-%sin (x)
1245 (declare (double-double-float x))
1246 (cond ((minusp x)
1247 (- (dd-%sin (- x))))
1248 ((< (abs x) (/ pi 4))
1249 (dd-%%sin x))
1250 (t
1251 ;; Argument reduction needed
1252 (multiple-value-bind (n reduced)
1253 (reduce-arg x)
1254 (case (logand n 3)
1255 (0 (dd-%%sin reduced))
1256 (1 (dd-%%cos reduced))
1257 (2 (- (dd-%%sin reduced)))
1258 (3 (- (dd-%%cos reduced))))))))
1259
1260 (defun dd-%cos (x)
1261 (declare (double-double-float x))
1262 (cond ((minusp x)
1263 (dd-%cos (- x)))
1264 ((< (abs x) (/ pi 4))
1265 (dd-%%cos x))
1266 (t
1267 ;; Argument reduction needed
1268 (multiple-value-bind (n reduced)
1269 (reduce-arg x)
1270 (case (logand n 3)
1271 (0 (dd-%%cos reduced))
1272 (1 (- (dd-%%sin reduced)))
1273 (2 (- (dd-%%cos reduced)))
1274 (3 (dd-%%sin reduced)))))))
1275
1276 (defun dd-%tan (x)
1277 (declare (double-double-float x))
1278 (cond ((minusp x)
1279 (- (dd-%tan (- x))))
1280 ((< (abs x) (/ pi 4))
1281 (dd-%%tan x))
1282 (t
1283 ;; Argument reduction needed
1284 (multiple-value-bind (n reduced)
1285 (reduce-arg x)
1286 (if (evenp n)
1287 (dd-%%tan reduced)
1288 (- (/ (dd-%%tan reduced))))))))
1289
1290 ;;; dd-%log2
1291 ;;; Base 2 logarithm.
1292
1293 (let ((P (make-array 13 :element-type 'double-double-float
1294 :initial-contents
1295 '(
1296 1.538612243596254322971797716843006400388w-6
1297 4.998469661968096229986658302195402690910w-1
1298 2.321125933898420063925789532045674660756w1
1299 4.114517881637811823002128927449878962058w2
1300 3.824952356185897735160588078446136783779w3
1301 2.128857716871515081352991964243375186031w4
1302 7.594356839258970405033155585486712125861w4
1303 1.797628303815655343403735250238293741397w5
1304 2.854829159639697837788887080758954924001w5
1305 3.007007295140399532324943111654767187848w5
1306 2.014652742082537582487669938141683759923w5
1307 7.771154681358524243729929227226708890930w4
1308 1.313572404063446165910279910527789794488w4
1309 )))
1310 (Q (make-array 12 :element-type 'double-double-float
1311 :initial-contents
1312 '(
1313 ;; 1.000000000000000000000000000000000000000w0
1314 4.839208193348159620282142911143429644326w1
1315 9.104928120962988414618126155557301584078w2
1316 9.147150349299596453976674231612674085381w3
1317 5.605842085972455027590989944010492125825w4
1318 2.248234257620569139969141618556349415120w5
1319 6.132189329546557743179177159925690841200w5
1320 1.158019977462989115839826904108208787040w6
1321 1.514882452993549494932585972882995548426w6
1322 1.347518538384329112529391120390701166528w6
1323 7.777690340007566932935753241556479363645w5
1324 2.626900195321832660448791748036714883242w5
1325 3.940717212190338497730839731583397586124w4
1326 )))
1327 ;; Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
1328 ;; where z = 2(x-1)/(x+1)
1329 ;; 1/sqrt(2) <= x < sqrt(2)
1330 ;; Theoretical peak relative error = 1.1e-35,
1331 ;; relative peak error spread 1.1e-9
1332 (R (make-array 6 :element-type 'double-double-float
1333 :initial-contents
1334 '(
1335 -8.828896441624934385266096344596648080902w-1
1336 8.057002716646055371965756206836056074715w1
1337 -2.024301798136027039250415126250455056397w3
1338 2.048819892795278657810231591630928516206w4
1339 -8.977257995689735303686582344659576526998w4
1340 1.418134209872192732479751274970992665513w5
1341 )))
1342 (S (make-array 6 :element-type 'double-double-float
1343 :initial-contents
1344 '(
1345 ;; 1.000000000000000000000000000000000000000w0 */
1346 -1.186359407982897997337150403816839480438w2
1347 3.998526750980007367835804959888064681098w3
1348 -5.748542087379434595104154610899551484314w4
1349 4.001557694070773974936904547424676279307w5
1350 -1.332535117259762928288745111081235577029w6
1351 1.701761051846631278975701529965589676574w6
1352 ))))
1353 (defun dd-%log2 (x)
1354 (declare (type double-double-float x)
1355 (optimize (speed 3) (space 0)))
1356 (multiple-value-bind (x e)
1357 (decode-float x)
1358 (declare (type double-double-float x)
1359 (type double-float-exponent e))
1360 (let ((z 0w0)
1361 (y 0w0))
1362 (declare (type double-double-float z y))
1363 (cond ((or (> e 2)
1364 (< e -2))
1365 (cond ((< x sqrt-1/2)
1366 ;; 2*(2*x-1)/(2*x+1)
1367 (decf e)
1368 (setf z (- x 0.5w0))
1369 (setf y (+ (* 0.5w0 z) 0.5w0)))
1370 (t
1371 ;; 2*(x-1)/(x+1)
1372 (setf z (- x 0.5w0))
1373 (decf z 0.5w0)
1374 (setf y (+ (* 0.5w0 z) 0.5w0))))
1375 (setf x (/ z y))
1376 (setf z (* x x))
1377 (setf y (* x (/ (* z (poly-eval z r))
1378 (poly-eval-1 z s)))))
1379 (t
1380 (cond ((< x sqrt-1/2)
1381 (decf e)
1382 (setf x (- (scale-float x 1) 1)))
1383 (t
1384 (decf x)))
1385 (setf z (* x x))
1386 (setf y (* x (/ (* z (poly-eval x p))
1387 (poly-eval-1 x q))))
1388 (decf y (scale-float z -1))))
1389 ;; Multiply log of fraction by log2(e) and base 2 exponent by 1
1390 ;;
1391 ;; This sequence of operations is critical
1392 (setf z (* y log2ea))
1393 (setf z (+ z (* x log2ea)))
1394 (setf z (+ z y))
1395 (setf z (+ z x))
1396 (setf z (+ z e))
1397 z))))
1398
1399 ;;; dd-%exp2
1400 ;;; 2^x
1401
1402 (let ((P (make-array 5 :element-type 'double-double-float
1403 :initial-contents
1404 '(
1405 1.587171580015525194694938306936721666031w2
1406 6.185032670011643762127954396427045467506w5
1407 5.677513871931844661829755443994214173883w8
1408 1.530625323728429161131811299626419117557w11
1409 9.079594442980146270952372234833529694788w12
1410 )))
1411 (Q (make-array 5 :element-type 'double-double-float
1412 :initial-contents
1413 '(
1414 ;; 1.000000000000000000000000000000000000000w0 */
1415 1.236602014442099053716561665053645270207w4
1416 2.186249607051644894762167991800811827835w7
1417 1.092141473886177435056423606755843616331w10
1418 1.490560994263653042761789432690793026977w12
1419 2.619817175234089411411070339065679229869w13
1420 ))))
1421 (defun dd-%exp2 (x)
1422 (declare (type double-double-float x)
1423 (optimize (speed 3) (space 0)))
1424 (when (>= x 1024w0)
1425 (return-from dd-%exp2
1426 (%make-double-double-float ext:double-float-positive-infinity
1427 ext:double-float-positive-infinity)))
1428 (when (<= x -1024w0)
1429 (return-from dd-%exp2 0w0))
1430 (multiple-value-bind (n x)
1431 (floor (the (double-double-float -1024w0 1024w0) x))
1432 (declare (type double-double-float x))
1433 (let* ((xx (* x x))
1434 (px (* x (poly-eval xx p))))
1435 (setf x (/ px (- (poly-eval-1 xx q) px)))
1436 (setf x (+ 1 (scale-float x 1)))
1437 (scale-float x n)))))
1438
1439 ;;; dd-%powil
1440 ;;; x^n where n is an integer
1441
1442 (defun dd-%powil (x nn)
1443 (declare (type double-double-float x)
1444 (fixnum nn)
1445 (optimize (speed 3) (space 0)))
1446 (when (zerop x)
1447 (return-from dd-%powil
1448 (cond ((zerop nn)
1449 1w0)
1450 ((minusp nn)
1451 (%make-double-double-float ext:double-float-positive-infinity
1452 ext:double-float-positive-infinity))
1453 (t
1454 0w0))))
1455 (when (zerop nn)
1456 (return-from dd-%powil 1w0))
1457
1458 (let ((asign 0)
1459 (sign 0)
1460 (n nn))
1461 (declare (type (integer -1 0) asign sign)
1462 (fixnum n))
1463 (when (minusp x)
1464 (setf asign -1)
1465 (setf x (- x)))
1466 (cond ((minusp nn)
1467 (setf sign -1)
1468 (setf n (- nn)))
1469 (t
1470 (setf sign 0)
1471 (setf n nn)))
1472 ;; Overflow detection
1473
1474 ;; Calculate approximate log of answer
1475 (multiple-value-bind (s lx)
1476 (decode-float x)
1477 (declare (type double-double-float s)
1478 (type double-float-exponent lx))
1479 (let ((e (* n (1- lx))))
1480 (cond ((or (zerop e)
1481 (> e 64)
1482 (< e -64))
1483 (setf s (/ (- s 7.0710678118654752w-1)
1484 (+ s 7.0710678118654752w-1)))
1485 (setf s (* (+ (- (* 2.9142135623730950w0 s)
1486 0.5w0)
1487 lx)
1488 nn log2e)))
1489 (t
1490 (setf s (* loge2 e))))
1491 (when (> s max-log)
1492 ;; Overflow. What to do?
1493 (error "Overflow"))
1494 (when (< s min-log)
1495 (return-from dd-%powil 0w0))
1496
1497 ;; Handle denormal answer
1498 (when (< s (- 2 max-log))
1499 (setf x (/ x))
1500 (setf sign 0))
1501
1502 ;; First bit of the power
1503 (let ((y 0w0))
1504 (declare (type double-double-float y))
1505 (cond ((zerop (logand n 1))
1506 (setf y 1w0)
1507 (setf asign 0))
1508 (t
1509 (setf y x)))
1510 (let ((w x))
1511 (declare (type double-double-float w))
1512 (setf n (ash n -1))
1513 (loop while (not (zerop n))
1514 do
1515 (setf w (* w w))
1516 (unless (zerop (logand n 1))
1517 (setf y (* y w)))
1518 (setf n (ash n -1))))
1519 ;; Odd power of a negative number
1520 (unless (zerop asign)
1521 (setf y (- y)))
1522
1523 (unless (zerop sign)
1524 (setf y (/ y)))
1525 y)))))
1526
1527
1528 ;;; dd-real-pow
1529 ;;; x^y, for x and y real, and real result.
1530
1531 (defun dd-real-pow (x y)
1532 (declare (type double-double-float x y)
1533 (optimize (speed 3) (space 0)))
1534 (let ((nflg 0)
1535 (w (floor y)))
1536 ;; nflg = 1 if x < 0 raised to integer power
1537 (when (and (= w y)
1538 (< (abs w) 32768))
1539 (return-from dd-real-pow (dd-%powil x w)))
1540
1541 (when (<= x 0)
1542 (cond ((zerop x)
1543 (if (zerop y)
1544 ;; 0^0
1545 (return-from dd-real-pow 1w0)
1546 ;; 0^y
1547 (return-from dd-real-pow 0w0)))
1548 (t
1549 (when (/= w y)
1550 ;; noninteger power of negative number
1551 (let ((p (the double-double-float (dd-real-pow (abs x) y)))
1552 (y*pi (* y dd-pi)))
1553 (return-from dd-real-pow (complex (* p (dd-%cos y*pi))
1554 (* p (dd-%sin y*pi))))))
1555
1556 ;; For negative x, find out if the integer exponent is odd or even.
1557 (let ((w (scale-float y -1)))
1558 (declare (type double-double-float w))
1559 (setf w (ffloor w))
1560 (setf w (scale-float w 1))
1561 (when (/= w y)
1562 (setf nflg 1))
1563 (setf x (abs x))))))
1564 (let ((z (dd-%log2 x)))
1565 (declare (type double-double-float z))
1566 (setf z (dd-%exp2 (* y z)))
1567 (unless (zerop nflg)
1568 (setf z (- z)))
1569 z)))
1570
1571 (defun dd-%pow (x y)
1572 (declare (type double-double-float x y))
1573 (dd-real-pow x y))
1574
1575
1576 ;; These are essentially the same as in irrat.lisp, but very slightly
1577 ;; modified to operate on double-double-floats.
1578 (defun dd-cssqs (z)
1579 ;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
1580 ;; result is r + i*k, where k is an integer.
1581
1582 ;; Save all FP flags
1583 (let ((x (float (realpart z) 1w0))
1584 (y (float (imagpart z) 1w0)))
1585 ;; Would this be better handled using an exception handler to
1586 ;; catch the overflow or underflow signal? For now, we turn all
1587 ;; traps off and look at the accrued exceptions to see if any
1588 ;; signal would have been raised.
1589 ;;
1590 ;; Actually, for double-double-floats, we should probably
1591 ;; explicitly check for overflow instead of disabling the traps.
1592 ;; Why? Because instead of overflow, double-double signals
1593 ;; invalid operation.
1594 (with-float-traps-masked (:underflow :overflow :invalid)
1595 (let ((rho (+ (square x) (square y))))
1596 (declare (optimize (speed 3) (space 0)))
1597 (cond ((and (or (float-nan-p rho)
1598 (float-infinity-p rho))
1599 (or (float-infinity-p (abs x))
1600 (float-infinity-p (abs y))))
1601 (values (%make-double-double-float ext:double-float-positive-infinity 0d0)
1602 0))
1603 ((let ((threshold #.(/ least-positive-double-float
1604 double-float-epsilon))
1605 (traps (ldb vm::float-sticky-bits
1606 (vm:floating-point-modes))))
1607 ;; Overflow raised or (underflow raised and rho <
1608 ;; lambda/eps)
1609 (or (not (zerop (logand vm:float-overflow-trap-bit traps)))
1610 (and (not (zerop (logand vm:float-underflow-trap-bit traps)))
1611 (< rho threshold))))
1612 ;; If we're here, neither x nor y are infinity and at
1613 ;; least one is non-zero.. Thus logb returns a nice
1614 ;; integer.
1615 (let ((k (- (logb-finite (max (abs x) (abs y))))))
1616 (values (+ (square (scalb x k))
1617 (square (scalb y k)))
1618 (- k))))
1619 (t
1620 (values rho 0)))))))
1621
1622 (defun dd-complex-sqrt (z)
1623 "Principle square root of Z
1624
1625 Z may be any number, but the result is always a complex."
1626 (declare (number z))
1627 (multiple-value-bind (rho k)
1628 (dd-cssqs z)
1629 (declare (type (or (member 0w0) (double-double-float 0w0)) rho)
1630 (type fixnum k))
1631 (let ((x (float (realpart z) 1.0w0))
1632 (y (float (imagpart z) 1.0w0))
1633 (eta 0w0)
1634 (nu 0w0))
1635 (declare (type double-double-float x y eta nu))
1636
1637 (locally
1638 ;; space 0 to get maybe-inline functions inlined.
1639 (declare (optimize (speed 3) (space 0)))
1640
1641 (if (not (float-nan-p x))
1642 (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
1643
1644 (cond ((oddp k)
1645 (setf k (ash k -1)))
1646 (t
1647 (setf k (1- (ash k -1)))
1648 (setf rho (+ rho rho))))
1649
1650 (setf rho (scalb (sqrt rho) k))
1651
1652 (setf eta rho)
1653 (setf nu y)
1654
1655 (when (/= rho 0d0)
1656 (when (not (float-infinity-p (abs nu)))
1657 (setf nu (/ (/ nu rho) 2d0)))
1658 (when (< x 0d0)
1659 (setf eta (abs nu))
1660 (setf nu (float-sign y rho))))
1661 (complex eta nu)))))
1662
1663 (defun dd-complex-log-scaled (z j)
1664 "Compute log(2^j*z).
1665
1666 This is for use with J /= 0 only when |z| is huge."
1667 (declare (number z)
1668 (fixnum j))
1669 ;; The constants t0, t1, t2 should be evaluated to machine
1670 ;; precision. In addition, Kahan says the accuracy of log1p
1671 ;; influences the choices of these constants but doesn't say how to
1672 ;; choose them. We'll just assume his choices matches our
1673 ;; implementation of log1p.
1674 (let ((t0 #.(/ 1 (sqrt 2.0w0)))
1675 (t1 1.2w0)
1676 (t2 3w0)
1677 (ln2 #.(log 2w0))
1678 (x (float (realpart z) 1.0w0))
1679 (y (float (imagpart z) 1.0w0)))
1680 (multiple-value-bind (rho k)
1681 (dd-cssqs z)
1682 (declare (optimize (speed 3)))
1683 (let ((beta (max (abs x) (abs y)))
1684 (theta (min (abs x) (abs y))))
1685 (complex (if (and (zerop k)
1686 (< t0 beta)
1687 (or (<= beta t1)
1688 (< rho t2)))
1689 (/ (dd-%log1p (+ (* (- beta 1.0d0)
1690 (+ beta 1.0d0))
1691 (* theta theta)))
1692 2d0)
1693 (+ (/ (log rho) 2d0)
1694 (* (+ k j) ln2)))
1695 (atan y x))))))
1696
1697 (defun dd-complex-log (z)
1698 "Log of Z = log |Z| + i * arg Z
1699
1700 Z may be any number, but the result is always a complex."
1701 (declare (number z))
1702 (dd-complex-log-scaled z 0))
1703
1704 ;; Let us note the following "strange" behavior. atanh 1.0d0 is
1705 ;; +infinity, but the following code returns approx 176 + i*pi/4. The
1706 ;; reason for the imaginary part is caused by the fact that arg i*y is
1707 ;; never 0 since we have positive and negative zeroes.
1708
1709 (defun dd-complex-atanh (z)
1710 "Compute atanh z = (log(1+z) - log(1-z))/2"
1711 (declare (number z))
1712 (cond ((realp z)
1713 ;; Look at the definition:
1714 ;;
1715 ;; atanh(z) = 1/2*(log(1+z)-log(1-z))
1716 ;;
1717 (cond ((> z 1)
1718 ;; Let x = z, x > 1. Then
1719 ;;
1720 ;; atanh(x) = 1/2*(log(1+x)-log(1-x))
1721 ;;
1722 ;; Only the term log(1-x) requires care since the
1723 ;; other term is purely real. The CLHS says atanh for
1724 ;; x > 1 is continuous with quadrant I. Assume x is
1725 ;; really x0 + i*eps, where eps > 0. Then
1726 ;;
1727 ;; log(1-x) = log(x0-1) - i*pi/2
1728 ;;
1729 ;; because arg(1-x) = arg(1-x0-i*eps) = -pi
1730 ;;
1731 ;; Thus
1732 ;;
1733 ;; atanh(x) = 1/2*log((x+1)/(x-1)) + i*pi/2
1734 ;; = 1/2*log(1+2/(x-1)) + i*pi/2
1735 (complex (* 0.5w0 (dd-%log1p (/ 2 (- z 1))))
1736 dd-pi/2))
1737 (t
1738 ;; As above, but z = -x, x > 1. Then
1739 ;;
1740 ;; atanh(z) = 1/2*(log(1-x)-log(1+x))
1741 ;;
1742 ;; And log(1-x) is the interesting term. The CLHS
1743 ;; says in this case atanh is continuous with quadrant
1744 ;; III. Let x = x0-i*eps. Then
1745 ;;
1746 ;; log(1-x) = log(x0-1) + i*pi/2
1747 ;;
1748 ;; because arg(1-x) = arg(1-x0-i*eps) = pi. Thus
1749 ;;
1750 ;; atanh(z) = 1/2*log((x-1)/(x+1)) - i*pi/2
1751 ;; = -1/2*log((x+1)/(x-1)) - i*pi/2
1752 (complex (* -0.5w0 (dd-%log1p (/ 2 (- (abs z) 1))))
1753 (- dd-pi/2)))))
1754 (t
1755 (flet ((careful-mul (a b)
1756 ;; Carefully multiply a and b, taking care to handle
1757 ;; signed zeroes. Only need to handle the case of b
1758 ;; being zero.
1759 (if (zerop b)
1760 (if (minusp (* (float-sign a) (float-sign b)))
1761 -0w0
1762 0w0)
1763 (* a b))))
1764 (let* ( ;; Constants
1765 (theta (/ (sqrt most-positive-double-float) 4.0w0))
1766 (rho (/ 4.0w0 (sqrt most-positive-double-float)))
1767 (half-pi dd-pi/2)
1768 (rp (float (realpart z) 1.0w0))
1769 (beta (float-sign rp 1.0w0))
1770 (x (* beta rp))
1771 (y (careful-mul beta (- (float (imagpart z) 1.0w0))))
1772 (eta 0.0w0)
1773 (nu 0.0w0))
1774 ;; Shouldn't need this declare.
1775 (declare (double-double-float x y))
1776 (locally
1777 (declare (optimize (speed 3)))
1778 (cond ((or (> x theta)
1779 (> (abs y) theta))
1780 ;; To avoid overflow...
1781 (setf nu (float-sign y half-pi))
1782 ;; eta is real part of 1/(x + iy). This is x/(x^2+y^2),
1783 ;; which can cause overflow. Arrange this computation so
1784 ;; that it won't overflow.
1785 (setf eta (let* ((x-bigger (> x (abs y)))
1786 (r (if x-bigger (/ y x) (/ x y)))
1787 (d (+ 1.0d0 (* r r))))
1788 (if x-bigger
1789 (/ (/ x) d)
1790 (/ (/ r y) d)))))
1791 ((= x 1.0w0)
1792 ;; Should this be changed so that if y is zero, eta is set
1793 ;; to +infinity instead of approx 176? In any case
1794 ;; tanh(176) is 1.0d0 within working precision.
1795 (let ((t1 (+ 4w0 (square y)))
1796 (t2 (+ (abs y) rho)))
1797 (setf eta (dd-%log (/ (sqrt (sqrt t1))
1798 (sqrt t2))))
1799 (setf nu (* 0.5d0
1800 (float-sign y
1801 (+ half-pi (dd-%atan (* 0.5d0 t2))))))))
1802 (t
1803 (let ((t1 (+ (abs y) rho)))
1804 ;; Normal case using log1p(x) = log(1 + x)
1805 (setf eta (* 0.25d0
1806 (dd-%log1p (/ (* 4.0d0 x)
1807 (+ (square (- 1.0d0 x))
1808 (square t1))))))
1809 (setf nu (* 0.5d0
1810 (dd-%atan2 (careful-mul 2.0d0 y)
1811 (- (* (- 1.0d0 x)
1812 (+ 1.0d0 x))
1813 (square t1))))))))
1814 (complex (* beta eta)
1815 (- (* beta nu)))))))))
1816
1817 (defun dd-complex-tanh (z)
1818 "Compute tanh z = sinh z / cosh z"
1819 (declare (number z))
1820 (let ((x (float (realpart z) 1.0w0))
1821 (y (float (imagpart z) 1.0w0)))
1822 (locally
1823 ;; space 0 to get maybe-inline functions inlined
1824 (declare (optimize (speed 3) (space 0)))
1825 (cond ((> (abs x)
1826 #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)
1827 ;; This is more accurate under linux.
1828 #+(or linux hpux) #.(/ (+ (%log 2.0d0)
1829 (%log most-positive-double-float)) 4d0))
1830 (complex (float-sign x)
1831 (float-sign y)))
1832 (t
1833 (let* ((tv (dd-%tan y))
1834 (beta (+ 1.0d0 (* tv tv)))
1835 (s (sinh x))
1836 (rho (sqrt (+ 1.0d0 (* s s)))))
1837 (if (float-infinity-p (abs tv))
1838 (complex (/ rho s)
1839 (/ tv))
1840 (let ((den (+ 1.0d0 (* beta s s))))
1841 (complex (/ (* beta rho s)
1842 den)
1843 (/ tv den))))))))))
1844
1845 ;; Kahan says we should only compute the parts needed. Thus, the
1846 ;; realpart's below should only compute the real part, not the whole
1847 ;; complex expression. Doing this can be important because we may get
1848 ;; spurious signals that occur in the part that we are not using.
1849 ;;
1850 ;; However, we take a pragmatic approach and just use the whole
1851 ;; expression.
1852
1853 ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
1854 ;; it's the conjugate of the square root or the square root of the
1855 ;; conjugate. This needs to be checked.
1856
1857 ;; I checked. It doesn't matter because (conjugate (sqrt z)) is the
1858 ;; same as (sqrt (conjugate z)) for all z. This follows because
1859 ;;
1860 ;; (conjugate (sqrt z)) = exp(0.5*log |z|)*exp(-0.5*j*arg z).
1861 ;;
1862 ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
1863 ;;
1864 ;; and these two expressions are equal if and only if arg conj z =
1865 ;; -arg z, which is clearly true for all z.
1866
1867 ;; NOTE: The rules of Common Lisp says that if you mix a real with a
1868 ;; complex, the real is converted to a complex before performing the
1869 ;; operation. However, Kahan says in this paper (pg 176):
1870 ;;
1871 ;; (iii) Careless handling can turn infinity or the sign of zero into
1872 ;; misinformation that subsequently disappears leaving behind
1873 ;; only a plausible but incorrect result. That is why compilers
1874 ;; must not transform z-1 into z-(1+i*0), as we have seen above,
1875 ;; nor -(-x-x^2) into (x+x^2), as we shall see below, lest a
1876 ;; subsequent logarithm or square root produce a non-zero
1877 ;; imaginary part whose sign is opposite to what was intended.
1878 ;;
1879 ;; The interesting examples are too long and complicated to reproduce
1880 ;; here. We refer the reader to his paper.
1881 ;;
1882 ;; The functions below are intended to handle the cases where a real
1883 ;; is mixed with a complex and we don't want CL complex contagion to
1884 ;; occur..
1885
1886 (declaim (inline 1+z 1-z z-1 z+1))
1887 (defun dd-1+z (z)
1888 (complex (+ 1 (realpart z)) (imagpart z)))
1889 (defun dd-1-z (z)
1890 (complex (- 1 (realpart z)) (- (imagpart z))))
1891 (defun dd-z-1 (z)
1892 (complex (- (realpart z) 1) (imagpart z)))
1893 (defun dd-z+1 (z)
1894 (complex (+ (realpart z) 1) (imagpart z)))
1895
1896 (defun dd-complex-acos (z)
1897 "Compute acos z = pi/2 - asin z
1898
1899 Z may be any number, but the result is always a complex."
1900 (declare (number z))
1901 (if (and (realp z) (> z 1))
1902 ;; acos is continuous in quadrant IV in this case.
1903 (complex-acos (complex z -0f0))
1904 (let ((sqrt-1+z (complex-sqrt (1+z z)))
1905 (sqrt-1-z (complex-sqrt (1-z z))))
1906 (cond ((zerop (realpart sqrt-1+z))
1907 ;; Same as below, but we compute atan ourselves (because we
1908 ;; have atan +/- infinity).
1909 (complex
1910 (if (minusp (float-sign (* (realpart sqrt-1-z)
1911 (realpart sqrt-1+z))))
1912 (- dd-pi)
1913 dd-pi)
1914 (asinh (imagpart (* (conjugate sqrt-1+z)
1915 sqrt-1-z)))))
1916 (t
1917 (complex (* 2 (atan (/ (realpart sqrt-1-z)
1918 (realpart sqrt-1+z))))
1919 (asinh (imagpart (* (conjugate sqrt-1+z)
1920 sqrt-1-z)))))))))
1921
1922 (defun dd-complex-acosh (z)
1923 "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
1924
1925 Z may be any number, but the result is always a complex."
1926 (declare (number z))
1927 (let* ((sqrt-z-1 (complex-sqrt (z-1 z)))
1928 (sqrt-z+1 (complex-sqrt (z+1 z))))
1929 ;; We need to handle the case where real part of sqrt-z+1 is zero,
1930 ;; because division by zero with double-double-floats doesn't
1931 ;; produce infinity.
1932 (cond ((zerop (realpart sqrt-z+1))
1933 ;; Same as below, but we compute atan ourselves (because we
1934 ;; have atan +/- infinity).
1935 (complex (asinh (realpart (* (conjugate sqrt-z-1)
1936 sqrt-z+1)))
1937 (if (minusp (float-sign (* (imagpart sqrt-z-1)
1938 (realpart sqrt-z+1))))
1939 (- dd-pi)
1940 dd-pi)))
1941 (t
1942 (complex (asinh (realpart (* (conjugate sqrt-z-1)
1943 sqrt-z+1)))
1944 (* 2 (atan (/ (imagpart sqrt-z-1)
1945 (realpart sqrt-z+1)))))))))
1946
1947
1948 (defun dd-complex-asin (z)
1949 "Compute asin z = asinh(i*z)/i
1950
1951 Z may be any number, but the result is always a complex."
1952 (declare (number z))
1953 (cond ((realp z)
1954 ;; Z is real, and |Z| > 1.
1955 ;;
1956 ;; Note that
1957 ;;
1958 ;; sin(pi/2+i*y) = cosh(y)
1959 ;;
1960 ;; and
1961 ;;
1962 ;; sin(-pi/2+i*y) = -cosh(y).
1963 ;;
1964 ;; Since cosh(y) >= 1 for all real y, we see that
1965 ;;
1966 ;; asin(z) = pi/2*sign(z) + i*acosh(abs(z))
1967 ;;
1968 (complex (if (minusp z)
1969 (- dd-pi/2)
1970 dd-pi/2)
1971 (dd-%acosh (abs z))))
1972 (t
1973 (let* ((sqrt-1-z (complex-sqrt (1-z z)))
1974 (sqrt-1+z (complex-sqrt (1+z z)))
1975 (den (realpart (* sqrt-1-z sqrt-1+z))))
1976 (cond ((zerop den)
1977 ;; Like below but we handle atan part ourselves.
1978 (complex (if (minusp (float-sign den))
1979 (- dd-pi/2)
1980 dd-pi/2)
1981 (asinh (imagpart (* (conjugate sqrt-1-z)
1982 sqrt-1+z)))))
1983 (t
1984 (with-float-traps-masked (:divide-by-zero)
1985 ;; We get a invalid operation here when z is real and |z| > 1.
1986 (complex (atan (/ (realpart z)
1987 den))
1988 (asinh (imagpart (* (conjugate sqrt-1-z)
1989 sqrt-1+z)))))))))))
1990
1991 (defun dd-complex-asinh (z)
1992 "Compute asinh z = log(z + sqrt(1 + z*z))
1993
1994 Z may be any number, but the result is always a complex."
1995 (declare (number z))
1996 ;; asinh z = -i * asin (i*z)
1997 (let* ((iz (complex (- (imagpart z)) (realpart z)))
1998 (result (complex-asin iz)))
1999 (complex (imagpart result)
2000 (- (realpart result)))))
2001
2002 (defun dd-complex-atan (z)
2003 "Compute atan z = atanh (i*z) / i
2004
2005 Z may be any number, but the result is always a complex."
2006 (declare (number z))
2007 ;; atan z = -i * atanh (i*z)
2008 (let* ((iz (complex (- (imagpart z)) (realpart z)))
2009 (result (complex-atanh iz)))
2010 (complex (imagpart result)
2011 (- (realpart result)))))
2012
2013 (defun dd-complex-tan (z)
2014 "Compute tan z = -i * tanh(i * z)
2015
2016 Z may be any number, but the result is always a complex."
2017 (declare (number z))
2018 ;; tan z = -i * tanh(i*z)
2019 (let* ((iz (complex (- (imagpart z)) (realpart z)))
2020 (result (complex-tanh iz)))
2021 (complex (imagpart result)
2022 (- (realpart result)))))

  ViewVC Help
Powered by ViewVC 1.1.5