Esta semana he publicado en Exercitium las soluciones de los siguientes problemas:
- 1. Primos con cubos
- 2. Sumas alternas de factoriales
- 3. Potencias perfectas
- 4. Sucesión de suma de cuadrados de los dígitos
- 5. La función indicatriz de Euler
A continuación se muestran las soluciones.
1. Primos con cubos
Un primo con cubo es un número primo p para el que existe algún entero positivo n tal que la expresión n^3 + n^2p es un cubo perfecto. Por ejemplo, 19 es un primo con cubo ya que 8^3 + 8^2×19 = 12^3.
Definir la sucesión
primosConCubos :: [Integer] |
tal que sus elementos son los primos con cubo. Por ejemplo,
λ> take 6 primosConCubos [7,19,37,61,127,271] λ> length (takeWhile (< 1000000) primosConCubos) 173 |
Soluciones
import Data.Numbers.Primes (isPrime) import Test.QuickCheck (NonNegative (NonNegative), maxSize, quickCheckWith, stdArgs) -- 1ª solución -- =========== primosConCubos1 :: [Integer] primosConCubos1 = [p | x <- [1..], n <- [1..x], (x^3 - n^3) `mod` (n^2) == 0, let p = (x^3 - n^3) `div` (n^2), isPrime p] -- 2ª solución -- =========== -- Para analizar la respuesta, en esta solución se calculan los pares -- (p,n) tales que p es un primo con cubo y n es un número positivo tal -- que n^3 + n^2p es un cubo primosConCubos2' :: [(Integer,Integer)] primosConCubos2' = [(p,n) | x <- [1..], n <- [1..x], (x^3 - n^3) `mod` (n^2) == 0, let p = (x^3 - n^3) `div` (n^2), isPrime p] -- El cálculo es -- λ> take 7 primosConCubos2 -- [(7,1),(19,8),(37,27),(61,64),(127,216),(271,729),(331,1000)] -- Se observa que la sucesión de los segundos elementos [1,8,27,64,...] -- es la de los cubos y que los primeros elementos se obtienen restando -- los segundos elementos consecutivos; es decir, -- 7 = 8 - 1 = 2^3 - 1^3 -- 19 = 27 - 8 = 3^3 - 2^3 -- 37 = 64 - 27 = 4^3 - 3^3 -- Continuando el patrón, -- 61 = 5^3 - 4^3 = 125 - 64 -- 91 = 6^3 - 5^3 = 216 - 125 -- 127 = 7^3 - 6^3 = 343 - 216 -- 271 = 10^3 - 9^3 = 1000 - 729 -- 331 = 11^3 -10^3 = 1331 - 1000 -- Por tanto, los primos con cubos son diferencias de dos cubos -- consecutivos; es decir, coinciden con los números cubanos del -- ejercicio anterior. A partir de la conjetura se obtienen las -- siguientes definiciones -- Basado en las anteriores observaciones se obtiene la siguiente -- definición primosConCubos2 :: [Integer] primosConCubos2 = filter isPrime [(x+1)^3 - x^3 | x <- [1..]] -- 3ª definición -- ============= primosConCubos3 :: [Integer] primosConCubos3 = filter isPrime diferenciasCubosConsecutivos diferenciasCubosConsecutivos :: [Integer] diferenciasCubosConsecutivos = zipWith (-) (tail cubos) cubos cubos :: [Integer] cubos = map (^3) [0..] -- 4ª solución -- =========== -- Simplificando la expresión (x+1)^3 - x^3 se obtiene 3*x^2+ 3*x + 1, -- con lo que la 3ª definición se reduce a primosConCubos4 :: [Integer] primosConCubos4 = [p | x <- [1..], let p = 3*x^2+ 3*x + 1, isPrime p] -- Comprobación de equivalencia -- ============================ -- La propiedad es prop_primosConCubos :: NonNegative Int -> Bool prop_primosConCubos (NonNegative n) = all (== primosConCubos1 !! n) [primosConCubos2 !! n, primosConCubos3 !! n, primosConCubos4 !! n] -- La comprobación es -- λ> quickCheckWith (stdArgs {maxSize=10}) prop_primosConCubos -- +++ OK, passed 100 tests. -- Comparación de eficiencia -- ========================= -- La comparación es -- λ> primosConCubos1 !! 6 -- 331 -- (1.15 secs, 1,105,612,336 bytes) -- λ> primosConCubos1 !! 7 -- 397 -- (1.96 secs, 1,909,369,592 bytes) -- λ> primosConCubos2 !! 7 -- 397 -- -- (0.01 secs, 648,840 bytes) -- λ> primosConCubos2 !! (10^3) -- 65580901 -- (0.53 secs, 1,726,837,688 bytes) -- λ> primosConCubos3 !! (10^3) -- 65580901 -- (0.49 secs, 1,724,258,632 bytes) -- λ> primosConCubos4 !! (10^3) -- 65580901 -- (0.47 secs, 1,724,833,992 bytes) -- Demostración de la conjetura -- ============================ -- Vamos a demostrar que los primos con cubos son diferencias de dos -- cubos consecutivos. -- -- Sea p un primo con cubo. Por la definición, existe un entero -- positivo n tal que n³ + n²p es un cubo. -- -- Lema 1: Los números n y p son coprimos (es decir, mcd(n,p) = 1). -- Dem.: En caso contrario, puesto que p es primo, existe un a tal que -- n = ap. Luego n³ + n²p = (a³+a²)p³ es un cubo y, por tanto, -- a³+a² es un cubo lo que es imposible ya que el siguiente cubo de -- a³ es a³+3a²+3a+1. -- -- Lema 2: Los números n² y n+p son coprimos. -- Dem.: Sea k = mcd(n^2,n+p). Por k divide n², luego k divide a n; -- además, k divide a n+p y (usando el lema 1 y el ser p primo), se -- tiene que k = 1. -- -- Puesto que n³+n²p = n²(n+p) es un cubo, usando el lema 2, se tiene -- que n² y n+p son cubos y, por serlo n², n también es un cubo. Es -- decir, existen enteros positivos x e y tales que n = x³ y -- n+p = y³. Por tanto, p = y³-x³. Sea k = y-x. Se tiene que k = 1 ya -- que -- p = y³-x³ -- = (n+k)³-n³ -- = 3k+3k²+k³ -- no es primo para k > 1. -- -- Por consiguiente, p = (x+1)³-x³. |
El código se encuentra en GitHub.
2. Sumas alternas de factoriales
Las primeras sumas alternas de los factoriales son números primos; en efecto,
3! - 2! + 1! = 5 4! - 3! + 2! - 1! = 19 5! - 4! + 3! - 2! + 1! = 101 6! - 5! + 4! - 3! + 2! - 1! = 619 7! - 6! + 5! - 4! + 3! - 2! + 1! = 4421 8! - 7! + 6! - 5! + 4! - 3! + 2! - 1! = 35899 |
son primos, pero
9! - 8! + 7! - 6! + 5! - 4! + 3! - 2! + 1! = 326981 |
no es primo.
Definir las funciones
sumaAlterna :: Integer -> Integer sumasAlternas :: [Integer] conSumaAlternaPrima :: [Integer] |
tales que
- (sumaAlterna n) es la suma alterna de los factoriales desde n hasta 1. Por ejemplo,
sumaAlterna 3 == 5
sumaAlterna 4 == 19
sumaAlterna 5 == 101
sumaAlterna 6 == 619
sumaAlterna 7 == 4421
sumaAlterna 8 == 35899
sumaAlterna 9 == 326981
sumaAlterna (5*10^4) `mod` (10^6) == 577019 |
- sumasAlternas es la sucesión de las sumas alternas de factoriales. Por ejemplo,
λ> take 10 sumasAlternas1
[0,1,1,5,19,101,619,4421,35899,326981] |
- conSumaAlternaPrima es la sucesión de los números cuya suma alterna de factoriales es prima. Por ejemplo,
λ> take 8 conSumaAlternaPrima
[3,4,5,6,7,8,10,15] |
Soluciones
import Data.List (genericTake) import Data.Numbers.Primes (isPrime) import Test.QuickCheck -- 1ª definición de sumaAlterna -- ============================ sumaAlterna1 :: Integer -> Integer sumaAlterna1 1 = 1 sumaAlterna1 n = factorial n - sumaAlterna1 (n-1) factorial :: Integer -> Integer factorial n = product [1..n] -- 2ª definición de sumaAlterna -- ============================ sumaAlterna2 :: Integer -> Integer sumaAlterna2 n = sum (genericTake n (zipWith (*) signos (tail factoriales))) where signos | odd n = concat (repeat [1,-1]) | otherwise = concat (repeat [-1,1]) -- factoriales es la lista de los factoriales. Por ejemplo, -- take 7 factoriales == [1,1,2,6,24,120,720] factoriales :: [Integer] factoriales = 1 : scanl1 (*) [1..] -- 3ª definición de sumaAlterna -- ============================ sumaAlterna3 :: Integer -> Integer sumaAlterna3 n = sum (genericTake n (zipWith (*) signos (tail factoriales))) where signos | odd n = cycle [1,-1] | otherwise = cycle [-1,1] -- 3ª definición de sumaAlterna -- ============================ sumaAlterna4 :: Integer -> Integer sumaAlterna4 n = foldl (flip (-)) 0 (scanl1 (*) [1..n]) -- Comprobación de equivalencia de sumaAlterna -- =========================================== -- La propiedad es prop_sumaAlterna :: Positive Integer -> Bool prop_sumaAlterna (Positive n) = all (== sumaAlterna1 n) [sumaAlterna2 n, sumaAlterna3 n, sumaAlterna4 n] -- La comprobación es -- λ> quickCheck prop_sumaAlterna -- +++ OK, passed 100 tests. -- Comparación de eficiencia de sumaAlterna -- ======================================== -- La comparación es -- λ> sumaAlterna1 4000 `mod` (10^6) -- 577019 -- (6.21 secs, 16,154,113,192 bytes) -- λ> sumaAlterna2 4000 `mod` (10^6) -- 577019 -- (0.01 secs, 24,844,664 bytes) -- -- λ> sumaAlterna2 (5*10^4) `mod` (10^6) -- 577019 -- (1.81 secs, 4,729,583,864 bytes) -- λ> sumaAlterna3 (5*10^4) `mod` (10^6) -- 577019 -- (0.89 secs, 4,725,983,928 bytes) -- λ> sumaAlterna4 (5*10^4) `mod` (10^6) -- 577019 -- (0.70 secs, 4,710,770,592 bytes) -- En lo que sigue se usa la 3ª definición sumaAlterna :: Integer -> Integer sumaAlterna = sumaAlterna3 -- 1ª definición de sumasAlternas -- ============================== sumasAlternas1 :: [Integer] sumasAlternas1 = map sumaAlterna [0..] -- 2ª definición de sumasAlternas -- ============================== sumasAlternas2 :: [Integer] sumasAlternas2 = 0 : zipWith (-) (tail factoriales) sumasAlternas2 -- 3ª definición de sumasAlternas -- ============================== sumasAlternas3 :: [Integer] sumasAlternas3 = scanl (flip (-)) 0 $ scanl1 (*) [1..] -- Comprobación de equivalencia -- ============================ -- La propiedad es prop_sumasAlternas :: NonNegative Int -> Bool prop_sumasAlternas (NonNegative n) = all (== sumasAlternas1 !! n) [sumasAlternas2 !! n, sumasAlternas3 !! n] -- La comprobación es -- λ> quickCheck prop_sumasAlternas -- +++ OK, passed 100 tests. -- Comparación de eficiencia -- ========================= -- La comparación es -- λ> length (show (sumasAlternas1 !! (5*10^4))) -- 213237 -- (4.90 secs, 4,731,620,600 bytes) -- λ> length (show (sumasAlternas2 !! (5*10^4))) -- 213237 -- (2.39 secs, 4,726,820,456 bytes) -- λ> length (show (sumasAlternas3 !! (5*10^4))) -- 213237 -- (1.78 secs, 4,726,820,280 bytes) -- 1ª definición de conSumaAlternaPrima -- ==================================== conSumaAlternaPrima1 :: [Integer] conSumaAlternaPrima1 = [n | n <- [0..], isPrime (sumaAlterna n)] -- 2ª definición de conSumaAlternaPrima -- ==================================== conSumaAlternaPrima2 :: [Integer] conSumaAlternaPrima2 = [x | (x,y) <- zip [0..] sumasAlternas2, isPrime y] -- 3ª definición de conSumaAlternaPrima -- ==================================== conSumaAlternaPrima3 :: [Integer] conSumaAlternaPrima3 = filter (isPrime . sumaAlterna) [0..] -- Comprobación de equivalencia de conSumaAlternaPrima -- =================================================== -- La propiedad es prop_conSumaAlternaPrima :: NonNegative Int -> Bool prop_conSumaAlternaPrima (NonNegative n) = all (== conSumaAlternaPrima1 !! n) [conSumaAlternaPrima2 !! n, conSumaAlternaPrima3 !! n] -- La comprobación es -- λ> quickCheckWith (stdArgs {maxSize=5}) prop_conSumaAlternaPrima -- +++ OK, passed 100 tests. |
El código se encuentra en GitHub.
3. Potencias perfectas
Un número natural n es una potencia perfecta si existen dos números naturales m > 1 y k > 1 tales que n = m^k. Las primeras potencias perfectas son
4 = 2², 8 = 2³, 9 = 3², 16 = 2⁴, 25 = 5², 27 = 3³, 32 = 2⁵, 36 = 6², 49 = 7², 64 = 2⁶, ... |
Definir la sucesión
potenciasPerfectas :: [Integer] |
cuyos términos son las potencias perfectas. Por ejemplo,
take 10 potenciasPerfectas == [4,8,9,16,25,27,32,36,49,64] potenciasPerfectas !! 3000 == 7778521 |
Definir el procedimiento
grafica :: Int -> IO () |
tal que (grafica n) es la representación gráfica de las n primeras potencias perfectas. Por ejemplo, para (grafica 30) dibuja
Soluciones
import Data.List (group) import Data.Numbers.Primes (primeFactors) import Graphics.Gnuplot.Simple (Attribute (Key, PNG), plotList) import Test.QuickCheck (NonNegative (NonNegative), quickCheck) -- 1ª solución -- =========== potenciasPerfectas1 :: [Integer] potenciasPerfectas1 = filter esPotenciaPerfecta [4..] -- (esPotenciaPerfecta x) se verifica si x es una potencia perfecta. Por -- ejemplo, -- esPotenciaPerfecta 36 == True -- esPotenciaPerfecta 72 == False esPotenciaPerfecta :: Integer -> Bool esPotenciaPerfecta = not . null. potenciasPerfectasDe -- (potenciasPerfectasDe x) es la lista de pares (a,b) tales que -- x = a^b. Por ejemplo, -- potenciasPerfectasDe 64 == [(2,6),(4,3),(8,2)] -- potenciasPerfectasDe 72 == [] potenciasPerfectasDe :: Integer -> [(Integer,Integer)] potenciasPerfectasDe n = [(m,k) | m <- takeWhile (\x -> x*x <= n) [2..] , k <- takeWhile (\x -> m^x <= n) [2..] , m^k == n] -- 2ª solución -- =========== potenciasPerfectas2 :: [Integer] potenciasPerfectas2 = [x | x <- [4..], esPotenciaPerfecta2 x] -- (esPotenciaPerfecta2 x) se verifica si x es una potencia perfecta. Por -- ejemplo, -- esPotenciaPerfecta2 36 == True -- esPotenciaPerfecta2 72 == False esPotenciaPerfecta2 :: Integer -> Bool esPotenciaPerfecta2 x = mcd (exponentes x) > 1 -- (exponentes x) es la lista de los exponentes de l factorización prima -- de x. Por ejemplos, -- exponentes 36 == [2,2] -- exponentes 72 == [3,2] exponentes :: Integer -> [Int] exponentes x = [length ys | ys <- group (primeFactors x)] -- (mcd xs) es el máximo común divisor de la lista xs. Por ejemplo, -- mcd [4,6,10] == 2 -- mcd [4,5,10] == 1 mcd :: [Int] -> Int mcd = foldl1 gcd -- 3ª solución -- =========== potenciasPerfectas3 :: [Integer] potenciasPerfectas3 = mezclaTodas potencias -- potencias es la lista las listas de potencias de todos los números -- mayores que 1 con exponentes mayores que 1. Por ejemplo, -- λ> map (take 3) (take 4 potencias) -- [[4,8,16],[9,27,81],[16,64,256],[25,125,625]] potencias:: [[Integer]] potencias = [[n^k | k <- [2..]] | n <- [2..]] -- (mezclaTodas xss) es la mezcla ordenada sin repeticiones de las -- listas ordenadas xss. Por ejemplo, -- take 7 (mezclaTodas potencias) == [4,8,9,16,25,27,32] mezclaTodas :: Ord a => [[a]] -> [a] mezclaTodas = foldr1 xmezcla where xmezcla (x:xs) ys = x : mezcla xs ys -- (mezcla xs ys) es la mezcla ordenada sin repeticiones de las -- listas ordenadas xs e ys. Por ejemplo, -- take 7 (mezcla [2,5..] [4,6..]) == [2,4,5,6,8,10,11] mezcla :: Ord a => [a] -> [a] -> [a] mezcla (x:xs) (y:ys) | x < y = x : mezcla xs (y:ys) | x == y = x : mezcla xs ys | x > y = y : mezcla (x:xs) ys -- Comprobación de equivalencia -- ============================ -- La propiedad es prop_potenciasPerfectas :: NonNegative Int -> Bool prop_potenciasPerfectas (NonNegative n) = all (== potenciasPerfectas1 !! n) [potenciasPerfectas2 !! n, potenciasPerfectas3 !! n] -- La comprobación es -- λ> quickCheck prop_potenciasPerfectas -- +++ OK, passed 100 tests. -- Comparación de eficiencia -- ========================= -- La comparación es -- λ> potenciasPerfectas1 !! 200 -- 28224 -- (10.56 secs, 8,434,647,368 bytes) -- λ> potenciasPerfectas2 !! 200 -- 28224 -- (0.36 secs, 825,040,416 bytes) -- λ> potenciasPerfectas3 !! 200 -- 28224 -- (0.05 secs, 7,474,280 bytes) -- -- λ> potenciasPerfectas2 !! 500 -- 191844 -- (4.16 secs, 9,899,367,112 bytes) -- λ> potenciasPerfectas3 !! 500 -- 191844 -- (0.09 secs, 51,275,464 bytes) -- En lo que sigue se usa la 3ª solución potenciasPerfectas :: [Integer] potenciasPerfectas = potenciasPerfectas3 -- Representación gráfica -- ====================== grafica :: Int -> IO () grafica n = plotList [ Key Nothing , PNG "Potencias_perfectas.png" ] (take n potenciasPerfectas) |
El código se encuentra en GitHub.
4. Sucesión de suma de cuadrados de los dígitos
Definir la función
sucSumaCuadradosDigitos :: Integer -> [Integer] |
tal que (sucSumaCuadradosDigitos n) es la sucesión cuyo primer término es n y los restantes se obtienen sumando los cuadrados de los dígitos de su término anterior. Por ejemplo,
λ> take 20 (sucSumaCuadradosDigitos1 2000) [2000,4,16,37,58,89,145,42,20,4,16,37,58,89,145,42,20,4,16,37] λ> take 20 (sucSumaCuadradosDigitos 1976) [1976,167,86,100,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1] λ> sucSumaCuadradosDigitos 2000 !! (10^9) 20 |
Soluciones
import Test.QuickCheck (Positive (Positive), NonNegative (NonNegative), quickCheck) -- 1ª solución -- =========== sucSumaCuadradosDigitos1 :: Integer -> [Integer] sucSumaCuadradosDigitos1 n = n : [sumaCuadradosDigitos x | x <- sucSumaCuadradosDigitos1 n] -- (sumaCuadradosDigitos n) es la suma de los cuadrados de los dígitos -- de n. Por ejemplo, -- sumaCuadradosDigitos 2016 == 41 sumaCuadradosDigitos :: Integer -> Integer sumaCuadradosDigitos n = sum (map (^2) (digitos n)) -- (digitos n) es la lista de los dígitos de n. Por ejemplo, -- digitos 325 == [3,2,5] digitos :: Integer -> [Integer] digitos n = [read [d] | d <- show n] -- 2ª solución -- =========== sucSumaCuadradosDigitos2 :: Integer -> [Integer] sucSumaCuadradosDigitos2 = iterate sumaCuadradosDigitos -- 3ª solución -- =========== -- A partir de los cálculos con las definiciones anteriores, se observa -- que para todo n (sucSumaCuadradosDigitos n) tiene una parte pura y -- otra periódica. Por ejemplo, para n = 2016, -- λ> take 20 (sucSumaCuadradosDigitos 2016) -- [2016,41,17,50,25,29,85,89,145,42,20,4,16,37,58,89,145,42,20,4] -- la parte pura es -- [2016,41,17,50,25,29,85] -- y la parte periódica es -- [89,145,42,20,4,16,37,58]) sucSumaCuadradosDigitos3 :: Integer -> [Integer] sucSumaCuadradosDigitos3 n = xs ++ cycle ys where (xs,ys) = sucCompactaSumaCuadradosDigitos n -- (sucCompactaSumaCuadradosDigitos n) es el par formado por la parte -- pura y la periódica de (sucSumaCuadradosDigitos n). Por ejemplo, -- λ> sucCompactaSumaCuadradosDigitos 2016 -- ([2016,41,17,50,25,29,85],[89,145,42,20,4,16,37,58]) -- λ> sucCompactaSumaCuadradosDigitos 1976 -- ([1976,167,86,100],[1]) sucCompactaSumaCuadradosDigitos :: Integer -> ([Integer],[Integer]) sucCompactaSumaCuadradosDigitos = partePuraPeriodica . sucSumaCuadradosDigitos1 -- (partePuraPeriodica xs) es el par formado por la parte pura y la -- periódica de xs. Por ejemplo, -- λ> partePuraPeriodica (sucSumaCuadradosDigitos 2016) -- ([2016,41,17,50,25,29,85],[89,145,42,20,4,16,37,58]) -- λ> partePuraPeriodica (sucSumaCuadradosDigitos 1976) -- ([1976,167,86,100],[1]) partePuraPeriodica :: [Integer] -> ([Integer],[Integer]) partePuraPeriodica = aux [] where aux as (b:bs) | b `elem` as = span (/=b) (reverse as) | otherwise = aux (b:as) bs -- 4ª solución -- =========== sucSumaCuadradosDigitos4 :: Integer -> [Integer] sucSumaCuadradosDigitos4 1 = repeat 1 sucSumaCuadradosDigitos4 89 = cycle [89,145,42,20,4,16,37,58] sucSumaCuadradosDigitos4 n = n : sucSumaCuadradosDigitos4 (sumaCuadradosDigitos n) -- Comprobación de equivalencia -- ============================ -- La propiedad es prop_sucSumaCuadradosDigitos :: Positive Integer -> NonNegative Int -> Bool prop_sucSumaCuadradosDigitos (Positive n) (NonNegative k) = all (== sucSumaCuadradosDigitos1 n !! k) [sucSumaCuadradosDigitos2 n !! k, sucSumaCuadradosDigitos3 n !! k, sucSumaCuadradosDigitos4 n !! k] -- La comprobación es -- λ> quickCheck prop_sucSumaCuadradosDigitos -- +++ OK, passed 100 tests. -- Comparación de eficiencia -- ========================= -- La comparación es -- λ> sucSumaCuadradosDigitos1 2000 !! (10^4) -- 20 -- (6.96 secs, 8,049,886,312 bytes) -- λ> sucSumaCuadradosDigitos2 2000 !! (10^4) -- 20 -- (0.08 secs, 91,024,688 bytes) -- λ> sucSumaCuadradosDigitos3 2000 !! (10^4) -- 20 -- (0.01 secs, 995,560 bytes) -- λ> sucSumaCuadradosDigitos4 2000 !! (10^4) -- 20 -- (0.02 secs, 587,040 bytes) -- -- λ> sucSumaCuadradosDigitos2 2000 !! (3*10^5) -- 20 -- (1.96 secs, 2,715,501,416 bytes) -- λ> sucSumaCuadradosDigitos3 2000 !! (3*10^5) -- 20 -- (0.02 secs, 995,872 bytes) -- λ> sucSumaCuadradosDigitos4 2000 !! (3*10^5) -- 20 -- (0.02 secs, 587,352 bytes) -- -- λ> sucSumaCuadradosDigitos3 2000 !! (10^9) -- 20 -- (2.85 secs, 996,016 bytes) -- λ> sucSumaCuadradosDigitos4 2000 !! (10^9) -- 20 -- (2.54 secs, 587,496 bytes) |
El código se encuentra en GitHub.
5. La función indicatriz de Euler
La indicatriz de Euler (también función φ de Euler) es una función importante en teoría de números. Si n es un entero positivo, entonces φ(n) se define como el número de enteros positivos menores o iguales a n y coprimos con n. Por ejemplo, φ(36) = 12 ya que los números menores o iguales a 36 y coprimos con 36 son doce: 1, 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, y 35.
Definir la función
phi :: Integer -> Integer |
tal que (phi n) es igual a φ(n). Por ejemplo,
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 |
Comprobar con QuickCheck que, para todo n > 0, φ(10^n) tiene n dígitos.
Soluciones
import Data.List (genericLength, group) import Data.Numbers.Primes (primeFactors) import Math.NumberTheory.ArithmeticFunctions (totient) import Test.QuickCheck (Positive (Positive), quickCheck) -- 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 -- 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. |
El código se encuentra en GitHub.