Enhancement to %HYPOT

Bug #888410 reported by Robert Smith on 2011-11-10
This bug affects 1 person
Affects Status Importance Assigned to Milestone

Bug Description

1. Lines 128-134 of `src/code/irrat.lisp` read

  ;; FIXME: libc hypot "computes the sqrt(x*x+y*y) without undue overflow or underflow"
  ;; ...we just do the stupid simple thing.
  (declaim (inline %hypot))
  (defun %hypot (x y)
    (sqrt (+ (* x x) (* y y)))))

For more accuracy, this can be replaced with

  ;; This is written in a peculiar way to avoid overflow. Note that in
  ;; sqrt(x^2 + y^2), either square or the sum can overflow.
  ;; Without loss of generality, assume |x| >= |y|. Then (y/x)^2 lies
  ;; in the unit interval [0,1], so 1 + (y/x)^2 lies in [1,2]. The
  ;; square root of that lies in [1,1.414] approximately.
  ;; Factoring x^2 out of sqrt(x^2 + y^2) gives us the expression
  ;; |x|sqrt(1 + (y/x)^2), which can only overflow if |x| is
  ;; sufficiently large.
  (declaim (inline %hypot))
  (defun %hypot (x y)
      ((zerop x) (abs y))
      ((zerop y) (abs x))
      ((= (abs x) (abs y)) (* (abs x) (sqrt 2)))
      (t (progn
           ;; Swap to make quotient in unit interval.
           (when (< (abs x) (abs y))
             (rotatef x y))

           ;; Compute the hypotenuse length.
           (let ((y/x (/ y x)))
             (* (abs x)
                (sqrt (1+ (* y/x y/x))))))))))

2. No real test cases.

3. SBCL 1.0.52

4. Darwin DiapentePomum.local 9.8.0 Darwin Kernel Version 9.8.0: Wed Jul 15 16:57:01 PDT 2009; root:xnu-1228.15.4~1/RELEASE_PPC Power Macintosh


description: updated
description: updated
description: updated
description: updated
description: updated
Paul Khuong (pvk) on 2011-11-11
Changed in sbcl:
status: New → Triaged
importance: Undecided → Low
assignee: nobody → Paul Khuong (pvk)
Paul Khuong (pvk) wrote :

Slightly modified and committed in 79c4a7f More numerically stable %hypot (ABS of complex floats) on win32.

Changed in sbcl:
status: Triaged → Fix Committed
assignee: Paul Khuong (pvk) → nobody
Changed in sbcl:
status: Fix Committed → Fix Released
To post a comment you must log in.
This report contains Public information  Edit
Everyone can see this information.

Other bug subscribers