```{-
-      ``Data/Random/Internal/Find''
-  Utilities for searching fractional domains.  Needs cleanup, testing,
-  and such.  Used for constructing generic ziggurats.
-}

module Data.Random.Internal.Find where

findMax :: (Fractional a, Ord a) => (a -> Bool) -> a
findMax p = negate (findMin (p.negate))

-- |Given an upward-closed predicate on an ordered Fractional type,
-- find the smallest value satisfying the predicate.
findMin :: (Fractional a, Ord a) => (a -> Bool) -> a
findMin = findMinFrom 0 1

-- |Given an upward-closed predicate on an ordered Fractional type,
-- find the smallest value satisfying the predicate.  Starts at the
-- specified point with the specified stepsize, performs an exponential
-- search out from there until it finds an interval bracketing the
-- change-point of the predicate, and then performs a bisection search
-- to isolate the change point.  Note that infinitely-divisible domains
-- such as 'Rational' cannot be searched by this function because it does
-- not terminate until it reaches a point where further subdivision of the
-- interval has no effect.
findMinFrom :: (Fractional a, Ord a) => a -> a -> (a -> Bool) -> a
findMinFrom z0 0 p = findMinFrom z0 1 p
findMinFrom z0 step1 p
| p z0      = descend (z0-step1) z0
| otherwise = fixZero (ascend z0 (z0+step1))
where
-- eliminate negative zero, which, in many domains, is technically
fixZero 0 = 0
fixZero z = z

-- preconditions:
-- not (p l)
-- 0 <= l < x
ascend l x
| p x       = bisect l x
| otherwise = ascend x \$! 2*x-z0

-- preconditions:
-- p h
-- x < h <= 0
descend x h
| p x       = (descend \$! 2*x-z0) x
| otherwise = bisect x h

-- preconditions:
-- not (p l)
-- p h
-- l <= h
bisect l h
| l /< h    = h
| l /< mid || mid /< h
= if p mid then mid else h
| p mid     = bisect l mid
| otherwise = bisect mid h
where
a /< b = not (a < b)
mid = (l+h)*0.5
```