teoremaMidy :: Integer -> Bool

teoremaMidy 200  ==  True

import Data.List (genericTake)
import Data.Numbers.Primes (primes, isPrime)
import Test.QuickCheck

teoremaMidy :: Integer -> Bool
teoremaMidy n =
    all conclusionMidy [(a,p) | p <- takeWhile (<n) (drop 3 primes),
                                a <- [1..p-1],
                                tienePeriodoDeLongitudPar (a,p)]

-- (tienePeriodoDeLongitudPar (a,p)) se verifica si el período de a/p
-- (con p un primo mayor que 5 y que a) es de longitud par. Por ejemplo,
--    tienePeriodoDeLongitudPar (1,7)   ==  True
--    tienePeriodoDeLongitudPar (1,37)  ==  False
tienePeriodoDeLongitudPar :: (Integer,Integer) -> Bool
tienePeriodoDeLongitudPar (a,p) =
    even (longitudPeriodo p)

-- (expansionDec (a,p)) es la expansión decimal de a/p. Por ejemplo,
--    take 10 (expansionDec (1,4))    ==  [0,2,5]
--    take 10 (expansionDec (1,7))    ==  [0,1,4,2,8,5,7,1,4,2]
--    take 12 (expansionDec (90,7))   ==  [12,8,5,7,1,4,2,8,5,7,1,4]
--    take 12 (expansionDec (23,14))  ==  [1,6,4,2,8,5,7,1,4,2,8,5]
expansionDec :: (Integer,Integer) -> [Integer]
expansionDec (x,y)
    | r == 0    = [q]
    | otherwise = q : expansionDec (r*10,y)
    where (q,r) = quotRem x y

-- (longitudPeriodo p) es la longitud del período de 1/p (y de
-- cualquier fracción irreducible de denominador p), donde p es un
-- número primo mayor que 5. Por ejemplo,
--    longitudPeriodo 7   ==  6
--    longitudPeriodo 37  ==  3
--    longitudPeriodo 83  ==  41
longitudPeriodo :: Integer -> Integer
longitudPeriodo p = head [k | k <- [1..], 10^k `mod` p == 1]

-- (periodo (a,p)) es el período de a/b (para a < p y p > 5 es primo).
-- Por ejemplo,
--    periodo (1,7)   ==  [1,4,2,8,5,7]
--    periodo (2,7)   ==  [2,8,5,7,1,4]
--    periodo (6,7)   ==  [8,5,7,1,4,2]
--    periodo (1,73)  ==  [0,1,3,6,9,8,6,3]
periodo :: (Integer,Integer) -> [Integer]
periodo (a,p) = genericTake t xs
    where (_:xs) = expansionDec (a,p)
          t      = longitudPeriodo p

-- (conclusionMidy (a,p)) se verifica si a/p verifica la conclusión del
-- teorema de Midy; es decir, podemos partir el período de a/p en dos
-- mitades, cuya suma es un número formado únicamente por nueves. Por
-- ejemplo,
--    conclusionMidy (1,7)   ==  True
--    conclusionMidy (1,37)  ==  False
conclusionMidy :: (Integer,Integer) -> Bool
conclusionMidy (a,p) =
    all (==9) (zipWith (+) ys zs)
    where xs      = periodo (a,p)
          (ys,zs) = splitAt (length xs `div` 2) xs

-- ---------------------------------------------------------------------
-- § Comprobación con QuickCheck                                      --
-- ---------------------------------------------------------------------

-- (condicionMidy (a,p)) se verifica si a/p cumple las condiciones del
-- teorema de Midy. Por ejemplo,
--    condicionMidy (1,7)   ==  True
--    condicionMidy (1,37)  ==  False
condicionMidy :: (Integer,Integer) -> Bool
condicionMidy (a,p) =
    a > 0 &&
    a < p &&
    isPrime p &&
    p > 5 &&
    tienePeriodoDeLongitudPar (a,p)

-- La propiedad es
prop_teoremaMidy :: (Integer,Integer) -> Property
prop_teoremaMidy (a,p) =
    condicionMidy (a,p) ==> conclusionMidy (a,p)

-- La comprobación es
--    λ> quickCheck prop_teoremaMidy
--    *** Gave up! Passed only 24 tests.

-- En la comprobación se observa que la mayoría de los ejemplos
-- generados no cumplen la condición del teorema de Midy y sólo 24 de
-- ellos la cumplen y también la conclusión.

-- Para aumentar el número de casos en lo que se cumpla la condición, se
-- puede aumentar el valor de maxSuccess; por ejemplo,
--    λ> quickCheckWith (stdArgs {maxSuccess=700}) prop_teoremaMidy
--    *** Gave up! Passed only 137 tests.

-- Otra forma de conseguirlo es definir un generador de números de Midy
-- (es decir, de pares (a,p) que cumplan la condición del teorema de
-- Midy). Por ejemplo,
--    λ> sample numeroMidy
--    (15,19)
--    (5,7)
--    (1,23)
--    (7,11)
numeroMidy :: Gen (Integer,Integer)
numeroMidy = suchThat arbitrary condicionMidy

-- La propiedad es
prop_teoremaMidy2 :: Property
prop_teoremaMidy2 = forAll numeroMidy conclusionMidy

-- La comprobación es
--    λ> quickCheck prop_teoremaMidy2
--    +++ OK, passed 100 tests.