/[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.5 - (hide annotations)
Wed Jul 19 15:02:00 2006 UTC (7 years, 8 months ago) by rtoy
Branch: MAIN
Changes since 1.4: +116 -7 lines
Implement accurate arg reduction for the trig functions sin, cos, tan,
using __kernel_rem_pio2.  Use the accurate reduction in dd-%sin,
dd-%cos, dd-%tan.  The original version were renamed dd-%%sin,
dd-%%cos, dd-%%tan because we still need them to evaluate the
functions for small args.

This also has the side effect that we don't signal overflow for
moderately large args (~1w50) anymore.
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.5 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/irrat-dd.lisp,v 1.5 2006/07/19 15:02: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     (if (< (abs x) (/ pi 4))
1250     (dd-%%sin x)
1251     ;; Argument reduction needed
1252     (multiple-value-bind (n reduced)
1253     (reduce-arg x)
1254     (case (logand n 3)
1255     (0 (dd-%%sin reduced))
1256     (1 (dd-%%cos reduced))
1257     (2 (- (dd-%%sin reduced)))
1258     (3 (- (dd-%%cos reduced)))))))
1259    
1260     (defun dd-%cos (x)
1261     (declare (double-double-float x))
1262     (if (< (abs x) (/ pi 4))
1263     (dd-%%cos x)
1264     ;; Argument reduction needed
1265     (multiple-value-bind (n reduced)
1266     (reduce-arg x)
1267     (case (logand n 3)
1268     (0 (dd-%%cos reduced))
1269     (1 (- (dd-%%sin reduced)))
1270     (2 (- (dd-%%cos reduced)))
1271     (3 (dd-%%sin reduced))))))
1272    
1273     (defun dd-%tan (x)
1274     (declare (double-double-float x))
1275     (if (< (abs x) (/ pi 4))
1276     (dd-%%tan x)
1277     ;; Argument reduction needed
1278     (multiple-value-bind (n reduced)
1279     (reduce-arg x)
1280     (if (evenp n)
1281     (dd-%%tan reduced)
1282     (- (/ (dd-%%tan reduced)))))))
1283    
1284 rtoy 1.2 ;;; dd-%log2
1285     ;;; Base 2 logarithm.
1286    
1287     (let ((P (make-array 13 :element-type 'double-double-float
1288     :initial-contents
1289     '(
1290     1.538612243596254322971797716843006400388w-6
1291     4.998469661968096229986658302195402690910w-1
1292     2.321125933898420063925789532045674660756w1
1293     4.114517881637811823002128927449878962058w2
1294     3.824952356185897735160588078446136783779w3
1295     2.128857716871515081352991964243375186031w4
1296     7.594356839258970405033155585486712125861w4
1297     1.797628303815655343403735250238293741397w5
1298     2.854829159639697837788887080758954924001w5
1299     3.007007295140399532324943111654767187848w5
1300     2.014652742082537582487669938141683759923w5
1301     7.771154681358524243729929227226708890930w4
1302     1.313572404063446165910279910527789794488w4
1303     )))
1304     (Q (make-array 12 :element-type 'double-double-float
1305     :initial-contents
1306     '(
1307     ;; 1.000000000000000000000000000000000000000w0
1308     4.839208193348159620282142911143429644326w1
1309     9.104928120962988414618126155557301584078w2
1310     9.147150349299596453976674231612674085381w3
1311     5.605842085972455027590989944010492125825w4
1312     2.248234257620569139969141618556349415120w5
1313     6.132189329546557743179177159925690841200w5
1314     1.158019977462989115839826904108208787040w6
1315     1.514882452993549494932585972882995548426w6
1316     1.347518538384329112529391120390701166528w6
1317     7.777690340007566932935753241556479363645w5
1318     2.626900195321832660448791748036714883242w5
1319     3.940717212190338497730839731583397586124w4
1320     )))
1321     ;; Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
1322     ;; where z = 2(x-1)/(x+1)
1323     ;; 1/sqrt(2) <= x < sqrt(2)
1324     ;; Theoretical peak relative error = 1.1e-35,
1325     ;; relative peak error spread 1.1e-9
1326     (R (make-array 6 :element-type 'double-double-float
1327     :initial-contents
1328     '(
1329     -8.828896441624934385266096344596648080902w-1
1330     8.057002716646055371965756206836056074715w1
1331     -2.024301798136027039250415126250455056397w3
1332     2.048819892795278657810231591630928516206w4
1333     -8.977257995689735303686582344659576526998w4
1334     1.418134209872192732479751274970992665513w5
1335     )))
1336     (S (make-array 6 :element-type 'double-double-float
1337     :initial-contents
1338     '(
1339     ;; 1.000000000000000000000000000000000000000w0 */
1340     -1.186359407982897997337150403816839480438w2
1341     3.998526750980007367835804959888064681098w3
1342     -5.748542087379434595104154610899551484314w4
1343     4.001557694070773974936904547424676279307w5
1344     -1.332535117259762928288745111081235577029w6
1345     1.701761051846631278975701529965589676574w6
1346     ))))
1347     (defun dd-%log2 (x)
1348     (declare (type double-double-float x)
1349 rtoy 1.3 (optimize (speed 3) (space 0)))
1350 rtoy 1.2 (multiple-value-bind (x e)
1351     (decode-float x)
1352     (declare (type double-double-float x)
1353     (type double-float-exponent e))
1354     (let ((z 0w0)
1355     (y 0w0))
1356     (declare (type double-double-float z y))
1357     (cond ((or (> e 2)
1358     (< e -2))
1359     (cond ((< x sqrt-1/2)
1360     ;; 2*(2*x-1)/(2*x+1)
1361     (decf e)
1362     (setf z (- x 0.5w0))
1363     (setf y (+ (* 0.5w0 z) 0.5w0)))
1364     (t
1365     ;; 2*(x-1)/(x+1)
1366     (setf z (- x 0.5w0))
1367     (decf z 0.5w0)
1368     (setf y (+ (* 0.5w0 z) 0.5w0))))
1369     (setf x (/ z y))
1370     (setf z (* x x))
1371     (setf y (* x (/ (* z (poly-eval z r))
1372     (poly-eval-1 z s)))))
1373     (t
1374     (cond ((< x sqrt-1/2)
1375     (decf e)
1376     (setf x (- (scale-float x 1) 1)))
1377     (t
1378     (decf x)))
1379     (setf z (* x x))
1380     (setf y (* x (/ (* z (poly-eval x p))
1381     (poly-eval-1 x q))))
1382     (decf y (scale-float z -1))))
1383     ;; Multiply log of fraction by log2(e) and base 2 exponent by 1
1384     ;;
1385     ;; This sequence of operations is critical
1386     (setf z (* y log2ea))
1387     (setf z (+ z (* x log2ea)))
1388     (setf z (+ z y))
1389     (setf z (+ z x))
1390     (setf z (+ z e))
1391     z))))
1392    
1393     ;;; dd-%exp2
1394     ;;; 2^x
1395    
1396     (let ((P (make-array 5 :element-type 'double-double-float
1397     :initial-contents
1398     '(
1399     1.587171580015525194694938306936721666031w2
1400     6.185032670011643762127954396427045467506w5
1401     5.677513871931844661829755443994214173883w8
1402     1.530625323728429161131811299626419117557w11
1403     9.079594442980146270952372234833529694788w12
1404     )))
1405     (Q (make-array 5 :element-type 'double-double-float
1406     :initial-contents
1407     '(
1408     ;; 1.000000000000000000000000000000000000000w0 */
1409     1.236602014442099053716561665053645270207w4
1410     2.186249607051644894762167991800811827835w7
1411     1.092141473886177435056423606755843616331w10
1412     1.490560994263653042761789432690793026977w12
1413     2.619817175234089411411070339065679229869w13
1414     ))))
1415     (defun dd-%exp2 (x)
1416     (declare (type double-double-float x)
1417 rtoy 1.3 (optimize (speed 3) (space 0)))
1418 rtoy 1.2 (when (>= x 1024w0)
1419     (return-from dd-%exp2
1420     (%make-double-double-float ext:double-float-positive-infinity
1421     ext:double-float-positive-infinity)))
1422     (when (<= x -1024w0)
1423     (return-from dd-%exp2 0w0))
1424     (multiple-value-bind (n x)
1425     (floor (the (double-double-float -1024w0 1024w0) x))
1426     (declare (type double-double-float x))
1427     (let* ((xx (* x x))
1428     (px (* x (poly-eval xx p))))
1429     (setf x (/ px (- (poly-eval-1 xx q) px)))
1430     (setf x (+ 1 (scale-float x 1)))
1431     (scale-float x n)))))
1432    
1433     ;;; dd-%powil
1434     ;;; x^n where n is an integer
1435    
1436     (defun dd-%powil (x nn)
1437     (declare (type double-double-float x)
1438     (fixnum nn)
1439 rtoy 1.3 (optimize (speed 3) (space 0)))
1440 rtoy 1.2 (when (zerop x)
1441     (return-from dd-%powil
1442     (cond ((zerop nn)
1443     1w0)
1444     ((minusp nn)
1445     (%make-double-double-float ext:double-float-positive-infinity
1446     ext:double-float-positive-infinity))
1447     (t
1448     0w0))))
1449     (when (zerop nn)
1450     (return-from dd-%powil 1w0))
1451    
1452     (let ((asign 0)
1453     (sign 0)
1454     (n nn))
1455     (declare (type (integer -1 0) asign sign)
1456     (fixnum n))
1457     (when (minusp x)
1458     (setf asign -1)
1459     (setf x (- x)))
1460     (cond ((minusp nn)
1461     (setf sign -1)
1462     (setf n (- nn)))
1463     (t
1464     (setf sign 0)
1465     (setf n nn)))
1466     ;; Overflow detection
1467    
1468     ;; Calculate approximate log of answer
1469     (multiple-value-bind (s lx)
1470     (decode-float x)
1471     (declare (type double-double-float s)
1472     (type double-float-exponent lx))
1473     (let ((e (* n (1- lx))))
1474     (cond ((or (zerop e)
1475     (> e 64)
1476     (< e -64))
1477     (setf s (/ (- s 7.0710678118654752w-1)
1478     (+ s 7.0710678118654752w-1)))
1479     (setf s (* (+ (- (* 2.9142135623730950w0 s)
1480     0.5w0)
1481     lx)
1482     nn log2e)))
1483     (t
1484     (setf s (* loge2 e))))
1485     (when (> s max-log)
1486     ;; Overflow. What to do?
1487     (error "Overflow"))
1488     (when (< s min-log)
1489     (return-from dd-%powil 0w0))
1490    
1491     ;; Handle denormal answer
1492     (when (< s (- 2 max-log))
1493     (setf x (/ x))
1494     (setf sign 0))
1495    
1496     ;; First bit of the power
1497     (let ((y 0w0))
1498     (declare (type double-double-float y))
1499     (cond ((zerop (logand n 1))
1500     (setf y 1w0)
1501     (setf asign 0))
1502     (t
1503     (setf y x)))
1504     (let ((w x))
1505     (declare (type double-double-float w))
1506     (setf n (ash n -1))
1507     (loop while (not (zerop n))
1508     do
1509     (setf w (* w w))
1510     (unless (zerop (logand n 1))
1511     (setf y (* y w)))
1512     (setf n (ash n -1))))
1513     ;; Odd power of a negative number
1514     (unless (zerop asign)
1515     (setf y (- y)))
1516    
1517     (unless (zerop sign)
1518     (setf y (/ y)))
1519     y)))))
1520    
1521    
1522     ;;; dd-real-pow
1523     ;;; x^y, for x and y real, and real result.
1524    
1525     (defun dd-real-pow (x y)
1526     (declare (type double-double-float x y)
1527 rtoy 1.3 (optimize (speed 3) (space 0)))
1528 rtoy 1.2 (let ((nflg 0)
1529     (w (floor y)))
1530     ;; nflg = 1 if x < 0 raised to integer power
1531     (when (and (= w y)
1532     (< (abs w) 32768))
1533     (return-from dd-real-pow (dd-%powil x w)))
1534    
1535     (when (<= x 0)
1536     (cond ((zerop x)
1537     (if (zerop y)
1538     ;; 0^0
1539     (return-from dd-real-pow 1w0)
1540     ;; 0^y
1541     (return-from dd-real-pow 0w0)))
1542     (t
1543     (when (/= w y)
1544     ;; noninteger power of negative number
1545     (let ((p (the double-double-float (dd-real-pow (abs x) y)))
1546     (y*pi (* y dd-pi)))
1547     (return-from dd-real-pow (complex (* p (dd-%cos y*pi))
1548     (* p (dd-%sin y*pi))))))
1549    
1550     ;; For negative x, find out if the integer exponent is odd or even.
1551     (let ((w (scale-float y -1)))
1552     (declare (type double-double-float w))
1553     (setf w (ffloor w))
1554     (setf w (scale-float w 1))
1555     (when (/= w y)
1556     (setf nflg 1))
1557     (setf x (abs x))))))
1558     (let ((z (dd-%log2 x)))
1559     (declare (type double-double-float z))
1560     (setf z (dd-%exp2 (* y z)))
1561     (unless (zerop nflg)
1562     (setf z (- z)))
1563     z)))
1564    
1565     (defun dd-%pow (x y)
1566     (declare (type double-double-float x y))
1567     (dd-real-pow x y))
1568    
1569    
1570     ;; These are essentially the same as in irrat.lisp, but very slightly
1571     ;; modified to operate on double-double-floats.
1572     (defun dd-cssqs (z)
1573     ;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
1574     ;; result is r + i*k, where k is an integer.
1575    
1576     ;; Save all FP flags
1577     (let ((x (float (realpart z) 1w0))
1578     (y (float (imagpart z) 1w0)))
1579     ;; Would this be better handled using an exception handler to
1580     ;; catch the overflow or underflow signal? For now, we turn all
1581     ;; traps off and look at the accrued exceptions to see if any
1582     ;; signal would have been raised.
1583 rtoy 1.5 ;;
1584     ;; Actually, for double-double-floats, we should probably
1585     ;; explicitly check for overflow instead of disabling the traps.
1586     ;; Why? Because instead of overflow, double-double signals
1587     ;; invalid operation.
1588     (with-float-traps-masked (:underflow :overflow :invalid)
1589 rtoy 1.2 (let ((rho (+ (square x) (square y))))
1590     (declare (optimize (speed 3) (space 0)))
1591     (cond ((and (or (float-nan-p rho)
1592     (float-infinity-p rho))
1593     (or (float-infinity-p (abs x))
1594     (float-infinity-p (abs y))))
1595 rtoy 1.4 (values (%make-double-double-float ext:double-float-positive-infinity 0d0)
1596 rtoy 1.2 0))
1597     ((let ((threshold #.(/ least-positive-double-float
1598     double-float-epsilon))
1599     (traps (ldb vm::float-sticky-bits
1600     (vm:floating-point-modes))))
1601     ;; Overflow raised or (underflow raised and rho <
1602     ;; lambda/eps)
1603     (or (not (zerop (logand vm:float-overflow-trap-bit traps)))
1604     (and (not (zerop (logand vm:float-underflow-trap-bit traps)))
1605     (< rho threshold))))
1606     ;; If we're here, neither x nor y are infinity and at
1607     ;; least one is non-zero.. Thus logb returns a nice
1608     ;; integer.
1609     (let ((k (- (logb-finite (max (abs x) (abs y))))))
1610     (values (+ (square (scalb x k))
1611     (square (scalb y k)))
1612     (- k))))
1613     (t
1614     (values rho 0)))))))
1615    
1616     (defun dd-complex-sqrt (z)
1617     "Principle square root of Z
1618    
1619     Z may be any number, but the result is always a complex."
1620     (declare (number z))
1621     (multiple-value-bind (rho k)
1622     (dd-cssqs z)
1623     (declare (type (or (member 0w0) (double-double-float 0w0)) rho)
1624     (type fixnum k))
1625     (let ((x (float (realpart z) 1.0w0))
1626     (y (float (imagpart z) 1.0w0))
1627     (eta 0w0)
1628     (nu 0w0))
1629     (declare (type double-double-float x y eta nu))
1630    
1631     (locally
1632     ;; space 0 to get maybe-inline functions inlined.
1633     (declare (optimize (speed 3) (space 0)))
1634    
1635     (if (not (float-nan-p x))
1636     (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
1637    
1638     (cond ((oddp k)
1639     (setf k (ash k -1)))
1640     (t
1641     (setf k (1- (ash k -1)))
1642     (setf rho (+ rho rho))))
1643    
1644     (setf rho (scalb (sqrt rho) k))
1645    
1646     (setf eta rho)
1647     (setf nu y)
1648    
1649     (when (/= rho 0d0)
1650     (when (not (float-infinity-p (abs nu)))
1651     (setf nu (/ (/ nu rho) 2d0)))
1652     (when (< x 0d0)
1653     (setf eta (abs nu))
1654     (setf nu (float-sign y rho))))
1655     (complex eta nu)))))
1656    
1657     (defun dd-complex-log-scaled (z j)
1658     "Compute log(2^j*z).
1659    
1660     This is for use with J /= 0 only when |z| is huge."
1661     (declare (number z)
1662     (fixnum j))
1663     ;; The constants t0, t1, t2 should be evaluated to machine
1664     ;; precision. In addition, Kahan says the accuracy of log1p
1665     ;; influences the choices of these constants but doesn't say how to
1666     ;; choose them. We'll just assume his choices matches our
1667     ;; implementation of log1p.
1668     (let ((t0 #.(/ 1 (sqrt 2.0w0)))
1669     (t1 1.2w0)
1670     (t2 3w0)
1671     (ln2 #.(log 2w0))
1672     (x (float (realpart z) 1.0w0))
1673     (y (float (imagpart z) 1.0w0)))
1674     (multiple-value-bind (rho k)
1675     (dd-cssqs z)
1676     (declare (optimize (speed 3)))
1677     (let ((beta (max (abs x) (abs y)))
1678     (theta (min (abs x) (abs y))))
1679     (complex (if (and (zerop k)
1680     (< t0 beta)
1681     (or (<= beta t1)
1682     (< rho t2)))
1683     (/ (dd-%log1p (+ (* (- beta 1.0d0)
1684     (+ beta 1.0d0))
1685     (* theta theta)))
1686     2d0)
1687     (+ (/ (log rho) 2d0)
1688     (* (+ k j) ln2)))
1689     (atan y x))))))
1690    
1691     (defun dd-complex-log (z)
1692     "Log of Z = log |Z| + i * arg Z
1693    
1694     Z may be any number, but the result is always a complex."
1695     (declare (number z))
1696     (dd-complex-log-scaled z 0))
1697    
1698     ;; Let us note the following "strange" behavior. atanh 1.0d0 is
1699     ;; +infinity, but the following code returns approx 176 + i*pi/4. The
1700     ;; reason for the imaginary part is caused by the fact that arg i*y is
1701     ;; never 0 since we have positive and negative zeroes.
1702    
1703     (defun dd-complex-atanh (z)
1704     "Compute atanh z = (log(1+z) - log(1-z))/2"
1705     (declare (number z))
1706     (if (and (realp z) (< z -1))
1707     ;; atanh is continuous in quadrant III in this case.
1708     (dd-complex-atanh (complex z -0f0))
1709     (let* ( ;; Constants
1710     (theta (/ (sqrt most-positive-double-float) 4.0w0))
1711     (rho (/ 4.0w0 (sqrt most-positive-double-float)))
1712     (half-pi dd-pi/2)
1713     (rp (float (realpart z) 1.0w0))
1714     (beta (float-sign rp 1.0w0))
1715     (x (* beta rp))
1716     (y (* beta (- (float (imagpart z) 1.0w0))))
1717     (eta 0.0w0)
1718     (nu 0.0w0))
1719     ;; Shouldn't need this declare.
1720     (declare (double-double-float x y))
1721     (locally
1722     (declare (optimize (speed 3)))
1723     (cond ((or (> x theta)
1724     (> (abs y) theta))
1725     ;; To avoid overflow...
1726     (setf nu (float-sign y half-pi))
1727     ;; eta is real part of 1/(x + iy). This is x/(x^2+y^2),
1728     ;; which can cause overflow. Arrange this computation so
1729     ;; that it won't overflow.
1730     (setf eta (let* ((x-bigger (> x (abs y)))
1731     (r (if x-bigger (/ y x) (/ x y)))
1732     (d (+ 1.0d0 (* r r))))
1733     (if x-bigger
1734     (/ (/ x) d)
1735     (/ (/ r y) d)))))
1736     ((= x 1.0w0)
1737     ;; Should this be changed so that if y is zero, eta is set
1738     ;; to +infinity instead of approx 176? In any case
1739     ;; tanh(176) is 1.0d0 within working precision.
1740     (let ((t1 (+ 4w0 (square y)))
1741     (t2 (+ (abs y) rho)))
1742     (setf eta (dd-%log (/ (sqrt (sqrt t1))
1743     (sqrt t2))))
1744     (setf nu (* 0.5d0
1745     (float-sign y
1746     (+ half-pi (dd-%atan (* 0.5d0 t2))))))))
1747     (t
1748     (let ((t1 (+ (abs y) rho)))
1749     ;; Normal case using log1p(x) = log(1 + x)
1750     (setf eta (* 0.25d0
1751     (dd-%log1p (/ (* 4.0d0 x)
1752     (+ (square (- 1.0d0 x))
1753     (square t1))))))
1754     (setf nu (* 0.5d0
1755     (dd-%atan2 (* 2.0d0 y)
1756     (- (* (- 1.0d0 x)
1757     (+ 1.0d0 x))
1758     (square t1))))))))
1759     (complex (* beta eta)
1760     (- (* beta nu)))))))
1761    
1762     (defun dd-complex-tanh (z)
1763     "Compute tanh z = sinh z / cosh z"
1764     (declare (number z))
1765     (let ((x (float (realpart z) 1.0w0))
1766     (y (float (imagpart z) 1.0w0)))
1767     (locally
1768     ;; space 0 to get maybe-inline functions inlined
1769     (declare (optimize (speed 3) (space 0)))
1770     (cond ((> (abs x)
1771     #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)
1772     ;; This is more accurate under linux.
1773     #+(or linux hpux) #.(/ (+ (%log 2.0d0)
1774     (%log most-positive-double-float)) 4d0))
1775     (complex (float-sign x)
1776     (float-sign y)))
1777     (t
1778     (let* ((tv (dd-%tan y))
1779     (beta (+ 1.0d0 (* tv tv)))
1780     (s (sinh x))
1781     (rho (sqrt (+ 1.0d0 (* s s)))))
1782     (if (float-infinity-p (abs tv))
1783     (complex (/ rho s)
1784     (/ tv))
1785     (let ((den (+ 1.0d0 (* beta s s))))
1786     (complex (/ (* beta rho s)
1787     den)
1788     (/ tv den))))))))))
1789    
1790     ;; Kahan says we should only compute the parts needed. Thus, the
1791     ;; realpart's below should only compute the real part, not the whole
1792     ;; complex expression. Doing this can be important because we may get
1793     ;; spurious signals that occur in the part that we are not using.
1794     ;;
1795     ;; However, we take a pragmatic approach and just use the whole
1796     ;; expression.
1797    
1798     ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
1799     ;; it's the conjugate of the square root or the square root of the
1800     ;; conjugate. This needs to be checked.
1801    
1802     ;; I checked. It doesn't matter because (conjugate (sqrt z)) is the
1803     ;; same as (sqrt (conjugate z)) for all z. This follows because
1804     ;;
1805     ;; (conjugate (sqrt z)) = exp(0.5*log |z|)*exp(-0.5*j*arg z).
1806     ;;
1807     ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
1808     ;;
1809     ;; and these two expressions are equal if and only if arg conj z =
1810     ;; -arg z, which is clearly true for all z.
1811    
1812     ;; NOTE: The rules of Common Lisp says that if you mix a real with a
1813     ;; complex, the real is converted to a complex before performing the
1814     ;; operation. However, Kahan says in this paper (pg 176):
1815     ;;
1816     ;; (iii) Careless handling can turn infinity or the sign of zero into
1817     ;; misinformation that subsequently disappears leaving behind
1818     ;; only a plausible but incorrect result. That is why compilers
1819     ;; must not transform z-1 into z-(1+i*0), as we have seen above,
1820     ;; nor -(-x-x^2) into (x+x^2), as we shall see below, lest a
1821     ;; subsequent logarithm or square root produce a non-zero
1822     ;; imaginary part whose sign is opposite to what was intended.
1823     ;;
1824     ;; The interesting examples are too long and complicated to reproduce
1825     ;; here. We refer the reader to his paper.
1826     ;;
1827     ;; The functions below are intended to handle the cases where a real
1828     ;; is mixed with a complex and we don't want CL complex contagion to
1829     ;; occur..
1830    
1831     (declaim (inline 1+z 1-z z-1 z+1))
1832     (defun dd-1+z (z)
1833     (complex (+ 1 (realpart z)) (imagpart z)))
1834     (defun dd-1-z (z)
1835     (complex (- 1 (realpart z)) (- (imagpart z))))
1836     (defun dd-z-1 (z)
1837     (complex (- (realpart z) 1) (imagpart z)))
1838     (defun dd-z+1 (z)
1839     (complex (+ (realpart z) 1) (imagpart z)))
1840    
1841     (defun dd-complex-acos (z)
1842     "Compute acos z = pi/2 - asin z
1843    
1844     Z may be any number, but the result is always a complex."
1845     (declare (number z))
1846     (if (and (realp z) (> z 1))
1847     ;; acos is continuous in quadrant IV in this case.
1848     (complex-acos (complex z -0f0))
1849     (let ((sqrt-1+z (complex-sqrt (1+z z)))
1850     (sqrt-1-z (complex-sqrt (1-z z))))
1851     (cond ((zerop (realpart sqrt-1+z))
1852     ;; Same as below, but we compute atan ourselves (because we
1853     ;; have atan +/- infinity).
1854     (complex
1855     (if (minusp (float-sign (* (realpart sqrt-1-z)
1856     (realpart sqrt-1+z))))
1857     (- dd-pi)
1858     dd-pi)
1859     (asinh (imagpart (* (conjugate sqrt-1+z)
1860     sqrt-1-z)))))
1861     (t
1862     (complex (* 2 (atan (/ (realpart sqrt-1-z)
1863     (realpart sqrt-1+z))))
1864     (asinh (imagpart (* (conjugate sqrt-1+z)
1865     sqrt-1-z)))))))))
1866    
1867     (defun dd-complex-acosh (z)
1868     "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
1869    
1870     Z may be any number, but the result is always a complex."
1871     (declare (number z))
1872     (let* ((sqrt-z-1 (complex-sqrt (z-1 z)))
1873     (sqrt-z+1 (complex-sqrt (z+1 z))))
1874     ;; We need to handle the case where real part of sqrt-z+1 is zero,
1875     ;; because division by zero with double-double-floats doesn't
1876     ;; produce infinity.
1877     (cond ((zerop (realpart sqrt-z+1))
1878     ;; Same as below, but we compute atan ourselves (because we
1879     ;; have atan +/- infinity).
1880     (complex (asinh (realpart (* (conjugate sqrt-z-1)
1881     sqrt-z+1)))
1882     (if (minusp (float-sign (* (imagpart sqrt-z-1)
1883     (realpart sqrt-z+1))))
1884     (- dd-pi)
1885     dd-pi)))
1886     (t
1887     (complex (asinh (realpart (* (conjugate sqrt-z-1)
1888     sqrt-z+1)))
1889     (* 2 (atan (/ (imagpart sqrt-z-1)
1890     (realpart sqrt-z+1)))))))))
1891    
1892    
1893     (defun dd-complex-asin (z)
1894     "Compute asin z = asinh(i*z)/i
1895    
1896     Z may be any number, but the result is always a complex."
1897     (declare (number z))
1898     (if (and (realp z) (> z 1))
1899     ;; asin is continuous in quadrant IV in this case.
1900     (dd-complex-asin (complex z -0f0))
1901     (let* ((sqrt-1-z (complex-sqrt (1-z z)))
1902     (sqrt-1+z (complex-sqrt (1+z z)))
1903     (den (realpart (* sqrt-1-z sqrt-1+z))))
1904     (cond ((zerop den)
1905     ;; Like below but we handle atan part ourselves.
1906     (complex (if (minusp (float-sign den))
1907     (- dd-pi/2)
1908     dd-pi/2)
1909     (asinh (imagpart (* (conjugate sqrt-1-z)
1910     sqrt-1+z)))))
1911     (t
1912     (with-float-traps-masked (:divide-by-zero)
1913     ;; We get a invalid operation here when z is real and |z| > 1.
1914     (complex (atan (/ (realpart z)
1915     (realpart (* sqrt-1-z sqrt-1+z))))
1916     (asinh (imagpart (* (conjugate sqrt-1-z)
1917     sqrt-1+z))))))))))
1918    
1919     (defun dd-complex-asinh (z)
1920     "Compute asinh z = log(z + sqrt(1 + z*z))
1921    
1922     Z may be any number, but the result is always a complex."
1923     (declare (number z))
1924     ;; asinh z = -i * asin (i*z)
1925     (let* ((iz (complex (- (imagpart z)) (realpart z)))
1926     (result (complex-asin iz)))
1927     (complex (imagpart result)
1928     (- (realpart result)))))
1929    
1930     (defun dd-complex-atan (z)
1931     "Compute atan z = atanh (i*z) / i
1932    
1933     Z may be any number, but the result is always a complex."
1934     (declare (number z))
1935     ;; atan z = -i * atanh (i*z)
1936     (let* ((iz (complex (- (imagpart z)) (realpart z)))
1937     (result (complex-atanh iz)))
1938     (complex (imagpart result)
1939     (- (realpart result)))))
1940    
1941     (defun dd-complex-tan (z)
1942     "Compute tan z = -i * tanh(i * z)
1943    
1944     Z may be any number, but the result is always a complex."
1945     (declare (number z))
1946     ;; tan z = -i * tanh(i*z)
1947     (let* ((iz (complex (- (imagpart z)) (realpart z)))
1948     (result (complex-tanh iz)))
1949     (complex (imagpart result)
1950     (- (realpart result)))))

  ViewVC Help
Powered by ViewVC 1.1.5