-----------------------------------------------------------------------------
-- |
-- 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 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.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 -> ite (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.
-- [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.
-- [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.
-- [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.
-- [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: negGCD                           Q.E.D.
-- [Proven] negGCD :: Ɐa ∷ Integer → Bool
zeroGCD :: TP (Proof (Forall "a" Integer -> SBool))
zeroGCD = lemma "negGCD" (\(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 we @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)

-- | \(a \mid |b| \iff a \mid b\)
--
-- A number divides another exactly when it also divides its absolute value. While this property
-- seems obvious, I was unable to get z3 to prove it. Even CVC5 needs a bit of help to guide it through
-- the case split on @b@.
--
-- ==== __Proof__
-- >>> runTP dvdAbs
-- 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
   l2r <- calcWith cvc5 "dvdAbs_l2r"
                   (\(Forall @"a" a) (Forall @"b" b) -> a `dvd` abs b .=> a `dvd` b) $
                   \a b -> [a `dvd` abs b]
                        |- cases [ b .<  0 ==> sTrue =: qed
                                 , b .>= 0 ==> sTrue =: qed
                                 ]

   r2l <- calcWith cvc5 "dvdAbs_r2l"
                   (\(Forall @"a" a) (Forall @"b" b) -> a `dvd` b .=> a `dvd` abs b) $
                   \a b -> [a `dvd` b]
                        |- cases [ b .<  0 ==> sTrue =: qed
                                 , b .>= 0 ==> sTrue =: qed
                                 ]

   lemma "dvdAbs"
         (\(Forall @"a" a) (Forall @"b" b) -> a `dvd` b .== a `dvd` abs b)
         [proofOf l2r, proofOf r2l]

-- | \(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                           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)
                                         ?? a .== d * a `sEDiv` d
                                         =: d `dvd` (k * d * a `sEDiv` d)
                                         =: qed
                              ]

-- | \(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, we now now m is even
                            =: 2 `sDivides` m

                            -- Give that divisor a name:
                            =: let n = m `sEDiv` 2

                            -- It follows that 2a = d(2n) = 2(dn)
                            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` 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.
-- [Proven] gcdDivides :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
gcdDivides :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
gcdDivides = do

   dAbs <- recall "dvdAbs" 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: gcdComm                          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.
-- [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" dvdAbs
   comm <- recall "gcdComm" 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.
-- [Proven] gcdCorrect :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
gcdCorrect :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
gcdCorrect = do
  divides <- recall "gcdDivides" gcdDivides
  maximal <- recall "gcdMaximal" 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: gcdNonNegative                   Q.E.D.
-- Lemma: gcdLargest
--   Step: 1                               Q.E.D.
--   Step: 2                               Q.E.D.
--   Result:                               Q.E.D.
-- [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"     gcdMaximal
   zero    <- recall "gcdZero"        gcdZero
   nn      <- recall "gcdNonNegative" 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)
               ?? zero  `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.
-- [Proven] gcdAdd :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
gcdAdd :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
gcdAdd = do

   dSum1   <- recall "dvdSum1"    dvdSum1
   dSum2   <- recall "dvdSum2"    dvdSum2
   divides <- recall "gcdDivides" gcdDivides
   largest <- recall "gcdLargest" 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: modEE                            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.
-- [Proven] gcdEvenEven :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
gcdEvenEven :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
gcdEvenEven = do

   modEE <- lemma "modEE"
                  (\(Forall @"a" a) (Forall @"b" b) -> b ./= 0 .=> (2 * a) `sEMod` (2 * b) .== 2 * (a `sEMod` b))
                  []

   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.
-- 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.
-- [Proven] gcdOddEven :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
gcdOddEven :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
gcdOddEven = do

   divides      <- recall "gcdDivides"     gcdDivides
   largest      <- recall "gcdLargest"     gcdLargest
   dMul         <- recall "dvdMul"         dvdMul
   dOddThenOdd  <- recall "dvdOddThenOdd"  dvdOddThenOdd
   dEvenWhenOdd <- recall "dvdEvenWhenOdd" 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 verision 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 -> ite (a .== b) a
                                        $ ite (a .== 0) b
                                        $ ite (b .== 0) a
                                        $ ite (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.
-- [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" commutative
   addG <- recall "gcdAdd"      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 -> ite (a .== 0)               b
                                        $ ite (b .== 0)               a
                                        $ ite (isEven a .&& isEven b) (2 * nGCDBin (a `sEDiv` 2) (b `sEDiv` 2))
                                        $ ite (isOdd  a .&& isEven b) (    nGCDBin a             (b `sEDiv` 2))
                                        $ ite (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.
-- 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.
-- [Proven] gcdBinEquiv :: Ɐa ∷ Integer → Ɐb ∷ Integer → Bool
gcdBinEquiv :: TP (Proof (Forall "a" Integer -> Forall "b" Integer -> SBool))
gcdBinEquiv = do
   gEvenEven <- recall "gcdEvenEven" gcdEvenEven
   gOddEven  <- recall "gcdOddEven"  gcdOddEven
   gAdd      <- recall "gcdAdd"      gcdAdd
   comm      <- recall "commutative" 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"    -}