From: Michael Orlitzky <michael@orlitzky.com>
Date: Mon, 1 Jul 2013 20:00:59 +0000 (-0400)
Subject: Add the unit_roundoff function.
X-Git-Url: http://gitweb.michael.orlitzky.com/?p=numerical-analysis.git;a=commitdiff_plain;h=a615055dd6d6a0bd59aa7ba03305388f786306cb;hp=04f56a8882bb0c574b603f8c3fed9481ea934f7f

Add the unit_roundoff function.
---

diff --git a/src/Misc.hs b/src/Misc.hs
index 24fd5b7..ddaad6e 100644
--- 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) ]