0 | module Data.Double.Bounded
  1 |
  2 | import Data.DPair
  3 |
  4 | import public Data.So
  5 |
  6 | %default total
  7 |
  8 | --- Auxiliary functions ---
  9 |
 10 | public export %inline
 11 | OR : Type -> Type -> Type
 12 | OR = Either
 13 |
 14 | export infixr 0 `OR`
 15 |
 16 | namespace DoubleUtils
 17 |
 18 |   public export %inline
 19 |   MaxDouble, MinDouble, PosInf, NegInf : Double
 20 |   MaxDouble = 1.79769e+308
 21 |   MinDouble = -MaxDouble
 22 |   PosInf = 1/0
 23 |   NegInf = -1/0
 24 |
 25 |   public export %inline
 26 |   min4 : Double -> Double -> Double -> Double -> Double
 27 |   min4 a b c d = a `min` (b `min` (c `min` d))
 28 |
 29 |   public export %inline
 30 |   min6 : Double -> Double -> Double -> Double -> Double -> Double -> Double
 31 |   min6 a b c d e f = a `min` (b `min` min4 c d e f)
 32 |
 33 |   public export %inline
 34 |   max4 : Double -> Double -> Double -> Double -> Double
 35 |   max4 a b c d = a `max` (b `max` (c `max` d))
 36 |
 37 |   public export %inline
 38 |   max6 : Double -> Double -> Double -> Double -> Double -> Double -> Double
 39 |   max6 a b c d e f = a `max` (b `max` max4 c d e f)
 40 |
 41 |   -- Returns zero if `l <= 0 <= u`, or `l` or `u` otherwise
 42 |   public export
 43 |   zormin : (l, u : Double) -> Double
 44 |   zormin l u = max l (0 `min` u)
 45 |
 46 |   public export %inline
 47 |   Finite : Double -> Type
 48 |   Finite x = So (x /= NegInf && x /= PosInf)
 49 |
 50 |   public export %inline
 51 |   NonZero : Double -> Type
 52 |   NonZero x = So (x /= 0)
 53 |
 54 | namespace DoubleProperties
 55 |
 56 |   export
 57 |   lteRefl : {0 x : Double} -> (0 _ : So $ x == x) => So $ x <= x
 58 |   lteRefl = believe_me Oh
 59 |
 60 |   export
 61 |   lteNotNaNL : {0 x, y : Double} -> (0 _ : So $ x <= y) => So $ x == x
 62 |   lteNotNaNL = believe_me Oh
 63 |
 64 |   export
 65 |   lteNotNaNR : {0 x, y : Double} -> (0 _ : So $ y <= x) => So $ x == x
 66 |   lteNotNaNR = believe_me Oh
 67 |
 68 |   export
 69 |   lteTrans : {0 a, b, c : Double} -> (0 _ : So $ a <= b) => (0 _ : So $ b <= c) => So $ a <= c
 70 |   lteTrans = believe_me Oh
 71 |
 72 |   export
 73 |   lteNegInf : {0 x : Double} -> (0 _ : So $ x == x) => So $ NegInf <= x
 74 |   lteNegInf = believe_me Oh
 75 |
 76 |   export
 77 |   ltePosInf : {0 x : Double} -> (0 _ : So $ x == x) => So $ x <= PosInf
 78 |   ltePosInf = believe_me Oh
 79 |
 80 |   export
 81 |   lteMin : {0 x : Double} -> (0 _ : So $ x == x) => (0 _ : Finite x) => So $ MinDouble <= x
 82 |   lteMin = believe_me Oh
 83 |
 84 |   export
 85 |   lteMax : {0 x : Double} -> (0 _ : So $ x == x) => (0 _ : Finite x) => So $ x <= MaxDouble
 86 |   lteMax = believe_me Oh
 87 |
 88 |   export
 89 |   lteFromLt : {0 x, y : Double} -> (0 _ : So $ x < y) => So $ x <= y
 90 |   lteFromLt = believe_me Oh
 91 |
 92 |   export
 93 |   lteRev : {0 x, y : Double} -> (0 _ : So $ x == x) => (0 _ : So $ y == y) => (0 _ : So $ not $ x <= y) => So $ y < x
 94 |   lteRev = believe_me Oh
 95 |
 96 | --- Type definitions ---
 97 |
 98 | ||| The type of double in the given bounds.
 99 | |||
100 | ||| Both bounds and the value can be infinite.
101 | ||| If either bounds is `NaN`, the type is uninhabited.
102 | ||| The internal double value cannot be `NaN` since to inhabit the constructor,
103 | ||| one needs to satisfy the bounds which is impossible with `NaN` value.
104 | public export
105 | data DoubleBetween : (lower, upper : Double) -> Type where
106 |   BoundedDouble :
107 |     {0 lower, upper : Double} ->
108 |     (x : Double) ->
109 |     (0 _ : So $ lower <= x) =>
110 |     (0 _ : So $ x <= upper) =>
111 |     DoubleBetween lower upper
112 |
113 | ||| Represents a bounded double which can hold any `Double` except the `NaN`.
114 | ||| Actually, an alias for `DoubleBetween` with infinite bounds.
115 | public export %inline
116 | SolidDouble : Type
117 | SolidDouble = DoubleBetween NegInf PosInf
118 |
119 | public export %inline
120 | FiniteDouble : Type
121 | FiniteDouble = DoubleBetween MinDouble MaxDouble
122 |
123 | --- Conversions ---
124 |
125 | public export
126 | maybeBoundedDouble : {lower, upper : Double} -> (x : Double) -> Maybe $ DoubleBetween lower upper
127 | maybeBoundedDouble x with (lower <= x) proof lx
128 |   _ | False = Nothing
129 |   _ | True with (x <= upper) proof xu
130 |     _ | False = Nothing
131 |     _ | True  = Just $ BoundedDouble x @{eqToSo lx} @{eqToSo xu}
132 |
133 | public export
134 | roughlyFit : {l, u : _} -> (0 _ : So $ l <= u) => (x : Double) -> DoubleBetween l u
135 | roughlyFit @{lu} x with (x <= u) proof xu
136 |   _ | False = BoundedDouble u @{lu} @{lteRefl @{lteNotNaNR @{lu}}}
137 |   _ | True with (l <= x) proof lx
138 |     _ | False = BoundedDouble l @{lteRefl @{lteNotNaNL @{lu}}}
139 |     _ | True  = BoundedDouble x @{eqToSo lx} @{eqToSo xu}
140 |
141 | public export
142 | Cast (DoubleBetween l u) Double where
143 |   cast $ BoundedDouble x = x
144 |
145 | public export %inline
146 | (.asDouble) : DoubleBetween l u -> Double
147 | (.asDouble) = cast
148 |
149 | public export
150 | weakenBounds : DoubleBetween l u -> (0 _ : So $ l' <= l) => (0 _ : So $ u <= u') => DoubleBetween l' u'
151 | weakenBounds $ BoundedDouble x = BoundedDouble x @{lteTrans {b=l}} @{lteTrans {b=u}}
152 |
153 | public export %inline
154 | relaxToSolid : DoubleBetween l u -> SolidDouble
155 | relaxToSolid x@(BoundedDouble _ @{lb} @{rb}) = weakenBounds x @{lteNegInf @{lteNotNaNL @{lb}}} @{ltePosInf @{lteNotNaNR @{rb}}}
156 |
157 | public export %inline
158 | relaxToFinite : (0 _ : Finite l) => (0 _ : Finite u) => DoubleBetween l u -> FiniteDouble
159 | relaxToFinite x@(BoundedDouble _ @{lb} @{rb}) = weakenBounds x @{lteMin @{lteNotNaNL @{lb}}} @{lteMax @{lteNotNaNR @{rb}}}
160 |
161 | public export
162 | strengthenBounds : {l', u' : _} -> DoubleBetween l u -> Maybe $ DoubleBetween l' u'
163 | strengthenBounds = maybeBoundedDouble . cast
164 |
165 | export
166 | Cast Nat (DoubleBetween 0 PosInf) where
167 |   cast = roughlyFit . cast
168 |
169 | export %inline
170 | fromNat : Nat -> DoubleBetween 0 PosInf
171 | fromNat = cast
172 |
173 | export
174 | Cast Integer SolidDouble where
175 |   cast = roughlyFit . cast
176 |
177 | --- Literals syntax ---
178 |
179 | namespace MinimalBounds
180 |
181 |   public export %inline
182 |   fromDouble : (x : Double) -> (0 _ : So $ x == x) => DoubleBetween x x
183 |   fromDouble x = BoundedDouble x @{lteRefl} @{lteRefl}
184 |
185 |   public export %inline
186 |   fromInteger : (x : Integer) -> (0 _ : So $ cast {to=Double} x == cast x) => DoubleBetween (cast x) (cast x)
187 |   fromInteger x = BoundedDouble (cast x) @{lteRefl} @{lteRefl}
188 |
189 | namespace KnownBounds
190 |
191 |   public export %inline
192 |   fromDouble : (x : Double) -> (0 _ : So $ lower <= x) => (0 _ : So $ x <= upper) => DoubleBetween lower upper
193 |   fromDouble x = BoundedDouble x
194 |
195 |   public export %inline
196 |   fromInteger : (x : Integer) -> (0 _ : So $ lower <= cast x) => (0 _ : So $ cast x <= upper) => DoubleBetween lower upper
197 |   fromInteger x = BoundedDouble $ cast x
198 |
199 | --- Comparison ---
200 |
201 | public export
202 | Eq (DoubleBetween l u) where
203 |   (==) = (==) `on` (.asDouble)
204 |   (/=) = (/=) `on` (.asDouble)
205 |
206 | public export
207 | Ord (DoubleBetween l u) where
208 |   compare = compare `on` (.asDouble)
209 |   (<)  = (<)  `on` (.asDouble)
210 |   (<=) = (<=) `on` (.asDouble)
211 |   (>)  = (>)  `on` (.asDouble)
212 |   (>=) = (>=) `on` (.asDouble)
213 |
214 | --- Printing ---
215 |
216 | public export
217 | Show (DoubleBetween l u) where
218 |   show $ BoundedDouble x = show x
219 |
220 | --- Important constants ---
221 |
222 | public export %inline
223 | NegInf : (0 _ : So $ u == u) => DoubleBetween NegInf u
224 | NegInf = BoundedDouble NegInf @{%search} @{lteNegInf}
225 |
226 | public export %inline
227 | PosInf : (0 _ : So $ l == l) => DoubleBetween l PosInf
228 | PosInf = BoundedDouble PosInf @{ltePosInf}
229 |
230 | public export %inline
231 | pi, Pi : DoubleBetween Prelude.pi Prelude.pi
232 | pi = MinimalBounds.fromDouble pi
233 | Pi = pi
234 |
235 | public export %inline
236 | euler, Euler : DoubleBetween Prelude.euler Prelude.euler
237 | euler = MinimalBounds.fromDouble euler
238 | Euler = euler
239 |
240 | --- Analysis operations ---
241 |
242 | export
243 | strengthenRight : (x, y : DoubleBetween l u) -> (0 _ : So $ x <= y) => Subset (DoubleBetween x.asDouble u) $ \y' => y'.asDouble = y.asDouble
244 | strengthenRight (BoundedDouble x) (BoundedDouble y) @{xy} = Element (BoundedDouble y) Refl
245 |
246 | export
247 | bisect : (m : DoubleBetween l u) -> DoubleBetween l u -> DoubleBetween l m.asDouble `OR` DoubleBetween m.asDouble u
248 | bisect (BoundedDouble m @{lm}) (BoundedDouble x @{lb}) with (x <= m) proof xm
249 |   _ | True  = Left $ BoundedDouble x @{%search} @{eqToSo xm}
250 |   _ | False = Right $ BoundedDouble x @{do
251 |                 let smx = eqToSo $ cong not xm
252 |                 let xx = lteNotNaNR @{lb}
253 |                 let mm = lteNotNaNR @{lm}
254 |                 lteFromLt @{lteRev}
255 |               }
256 |
257 | export
258 | trisect : (m1, m2 : DoubleBetween l u) ->
259 |           (0 _ : So $ m1 <= m2) =>
260 |           DoubleBetween l u ->
261 |           DoubleBetween l m1.asDouble `OR` DoubleBetween m1.asDouble m2.asDouble `OR` DoubleBetween m2.asDouble u
262 | trisect m1 m2 x = case bisect m1 x of
263 |                     Left l  => Left l
264 |                     Right r => Right $ do
265 |                       let Element m2' prf = strengthenRight m1 m2
266 |                       rewrite sym prf
267 |                       bisect m2' r
268 |
269 | --- Basic arithmetics ---
270 |
271 | export
272 | (+) : DoubleBetween l u ->
273 |       DoubleBetween l' u' ->
274 |       (0 _ : Finite l `OR` Finite l' `OR` So (l == l')) =>
275 |       (0 _ : Finite u `OR` Finite u' `OR` So (u == u')) =>
276 |       DoubleBetween (l+l') (u+u')
277 | BoundedDouble x + BoundedDouble y = BoundedDouble (x + y) @{believe_me Oh} @{believe_me Oh}
278 |
279 | export
280 | (-) : DoubleBetween l u ->
281 |       DoubleBetween l' u' ->
282 |       (0 _ : Finite l `OR` Finite u' `OR` So (l /= u')) =>
283 |       (0 _ : Finite u `OR` Finite l' `OR` So (u /= l')) =>
284 |       DoubleBetween (l-u') (u-l')
285 | BoundedDouble x - BoundedDouble y = BoundedDouble (x - y) @{believe_me Oh} @{believe_me Oh}
286 |
287 | export
288 | (*) : DoubleBetween l u ->
289 |       DoubleBetween l' u' ->
290 |       (0 _ : Finite l `OR` NonZero l') =>
291 |       (0 _ : Finite l `OR` NonZero u') =>
292 |       (0 _ : Finite u `OR` NonZero l') =>
293 |       (0 _ : Finite u `OR` NonZero u') =>
294 |       (0 _ : Finite l' `OR` NonZero l) =>
295 |       (0 _ : Finite l' `OR` NonZero u) =>
296 |       (0 _ : Finite u' `OR` NonZero l) =>
297 |       (0 _ : Finite u' `OR` NonZero u) =>
298 |       DoubleBetween (min4 (l*l') (l*u') (u*l') (u*u')) (max4 (l*l') (l*u') (u*l') (u*u'))
299 | BoundedDouble x * BoundedDouble y = BoundedDouble (x * y) @{believe_me Oh} @{believe_me Oh}
300 |
301 | export
302 | (/) : {l, u, l', u' : _} ->
303 |       (num : DoubleBetween l u) ->
304 |       (den : DoubleBetween l' u') ->
305 |       (0 _ : So (0 < l') `OR` So (u' < 0) `OR` So (l' < 0 && 0 < u' && l /= 0 && u /= 0)) =>
306 |       (0 _ : So (0 < l) `OR` So (u < 0) `OR` So (0 < l') `OR` So (u' < 0) `OR` (NonZero num.asDouble, NonZero den.asDouble)) =>
307 |       (0 _ : Finite l `OR` Finite l') =>
308 |       (0 _ : Finite l `OR` Finite u') =>
309 |       (0 _ : Finite u `OR` Finite l') =>
310 |       (0 _ : Finite u `OR` Finite u') =>
311 |       (0 _ : NonZero l `OR` NonZero l') =>
312 |       (0 _ : NonZero l `OR` NonZero u') =>
313 |       (0 _ : NonZero u `OR` NonZero l') =>
314 |       (0 _ : NonZero u `OR` NonZero u') =>
315 |       DoubleBetween
316 |         (min6 (l/l') (l/u') (u/l') (u/u') (l/zormin l' u') (u/zormin l' u'))
317 |         (max6 (l/l') (l/u') (u/l') (u/u') (l/zormin l' u') (u/zormin l' u'))
318 | BoundedDouble x / BoundedDouble y = roughlyFit @{believe_me Oh} (x / y)
319 |
320 | export
321 | negate : DoubleBetween l u -> DoubleBetween (-u) (-l)
322 | negate x = BoundedDouble (negate x.asDouble) @{believe_me Oh} @{believe_me Oh}
323 |
324 | --- Math functions ---
325 |
326 | export
327 | sqrt : (0 _ : So $ 0 <= l) => DoubleBetween l u -> DoubleBetween (sqrt l) (sqrt u)
328 | sqrt x = BoundedDouble (sqrt x.asDouble) @{believe_me Oh} @{believe_me Oh}
329 |
330 | export
331 | sqrtRelaxed : (0 _ : So $ 0 <= l) => DoubleBetween l u -> DoubleBetween 0 PosInf
332 | sqrtRelaxed x = BoundedDouble (sqrt x.asDouble) @{believe_me Oh} @{believe_me Oh}
333 |
334 | export
335 | log : (0 _ : So $ 0 <= l) => DoubleBetween l u -> DoubleBetween (log l) (log u)
336 | log x = BoundedDouble (log x.asDouble) @{believe_me Oh} @{believe_me Oh}
337 |
338 | export
339 | exp : DoubleBetween l u -> DoubleBetween (exp l) (exp u)
340 | exp x = BoundedDouble (exp x.asDouble) @{believe_me Oh} @{believe_me Oh}
341 |
342 | export
343 | expRelaxed : DoubleBetween l u -> DoubleBetween 0 PosInf
344 | expRelaxed x = BoundedDouble (exp x.asDouble) @{believe_me Oh} @{believe_me Oh}
345 |
346 | export
347 | pow2 : DoubleBetween l u -> DoubleBetween (pow 2 l) (pow 2 u)
348 | pow2 x = BoundedDouble (pow 2 x.asDouble) @{believe_me Oh} @{believe_me Oh}
349 |
350 | export
351 | sin : FiniteDouble -> DoubleBetween (-1) 1
352 | sin x = BoundedDouble (sin x.asDouble) @{believe_me Oh} @{believe_me Oh}
353 |
354 | export
355 | cos : FiniteDouble -> DoubleBetween (-1) 1
356 | cos x = BoundedDouble (cos x.asDouble) @{believe_me Oh} @{believe_me Oh}
357 |
358 | export
359 | asin : DoubleBetween (-1) 1 -> DoubleBetween (-Prelude.pi/2) (Prelude.pi/2)
360 | asin x = BoundedDouble (asin x.asDouble) @{believe_me Oh} @{believe_me Oh}
361 |
362 | export
363 | acos : DoubleBetween (-1) 1 -> DoubleBetween 0 Prelude.pi
364 | acos x = BoundedDouble (acos x.asDouble) @{believe_me Oh} @{believe_me Oh}
365 |