Changeset 14673


Ignore:
Timestamp:
Mar 11, 2011, 7:30:08 PM (9 years ago)
Author:
rme
Message:

When computing the absolute value of a complex number,
be a little more careful about it.

hypot in libm is even more careful, and one could
certainly make the argument that we should just arrange
to ff-call that.

This new version conses more than what we were
previously doing, but seems to be about the same speed.
If this ends up being a performance issue, writing a
LAP version probably would not be too tough.

Fixes ticket:830.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/source/level-0/l0-numbers.lisp

    r14119 r14673  
    915915         (sqrt (+ (* rx rx) (* ix ix))))
    916916        (short-float
    917          (%short-float (%hypot (%double-float rx)
    918                                (%double-float ix))))
     917         (%short-float (%double-float-hypot (%double-float rx)
     918                                            (%double-float ix))))
    919919        (double-float
    920          (%hypot rx ix)))))))
     920         (%double-float-hypot rx ix)))))))
    921921
    922922
     
    20022002  (/ 1 n))
    20032003
    2004 ; x & y must both be double floats
    2005 (defun %hypot (x y)
    2006   (with-stack-double-floats ((x**2) (y**2))
    2007     (let ((res**2 x**2))
    2008       (%double-float*-2! x x x**2)
    2009       (%double-float*-2! y y y**2)
    2010       (%double-float+-2! x**2 y**2 res**2)
    2011       (fsqrt res**2))))
    2012 
    2013 
     2004;; Compute (sqrt (+ (* x x) (* y y))), but
     2005;; try to be a little more careful about it.
     2006;; Both x and y must be double-floats.
     2007(defun %double-float-hypot (x y)
     2008  (with-stack-double-floats ((a) (b) (c))
     2009    (%%double-float-abs! x a)
     2010    (%%double-float-abs! y b)
     2011    (when (> a b)
     2012      (psetq a b b a))
     2013    (if (= b 0d0)
     2014      0d0
     2015      (progn
     2016        (%double-float/-2! a b c)
     2017        (* b (fsqrt (+ 1d0 (* c c))))))))
     2018                                       
Note: See TracChangeset for help on using the changeset viewer.