PFH: La semana en Exercitium (22 de julio de 2022)
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
|
1 |
primosConCubos :: [Integer] |
tal que sus elementos son los primos con cubo. Por ejemplo,
|
1 2 3 4 |
λ> take 6 primosConCubos [7,19,37,61,127,271] λ> length (takeWhile (< 1000000) primosConCubos) 173 |
Soluciones
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 |
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,
|
1 2 3 4 5 6 |
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
|
1 |
9! - 8! + 7! - 6! + 5! - 4! + 3! - 2! + 1! = 326981 |
no es primo.
Definir las funciones
|
1 2 3 |
sumaAlterna :: Integer -> Integer sumasAlternas :: [Integer] conSumaAlternaPrima :: [Integer] |
tales que
- (sumaAlterna n) es la suma alterna de los factoriales desde n hasta 1. Por ejemplo,
|
1 2 3 4 5 6 7 8 |
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,
|
1 2 |
λ> 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,
|
1 2 |
λ> take 8 conSumaAlternaPrima [3,4,5,6,7,8,10,15] |
Soluciones
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 |
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
|
1 2 |
4 = 2², 8 = 2³, 9 = 3², 16 = 2⁴, 25 = 5², 27 = 3³, 32 = 2⁵, 36 = 6², 49 = 7², 64 = 2⁶, ... |
Definir la sucesión
|
1 |
potenciasPerfectas :: [Integer] |
cuyos términos son las potencias perfectas. Por ejemplo,
|
1 2 |
take 10 potenciasPerfectas == [4,8,9,16,25,27,32,36,49,64] potenciasPerfectas !! 3000 == 7778521 |
Definir el procedimiento
|
1 |
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
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 |
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
|
1 |
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,
|
1 2 3 4 5 6 |
λ> 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
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 |
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
|
1 |
phi :: Integer -> Integer |
tal que (phi n) es igual a φ(n). Por ejemplo,
|
1 2 3 4 |
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
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 |
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.