/[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 - (hide annotations)
Tue Apr 20 17:57:44 2010 UTC (4 years 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 rtoy 1.2 ;;; -*- 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 rtoy 1.20 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/irrat-dd.lisp,v 1.20 2010/04/20 17:57:44 rtoy Rel $")
9 rtoy 1.2 ;;;
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 rtoy 1.19 (intl:textdomain "cmucl")
23    
24 rtoy 1.2 ;;;; Random constants, utility functions.
25    
26     (defconstant max-log
27     7.0978271289338399678773454114191w2
28 rtoy 1.19 _N"log(most-positive-double-double-float)")
29 rtoy 1.2
30     (defconstant min-log
31     -7.4444007192138126231410729844608w2
32 rtoy 1.19 _N"log(least-positive-double-double-float")
33 rtoy 1.2
34    
35     (defconstant loge2
36     0.6931471805599453094172321214581765680755001w0
37 rtoy 1.19 _N"log(2)")
38 rtoy 1.2
39     (defconstant log2e
40     1.442695040888963407359924681001892137426646w0
41 rtoy 1.19 _N"Log base 2 of e")
42 rtoy 1.2
43     (defconstant log2ea
44     4.4269504088896340735992468100189213742664595w-1
45 rtoy 1.19 _N"log2(e)-1")
46 rtoy 1.2
47     (defconstant dd-pi
48     3.141592653589793238462643383279502884197169w0
49 rtoy 1.19 _N"Pi")
50 rtoy 1.2
51     (defconstant dd-pi/2
52     1.570796326794896619231321691639751442098585w0
53 rtoy 1.19 _N"Pi/2")
54 rtoy 1.2
55     (defconstant dd-pi/4
56     0.7853981633974483096156608458198757210492923w0
57 rtoy 1.19 _N"Pi/4")
58 rtoy 1.2
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 rtoy 1.19 _N"Sqrt(1/2)")
70 rtoy 1.2
71     ;; Evaluate polynomial
72 rtoy 1.18 (declaim (maybe-inline poly-eval poly-eval-1))
73 rtoy 1.2 (defun poly-eval (x coef)
74     (declare (type double-double-float x)
75 rtoy 1.3 (type (simple-array double-double-float (*)) coef)
76 rtoy 1.16 (optimize (speed 3) (space 0) (inhibit-warnings 3)))
77 rtoy 1.2 ;; 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 rtoy 1.3 (type (simple-array double-double-float (*)) coef)
91 rtoy 1.16 (optimize (speed 3) (space 0) (inhibit-warnings 3)))
92 rtoy 1.2 ;; 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 rtoy 1.7 -8.802340681794263968892934703309274564037w1
135     3.697714952261803935521187272204485251835w3
136     -9.615511549171441430850103489315371768998w4
137     1.682912729190313538934190635536631941751w6
138     -2.019684072836541751428967854947019415698w7
139     1.615869009634292424463780387327037251069w8
140 rtoy 1.2 -7.848989743695296475743081255027098295771w8
141 rtoy 1.7 1.766112549341972444333352727998584753865w9)))
142 rtoy 1.2 ;; 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 rtoy 1.19 _N"exp(x) - 1"
150 rtoy 1.2 (declare (type double-double-float x)
151 rtoy 1.17 (optimize (speed 3) (space 0)
152     (inhibit-warnings 3)))
153 rtoy 1.2 ;; 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 rtoy 1.7 (k (floor (+ 0.5w0 (/ (the (double-double-float * 710w0) x) xx))))
171 rtoy 1.2 (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 rtoy 1.7 (setf px (* x (poly-eval x p)))
180 rtoy 1.2 (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 rtoy 1.17 (optimize (speed 3) (space 0)
250     (inhibit-warnings 3)))
251 rtoy 1.2 (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 rtoy 1.17 (optimize (speed 3) (space 0)
363     (inhibit-warnings 3)))
364 rtoy 1.2 ;; 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 rtoy 1.3 (declare (type double-double-float xm1)
487     (optimize (space 0)))
488 rtoy 1.2 (let ((x (+ xm1 1))
489     (z 0w0)
490     (y 0w0))
491     (declare (type double-double-float x y z)
492 rtoy 1.17 (optimize (speed 3)
493     (inhibit-warnings 3)))
494 rtoy 1.2 (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 rtoy 1.17 (optimize (speed 3) (space 0)
556     (inhibit-warnings 3)))
557 rtoy 1.2 (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 rtoy 1.3 (declare (type double-double-float x)
591 rtoy 1.17 (optimize (speed 3) (space 0)
592     (inhibit-warnings 3)))
593 rtoy 1.2 (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 rtoy 1.17 (optimize (speed 3) (space 0)
615     (inhibit-warnings 3)))
616 rtoy 1.2 (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 rtoy 1.17 (optimize (speed 3) (space 0)
666     (inhibit-warnings 3)))
667 rtoy 1.2 (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 rtoy 1.17 (optimize (speed 3) (space 0)
707     (inhibit-warnings 3)))
708 rtoy 1.2 (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 rtoy 1.16 (declare (type (double-double-float 1w0) x)
749 rtoy 1.17 (optimize (speed 3) (space 0)
750     (inhibit-warnings 3)))
751 rtoy 1.2 (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 rtoy 1.16
763 rtoy 1.2 (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 rtoy 1.16 (declare (type (double-double-float -1w0 1w0) x)
792 rtoy 1.17 (optimize (speed 3) (space 0)
793     (inhibit-warnings 3)))
794 rtoy 1.2 (cond ((minusp x)
795     (- (the double-double-float (dd-%asin (- x)))))
796     #+nil
797     ((< x 1w-8)
798     x)
799     (t
800 rtoy 1.16 (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 rtoy 1.2 (declare (type double-double-float z zz))
808 rtoy 1.16
809 rtoy 1.2 (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 rtoy 1.16 (declare (type (double-double-float -1w0 1w0) x)
819 rtoy 1.17 (optimize (speed 3) (space 0)
820     (inhibit-warnings 3)))
821 rtoy 1.2 (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 rtoy 1.17 (optimize (speed 3) (space 0)
862     (inhibit-warnings 3)))
863 rtoy 1.2 (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 rtoy 1.17 (optimize (speed 3) (space 0)
889     (inhibit-warnings 3)))
890 rtoy 1.2 (let ((code 0)
891 rtoy 1.10 (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 rtoy 1.2 (setf code 2))
896 rtoy 1.10 (when neg-y
897 rtoy 1.2 (setf code (logior code 1)))
898 rtoy 1.10 ;; 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 rtoy 1.2
944 rtoy 1.10 (+ w (dd-%atan (/ y x))))))))
945 rtoy 1.2
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 rtoy 1.5 (defun dd-%%sin (x)
1032 rtoy 1.3 (declare (type double-double-float x)
1033 rtoy 1.17 (optimize (speed 2) (space 0)
1034     (inhibit-warnings 3)))
1035 rtoy 1.2 (when (minusp x)
1036 rtoy 1.5 (return-from dd-%%sin (- (the double-double-float (dd-%%sin (- x))))))
1037 rtoy 1.2 ;; 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 rtoy 1.16 (declare (type (integer -1 1) sign))
1047 rtoy 1.2 (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 rtoy 1.5 (defun dd-%%cos (x)
1071 rtoy 1.2 (declare (type double-double-float x)
1072 rtoy 1.17 (optimize (speed 2) (space 0)
1073     (inhibit-warnings 3)))
1074 rtoy 1.2 (when (minusp x)
1075 rtoy 1.5 (return-from dd-%%cos (dd-%%cos (- x))))
1076 rtoy 1.2 ;; 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 rtoy 1.16 (optimize (speed 2) (space 0)))
1138 rtoy 1.2 (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 rtoy 1.5 (defun dd-%%tan (x)
1181 rtoy 1.2 (declare (type double-double-float x))
1182     (dd-tancot x nil))
1183    
1184 rtoy 1.5 (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 rtoy 1.19 _N"396 (hex) digits of 2/pi")
1212 rtoy 1.5
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 rtoy 1.17 (optimize (speed 3) (space 0)
1224     (inhibit-warnings 3)))
1225 rtoy 1.5 (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 rtoy 1.17 (optimize (speed 3)
1238     (inhibit-warnings 3)))
1239 rtoy 1.5 (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 rtoy 1.16 (declaim (ftype (function (double-double-float) double-double-float)
1254     dd-%sin))
1255 rtoy 1.5 (defun dd-%sin (x)
1256     (declare (double-double-float x))
1257 rtoy 1.13 (cond ((minusp (float-sign x))
1258 rtoy 1.6 (- (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 rtoy 1.16 (ecase (logand n 3)
1266 rtoy 1.6 (0 (dd-%%sin reduced))
1267     (1 (dd-%%cos reduced))
1268     (2 (- (dd-%%sin reduced)))
1269     (3 (- (dd-%%cos reduced))))))))
1270 rtoy 1.5
1271 rtoy 1.16 (declaim (ftype (function (double-double-float) double-double-float)
1272     dd-%cos))
1273 rtoy 1.5 (defun dd-%cos (x)
1274     (declare (double-double-float x))
1275 rtoy 1.6 (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 rtoy 1.16 (ecase (logand n 3)
1284 rtoy 1.6 (0 (dd-%%cos reduced))
1285     (1 (- (dd-%%sin reduced)))
1286     (2 (- (dd-%%cos reduced)))
1287     (3 (dd-%%sin reduced)))))))
1288 rtoy 1.5
1289 rtoy 1.16 (declaim (ftype (function (double-double-float) double-double-float)
1290     dd-%tan))
1291 rtoy 1.5 (defun dd-%tan (x)
1292     (declare (double-double-float x))
1293 rtoy 1.13 (cond ((minusp (float-sign x))
1294 rtoy 1.6 (- (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 rtoy 1.5
1305 rtoy 1.2 ;;; 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 rtoy 1.17 (optimize (speed 3) (space 0)
1371     (inhibit-warnings 3)))
1372 rtoy 1.2 (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 rtoy 1.17 (optimize (speed 3) (space 0)
1440     (inhibit-warnings 3)))
1441 rtoy 1.2 (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 rtoy 1.17 (optimize (speed 3) (space 0)
1463     (inhibit-warnings 3)))
1464 rtoy 1.2 (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 rtoy 1.20 (error (intl:gettext "Overflow")))
1512 rtoy 1.2 (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 rtoy 1.17 (optimize (speed 3) (space 0)
1552     (inhibit-warnings 3)))
1553 rtoy 1.2 (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 rtoy 1.5 ;;
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 rtoy 1.2 (let ((rho (+ (square x) (square y))))
1615 rtoy 1.17 (declare (optimize (speed 3) (space 0)
1616     (inhibit-warnings 3)))
1617 rtoy 1.2 (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 rtoy 1.4 (values (%make-double-double-float ext:double-float-positive-infinity 0d0)
1622 rtoy 1.2 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 rtoy 1.19 _N"Principle square root of Z
1644 rtoy 1.2
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 rtoy 1.17 (declare (optimize (speed 3) (space 0)
1660     (inhibit-warnings 3)))
1661 rtoy 1.2
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 rtoy 1.19 _N"Compute log(2^j*z).
1686 rtoy 1.2
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 rtoy 1.17 (declare (optimize (speed 3)
1704     (inhibit-warnings 3)))
1705 rtoy 1.2 (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 rtoy 1.19 _N"Log of Z = log |Z| + i * arg Z
1721 rtoy 1.2
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 rtoy 1.19 _N"Compute atanh z = (log(1+z) - log(1-z))/2"
1733 rtoy 1.2 (declare (number z))
1734 rtoy 1.12 (cond ((and (realp z) (< z -1))
1735     ;; ATANH is continuous with quadrant III in this case.
1736     (dd-complex-atanh (complex z -0d0)))
1737 rtoy 1.9 (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 rtoy 1.14 (y (* beta (- (float (imagpart z) 1.0w0))))
1746 rtoy 1.9 (eta 0.0w0)
1747     (nu 0.0w0))
1748     ;; Shouldn't need this declare.
1749     (declare (double-double-float x y))
1750     (locally
1751 rtoy 1.17 (declare (optimize (speed 3)
1752     (inhibit-warnings 3)))
1753 rtoy 1.9 (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 rtoy 1.14 (dd-%atan2 (* 2.0d0 y)
1786 rtoy 1.9 (- (* (- 1.0d0 x)
1787     (+ 1.0d0 x))
1788     (square t1))))))))
1789     (complex (* beta eta)
1790 rtoy 1.15 (- (* beta nu))))))))
1791 rtoy 1.2
1792     (defun dd-complex-tanh (z)
1793 rtoy 1.19 _N"Compute tanh z = sinh z / cosh z"
1794 rtoy 1.2 (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 rtoy 1.17 (declare (optimize (speed 3) (space 0)
1800     (inhibit-warnings 3)))
1801 rtoy 1.2 (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 rtoy 1.19 _N"Compute acos z = pi/2 - asin z
1874 rtoy 1.2
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 rtoy 1.19 _N"Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
1900 rtoy 1.2
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 rtoy 1.19 _N"Compute asin z = asinh(i*z)/i
1926 rtoy 1.2
1927     Z may be any number, but the result is always a complex."
1928     (declare (number z))
1929 rtoy 1.11 (cond ((and (realp z) (> z 1))
1930     (dd-complex-asin (complex z -0w0)))
1931 rtoy 1.8 (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 rtoy 1.11 ;; 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 rtoy 1.8 (- 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 rtoy 1.2
1953     (defun dd-complex-asinh (z)
1954 rtoy 1.19 _N"Compute asinh z = log(z + sqrt(1 + z*z))
1955 rtoy 1.2
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 rtoy 1.19 _N"Compute atan z = atanh (i*z) / i
1966 rtoy 1.2
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 rtoy 1.19 _N"Compute tan z = -i * tanh(i * z)
1977 rtoy 1.2
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