Bezout.hs

{-# LANGUAGE BangPatterns, ScopedTypeVariables #-}

import Test.QuickCheck
import qualified Test.QuickCheck.Monadic as QM
import qualified Z3.Monad as Z3

{-
a x + b y = n d where d = gcd(a,b) の解は必ず
(x, y) = (n x0 + b/d k, n y0 - a/d k)
の形で書けることを示したい。

quickCheck prop_bezout2 する限りはそういうのはなさそう。
-}
prop_bezout2 = QM.monadicIO $ do
  NonZero (a :: Integer) <- QM.pick arbitrary
  NonZero (b :: Integer) <- QM.pick arbitrary
  let (d, x0, y0) = exgcd a b
  QM.monitor (counterexample (show (d,x0,y0)))
  (n :: Integer) <- QM.pick arbitrary

  ret <- QM.run $ Z3.evalZ3 $ do
    a_ <- Z3.mkIntNum a
    b_ <- Z3.mkIntNum b
    n_ <- Z3.mkIntNum n
    x <- Z3.mkIntVar =<< Z3.mkStringSymbol "x"
    y <- Z3.mkIntVar =<< Z3.mkStringSymbol "y"
    e <- do
      ax <- Z3.mkMul [a_, x]
      by <- Z3.mkMul [b_, y]
      Z3.mkAdd [ax, by]
    Z3.solverAssertCnstr =<< Z3.mkEq e =<< Z3.mkIntNum (n * d)

    x0_ <- Z3.mkIntNum x0
    y0_ <- Z3.mkIntNum y0
    tmp1 <- Z3.mkMul [n_, x0_] >>= \nx0 -> Z3.mkSub [x, nx0]
    tmp2 <- Z3.mkMul [n_, y0_] >>= \ny0 -> Z3.mkSub [y, ny0]
    zero <- Z3.mkIntNum 0
    cond1 <- Z3.mkEq zero =<< Z3.mkMod tmp1 =<< Z3.mkIntNum (b `div` d)
    cond2 <- Z3.mkEq zero =<< Z3.mkMod tmp2 =<< Z3.mkIntNum (a `div` d)
    k1 <- Z3.mkDiv tmp1 =<< Z3.mkIntNum (b `div` d)
    k2 <- Z3.mkUnaryMinus =<< Z3.mkDiv tmp2 =<< Z3.mkIntNum (a `div` d)
    cond3 <- Z3.mkEq k1 k2
    Z3.solverAssertCnstr =<< Z3.mkNot =<< Z3.mkAnd [cond1, cond2, cond3]

    ret <- Z3.solverCheck
    case ret of
      Z3.Sat -> do
        model <- Z3.solverGetModel
        Just xVal <- Z3.evalInt model x
        Just yVal <- Z3.evalInt model y
        return $ Just (xVal, yVal)
      _ -> return Nothing

  case ret of
    Nothing -> return ()
    Just (x,y) -> do
      QM.monitor (counterexample (show (x,y)))
      QM.assert False


exgcd :: Integral a => a -> a -> (a, a, a)
exgcd f1 f2 = f $ go f1 f2 1 0 0 1
  where
    go !r0 !r1 !s0 !s1 !t0 !t1
      | r1 == 0   = (r0, s0, t0)
      | otherwise = go r1 r2 s1 s2 t1 t2
      where
        (q, r2) = r0 `divMod` r1
        s2 = s0 - q*s1
        t2 = t0 - q*t1
    f (g,u,v)
      | g < 0 = (-g, -u, -v)
      | otherwise = (g,u,v)