sbv-14.4: Documentation/SBV/Examples/TP/GCD.hs
-----------------------------------------------------------------------------
-- |
-- Module : Documentation.SBV.Examples.TP.GCD
-- Copyright : (c) Levent Erkok
-- License : BSD3
-- Maintainer: erkokl@gmail.com
-- Stability : experimental
--
-- We define three different versions of the GCD algorithm: (1) Regular
-- version using the modulus operator, (2) the more basic version using
-- subtraction, and (3) the so called binary GCD. We prove that the modulus
-- based algorithm correct, i.e., that it calculates the greatest-common-divisor
-- of its arguments. We then prove that the other two variants are equivalent
-- to this version, thus establishing their correctness as well.
-----------------------------------------------------------------------------
{-# LANGUAGE CPP #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE QuasiQuotes #-}
{-# LANGUAGE TypeAbstractions #-}
{-# LANGUAGE TypeApplications #-}
{-# OPTIONS_GHC -Wall -Werror #-}
module Documentation.SBV.Examples.TP.GCD where
import Prelude hiding (gcd)
import Data.SBV
import Data.SBV.TP
import Data.SBV.Tuple
#ifdef DOCTEST
-- $setup
-- >>> import Data.SBV
-- >>> import Data.SBV.TP
#endif
-- * Calculating GCD
-- | @nGCD@ is the version of GCD that works on non-negative integers.
--
-- Ideally, we should make this function local to @gcd@, but then we can't refer to it explicitly in our proofs.
--
-- Note on maximality: Note that, by definition @gcd 0 0 = 0@. Since any number divides @0@,
-- there is no greatest common divisor for the pair @(0, 0)@. So, maximality here is meant
-- to be in terms of divisibility. That is, any divisor of @a@ and @b@ will also divide their @gcd@.
nGCD :: SInteger -> SInteger -> SInteger
nGCD = smtFunction "nGCD" $ \a b -> [sCase| b of
_ | b .== 0 -> a
_ -> nGCD b (a `sEMod` b)
|]
-- | Generalized GCD, working for all integers. We simply call @nGCD@ with the absolute value of the arguments.
gcd :: SInteger -> SInteger -> SInteger
gcd a b = nGCD (abs a) (abs b)
-- * Basic properties
-- | \(\gcd\, a\ b \geq 0\)
--
-- ==== __Proof__
-- >>> runTP gcdNonNegative
-- Inductive lemma (strong): nonNegativeNGCD
-- Step: Measure is non-negative Q.E.D.
-- Step: 1 (2 way case split)
-- Step: 1.1 Q.E.D.
-- Step: 1.2.1 Q.E.D.
-- Step: 1.2.2 Q.E.D.
-- Step: 1.Completeness Q.E.D.
-- Result: Q.E.D.
-- Lemma: nonNegative Q.E.D.
-- Functions proven terminating: nGCD
-- [Proven] nonNegative :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
gcdNonNegative :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
gcdNonNegative = do
-- We first prove over nGCD, using strong induction with the measure @a+b@.
nn <- sInduct "nonNegativeNGCD"
(\(Forall a) (Forall b) -> a .>= 0 .&& b .>= 0 .=> nGCD a b .>= 0)
(\_a b -> b, []) $
\ih a b -> [a .>= 0, b .>= 0]
|- cases [ b .== 0 ==> trivial
, b ./= 0 ==> nGCD a b .>= 0
=: nGCD b (a `sEMod` b) .>= 0
?? ih `at` (Inst @"a" b, Inst @"b" (a `sEMod` b))
=: sTrue
=: qed
]
lemma "nonNegative"
(\(Forall a) (Forall b) -> gcd a b .>= 0)
[proofOf nn]
-- | \(\gcd\, a\ b=0\implies a=0\land b=0\)
--
-- ==== __Proof__
-- >>> runTP gcdZero
-- Inductive lemma (strong): nGCDZero
-- Step: Measure is non-negative Q.E.D.
-- Step: 1 (2 way case split)
-- Step: 1.1 Q.E.D.
-- Step: 1.2.1 Q.E.D.
-- Step: 1.2.2 Q.E.D.
-- Step: 1.Completeness Q.E.D.
-- Result: Q.E.D.
-- Lemma: gcdZero Q.E.D.
-- Functions proven terminating: nGCD
-- [Proven] gcdZero :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
gcdZero :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
gcdZero = do
-- First prove over nGCD:
nGCDZero <-
sInduct "nGCDZero"
(\(Forall @"a" a) (Forall @"b" b) -> a .>= 0 .&& b .>= 0 .&& nGCD a b .== 0 .=> a .== 0 .&& b .== 0)
(\_a b -> b, []) $
\ih a b -> [a .>= 0, b .>= 0]
|- (nGCD a b .== 0 .=> a .== 0 .&& b .== 0)
=: cases [ b .== 0 ==> trivial
, b .> 0 ==> (nGCD b (a `sEMod` b) .== 0 .=> a .== 0 .&& b .== 0)
?? ih `at` (Inst @"a" b, Inst @"b" (a `sEMod` b))
=: sTrue
=: qed
]
lemma "gcdZero"
(\(Forall @"a" a) (Forall @"b" b) -> gcd a b .== 0 .=> a .== 0 .&& b .== 0)
[proofOf nGCDZero]
-- | \(\gcd\, a\ b=\gcd\, b\ a\)
--
-- ==== __Proof__
-- >>> runTP commutative
-- Lemma: nGCDCommutative
-- Step: 1 Q.E.D.
-- Result: Q.E.D.
-- Lemma: commutative
-- Step: 1 Q.E.D.
-- Step: 2 Q.E.D.
-- Result: Q.E.D.
-- Functions proven terminating: nGCD
-- [Proven] commutative :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
commutative :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
commutative = do
-- First prove over nGCD. Simple enough proof, but quantifiers and recursive functions
-- cause z3 to diverge. So, we have to explicitly write it out.
nGCDComm <-
calc "nGCDCommutative"
(\(Forall @"a" a) (Forall @"b" b) -> a .>= 0 .&& b .>= 0 .=> nGCD a b .== nGCD b a) $
\a b -> [a .>= 0, b .>= 0]
|- nGCD a b
=: nGCD b a
=: qed
-- It's unfortunate we have to spell this out explicitly, a simple lemma call
-- that uses the above proof doesn't converge.
calc "commutative"
(\(Forall a) (Forall b) -> gcd a b .== gcd b a) $
\a b -> [] |- gcd a b
=: nGCD (abs a) (abs b)
?? nGCDComm `at` (Inst @"a" (abs a), Inst @"b" (abs b))
=: gcd b a
=: qed
-- | \(\gcd\,(-a)\,b = \gcd\,a\,b = \gcd\,a\,(-b)\)
--
-- ==== __Proof__
-- >>> runTP negGCD
-- Lemma: negGCD Q.E.D.
-- Functions proven terminating: nGCD
-- [Proven] negGCD :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
negGCD :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
negGCD = lemma "negGCD" (\(Forall a) (Forall b) -> let g = gcd a b in gcd (-a) b .== g .&& g .== gcd a (-b)) []
-- | \( \gcd\,a\,0 = \gcd\,0\,a = |a| \land \gcd\,0\,0 = 0\)
--
-- ==== __Proof__
-- >>> runTP zeroGCD
-- Lemma: zeroGCD Q.E.D.
-- Functions proven terminating: nGCD
-- [Proven] zeroGCD :: Ɐa ∷ Integer → Bool
zeroGCD :: TP (Proof (Forall "a" Integer -> SBool))
zeroGCD = lemma "zeroGCD" (\(Forall a) -> gcd a 0 .== gcd 0 a .&& gcd 0 a .== abs a .&& gcd 0 0 .== 0) []
-- * Even and odd
-- | Is the given integer even?
isEven :: SInteger -> SBool
isEven = (2 `sDivides`)
-- | Is the given integer odd?
isOdd :: SInteger -> SBool
isOdd = sNot . isEven
-- * Divisibility
-- | Divides relation. By definition @0@ only divides @0@. (But every number divides @0@).
dvd :: SInteger -> SInteger -> SBool
a `dvd` b = ite (a .== 0) (b .== 0) (b `sEMod` a .== 0)
-- | \(d \mid a \implies d \mid ka\)
--
-- ==== __Proof__
-- >>> runTP dvdMul
-- Lemma: dvdMul
-- Step: 1 (2 way case split)
-- Step: 1.1 Q.E.D.
-- Step: 1.2.1 Q.E.D.
-- Step: 1.2.2 Q.E.D.
-- Step: 1.Completeness Q.E.D.
-- Result: Q.E.D.
-- [Proven] dvdMul :: Ɐd ∷ Integer → Ɐa ∷ Integer → Ɐk ∷ Integer → Bool
dvdMul :: TP (Proof (Forall "d" Integer -> Forall "a" Integer -> Forall "k" Integer -> SBool))
dvdMul = calc "dvdMul"
(\(Forall d) (Forall a) (Forall k) -> d `dvd` a .=> d `dvd` (k*a)) $
\d a k -> [d `dvd` a]
|- cases [ d .== 0 ==> d `dvd` (k*a)
?? a .== 0
=: sTrue
=: qed
, d ./= 0 ==> d `dvd` (k*a)
=: (k*a) `sEMod` d .== 0
?? a .== d * a `sEDiv` d
?? k * a .== d * (k * a `sEDiv` d)
?? (d * (k * a `sEDiv` d)) `sEMod` d .== 0
=: sTrue
=: qed
]
-- | \(a \mid |b| \iff a \mid b\)
--
-- A number divides another exactly when it also divides its absolute value. This follows
-- from 'dvdMul', as both directions are an instance of multiplying by @-1@.
--
-- ==== __Proof__
-- >>> runTP dvdAbs
-- Lemma: dvdMul Q.E.D.
-- Lemma: dvdAbs_l2r
-- Step: 1 (2 way case split)
-- Step: 1.1 Q.E.D.
-- Step: 1.2 Q.E.D.
-- Step: 1.Completeness Q.E.D.
-- Result: Q.E.D.
-- Lemma: dvdAbs_r2l
-- Step: 1 (2 way case split)
-- Step: 1.1 Q.E.D.
-- Step: 1.2 Q.E.D.
-- Step: 1.Completeness Q.E.D.
-- Result: Q.E.D.
-- Lemma: dvdAbs Q.E.D.
-- [Proven] dvdAbs :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
dvdAbs :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
dvdAbs = do
dM <- recall dvdMul
l2r <- calc "dvdAbs_l2r"
(\(Forall @"a" a) (Forall @"b" b) -> a `dvd` abs b .=> a `dvd` b) $
\a b -> [a `dvd` abs b]
|- cases [ b .>= 0 ==> a `dvd` b
=: sTrue
=: qed
, b .< 0 ==> a `dvd` b
?? dM `at` (Inst @"d" a, Inst @"a" (abs b), Inst @"k" (-1))
=: sTrue
=: qed
]
r2l <- calc "dvdAbs_r2l"
(\(Forall @"a" a) (Forall @"b" b) -> a `dvd` b .=> a `dvd` abs b) $
\a b -> [a `dvd` b]
|- cases [ b .>= 0 ==> a `dvd` abs b
=: sTrue
=: qed
, b .< 0 ==> a `dvd` abs b
?? dM `at` (Inst @"d" a, Inst @"a" b, Inst @"k" (-1))
=: sTrue
=: qed
]
lemma "dvdAbs"
(\(Forall @"a" a) (Forall @"b" b) -> a `dvd` b .== a `dvd` abs b)
[proofOf l2r, proofOf r2l]
-- | \(d \mid (2a + 1) \implies \mathrm{isOdd}(d)\)
--
-- ==== __Proof__
-- >>> runTP dvdOddThenOdd
-- Lemma: dvdOddThenOdd
-- Step: 1 (2 way case split)
-- Step: 1.1 Q.E.D.
-- Step: 1.2.1 Q.E.D.
-- Step: 1.2.2 Q.E.D.
-- Step: 1.Completeness Q.E.D.
-- Result: Q.E.D.
-- [Proven] dvdOddThenOdd :: Ɐd ∷ Integer → Ɐa ∷ Integer → Bool
dvdOddThenOdd :: TP (Proof (Forall "d" Integer -> Forall "a" Integer -> SBool))
dvdOddThenOdd = calc "dvdOddThenOdd"
(\(Forall d) (Forall a) -> d `dvd` (2*a+1) .=> isOdd d) $
\d a -> [d `dvd` (2*a+1)]
|- cases [ isOdd d ==> trivial
, isEven d ==> (2 * (d `sEDiv` 2)) `dvd` (2*a+1)
=: 2 `dvd` (2*a+1)
=: contradiction
]
-- | \(\mathrm{isOdd}(d) \land d \mid 2a \implies d \mid a\)
--
-- ==== __Proof__
-- >>> runTP dvdEvenWhenOdd
-- Lemma: dvdEvenWhenOdd
-- Step: 1 Q.E.D.
-- Step: 2 Q.E.D.
-- Step: 3 Q.E.D.
-- Step: 4 Q.E.D.
-- Step: 5 Q.E.D.
-- Step: 6 Q.E.D.
-- Step: 7 Q.E.D.
-- Result: Q.E.D.
-- [Proven] dvdEvenWhenOdd :: Ɐd ∷ Integer → Ɐa ∷ Integer → Bool
dvdEvenWhenOdd :: TP (Proof (Forall "d" Integer -> Forall "a" Integer -> SBool))
dvdEvenWhenOdd = calc "dvdEvenWhenOdd"
(\(Forall d) (Forall a) -> isOdd d .&& d `dvd` (2*a) .=> d `dvd` a) $
\d a -> [isOdd d, d `dvd` (2*a)]
|- let t = (d - 1) `sEDiv` 2
m = (2*a) `sEDiv` d
in sTrue
-- Observe that d = 2t+1 and 2a = dm
=: d .== 2*t + 1 .&& 2*a .== d*m
-- So, 2a == (2t+1)m holds
=: 2*a .== (2*t+1) * m
-- Arithmetic gives us
=: 2*a .== 2*t*m + m .&& 2*(a-t*m) .== m
-- So m = 2*(a-t*m), i.e., m is even
=: m .== 2 * (a - t*m)
-- Let n = a - t*m, so m = 2n. It follows that 2a = d(2n) = 2(dn)
=: let n = a - t*m
in 2*a .== d * (2 * n) .&& 2 * a .== 2 * (d * n)
-- From which we can conclude a = dn
=: a .== d * n
-- Thus we can deduce d must divide a
?? d `dvd` (d * n)
=: d `dvd` a
-- Done!
=: qed
-- | \(d \mid a \land d \mid b \implies d \mid (a + b)\)
--
-- ==== __Proof__
-- >>> runTP dvdSum1
-- Lemma: dvdSum1
-- Step: 1 (2 way case split)
-- Step: 1.1 Q.E.D.
-- Step: 1.2.1 Q.E.D.
-- Step: 1.2.2 Q.E.D.
-- Step: 1.2.3 Q.E.D.
-- Step: 1.Completeness Q.E.D.
-- Result: Q.E.D.
-- [Proven] dvdSum1 :: Ɐd ∷ Integer → Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
dvdSum1 :: TP (Proof (Forall "d" Integer -> Forall "a" Integer -> Forall "b" Integer -> SBool))
dvdSum1 =
calc "dvdSum1"
(\(Forall d) (Forall a) (Forall b) -> d `dvd` a .&& d `dvd` b .=> d `dvd` (a + b)) $
\d a b -> [d `dvd` a .&& d `dvd` b]
|- cases [ a .== 0 .|| b .== 0 ==> trivial
, a ./= 0 .&& b ./= 0 ==> d `dvd` (a + b)
=: d `dvd` (a `sEDiv` d * d + b `sEDiv` d * d)
=: d `dvd` (d * (a `sEDiv` d + b `sEDiv` d))
=: sTrue
=: qed
]
-- | \(d \mid (a + b) \land d \mid b \implies d \mid a \)
--
-- ==== __Proof__
-- >>> runTP dvdSum2
-- Lemma: dvdSum2
-- Step: 1 (2 way case split)
-- Step: 1.1 Q.E.D.
-- Step: 1.2.1 Q.E.D.
-- Step: 1.2.2 Q.E.D.
-- Step: 1.2.3 Q.E.D.
-- Step: 1.Completeness Q.E.D.
-- Result: Q.E.D.
-- [Proven] dvdSum2 :: Ɐd ∷ Integer → Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
dvdSum2 :: TP (Proof (Forall "d" Integer -> Forall "a" Integer -> Forall "b" Integer -> SBool))
dvdSum2 =
calc "dvdSum2"
(\(Forall d) (Forall a) (Forall b) -> d `dvd` (a + b) .&& d `dvd` b .=> d `dvd` a) $
\d a b -> [d `dvd` (a + b) .&& d `dvd` b]
|- cases [ d .== 0 ==> trivial
, d ./= 0 ==> let k1 = (a + b) `sEDiv` d
k2 = b `sEDiv` d
in a `sEDiv` d
=: (a + b - b) `sEDiv` d
=: (k1 * d - k2 * d) `sEDiv` d
=: (k1 - k2) * d `sEDiv` d
=: qed
]
-- * Correctness of GCD
-- | \(\gcd\,a\,b \mid a \land \gcd\,a\,b \mid b\)
--
-- GCD of two numbers divide these numbers. This is part one of the proof, where we are
-- not concerned with maximality. Our goal is to show that the calculated gcd divides both inputs.
--
-- ==== __Proof__
-- >>> runTP gcdDivides
-- Lemma: dvdAbs Q.E.D.
-- Lemma: helper
-- Step: 1 Q.E.D.
-- Result: Q.E.D.
-- Inductive lemma (strong): dvdNGCD
-- Step: Measure is non-negative Q.E.D.
-- Step: 1 (2 way case split)
-- Step: 1.1 Q.E.D.
-- Step: 1.2.1 Q.E.D.
-- Step: 1.2.2 Q.E.D.
-- Step: 1.Completeness Q.E.D.
-- Result: Q.E.D.
-- Lemma: gcdDivides Q.E.D.
-- Functions proven terminating: nGCD
-- [Proven] gcdDivides :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
gcdDivides :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
gcdDivides = do
dAbs <- recall dvdAbs
-- Helper about divisibility. If x|b and x| a%b, then x|a.
helper <- calc "helper"
(\(Forall @"a" a) (Forall @"b" b) (Forall @"x" x) ->
b ./= 0 .&& x `dvd` b .&& x `dvd` (a `sEMod` b)
.=> -----------------------------------------------
x `dvd` a
) $
\a b x -> [b ./= 0, x `dvd` b, x `dvd` (a `sEMod` b)]
|- x `dvd` a
?? a `sEDiv` x .== (a `sEDiv` b) * (b `sEDiv` x) + (a `sEMod` b) `sEDiv` x
=: sTrue
=: qed
-- Use strong induction to prove divisibility over non-negative numbers.
dNGCD <- sInduct "dvdNGCD"
(\(Forall @"a" a) (Forall @"b" b) -> a .>= 0 .&& b .>= 0 .=> nGCD a b `dvd` a .&& nGCD a b `dvd` b)
(\_a b -> b, []) $
\ih a b -> [a .>= 0, b .>= 0]
|- let g = nGCD a b
in g `dvd` a .&& g `dvd` b
=: cases [ b .== 0 ==> trivial
, b .> 0 ==> let g' = nGCD b (a `sEMod` b)
in g' `dvd` a .&& g' `dvd` b
?? ih `at` (Inst @"a" b, Inst @"b" (a `sEMod` b))
?? helper
=: sTrue
=: qed
]
-- Now generalize to arbitrary integers.
lemma"gcdDivides"
(\(Forall a) (Forall b) -> gcd a b `dvd` a .&& gcd a b `dvd` b)
[proofOf dAbs, proofOf dNGCD]
-- | \(x \mid a \land x \mid b \implies x \mid \gcd\,a\,b\)
--
-- Maximality. Any divisor of the inputs divides the GCD.
--
-- ==== __Proof__
-- >>> runTP gcdMaximal
-- Lemma: dvdAbs Q.E.D.
-- Lemma: commutative Q.E.D.
-- Lemma: eDiv Q.E.D.
-- Lemma: helper
-- Step: 1 (x `dvd` a && x `dvd` b) Q.E.D.
-- Step: 2 Q.E.D.
-- Step: 3 Q.E.D.
-- Result: Q.E.D.
-- Inductive lemma (strong): mNGCD
-- Step: Measure is non-negative Q.E.D.
-- Step: 1 (2 way case split)
-- Step: 1.1 Q.E.D.
-- Step: 1.2.1 Q.E.D.
-- Step: 1.2.2 Q.E.D.
-- Step: 1.Completeness Q.E.D.
-- Result: Q.E.D.
-- Lemma: gcdMaximal
-- Step: 1 (2 way case split)
-- Step: 1.1.1 Q.E.D.
-- Step: 1.1.2 Q.E.D.
-- Step: 1.2.1 Q.E.D.
-- Step: 1.2.2 Q.E.D.
-- Step: 1.2.3 Q.E.D.
-- Step: 1.2.4 Q.E.D.
-- Step: 1.Completeness Q.E.D.
-- Result: Q.E.D.
-- Functions proven terminating: nGCD
-- [Proven] gcdMaximal :: Ɐa ∷ Integer → Ɐb ∷ Integer → Ɐx ∷ Integer → Bool
gcdMaximal :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> Forall "x" Integer -> SBool))
gcdMaximal = do
dAbs <- recall dvdAbs
comm <- recall commutative
eDiv <- lemma "eDiv"
(\(Forall @"x" x) (Forall @"y" y) -> y ./= 0 .=> x .== (x `sEDiv` y) * y + x `sEMod` y)
[]
-- Helper: If x|a, x|b then x|a%b.
helper <- calc "helper"
(\(Forall @"a" a) (Forall @"b" b) (Forall @"x" x) ->
x ./= 0 .&& b ./= 0 .&& x `dvd` a .&& x `dvd` b
.=> -----------------------------------------------
x `dvd` (a `sEMod` b)
) $
\a b x -> [x ./= 0, b ./= 0, x `dvd` a, x `dvd` b]
|- x `dvd` (a `sEMod` b)
?? "x `dvd` a && x `dvd` b"
=: let k1 = a `sDiv` x
k2 = b `sDiv` x
in x `dvd` ((k1*x) `sEMod` (k2*x))
?? eDiv `at` (Inst @"x" (k1*x), Inst @"y" (k2*x))
=: x `dvd` ((k1*x) - ((k1*x) `sEDiv` (k2*x)) * (k2*x))
=: sTrue
=: qed
-- Now prove maximality for non-negative integers:
mNGCD <- sInduct "mNGCD"
(\(Forall @"a" a) (Forall @"b" b) (Forall @"x" x) ->
a .>= 0 .&& b .>= 0 .&& x `dvd` a .&& x `dvd` b .=> x `dvd` nGCD a b)
(\_a b _x -> b, []) $
\ih a b x -> let g = nGCD a b
in [a .>= 0, b .>= 0, x `dvd` a .&& x `dvd` b]
|- x `dvd` g
=: cases [ b .== 0 ==> trivial
, b .> 0 ==> x `dvd` nGCD b (a `sEMod` b)
?? ih `at` (Inst @"a" b, Inst @"b" (a `sEMod` b), Inst @"x" x)
?? helper
=: sTrue
=: qed
]
-- Generalize to arbitrary integers:
calc "gcdMaximal"
(\(Forall @"a" a) (Forall @"b" b) (Forall @"x" x) -> x `dvd` a .&& x `dvd` b .=> x `dvd` gcd a b) $
\a b x -> [x `dvd` a, x `dvd` b]
|- x `dvd` gcd a b
=: cases [ abs a .>= abs b ==> x `dvd` nGCD (abs a) (abs b)
?? mNGCD `at` (Inst @"a" (abs a), Inst @"b" (abs b), Inst @"x" x)
?? dAbs `at` (Inst @"a" x, Inst @"b" a)
?? dAbs `at` (Inst @"a" x, Inst @"b" b)
=: sTrue
=: qed
, abs a .< abs b ==> x `dvd` gcd a b
?? comm `at` (Inst @"a" a, Inst @"b" b)
=: x `dvd` gcd b a
=: x `dvd` nGCD (abs b) (abs a)
?? mNGCD `at` (Inst @"a" (abs b), Inst @"b" (abs a), Inst @"x" x)
?? dAbs `at` (Inst @"a" x, Inst @"b" a)
?? dAbs `at` (Inst @"a" x, Inst @"b" b)
=: sTrue
=: qed
]
-- | \(\gcd\,a\,b \mid a \land \gcd\,a\,b \mid b \land (x \mid a \land x \mid b \implies x \mid \gcd\,a\,b)\)
--
-- Putting it all together: GCD divides both arguments, and its maximal.
--
-- ==== __Proof__
-- >>> runTP gcdCorrect
-- Lemma: gcdDivides Q.E.D.
-- Lemma: gcdMaximal Q.E.D.
-- Lemma: gcdCorrect
-- Step: 1 Q.E.D.
-- Step: 2 Q.E.D.
-- Result: Q.E.D.
-- Functions proven terminating: nGCD
-- [Proven] gcdCorrect :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
gcdCorrect :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
gcdCorrect = do
divides <- recall gcdDivides
maximal <- recall gcdMaximal
calc "gcdCorrect"
(\(Forall a) (Forall b) ->
let g = gcd a b
in g `dvd` a
.&& g `dvd` b
.&& quantifiedBool (\(Forall x) -> x `dvd` a .&& x `dvd` b .=> x `dvd` g)
) $
\a b -> []
|- let g = gcd a b
m = quantifiedBool (\(Forall x) -> x `dvd` a .&& x `dvd` b .=> x `dvd` g)
in g `dvd` a .&& g `dvd` b .&& m
?? divides `at` (Inst @"a" a, Inst @"b" b)
=: m
?? maximal
=: sTrue
=: qed
-- | \(\bigl((a \neq 0 \lor b \neq 0) \land x \mid a \land x \mid b \bigr) \implies x \leq \gcd\,a\,b\)
--
-- Additionally prove that GCD is really maximum, i.e., it is the largest in the regular sense. Note
-- that we have to make an exception for @gcd 0 0@ since by definition the GCD is @0@, which is clearly
-- not the largest divisor of @0@ and @0@. (Since any number is a GCD for the pair @(0, 0)@, there is
-- no maximum.)
--
-- ==== __Proof__
-- >>> runTP gcdLargest
-- Lemma: gcdMaximal Q.E.D.
-- Lemma: gcdZero Q.E.D.
-- Lemma: nonNegative Q.E.D.
-- Lemma: gcdLargest
-- Step: 1 Q.E.D.
-- Step: 2 Q.E.D.
-- Result: Q.E.D.
-- Functions proven terminating: nGCD
-- [Proven] gcdLargest :: Ɐa ∷ Integer → Ɐb ∷ Integer → Ɐx ∷ Integer → Bool
gcdLargest :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> Forall "x" Integer -> SBool))
gcdLargest = do
maximal <- recall gcdMaximal
gcdZ <- recall gcdZero
nn <- recall gcdNonNegative
calc "gcdLargest"
(\(Forall a) (Forall b) (Forall x) -> (a ./= 0 .|| b ./= 0) .&& x `dvd` a .&& x `dvd` b .=> x .<= gcd a b) $
\a b x -> [(a ./= 0 .|| b ./= 0) .&& x `dvd` a, x `dvd` b]
|- x .<= gcd a b
?? maximal `at` (Inst @"a" a, Inst @"b" b, Inst @"x" x)
=: (x `dvd` gcd a b .=> x .<= gcd a b)
?? gcdZ `at` (Inst @"a" a, Inst @"b" b)
?? nn `at` (Inst @"a" a, Inst @"b" b)
=: sTrue
=: qed
-- * Other GCD Facts
-- | \(\gcd\, a\, b = \gcd\, (a + b)\, b\)
--
-- ==== __Proof__
-- >>> runTP gcdAdd
-- Lemma: dvdSum1 Q.E.D.
-- Lemma: dvdSum2 Q.E.D.
-- Lemma: gcdDivides Q.E.D.
-- Lemma: gcdLargest Q.E.D.
-- Lemma: gcdAdd
-- Step: 1 Q.E.D.
-- Step: 2 Q.E.D.
-- Step: 3 Q.E.D.
-- Step: 4 Q.E.D.
-- Step: 5 Q.E.D.
-- Step: 6 Q.E.D.
-- Step: 7 Q.E.D.
-- Result: Q.E.D.
-- Functions proven terminating: nGCD
-- [Proven] gcdAdd :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
gcdAdd :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
gcdAdd = do
dSum1 <- recall dvdSum1
dSum2 <- recall dvdSum2
divides <- recall gcdDivides
largest <- recall gcdLargest
calc "gcdAdd"
(\(Forall @"a" a) (Forall @"b" b) -> gcd a b .== gcd (a + b) b) $
\a b -> [] |-> let g1 = gcd a b
g2 = gcd (a + b) b
in sTrue
-- First use the divides property to conclude that g1 divides a and b
?? divides `at` (Inst @"a" a, Inst @"b" b)
=: g1 `dvd` a .&& g1 `dvd` b
-- Same for g2 for a+b and b
?? divides `at` (Inst @"a" (a + b), Inst @"b" b)
=: g2 `dvd` (a+b) .&& g2 `dvd` b
-- Use dSum1 to show g1 divides a+b
?? dSum1 `at` (Inst @"d" g1, Inst @"a" a, Inst @"b" b)
=: g1 `dvd` (a+b)
-- Similarly, use dSum2 to show g2 divides a
?? dSum2 `at` (Inst @"d" g2, Inst @"a" a, Inst @"b" b)
=: g2 `dvd` a
-- Now use largest to show g1 >= g2
?? largest `at` (Inst @"a" a, Inst @"b" b, Inst @"x" g2)
=: g1 .>= g2
-- But again via largest, we can show g2 >= g1
?? largest `at` (Inst @"a" (a+b), Inst @"b" b, Inst @"x" g1)
=: g2 .>= g1
-- Finally conclude g1 = g2, since both are greater-than-equal to each other:
=: g1 .== g2
=: qed
-- | \(\gcd\, (2a)\, (2b) = 2 (\gcd\,a\, b)\)
--
-- ==== __Proof__
-- >>> runTP gcdEvenEven
-- Lemma: red2 Q.E.D.
-- Lemma: modEE
-- Step: 1 Q.E.D.
-- Step: 2 Q.E.D.
-- Step: 3 Q.E.D.
-- Result: Q.E.D.
-- Inductive lemma (strong): nGCDEvenEven
-- Step: Measure is non-negative Q.E.D.
-- Step: 1 (2 way case split)
-- Step: 1.1 Q.E.D.
-- Step: 1.2.1 Q.E.D.
-- Step: 1.2.2 Q.E.D.
-- Step: 1.2.3 Q.E.D.
-- Step: 1.2.4 Q.E.D.
-- Step: 1.Completeness Q.E.D.
-- Result: Q.E.D.
-- Lemma: gcdEvenEven
-- Step: 1 Q.E.D.
-- Step: 2 Q.E.D.
-- Step: 3 Q.E.D.
-- Step: 4 Q.E.D.
-- Result: Q.E.D.
-- Functions proven terminating: nGCD
-- [Proven] gcdEvenEven :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
gcdEvenEven :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
gcdEvenEven = do
red2 <- lemmaWith z3 "red2"
(\(Forall @"a" a) (Forall @"b" b) -> b ./= 0 .=> (2*a) `sEDiv` (2*b) .== a `sEDiv` b)
[]
modEE <- calcWith cvc5 "modEE"
(\(Forall @"a" a) (Forall @"b" b) -> b ./= 0 .=> (2*a) `sEMod` (2*b) .== 2 * (a `sEMod` b)) $
\a b -> [b ./= 0]
|- (2*a) `sEMod` (2*b)
?? red2 `at` (Inst @"a" a, Inst @"b" b)
=: 2*a - 2*b * (a `sEDiv` b)
=: 2 * (a - b * (a `sEDiv` b))
=: 2 * (a `sEMod` b)
=: qed
nGCDEvenEven <- sInduct "nGCDEvenEven"
(\(Forall @"a" a) (Forall @"b" b) -> a .>= 0 .&& b .>= 0 .=> nGCD (2*a) (2*b) .== 2 * nGCD a b)
(\_a b -> b, []) $
\ih a b -> [a .>= 0, b .>= 0]
|- nGCD (2*a) (2*b)
=: cases [ b .== 0 ==> trivial
, b ./= 0 ==> nGCD (2 * a) (2 * b)
=: nGCD (2 * b) ((2 * a) `sEMod` (2 * b))
?? modEE `at` (Inst @"a" a, Inst @"b" b)
=: nGCD (2 * b) (2 * (a `sEMod` b))
?? ih
=: 2 * nGCD a b
=: qed
]
calc "gcdEvenEven"
(\(Forall a) (Forall b) -> gcd (2*a) (2*b) .== 2 * gcd a b) $
\a b -> [] |- gcd (2*a) (2*b)
=: nGCD (abs (2*a)) (abs (2*b))
=: nGCD (2 * abs a) (2 * abs b)
?? nGCDEvenEven `at` (Inst @"a" (abs a), Inst @"b" (abs b))
=: 2 * nGCD (abs a) (abs b)
=: 2 * gcd a b
=: qed
-- | \(\gcd\, (2a+1)\, (2b) = \gcd\,(2a+1)\, b\)
--
-- ==== __Proof__
-- >>> runTP gcdOddEven
-- Lemma: gcdDivides Q.E.D.
-- Lemma: gcdLargest Q.E.D.
-- Lemma: dvdMul Q.E.D. [Cached]
-- Lemma: dvdOddThenOdd Q.E.D.
-- Lemma: dvdEvenWhenOdd Q.E.D.
-- Lemma: gcdOddEven
-- Step: 1 Q.E.D.
-- Step: 2 Q.E.D.
-- Step: 3 Q.E.D.
-- Step: 4 Q.E.D.
-- Step: 5 Q.E.D.
-- Step: 6 Q.E.D.
-- Step: 7 Q.E.D.
-- Step: 8 Q.E.D.
-- Result: Q.E.D.
-- Functions proven terminating: nGCD
-- [Proven] gcdOddEven :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
gcdOddEven :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
gcdOddEven = do
divides <- recall gcdDivides
largest <- recall gcdLargest
dMul <- recall dvdMul
dOddThenOdd <- recall dvdOddThenOdd
dEvenWhenOdd <- recall dvdEvenWhenOdd
calc "gcdOddEven"
(\(Forall a) (Forall b) -> gcd (2*a+1) (2*b) .== gcd (2*a+1) b) $
\a b -> [] |-> let g1 = gcd (2*a+1) (2*b)
g2 = gcd (2*a+1) b
in sTrue
-- First use the divides property to conclude that g1 divides both 2*a+1 and 2*b
?? divides `at` (Inst @"a" (2*a+1), Inst @"b" (2*b))
=: g1 `dvd` (2*a+1) .&& g1 `dvd` (2*b)
-- Same for g2, for 2*a+1 and b
?? divides `at` (Inst @"a" (2*a+1), Inst @"b" b)
=: g2 `dvd` (2*a+1) .&& g2 `dvd` b
-- By arithmetic, g2 divides 2*b
?? dMul `at` (Inst @"d" g2, Inst @"a" b, Inst @"k" 2)
=: g2 `dvd` (2*b)
-- Observe that g1 must be odd
?? dOddThenOdd `at` (Inst @"d" g1, Inst @"a" a)
=: isOdd g1
-- Conclude that g1 must divide b
?? dEvenWhenOdd `at` (Inst @"d" g1, Inst @"a" b)
=: g1 `dvd` b
-- Now use largest to show g1 >= g2
?? largest `at` (Inst @"a" (2*a+1), Inst @"b" (2*b), Inst @"x" g2)
=: g1 .>= g2
-- But again via largest, we can show g2 >= g1
?? largest `at` (Inst @"a" (2*a+1), Inst @"b" b, Inst @"x" g1)
=: g2 .>= g1
-- Finally conclude g1 = g2 since both are greater-than-equal to each other:
=: g1 .== g2
=: qed
-- * GCD via subtraction
-- | @nGCDSub@ is the original version of Euclid, which uses subtraction instead of modulus. This is the version that
-- works on non-negative numbers. It has the precondition that @a >= b >= 0@, and maintains this invariant in each
-- recursive call.
nGCDSub :: SInteger -> SInteger -> SInteger
nGCDSub = smtFunction "nGCDSub"
$ \a b -> [sCase| a of
_ | a .== b -> a
_ | a .<= 0 -> b
_ | b .<= 0 -> a
_ | a .> b -> nGCDSub (a - b) b
_ -> nGCDSub a (b - a)
|]
-- | Generalized version of subtraction based GCD, working over all integers.
gcdSub :: SInteger -> SInteger -> SInteger
gcdSub a b = nGCDSub (abs a) (abs b)
-- | \(\mathrm{gcdSub}\, a\, b = \gcd\, a\, b\)
--
-- Instead of proving @gcdSub@ correct, we'll simply show that it is equivalent to @gcd@, hence it has
-- all the properties we already established.
--
-- ==== __Proof__
-- >>> runTP gcdSubEquiv
-- Lemma: commutative Q.E.D.
-- Lemma: gcdAdd Q.E.D.
-- Inductive lemma (strong): nGCDSubEquiv
-- Step: Measure is non-negative Q.E.D.
-- Step: 1 (5 way case split)
-- Step: 1.1 Q.E.D.
-- Step: 1.2 Q.E.D.
-- Step: 1.3 Q.E.D.
-- Step: 1.4.1 Q.E.D.
-- Step: 1.4.2 Q.E.D.
-- Step: 1.4.3 Q.E.D.
-- Step: 1.5.1 Q.E.D.
-- Step: 1.5.2 Q.E.D.
-- Step: 1.5.3 Q.E.D.
-- Step: 1.5.4 Q.E.D.
-- Step: 1.5.5 Q.E.D.
-- Step: 1.Completeness Q.E.D.
-- Result: Q.E.D.
-- Lemma: gcdSubEquiv
-- Step: 1 Q.E.D.
-- Step: 2 Q.E.D.
-- Step: 3 Q.E.D.
-- Result: Q.E.D.
-- Functions proven terminating: nGCD, nGCDSub
-- [Proven] gcdSubEquiv :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
gcdSubEquiv :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
gcdSubEquiv = do
-- We'll be using the commutativity of GCD and the gcdAdd property
comm <- recall commutative
addG <- recall gcdAdd
-- First prove over the non-negative numbers:
nEq <- sInduct "nGCDSubEquiv"
(\(Forall @"a" a) (Forall @"b" b) -> a .>= 0 .&& b .>= 0 .=> nGCDSub a b .== nGCD a b)
(\a b -> a + b, []) $
\ih a b -> [a .>= 0, b .>= 0]
|- nGCDSub a b
=: cases [ a .== b ==> nGCD a b =: qed
, a .== 0 ==> nGCD a b =: qed
, b .== 0 ==> nGCD a b =: qed
, a .> b .&& b ./= 0 ==> nGCDSub (a - b) b
?? ih
=: nGCD (a - b) b
?? addG `at` (Inst @"a" (a - b), Inst @"b" b)
=: nGCD a b
=: qed
, a .< b .&& a ./= 0 ==> nGCDSub a (b - a)
?? ih
=: nGCD a (b - a)
?? comm
=: nGCD (b - a) a
?? addG `at` (Inst @"a" (b - a), Inst @"b" a)
=: nGCD b a
?? comm
=: nGCD a b
=: qed
]
-- Now prove over all integers
calcWith cvc5 "gcdSubEquiv"
(\(Forall a) (Forall b) -> gcd a b .== gcdSub a b) $
\a b -> [] |- gcd a b
=: nGCD (abs a) (abs b)
?? nEq `at` (Inst @"a" (abs a), Inst @"b" (abs b))
=: nGCDSub (abs a) (abs b)
=: gcdSub a b
=: qed
-- * Binary GCD
-- | @nGCDBin@ is the binary GCD algorithm that works on non-negative numbers.
nGCDBin :: SInteger -> SInteger -> SInteger
nGCDBin = smtFunction "nGCDBin"
$ \a b -> [sCase| a of
_ | a .<= 0 -> b
_ | b .<= 0 -> a
_ | isEven a .&& isEven b -> 2 * nGCDBin (a `sEDiv` 2) (b `sEDiv` 2)
_ | isOdd a .&& isEven b -> nGCDBin a (b `sEDiv` 2)
_ | a .<= b -> nGCDBin a (b - a)
_ -> nGCDBin (a - b) b
|]
-- | Generalized version that works on arbitrary integers.
gcdBin :: SInteger -> SInteger -> SInteger
gcdBin a b = nGCDBin (abs a) (abs b)
-- | \(\mathrm{gcdBin}\, a\, b = \gcd\, a\, b\)
--
-- Instead of proving @gcdBin@ correct, we'll simply show that it is equivalent to @gcd@, hence it has
-- all the properties we already established.
--
-- ==== __Proof__
-- >>> runTP gcdBinEquiv
-- Lemma: gcdEvenEven Q.E.D.
-- Lemma: gcdOddEven Q.E.D.
-- Lemma: gcdAdd Q.E.D.
-- Lemma: commutative Q.E.D. [Cached]
-- Inductive lemma (strong): nGCDBinEquiv
-- Step: Measure is non-negative Q.E.D.
-- Step: 1 (5 way case split)
-- Step: 1.1 Q.E.D.
-- Step: 1.2 Q.E.D.
-- Step: 1.3.1 Q.E.D.
-- Step: 1.3.2 Q.E.D.
-- Step: 1.3.3 Q.E.D.
-- Step: 1.4.1 Q.E.D.
-- Step: 1.4.2 Q.E.D.
-- Step: 1.4.3 Q.E.D.
-- Step: 1.5 (3 way case split)
-- Step: 1.5.1 Q.E.D.
-- Step: 1.5.2.1 Q.E.D.
-- Step: 1.5.2.2 Q.E.D.
-- Step: 1.5.2.3 Q.E.D.
-- Step: 1.5.2.4 Q.E.D.
-- Step: 1.5.2.5 Q.E.D.
-- Step: 1.5.2.6 Q.E.D.
-- Step: 1.5.3.1 Q.E.D.
-- Step: 1.5.3.2 Q.E.D.
-- Step: 1.5.3.3 Q.E.D.
-- Step: 1.5.3.4 Q.E.D.
-- Step: 1.5.Completeness Q.E.D.
-- Step: 1.Completeness Q.E.D.
-- Result: Q.E.D.
-- Lemma: gcdBinEquiv
-- Step: 1 Q.E.D.
-- Step: 2 Q.E.D.
-- Step: 3 Q.E.D.
-- Result: Q.E.D.
-- Functions proven terminating: nGCD, nGCDBin
-- [Proven] gcdBinEquiv :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
gcdBinEquiv :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
gcdBinEquiv = do
gEvenEven <- recallWith cvc5 gcdEvenEven
gOddEven <- recall gcdOddEven
gAdd <- recall gcdAdd
comm <- recall commutative
-- First prove over the non-negative numbers:
nEq <- sInduct "nGCDBinEquiv"
(\(Forall @"a" a) (Forall @"b" b) -> a .>= 0 .&& b .>= 0 .=> nGCDBin a b .== nGCD a b)
(\a b -> tuple (a, b), []) $
\ih a b -> [a .>= 0, b .>= 0]
|- nGCDBin a b
=: cases [ a .== 0 ==> trivial
, b .== 0 ==> trivial
, isEven a .&& isEven b ==> 2 * nGCDBin (a `sEDiv` 2) (b `sEDiv` 2)
?? ih `at` (Inst @"a" (a `sEDiv` 2), Inst @"b" (b `sEDiv` 2))
=: 2 * nGCD (a `sEDiv` 2) (b `sEDiv` 2)
?? a .== 2 * a `sEDiv` 2
?? b .== 2 * b `sEDiv` 2
?? gEvenEven `at` (Inst @"a" (a `sEDiv` 2), Inst @"b" (b `sEDiv` 2))
=: nGCD a b
=: qed
, isOdd a .&& isEven b ==> nGCDBin a (b `sEDiv` 2)
?? ih `at` (Inst @"a" a, Inst @"b" (b `sEDiv` 2))
=: nGCD a (b `sEDiv` 2)
?? a .== 2 * ((a-1) `sEDiv` 2) + 1
?? b .== 2 * b `sEDiv` 2
?? gOddEven `at` (Inst @"a" ((a-1) `sEDiv` 2), Inst @"b" (b `sEDiv` 2))
=: nGCD a b
=: qed
, isOdd b ==> cases [ a .== 0 ==> trivial
, a ./= 0 .&& a .<= b ==> nGCDBin a b
=: nGCDBin a (b - a)
?? ih `at` (Inst @"a" a, Inst @"b" (b - a))
=: nGCD a (b - a)
?? comm `at` (Inst @"a" a, Inst @"b" (b - a))
=: nGCD (b - a) a
?? gAdd `at` (Inst @"a" (b - a), Inst @"b" a)
=: nGCD b a
?? comm `at` (Inst @"a" b, Inst @"b" a)
=: nGCD a b
=: qed
, a .> b ==> nGCDBin a b
=: nGCDBin (a - b) b
?? ih `at` (Inst @"a" (a - b), Inst @"b" b)
=: nGCD (a - b) b
?? gAdd `at` (Inst @"a" a, Inst @"b" (-b))
=: nGCD a b
=: qed
]
]
-- Now prove over all integers
calcWith cvc5 "gcdBinEquiv"
(\(Forall a) (Forall b) -> gcd a b .== gcdBin a b) $
\a b -> [] |- gcd a b
=: nGCD (abs a) (abs b)
?? nEq `at` (Inst @"a" (abs a), Inst @"b" (abs b))
=: nGCDBin (abs a) (abs b)
=: gcdBin a b
=: qed
{- HLint ignore gcdSubEquiv "Avoid lambda" -}
{- HLint ignore gcdBinEquiv "Use curry" -}