Add the unit_roundoff function.

diff --git a/src/Misc.hs b/src/Misc.hs
index 24fd5b756febd19bd512cf9a7b06682f2aba9f7d..ddaad6e3ba7061f9a3ee2c119998058b3d298201 100644 (file)
--- a/src/Misc.hs
+++ b/src/Misc.hs
@@ -1,4 +1,5 @@
 {-# LANGUAGE RebindableSyntax #-}
+{-# LANGUAGE ScopedTypeVariables #-}
 
 -- | Stuff for which I'm too lazy to come up with a decent name.
 module Misc
@@ -6,6 +7,7 @@ where
 
 import NumericPrelude
 import Algebra.Field
+import Algebra.RealRing
 import Algebra.ToInteger
 
 -- | Partition the interval [@a@, @b@] into @n@ subintervals, which we
@@ -35,3 +37,13 @@ partition n a b
                  let xj = a + (k'+1)*h ]
     where
       h = (b-a)/(fromIntegral $ toInteger n)
+
+
+-- | Compute the unit roundoff (machine epsilon) for this machine. We
+--   find the largest number epsilon such that 1+epsilon <= 1. If you
+--   request anything other than a Float or Double from this, expect
+--   to wait a while.
+--
+unit_roundoff :: forall a. (Algebra.RealRing.C a, Algebra.Field.C a) => a
+unit_roundoff =
+  head [ 1/2^(k-1) | k <- [0..], 1 + 1/(2^k) <= (1::a) ]