Skip to content

Instantly share code, notes, and snippets.

@rrnewton
Created May 14, 2012 05:09
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save rrnewton/2691877 to your computer and use it in GitHub Desktop.
Save rrnewton/2691877 to your computer and use it in GitHub Desktop.
Better Haskell magnitude
import Data.Complex
-- This should be the biggest representable double:
big :: Double
big = 1.7976931348623157e308
-- This uses OCaml's approach at overflow-avoidance:
mag :: Complex Double -> Double
mag (x:+y) =
if r == 0 then i
else if i == 0 then r
else if r >= i then let q = i / r in r * sqrt(1 + q*q)
else let q = r / i in i * sqrt(1 + q*q)
where r = abs x
i = abs y
t0 = mag (big :+ big) -- Inf
t1 = mag (0 :+ big)
t2 = mag (big :+ 0)
t3 = magnitude (big :+ big) -- Inf (has to be)
biggish = big / 2
t4 = magnitude (biggish :+ biggish)
t5 = mag (biggish :+ biggish)
t6 = sqrt (biggish * biggish + biggish * biggish)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment