| blob | ed57efd9cb03057c6228bff26e6ef97cb748adc8 |
1 ;;;; -*- Mode: lisp -*-
2 ;;;;
3 ;;;; Copyright (c) 2007 Raymond Toy
4 ;;;;
5 ;;;; Permission is hereby granted, free of charge, to any person
6 ;;;; obtaining a copy of this software and associated documentation
7 ;;;; files (the "Software"), to deal in the Software without
8 ;;;; restriction, including without limitation the rights to use,
9 ;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
10 ;;;; copies of the Software, and to permit persons to whom the
11 ;;;; Software is furnished to do so, subject to the following
12 ;;;; conditions:
13 ;;;;
14 ;;;; The above copyright notice and this permission notice shall be
15 ;;;; included in all copies or substantial portions of the Software.
16 ;;;;
17 ;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
18 ;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
19 ;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
20 ;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
21 ;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
22 ;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
23 ;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
24 ;;;; OTHER DEALINGS IN THE SOFTWARE.
26 ;; Most of this code taken from CMUCL and slightly modified to support
27 ;; QD-COMPLEX.
134 #+cmu
136 ;; For now, convert B into a qd-complex and use that.
143 #-cmu
145 ;; For now, convert B into a qd-complex and use that.
194 ;; This can be simplified since X is real.
219 number)
226 number)
234 "Principle square root of Z
236 Z may be any number, but the result is always a complex."
249 ;; space 0 to get maybe-inline functions inlined.
273 "Compute log(2^j*z).
275 This is for use with J /= 0 only when |z| is huge."
278 ;; The constants t0, t1, t2 should be evaluated to machine
279 ;; precision. In addition, Kahan says the accuracy of log1p
280 ;; influences the choices of these constants but doesn't say how to
281 ;; choose them. We'll just assume his choices matches our
282 ;; implementation of log1p.
307 "Log of Z = log |Z| + i * arg Z
309 Z may be any number, but the result is always a complex."
314 ;; Let us note the following "strange" behavior. atanh 1.0d0 is
315 ;; +infinity, but the following code returns approx 176 + i*pi/4. The
316 ;; reason for the imaginary part is caused by the fact that arg i*y is
317 ;; never 0 since we have positive and negative zeroes.
319 ;;
320 ;; atanh(z) = (log(1+z) - log(1-z))/2
321 ;;
322 ;; The branch cut is on the real axis for |x| >= 1. For x =< -1,
323 ;; atanh is continuous with quadrant III; for x >= 1, continuous with
324 ;; quadrant I.
326 "Compute atanh z = (log(1+z) - log(1-z))/2"
332 ;; Carefully multiply a and b, taking care to handle
333 ;; signed zeroes. Only need to handle the case of b
334 ;; being zero.
337 #q-0
350 ;; Shouldn't need this declare.
356 ;; To avoid overflow...
358 ;; eta is real part of 1/(x + iy). This is x/(x^2+y^2),
359 ;; which can cause overflow. Arrange this computation so
360 ;; that it won't overflow.
368 ;; Should this be changed so that if y is zero, eta is set
369 ;; to +infinity instead of approx 176? In any case
370 ;; tanh(176) is 1.0d0 within working precision.
378 ))
381 ;; Normal case using log1p(x) = log(1 + x)
395 ;; tanh(z) = sinh(z)/cosh(z)
396 ;;
398 "Compute tanh z = sinh z / cosh z"
403 ;; space 0 to get maybe-inline functions inlined
408 ;; The threshold above is
409 ;; asinh(most-positive-double-float)/4, but many Lisps
410 ;; cannot actually compute that. Hence use the
411 ;; (accurate) approximation
412 ;; asinh(most-positive-double-float) =
413 ;; log(most-positive-double-float) + log(2)
417 ;; With quad-double's it happens that tan(pi/2) will
418 ;; actually produce a division by zero error. We need to
419 ;; handle that case carefully.
427 den)
430 ;; This means tan(y) produced some error. We'll
431 ;; assume it's an overflow error because y is
432 ;; pi/2 + 2*k*pi. But we need a value for tv to
433 ;; compute (/ tv). This would be a signed-zero.
434 ;; For now, just return +0.
438 ;; Kahan says we should only compute the parts needed. Thus, the
439 ;; realpart's below should only compute the real part, not the whole
440 ;; complex expression. Doing this can be important because we may get
441 ;; spurious signals that occur in the part that we are not using.
442 ;;
443 ;; However, we take a pragmatic approach and just use the whole
444 ;; expression.
446 ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
447 ;; it's the conjugate of the square root or the square root of the
448 ;; conjugate. This needs to be checked.
450 ;; I checked. It doesn't matter because (conjugate (sqrt z)) is the
451 ;; same as (sqrt (conjugate z)) for all z. This follows because
452 ;;
453 ;; (conjugate (sqrt z)) = exp(0.5*log |z|)*exp(-0.5*j*arg z).
454 ;;
455 ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
456 ;;
457 ;; and these two expressions are equal if and only if arg conj z =
458 ;; -arg z, which is clearly true for all z.
460 ;; NOTE: The rules of Common Lisp says that if you mix a real with a
461 ;; complex, the real is converted to a complex before performing the
462 ;; operation. However, Kahan says in this paper (pg 176):
463 ;;
464 ;; (iii) Careless handling can turn infinity or the sign of zero into
465 ;; misinformation that subsequently disappears leaving behind
466 ;; only a plausible but incorrect result. That is why compilers
467 ;; must not transform z-1 into z-(1+i*0), as we have seen above,
468 ;; nor -(-x-x^2) into (x+x^2), as we shall see below, lest a
469 ;; subsequent logarithm or square root produce a non-zero
470 ;; imaginary part whose sign is opposite to what was intended.
471 ;;
472 ;; The interesting examples are too long and complicated to reproduce
473 ;; here. We refer the reader to his paper.
474 ;;
475 ;; The functions below are intended to handle the cases where a real
476 ;; is mixed with a complex and we don't want CL complex contagion to
477 ;; occur..
490 "Compute acos z = pi/2 - asin z
492 Z may be any number, but the result is always a complex."
495 ;; acos is continuous in quadrant IV in this case.
500 ;; Same as below, but we compute atan ourselves (because we
501 ;; have atan +/- infinity).
508 sqrt-1-z)))))
513 sqrt-1-z)))))))))
516 "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
518 Z may be any number, but the result is always a complex."
522 ;; We need to handle the case where real part of sqrt-z+1 is zero,
523 ;; because division by zero with double-double-floats doesn't
524 ;; produce infinity.
526 ;; Same as below, but we compute atan ourselves (because we
527 ;; have atan +/- infinity).
540 ;; asin(z) = asinh(i*z)/i
541 ;; = -i log(i*z + sqrt(1-z^2))
543 "Compute asin z = asinh(i*z)/i
545 Z may be any number, but the result is always a complex."
554 ;; Like below but we handle atan part ourselves.
555 ;; Must be sure to take into account the sign of
556 ;; (realpart z) and den!
564 ;; We get a invalid operation here when z is real and |z| > 1.
571 "Compute asinh z = log(z + sqrt(1 + z*z))
573 Z may be any number, but the result is always a complex."
575 ;; asinh z = -i * asin (i*z)
582 "Compute atan z = atanh (i*z) / i
584 Z may be any number, but the result is always a complex."
586 ;; atan z = -i * atanh (i*z)
593 "Compute tan z = -i * tanh(i * z)
595 Z may be any number, but the result is always a complex."
597 ;; tan z = -i * tanh(i*z)
602 \f