i

i

i



x i+1 = (l i +u i )/2

(l i+1 , u i+1 ) = (l i , x i+1 ) if f(x i+1 )>0

(l i+1 , u i+1 ) = (x i+1 , u i ) otherwise



0

0

i+1

i

i

i

i

i

> import Ratio > infixl 6 .* > ivt0 :: (Float -> Float) -> Float > ivt0 f = ivt' 0 1 where > ivt' l u = > let z = 0.5*(l+u) > in if z==l || z==u || f z == 0 > then z > else if f z < 0 > then ivt' z u > else ivt' l z

> f0, g0 :: Float -> Float > f0 x = 3*x-1 > g0 x = 2*x-1

ivt0 f0

ivt0 f1

f (ivt f) == 0

Float

f x == 0

Float

Float

i

i

> data R = R { runR :: Integer -> Rational }

> instance Show R where > show (R f) = show (fromRational (f 30) :: Float)

Float

> instance Fractional R where > fromRational x = R $ const $ fromRational x

Eq

Num

> instance Eq R where > R f == R g = eq' 2 0 f g where > eq' delta i f g = > let fi = f i > gi = g i > in if fi<gi-delta > then False > else if fi>gi+delta > then False > else eq' (delta/2) (i+1) f g

f i

g i

f

g

Ord

compare

> instance Ord R where > compare (R f) (R g) = compare' 2 0 f g where > compare' delta i f g = > let fi = f i > gi = g i > in if fi<gi-delta > then LT > else if fi>gi+delta > then GT > else compare' (delta/2) (i+1) f g

min

max

> min (R f) (R g) = R $ \i -> min (f i) (g i) > max (R f) (R g) = R $ \i -> max (f i) (g i)

> instance Num R where > fromInteger n = R $ const $ fromInteger n > R f + R g = R $ \i -> let j = i+1 in f j + g j > negate (R f) = R $ negate . f

> x .* R f = R $ mul x x f where > mul x x' f i | abs x' < 1 = x * f i > | otherwise = mul x (x'/2) f (i+1)

ivt1

> ivt1 :: (R -> R) -> R > ivt1 f = R $ \i -> ivt' 0 1 i where > ivt' :: Rational -> Rational -> Integer -> Rational > ivt' x y i = > let z = (1%2)*(x+y) > in if i==0 > then z > else if f (fromRational z) < 0 > then ivt' z y (i-1) > else ivt' x z (i-1)

> f, g :: R -> R > f x = 3.*x-1 > g x = 2.*x-1

R -> R

ivt1 f

ivt1 g

ivt1

z

f z

f

R -> R

f

f

i

f

R -> R

ivt

ivt

g

g

Rational -> Rational -> Rational

g x y

x

y

ivt2

> ivt2 :: (Rational -> Rational -> Rational) -> (R -> R) -> R > ivt2 g f = R $ \i -> > let delta = (1%2)^i > ivt' :: Rational -> Rational -> Rational > ivt' x y = > let l = (1%3)*(2*x+y) > u = (1%3)*(x+2*y) > z = g l u > in if u-l < delta > then z > else if f (fromRational z) < 0 > then ivt' z y > else ivt' x z > in ivt' 0 1

pickNonZero

True

False

> pickNonZero (R f) (R g) = pickNonZero' 0 1 where > pickNonZero' i e = let fi = f i > gi = g i > in if fi > e || fi < -e > then True > else if gi > e || gi < -e > then False > else pickNonZero' (i+1) ((1%2)*e)

par

f

> monotonicPicker :: (R -> R) -> Rational -> Rational -> Rational > monotonicPicker f x y = if pickNonZero > (f (fromRational x)) (f (fromRational y)) > then x > else y > ivt3 :: (R -> R) -> R > ivt3 f = ivt2 (monotonicPicker f) f

monotonicPicker

ivt3 g











The Intermediate Value Theorem , as it is usually told, tells us that if we have a continuous real-valued function f on the closed interval [0,1], with f(0)<0 and f(1)>0 then there is a point x in (0,1) such that f(x)=0. Here's a sketch of a proof:Define three sequences l, u, xrecursively as follows:Set l=0, u=1.The xare bracketed by the land u, and |u-l|=2so the xform a Cauchy sequence with some limit x. Because f is continuous, we get f(x)=0.This proof not only shows that we have an intermediate value x where f crosses the x-axis, it also seems to describe a procedure for computing x. Let's use this strategy to write a numerical method.Here are some simple functions we can use to test it:andgive the expected results.Of course, it might not actually be the case that. We're dealing with the typehere and there might not actually be a solution toto the precision of a. And it's not clear that the word continuous means anything here anyway because the limited precision of aprecludes us from constructing the deltas and epsilons the definition of continuity demands.So can we do better? Can we implement a function to compute intermediate values for exact real arithmetic? I recently talked about computable real numbers, and they are what I intend to use here. I'll represent the real number x using a sequence of rationals, x, such that |x-x|<2. Here's a suitable type:(I'll only be calling it for non-negative arguments, but alas, Haskell has no natural number type.)We'll need to display our numbers. Here's something cheap and cheerful that's good enough for this article:It displays athat is formed from a rational approximation that is within 2of our value.We can inject the rationals into the computable reals by using constant sequences:An important point here is that two different sequences might define the same real. So when we construct functions mapping reals to reals we must be sure that they take different representations of the same real to representations of the same real. I'll call these 'well-defined'.Now our problems begin. Haskell wants us to implementbefore, but equality of computable reals is not decidable . Still, it is semi-decidable in that if two reals differ, we can detect this. Here's a partial equality function:The idea is that ifandcan be separated by more than 2*2then the reals thatandrepresent can't be equal.We have to implementtoo. If one computable real is bigger than another then they must differ by more than some power of 1/2. This means we can eventually distinguish between them. As a result,will always terminate if we compare distinct numbers:We can defineandstraightforwardly:And now we can define addition. We use the fact that if |x-x'|here And you'll find thatgives the result you expect. But sadly,never returns. What went wrong?After one iteration,setsto the value we want. This seems like a particularly easy case. But the algorithm immediately fails because it then proceeds to comparewith zero. We already know that our equality test must fail for this case. How do we fix our algorithm?Consider the actions of any higher order function that operates on a continuous functiontype. It will evaluateat various real values. The return values fromwill then be sampled at various values of. Each sample corresponds to a bracket around a value of. Consider a higher order function acting on the function below:If our higher order function terminates, it can only compute a finite number of these brackets. So I can choose an ε small enough that I could nudge the graph up or down by ε without invalidating the brackets. I could provide the higher order function with a different functionthat responds in exactly the same way to the finite number of samples.Now consider the case of anfunction. It will respond the same way to two functions differing only by ε. The function in the graph above has been designed so that it is flat from 0.4 to 0.6. Nudging the function above by ε>0 up and down causes the crossing of the x-axis to move from x<0.4 to x>0.6. In other words, no algorithm could possibly return a valid intermediate value for all possible inputs because I can always contrive a pair of functions, to whichwould respond identically, and yet which have widely separated points where they cross the x-axis. This shows that the problem of determining a crossing point of a continuous function that starts off less than zero, and ends up greater than zero, is uncomputable. We can do it for some functions. But no function can solve this problem for all continuous inputs. What's really weird about this is that it fails precisely when the problem seems easiest - when there are lots of zeroes!But we don't need to give up. If we could ensure that at each bisection we could pick a point that wasn't exactly a zero, then we'd be fine. But how can we do this? Well here's a cheating way: we don't bother. We let that problem be the caller's responsibility. Instead of bisecting at each stage, we'll trisect. We'll pick a point that isn't a zero in the middle third by using a function pased in by the caller. Call this functionis of typeand the specification is thatmust be a rational number betweenandthat isn't a zero. Here's our newfunction:How can we make use of this? Well suppose we want to find where a polynomial, f, crosses the x-axis. If the polynomial has degree n, it can only cross at n points. So if we pick n+1 values, f must evaluate to a non-zero value at one of them. We can examine the n+1 values of f(x) in parallel to find one that is non-zero. This must terminate. Here's an implementation for the case n+1=2. If we pass two real values toit willoraccording to whether the first or second argument is non-zero. It must only be called with a pair of numbers where one is non-zero:Note how it examines the two values in parallel. If it simply tested the two numbers in turn against zero it would fail if the first were zero. (Interesting side note: we could write this using.)Suppose the functionis strictly monotonically increasing. Then we know that if y>x, then f(y)>f(x) so we can't have that both f(x) and f(y) are zero. We can use this to write a suitable helper function to write an intermediate value function that is 100% guaranteed to work for all monotonically increasing functions (assuming we don't run out of memory):By varyingwe can find crossing points for all manner of continuous function. Try. But no matter how hard we try, we can never write one that works for all continuous functions.So why did the Intermediate Value Theorem, despite its appearance of constructing an intermediate value, fail to yield an algorithm? For an answer to that, you'll have to wait until my next article. The surprising thing is that it's not a problem of analysis, but a problem of logic.I have borrowed liberally from Paul Taylor , especially from here . However, I've chosen to phrase everything in terms of classical reasoning about Haskell programs. Any errors are, of course, mine. It was probably Andrej Bauer 's awesome blog that got me interested in this stuff.I want to stress that I have made absolutely no mention of intuitionism or constructivism in what I've said above. I have simply used ordinary classical mathematics to reason about some Haskell programs. I'm leaving that for the next article.