TwitterProblem.hs

{-# LANGUAGE TypeApplications #-}
module TwitterProblem where

import Data.SBV

--------------------------------------------------------------------------------
-- Complex Numbers, constructed from Algebraic Reals

data AlgComplex a = AlgComplex { real :: a, imaginary :: a}

instance (Floating a) => Num (AlgComplex a) where
    fromInteger i = AlgComplex (fromInteger i) 0
    (AlgComplex r0 i0) + (AlgComplex r1 i1) = AlgComplex (r0 + r1) (i0 + i1)
    (AlgComplex r0 i0) - (AlgComplex r1 i1) = AlgComplex (r0 - r1) (i0 - i1)
    (AlgComplex r0 i0) * (AlgComplex r1 i1) = AlgComplex (r0*r1 - i0*i1) (r0*i1 + r1*i0)
    abs (AlgComplex r i) = AlgComplex (sqrt (r^2 + i^2)) 0
    signum (AlgComplex r i) = AlgComplex (signum (sqrt (r^2 + i^2))) 0

instance (Floating a) => Fractional (AlgComplex a) where
    fromRational r = AlgComplex (fromRational r) 0
    (AlgComplex r0 i0) / (AlgComplex r1 i1) =
        AlgComplex ((r0 * r1 + i0 * i1)/d) ((i0*r1 - r0*i1)/d)
      where
        d = r1^2 + i1^2

--------------------------------------------------------------------------------
-- SMT Solving

solution :: IO (Maybe (AlgReal, AlgReal, AlgReal, AlgReal))
solution = extractModel `fmap` satWith z3 problem
    where
      problem = do
        setLogic QF_NRA
        xr <- free @AlgReal "xr"
        xi <- free @AlgReal "xi"
        yr <- free @AlgReal "yr"
        yi <- free @AlgReal "yi"
        let x = AlgComplex xr xi
        let y = AlgComplex yr yi
        pure $ sAnd
          [ (real $ x + y) .== 1
          , (imaginary $ x + y) .== 0
          , (real $ (1/x) + (1/y)) .== 1
          , (imaginary $ (1/x) + (1/y)) .== 0
          , xr ./= 0 .|| xi ./= 0
          , yr ./= 0 .|| yi ./= 0
          ]

-- > Just (0.5,root(2, 4x^2 = 3) = 0.8660254037844386...,0.5,root(1, 4x^2 = 3) = -0.8660254037844386...)