La función indicatriz de Euler
Publicado el 24 de enero de 2024 por José A. Alonso

phi :: Integer -> Integer

phi 36                          ==  12
map phi [10..20]                ==  [4,10,4,12,6,8,8,16,6,18,8]
phi (3^10^5) `mod` (10^9)       ==  681333334
length (show (phi (10^(10^5)))) ==  100000

module La_funcion_indicatriz_de_Euler where

import Data.List (genericLength, group)
import Data.Numbers.Primes (primeFactors)
import Math.NumberTheory.ArithmeticFunctions (totient)
import Test.QuickCheck (Positive (Positive), quickCheck)
import Test.Hspec (Spec, describe, hspec, it, shouldBe)

-- 1ª solución
-- ===========

phi1 :: Integer -> Integer
phi1 n = genericLength [x | x <- [1..n], gcd x n == 1]

-- 2ª solución
-- ===========

phi2 :: Integer -> Integer
phi2 n = product [(p-1)*p^(e-1) | (p,e) <- factorizacion n]

factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
  [(head xs,genericLength xs) | xs <- group (primeFactors n)]

-- 3ª solución
-- =============

phi3 :: Integer -> Integer
phi3 n =
  product [(x-1) * product xs | (x:xs) <- group (primeFactors n)]

-- 4ª solución
-- =============

phi4 :: Integer -> Integer
phi4 = totient

-- Verificación
-- ============

verifica :: IO ()
verifica = hspec spec

specG :: (Integer -> Integer) -> Spec
specG phi = do
  it "e1" $
    phi 36 `shouldBe` 12
  it "e2" $
    map phi [10..20] `shouldBe` [4,10,4,12,6,8,8,16,6,18,8]

spec :: Spec
spec = do
  describe "def. 1" $ specG phi1
  describe "def. 2" $ specG phi2
  describe "def. 3" $ specG phi3
  describe "def. 4" $ specG phi4

-- La verificación es
--    λ> verifica
--
--    8 examples, 0 failures

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_phi :: Positive Integer -> Bool
prop_phi (Positive n) =
  all (== phi1 n)
      [phi2 n,
       phi3 n,
       phi4 n]

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

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> phi1 (2*10^6)
--    800000
--    (2.49 secs, 2,117,853,856 bytes)
--    λ> phi2 (2*10^6)
--    800000
--    (0.02 secs, 565,664 bytes)
--
--    λ> length (show (phi2 (10^100000)))
--    100000
--    (2.80 secs, 5,110,043,208 bytes)
--    λ> length (show (phi3 (10^100000)))
--    100000
--    (4.81 secs, 7,249,353,896 bytes)
--    λ> length (show (phi4 (10^100000)))
--    100000
--    (0.78 secs, 1,467,573,768 bytes)

-- Verificación de la propiedad
-- ============================

-- La propiedad es
prop_phi2 :: Positive Integer -> Bool
prop_phi2 (Positive n) =
  genericLength (show (phi4 (10^n))) == n

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

from functools import reduce
from math import gcd
from operator import mul
from sys import set_int_max_str_digits
from timeit import Timer, default_timer

from hypothesis import given
from hypothesis import strategies as st
from sympy import factorint, totient

set_int_max_str_digits(10**6)

# 1ª solución
# ===========

def phi1(n: int) -> int:
    return len([x for x in range(1, n+1) if gcd(x, n) == 1])

# 2ª solución
# ===========

def producto(xs: list[int]) -> int:
    return reduce(mul, xs, 1)

def phi2(n: int) -> int:
    factores = factorint(n)
    return producto([(p-1)*p**(e-1) for p, e in factores.items()])

# 3ª solución
# =============

def phi3(n: int) -> int:
    return totient(n)

# Verificación
# ============

def test_phi() -> None:
    for phi in [phi1, phi2, phi3]:
        assert phi(36) == 12
        assert list(map(phi, range(10, 21))) == [4,10,4,12,6,8,8,16,6,18,8]
    print("Verificado")

# La verificación es
#    >>> test_phi()
#    Verificado

# Comprobación de equivalencia
# ============================

# La propiedad es
@given(st.integers(min_value=1, max_value=1000))
def test_phi_equiv(n: int) -> None:
    r = phi1(n)
    assert phi2(n) == r
    assert phi3(n) == r

# La comprobación es
#    >>> test_phi_equiv()
#    >>>

# Comparación de eficiencia
# =========================

def tiempo(e: str) -> None:
    """Tiempo (en segundos) de evaluar la expresión e."""
    t = Timer(e, "", default_timer, globals()).timeit(1)
    print(f"{t:0.2f} segundos")

# La comparación es
#    >>> tiempo('phi1(9*10**6)')
#    2.09 segundos
#    >>> tiempo('phi2(9*10**6)')
#    0.00 segundos
#    >>> tiempo('phi3(9*10**6)')
#    0.00 segundos
#
#    >>> tiempo('phi2(10**1000000)')
#    3.55 segundos
#    >>> tiempo('phi3(10**1000000)')
#    3.37 segundos

# Verificación de la propiedad
# ============================

# La propiedad es
@given(st.integers(min_value=1, max_value=1000))
def test_phi_prop(n: int) -> None:
    assert len(str(phi2(10**n))) == n

# La comprobación es
#    >>> test_phi_prop()
#    >>>

# Comprobación de todas las propiedades
# =====================================

# La comprobación es
#    src> poetry run pytest -v La_funcion_indicatriz_de_Euler.py
#    ===== test session starts =====
#       test_phi PASSED
#       test_phi_equiv PASSED
#       test_phi_prop PASSED
#    ===== passed in 0.55s =====