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

  ViewVC Help
Powered by ViewVC 1.1.5