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

  ViewVC Help
Powered by ViewVC 1.1.5