1----------------------------------------------------------------------------- 2-- | 3-- Module : Data.SBV.Core.AlgReals 4-- Copyright : (c) Levent Erkok 5-- License : BSD3 6-- Maintainer: erkokl@gmail.com 7-- Stability : experimental 8-- 9-- Algebraic reals in Haskell. 10----------------------------------------------------------------------------- 11 12{-# LANGUAGE FlexibleInstances #-} 13 14{-# OPTIONS_GHC -Wall -Werror -fno-warn-orphans #-} 15 16module Data.SBV.Core.AlgReals ( 17 AlgReal(..) 18 , AlgRealPoly(..) 19 , RationalCV(..) 20 , RealPoint(..), realPoint 21 , mkPolyReal 22 , algRealToSMTLib2 23 , algRealToHaskell 24 , algRealToRational 25 , mergeAlgReals 26 , isExactRational 27 , algRealStructuralEqual 28 , algRealStructuralCompare 29 ) 30 where 31 32import Data.Char (isDigit) 33 34import Data.List (sortBy, isPrefixOf, partition) 35import Data.Ratio ((%), numerator, denominator) 36import Data.Function (on) 37import System.Random 38import Test.QuickCheck (Arbitrary(..)) 39 40import Numeric (readSigned, readFloat) 41 42-- | Is the endpoint included in the interval? 43data RealPoint a = OpenPoint a -- ^ open: i.e., doesn't include the point 44 | ClosedPoint a -- ^ closed: i.e., includes the point 45 deriving (Show, Eq, Ord) 46 47-- | Extract the point associated with the open-closed point 48realPoint :: RealPoint a -> a 49realPoint (OpenPoint a) = a 50realPoint (ClosedPoint a) = a 51 52-- | Algebraic reals. Note that the representation is left abstract. We represent 53-- rational results explicitly, while the roots-of-polynomials are represented 54-- implicitly by their defining equation 55data AlgReal = AlgRational Bool Rational -- ^ bool says it's exact (i.e., SMT-solver did not return it with ? at the end.) 56 | AlgPolyRoot (Integer, AlgRealPoly) (Maybe String) -- ^ which root of this polynomial and an approximate decimal representation with given precision, if available 57 | AlgInterval (RealPoint Rational) (RealPoint Rational) -- ^ interval, with low and high bounds 58 59-- | Check whether a given argument is an exact rational 60isExactRational :: AlgReal -> Bool 61isExactRational (AlgRational True _) = True 62isExactRational _ = False 63 64-- | A univariate polynomial, represented simply as a 65-- coefficient list. For instance, "5x^3 + 2x - 5" is 66-- represented as [(5, 3), (2, 1), (-5, 0)] 67newtype AlgRealPoly = AlgRealPoly [(Integer, Integer)] 68 deriving (Eq, Ord) 69 70-- | Construct a poly-root real with a given approximate value (either as a decimal, or polynomial-root) 71mkPolyReal :: Either (Bool, String) (Integer, [(Integer, Integer)]) -> AlgReal 72mkPolyReal (Left (exact, str)) 73 = case (str, break (== '.') str) of 74 ("", _) -> AlgRational exact 0 75 (_, (x, '.':y)) | all isDigit (x ++ y) -> AlgRational exact (read (x++y) % (10 ^ length y)) 76 (_, (x, "")) | all isDigit x -> AlgRational exact (read x % 1) 77 _ -> error $ unlines [ "*** Data.SBV.mkPolyReal: Unable to read a number from:" 78 , "***" 79 , "*** " ++ str 80 , "***" 81 , "*** Please report this as a bug." 82 ] 83mkPolyReal (Right (k, coeffs)) 84 = AlgPolyRoot (k, AlgRealPoly (normalize coeffs)) Nothing 85 where normalize :: [(Integer, Integer)] -> [(Integer, Integer)] 86 normalize = merge . sortBy (flip compare `on` snd) 87 merge [] = [] 88 merge [x] = [x] 89 merge ((a, b):r@((c, d):xs)) 90 | b == d = merge ((a+c, b):xs) 91 | True = (a, b) : merge r 92 93instance Show AlgRealPoly where 94 show (AlgRealPoly xs) = chkEmpty (join (concat [term p | p@(_, x) <- xs, x /= 0])) ++ " = " ++ show c 95 where c = -1 * head ([k | (k, 0) <- xs] ++ [0]) 96 term ( 0, _) = [] 97 term ( 1, 1) = [ "x"] 98 term ( 1, p) = [ "x^" ++ show p] 99 term (-1, 1) = ["-x"] 100 term (-1, p) = ["-x^" ++ show p] 101 term (k, 1) = [show k ++ "x"] 102 term (k, p) = [show k ++ "x^" ++ show p] 103 join [] = "" 104 join (k:ks) = k ++ s ++ join ks 105 where s = case ks of 106 [] -> "" 107 (y:_) | "-" `isPrefixOf` y -> "" 108 | "+" `isPrefixOf` y -> "" 109 | True -> "+" 110 chkEmpty s = if null s then "0" else s 111 112instance Show AlgReal where 113 show (AlgRational exact a) = showRat exact a 114 show (AlgPolyRoot (i, p) mbApprox) = "root(" ++ show i ++ ", " ++ show p ++ ")" ++ maybe "" app mbApprox 115 where app v | last v == '?' = " = " ++ init v ++ "..." 116 | True = " = " ++ v 117 show (AlgInterval a b) = case (a, b) of 118 (OpenPoint l, OpenPoint h) -> "(" ++ show l ++ ", " ++ show h ++ ")" 119 (OpenPoint l, ClosedPoint h) -> "(" ++ show l ++ ", " ++ show h ++ "]" 120 (ClosedPoint l, OpenPoint h) -> "[" ++ show l ++ ", " ++ show h ++ ")" 121 (ClosedPoint l, ClosedPoint h) -> "[" ++ show l ++ ", " ++ show h ++ "]" 122 123-- lift unary op through an exact rational, otherwise bail 124lift1 :: String -> (Rational -> Rational) -> AlgReal -> AlgReal 125lift1 _ o (AlgRational e a) = AlgRational e (o a) 126lift1 nm _ a = error $ "AlgReal." ++ nm ++ ": unsupported argument: " ++ show a 127 128-- lift binary op through exact rationals, otherwise bail 129lift2 :: String -> (Rational -> Rational -> Rational) -> AlgReal -> AlgReal -> AlgReal 130lift2 _ o (AlgRational True a) (AlgRational True b) = AlgRational True (a `o` b) 131lift2 nm _ a b = error $ "AlgReal." ++ nm ++ ": unsupported arguments: " ++ show (a, b) 132 133-- The idea in the instances below is that we will fully support operations 134-- on "AlgRational" AlgReals, but leave everything else undefined. When we are 135-- on the Haskell side, the AlgReal's are *not* reachable. They only represent 136-- return values from SMT solvers, which we should *not* need to manipulate. 137instance Eq AlgReal where 138 AlgRational True a == AlgRational True b = a == b 139 a == b = error $ "AlgReal.==: unsupported arguments: " ++ show (a, b) 140 141instance Ord AlgReal where 142 AlgRational True a `compare` AlgRational True b = a `compare` b 143 a `compare` b = error $ "AlgReal.compare: unsupported arguments: " ++ show (a, b) 144 145-- | Structural equality for AlgReal; used when constants are Map keys 146algRealStructuralEqual :: AlgReal -> AlgReal -> Bool 147AlgRational a b `algRealStructuralEqual` AlgRational c d = (a, b) == (c, d) 148AlgPolyRoot a b `algRealStructuralEqual` AlgPolyRoot c d = (a, b) == (c, d) 149_ `algRealStructuralEqual` _ = False 150 151-- | Structural comparisons for AlgReal; used when constants are Map keys 152algRealStructuralCompare :: AlgReal -> AlgReal -> Ordering 153AlgRational a b `algRealStructuralCompare` AlgRational c d = (a, b) `compare` (c, d) 154AlgRational _ _ `algRealStructuralCompare` AlgPolyRoot _ _ = LT 155AlgRational _ _ `algRealStructuralCompare` AlgInterval _ _ = LT 156AlgPolyRoot _ _ `algRealStructuralCompare` AlgRational _ _ = GT 157AlgPolyRoot a b `algRealStructuralCompare` AlgPolyRoot c d = (a, b) `compare` (c, d) 158AlgPolyRoot _ _ `algRealStructuralCompare` AlgInterval _ _ = LT 159AlgInterval _ _ `algRealStructuralCompare` AlgRational _ _ = GT 160AlgInterval _ _ `algRealStructuralCompare` AlgPolyRoot _ _ = GT 161AlgInterval a b `algRealStructuralCompare` AlgInterval c d = (a, b) `compare` (c, d) 162 163instance Num AlgReal where 164 (+) = lift2 "+" (+) 165 (*) = lift2 "*" (*) 166 (-) = lift2 "-" (-) 167 negate = lift1 "negate" negate 168 abs = lift1 "abs" abs 169 signum = lift1 "signum" signum 170 fromInteger = AlgRational True . fromInteger 171 172-- | NB: Following the other types we have, we require `a/0` to be `0` for all a. 173instance Fractional AlgReal where 174 (AlgRational True _) / (AlgRational True b) | b == 0 = 0 175 a / b = lift2 "/" (/) a b 176 fromRational = AlgRational True 177 178instance Real AlgReal where 179 toRational (AlgRational True v) = v 180 toRational x = error $ "AlgReal.toRational: Argument cannot be represented as a rational value: " ++ algRealToHaskell x 181 182instance Random Rational where 183 random g = (a % b', g'') 184 where (a, g') = random g 185 (b, g'') = random g' 186 b' = if 0 < b then b else 1 - b -- ensures 0 < b 187 188 randomR (l, h) g = (r * d + l, g'') 189 where (b, g') = random g 190 b' = if 0 < b then b else 1 - b -- ensures 0 < b 191 (a, g'') = randomR (0, b') g' 192 193 r = a % b' 194 d = h - l 195 196instance Random AlgReal where 197 random g = let (a, g') = random g in (AlgRational True a, g') 198 randomR (AlgRational True l, AlgRational True h) g = let (a, g') = randomR (l, h) g in (AlgRational True a, g') 199 randomR lh _ = error $ "AlgReal.randomR: unsupported bounds: " ++ show lh 200 201-- | Render an 'AlgReal' as an SMTLib2 value. Only supports rationals for the time being. 202algRealToSMTLib2 :: AlgReal -> String 203algRealToSMTLib2 (AlgRational True r) 204 | m == 0 = "0.0" 205 | m < 0 = "(- (/ " ++ show (abs m) ++ ".0 " ++ show n ++ ".0))" 206 | True = "(/ " ++ show m ++ ".0 " ++ show n ++ ".0)" 207 where (m, n) = (numerator r, denominator r) 208algRealToSMTLib2 r@(AlgRational False _) 209 = error $ "SBV: Unexpected inexact rational to be converted to SMTLib2: " ++ show r 210algRealToSMTLib2 (AlgPolyRoot (i, AlgRealPoly xs) _) = "(root-obj (+ " ++ unwords (concatMap term xs) ++ ") " ++ show i ++ ")" 211 where term (0, _) = [] 212 term (k, 0) = [coeff k] 213 term (1, 1) = ["x"] 214 term (1, p) = ["(^ x " ++ show p ++ ")"] 215 term (k, 1) = ["(* " ++ coeff k ++ " x)"] 216 term (k, p) = ["(* " ++ coeff k ++ " (^ x " ++ show p ++ "))"] 217 coeff n | n < 0 = "(- " ++ show (abs n) ++ ")" 218 | True = show n 219algRealToSMTLib2 r@AlgInterval{} 220 = error $ "SBV: Unexpected inexact rational to be converted to SMTLib2: " ++ show r 221 222-- | Render an 'AlgReal' as a Haskell value. Only supports rationals, since there is no corresponding 223-- standard Haskell type that can represent root-of-polynomial variety. 224algRealToHaskell :: AlgReal -> String 225algRealToHaskell (AlgRational True r) = "((" ++ show r ++ ") :: Rational)" 226algRealToHaskell r = error $ unlines [ "" 227 , "SBV.algRealToHaskell: Unsupported argument:" 228 , "" 229 , " " ++ show r 230 , "" 231 , "represents an irrational number, and cannot be converted to a Haskell value." 232 ] 233 234-- | Conversion from internal rationals to Haskell values 235data RationalCV = RatIrreducible AlgReal -- ^ Root of a polynomial, cannot be reduced 236 | RatExact Rational -- ^ An exact rational 237 | RatApprox Rational -- ^ An approximated value 238 | RatInterval (RealPoint Rational) (RealPoint Rational) -- ^ Interval. Can be open/closed on both ends. 239 deriving Show 240 241-- | Convert an 'AlgReal' to a 'Rational'. If the 'AlgReal' is exact, then you get a 'Left' value. Otherwise, 242-- you get a 'Right' value which is simply an approximation. 243algRealToRational :: AlgReal -> RationalCV 244algRealToRational a = case a of 245 AlgRational True r -> RatExact r 246 AlgRational False r -> RatExact r 247 AlgPolyRoot _ Nothing -> RatIrreducible a 248 AlgPolyRoot _ (Just s) -> let trimmed = case reverse s of 249 '.':'.':'.':rest -> reverse rest 250 _ -> s 251 in case readSigned readFloat trimmed of 252 [(v, "")] -> RatApprox v 253 _ -> bad "represents a value that cannot be converted to a rational" 254 AlgInterval lo hi -> RatInterval lo hi 255 where bad w = error $ unlines [ "" 256 , "SBV.algRealToRational: Unsupported argument:" 257 , "" 258 , " " ++ show a 259 , "" 260 , w 261 ] 262 263-- Try to show a rational precisely if we can, with finite number of 264-- digits. Otherwise, show it as a rational value. 265showRat :: Bool -> Rational -> String 266showRat exact r = p $ case f25 (denominator r) [] of 267 Nothing -> show r -- bail out, not precisely representable with finite digits 268 Just (noOfZeros, num) -> let present = length num 269 in neg $ case noOfZeros `compare` present of 270 LT -> let (b, a) = splitAt (present - noOfZeros) num in b ++ "." ++ if null a then "0" else a 271 EQ -> "0." ++ num 272 GT -> "0." ++ replicate (noOfZeros - present) '0' ++ num 273 where p = if exact then id else (++ "...") 274 neg = if r < 0 then ('-':) else id 275 -- factor a number in 2's and 5's if possible 276 -- If so, it'll return the number of digits after the zero 277 -- to reach the next power of 10, and the numerator value scaled 278 -- appropriately and shown as a string 279 f25 :: Integer -> [Integer] -> Maybe (Int, String) 280 f25 1 sofar = let (ts, fs) = partition (== 2) sofar 281 [lts, lfs] = map length [ts, fs] 282 noOfZeros = lts `max` lfs 283 in Just (noOfZeros, show (abs (numerator r) * factor ts fs)) 284 f25 v sofar = let (q2, r2) = v `quotRem` 2 285 (q5, r5) = v `quotRem` 5 286 in case (r2, r5) of 287 (0, _) -> f25 q2 (2 : sofar) 288 (_, 0) -> f25 q5 (5 : sofar) 289 _ -> Nothing 290 -- compute the next power of 10 we need to get to 291 factor [] fs = product [2 | _ <- fs] 292 factor ts [] = product [5 | _ <- ts] 293 factor (_:ts) (_:fs) = factor ts fs 294 295-- | Merge the representation of two algebraic reals, one assumed to be 296-- in polynomial form, the other in decimal. Arguments can be the same 297-- kind, so long as they are both rationals and equivalent; if not there 298-- must be one that is precise. It's an error to pass anything 299-- else to this function! (Used in reconstructing SMT counter-example values with reals). 300mergeAlgReals :: String -> AlgReal -> AlgReal -> AlgReal 301mergeAlgReals _ f@(AlgRational exact r) (AlgPolyRoot kp Nothing) 302 | exact = f 303 | True = AlgPolyRoot kp (Just (showRat False r)) 304mergeAlgReals _ (AlgPolyRoot kp Nothing) f@(AlgRational exact r) 305 | exact = f 306 | True = AlgPolyRoot kp (Just (showRat False r)) 307mergeAlgReals _ f@(AlgRational e1 r1) s@(AlgRational e2 r2) 308 | (e1, r1) == (e2, r2) = f 309 | e1 = f 310 | e2 = s 311mergeAlgReals m _ _ = error m 312 313-- Quickcheck instance 314instance Arbitrary AlgReal where 315 arbitrary = AlgRational True `fmap` arbitrary 316