PFH: La semana en Exercitium (25 de junio de 2022)
Esta semana he publicado en Exercitium las soluciones de los siguientes problemas:
- 1. Menor número con una cantidad dada de divisores
- 2. Cálculo aproximado de integrales definidas
- 3. Cálculo de la suma 11! + 22! + 33! + … + nn!
- 4. Números para los que mcm(1,2,…n-1) = mcm(1,2,…,n)
- 5. Método de bisección para aproximar raíces de funciones
A continuación se muestran las soluciones.
1. Menor número con una cantidad dada de divisores
El menor número con 2 divisores es el 2, ya que tiene 2 divisores (el 1 y el 2) y el anterior al 2 (el 1) sólo tiene 1 divisor (el 1).
El menor número con 4 divisores es el 6, ya que tiene 4 divisores (el 1, 2, 3 y 6) y sus anteriores (el 1, 2, 3, 4 y 5) tienen menos de 4 divisores (tienen 1, 1, 1, 3 y 1, respectivamente).
El menor número con 8 divisores es el 24, ya que tiene 8 divisores (el 1, 2, 3, 4, 6, 8, 12 y 24) y sus anteriores (del 1 al 23) tienen menos de 8 divisores.
El menor número con 16 divisores es el 120, ya que tiene 16 divisores (el 1, 2, 3, 4, 5, 6, 8, 10, 12, 15, 20, 24, 30, 40, 60 y 120) y sus anteriores (del 1 al 119) tienen menos de 16 divisores.
Definir la función
|
1 |
menor :: Integer -> Integer |
tal que (menor n) es el menor número con 2^n divisores. Por ejemplo,
|
1 2 3 4 5 |
menor 1 == 2 menor 2 == 6 menor 3 == 24 menor 4 == 120 length (show (menor (4*10^4))) == 207945 |
Comprobar con QuickCheck que, para todo k >=0, (menor (2^k)) es un divisor de (menor (2^(k+1))).
Nota: Este ejercicio está basado en el problema N1 de la Olimpíada Internacional de Matemáticas (IMO) del 2011.
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 |
import Data.List (genericLength, genericTake, group) import Data.Numbers.Primes (primeFactors, primes) import Test.QuickCheck (Positive (Positive), maxSize, stdArgs, quickCheck, quickCheckWith) -- 1ª solución -- =========== menor1 :: Integer -> Integer menor1 n = head [x | x <- [1..], numeroDivisores x == 2^n] -- (numeroDivisores n) es el número de divisores de n. Por ejemplo, -- numeroDivisores 12 == 6 numeroDivisores :: Integer -> Integer numeroDivisores = genericLength . divisores -- (divisores x) es la lista de los divisores de x. Por ejemplo, -- divisores 12 == [1,3,2,6,4,12] -- divisores 25 == [1,5,25] divisores :: Integer -> [Integer] divisores n = [x | x <- [1..n], n `rem` x == 0] -- 2ª solución -- =========== menor2 :: Integer -> Integer menor2 n = head [x | x <- [1..], numeroDivisores2 x == 2^n] numeroDivisores2 :: Integer -> Integer numeroDivisores2 = product . map ((+1) . genericLength) . group . primeFactors -- 3ª solución -- =========== menor3 :: Integer -> Integer menor3 n = product (genericTake n potencias) -- potencias es la sucesión de las potencias de la forma p^(2^k), -- donde p es un número primo y k es un número natural, ordenadas de -- menor a mayor. Por ejemplo, -- take 14 potencias == [2,3,4,5,7,9,11,13,16,17,19,23,25,29] potencias :: [Integer] potencias = 2 : mezcla (tail primes) (map (^2) potencias) -- (mezcla xs ys) es la lista obtenida mezclando las dos listas xs e ys, -- que se suponen ordenadas y disjuntas. Por ejemplo, -- λ> take 15 (mezcla [2^n | n <- [1..]] [3^n | n <- [1..]]) -- [2,3,4,8,9,16,27,32,64,81,128,243,256,512,729] mezcla :: Ord a => [a] -> [a] -> [a] mezcla (x:xs) (y:ys) | x < y = x : mezcla xs (y:ys) | x > y = y : mezcla (x:xs) ys -- Comprobación de equivalencia -- ============================ -- La propiedad es prop_menor :: Positive Integer -> Bool prop_menor (Positive n) = all (== menor1 n) [menor2 n, menor3 n] -- La comprobación es -- λ> quickCheckWith (stdArgs {maxSize=5}) prop_menor -- +++ OK, passed 100 tests. -- Comparación de eficiencia -- ========================= -- La comparación es -- λ> menor 8 -- 1081080 -- (47.69 secs, 94,764,856,352 bytes) -- λ> menor2 8 -- 1081080 -- (36.17 secs, 94,764,856,368 bytes) -- λ> menor3 8 -- 1081080 -- (0.00 secs, 116,960 bytes) -- Definición de menor -- =================== -- En lo que sigue, usaremos menor3 como menor. menor :: Integer -> Integer menor = menor3 -- Propiedad -- ========= -- La propiedad es prop_menor_divide :: Positive Integer -> Bool prop_menor_divide (Positive n) = menor (n+1) `mod` menor n == 0 -- La comprobación es -- λ> quickCheck prop_menor_divide -- +++ OK, passed 100 tests. |
El código se encuentra en GitHub.
2. Cálculo aproximado de integrales definidas
La integral definida de una función f entre los límites a y b puede calcularse mediante la regla del rectángulo usando la fórmula
|
1 |
h * (f(a+h/2) + f(a+h+h/2) + f(a+2h+h/2) + ... + f(a+n*h+h/2)) |
con a+n*h+h/2 <= b < a+(n+1)*h+h/2 y usando valores pequeños para h.
Definir la función
|
1 |
integral :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a |
tal que (integral a b f h) es el valor de dicha expresión. Por ejemplo, el cálculo de la integral de f(x) = x^3 entre 0 y 1, con paso 0.01, es
|
1 |
integral 0 1 (^3) 0.01 == 0.24998750000000042 |
Otros ejemplos son
|
1 2 3 4 |
integral 0 1 (^4) 0.01 == 0.19998333362500048 integral 0 1 (\x -> 3*x^2 + 4*x^3) 0.01 == 1.9999250000000026 log 2 - integral 1 2 (\x -> 1/x) 0.01 == 3.124931644782336e-6 pi - 4 * integral 0 1 (\x -> 1/(x^2+1)) 0.01 == -8.333333331389525e-6 |
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 |
import Test.QuickCheck.HigherOrder (quickCheck') import Test.QuickCheck (Property, (==>), quickCheck) -- 1ª solución -- =========== integral1 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a integral1 a b f h | a+h/2 > b = 0 | otherwise = h * f (a+h/2) + integral1 (a+h) b f h -- 2ª solución -- =========== integral2 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a integral2 a b f h = aux a where aux x | x+h/2 > b = 0 | otherwise = h * f (x+h/2) + aux (x+h) -- 3ª solución -- =========== integral3 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a integral3 a b f h = h * suma (a+h/2) b (+h) f -- (suma a b s f) es l valor de -- f(a) + f(s(a)) + f(s(s(a)) + ... + f(s(...(s(a))...)) -- hasta que s(s(...(s(a))...)) > b. Por ejemplo, -- suma 2 5 (1+) (^3) == 224 suma :: (Ord t, Num a) => t -> t -> (t -> t) -> (t -> a) -> a suma a b s f = sum [f x | x <- sucesion a b s] -- (sucesion x y s) es la lista -- [a, s(a), s(s(a), ..., s(...(s(a))...)] -- hasta que s(s(...(s(a))...)) > b. Por ejemplo, -- sucesion 3 20 (+2) == [3,5,7,9,11,13,15,17,19] sucesion :: Ord a => a -> a -> (a -> a) -> [a] sucesion a b s = takeWhile (<=b) (iterate s a) -- Comprobación de equivalencia -- ============================ -- La propiedad es prop_integral :: Int -> Int -> (Int -> Int) -> Int -> Property prop_integral a b f h = a < b && h > 0 ==> all (=~ integral1 a' b' f' h') [integral2 a' b' f' h', integral3 a' b' f' h'] where a' = fromIntegral a b' = fromIntegral b h' = fromIntegral h f' = fromIntegral . f. round x =~ y = abs (x - y) < 0.001 -- La comprobación es -- λ> quickCheck' prop_integral -- +++ OK, passed 100 tests; 385 discarded. -- Comparación de eficiencia -- ========================= -- La comparación es -- λ> integral1 0 10 (^3) 0.00001 -- 2499.999999881125 -- (2.63 secs, 1,491,006,744 bytes) -- λ> integral2 0 10 (^3) 0.00001 -- 2499.999999881125 -- (1.93 secs, 1,419,006,696 bytes) -- λ> integral3 0 10 (^3) 0.00001 -- 2499.9999998811422 -- (1.28 secs, 817,772,216 bytes) |
El código se encuentra en GitHub.
3. Cálculo de la suma 11! + 22! + 33! + … + nn!
Definir la función
|
1 |
suma :: Integer -> Integer |
tal que (suma n) es la suma 1·1! + 2·2! + 3·3! + ... + n·n!. Por ejemplo,
|
1 2 3 4 5 6 |
suma 1 == 1 suma 2 == 5 suma 3 == 23 suma 4 == 119 suma 5 == 719 take 9 (show (suma 70000)) == "823780458" |
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 |
import Test.QuickCheck (Positive (Positive), quickCheck) -- 1ª solución -- =========== suma1 :: Integer -> Integer suma1 n = sum [k * factorial k | k <- [1..n]] factorial :: Integer -> Integer factorial n = product [1..n] -- 2ª solución -- =========== suma2 :: Integer -> Integer suma2 n = sum (zipWith (*) [1..n] factoriales) factoriales :: [Integer] factoriales = scanl (*) 1 [2..] -- 3ª solución -- =========== -- Basada en los siguientes cálculos -- λ> [suma1 n | n <- [0..10]] -- [0,1,5,23,119,719,5039,40319,362879,3628799,39916799] -- λ> [factorial n | n <- [0..10]] -- [1,1,2,6,24,120,720,5040,40320,362880,3628800] -- λ> [factorial n | n <- [1..11]] -- [1,2,6,24,120,720,5040,40320,362880,3628800,39916800] -- λ> [factorial n - 1 | n <- [1..11]] -- [0,1,5,23,119,719,5039,40319,362879,3628799,39916799] suma3 :: Integer -> Integer suma3 n = factorial (n+1) - 1 -- Comprobación de equivalencia -- ============================ -- La propiedad es prop_suma :: Positive Integer -> Bool prop_suma (Positive n) = all (== suma1 n) [suma2 n, suma3 n] -- La comprobación es -- λ> quickCheck prop_suma -- +++ OK, passed 100 tests. -- Comparación de eficiencia -- ========================= -- La comparación es -- λ> take 5 (show (suma1 4000)) -- "73170" -- (5.04 secs, 16,225,195,448 bytes) -- λ> take 5 (show (suma2 4000)) -- "73170" -- (0.08 secs, 35,862,152 bytes) -- λ> take 5 (show (suma3 4000)) -- "73170" -- (0.01 secs, 12,896,968 bytes) -- -- -- λ> take 5 (show (suma2 40000)) -- "83669" -- (1.82 secs, 4,549,612,264 bytes) -- λ> take 5 (show (suma3 40000)) -- "83669" -- (0.24 secs, 1,620,976,984 bytes) |
El código se encuentra en GitHub.
4. Números para los que mcm(1,2,…n-1) = mcm(1,2,…,n)
Un número n es especial si mcm(1,2,…,n-1) = mcm(1,2,…,n). Por ejemplo, el 6 es especial ya que
|
1 |
mcm(1,2,3,4,5) = 60 = mcm(1,2,3,4,5,6) |
Definir la sucesión
|
1 |
especiales :: [Integer] |
cuyos términos son los números especiales. Por ejemplo,
|
1 2 3 4 5 |
take 10 especiales == [1,6,10,12,14,15,18,20,21,22] especiales !! 50 == 84 especiales !! 500 == 638 especiales !! 5000 == 5806 especiales !! 50000 == 55746 |
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 |
import Test.QuickCheck (NonNegative(NonNegative), quickCheck) -- 1ª solución -- =========== especiales1 :: [Integer] especiales1 = filter especial1 [1..] especial1 :: Integer -> Bool especial1 n = mcm1 [1..n-1] == mcm1 [1..n] mcm1 :: [Integer] -> Integer mcm1 [] = 1 mcm1 (x:xs) = lcm x (mcm1 xs) -- 2ª solución -- =========== especiales2 :: [Integer] especiales2 = filter especial2 [1..] especial2 :: Integer -> Bool especial2 n = mcm2 [1..n-1] == mcm2 [1..n] mcm2 :: [Integer] -> Integer mcm2 = foldr lcm 1 -- 3ª solución -- =========== especiales3 :: [Integer] especiales3 = [n | ((n,x),(_,y)) <- zip mcms (tail mcms) , x == y] mcms :: [(Integer,Integer)] mcms = zip [1..] (scanl lcm 1 [1..]) -- Comprobación de equivalencia -- ============================ -- La propiedad es prop_especiales :: NonNegative Int -> Bool prop_especiales (NonNegative n) = all (== especiales1 !! n) [especiales2 !! n, especiales3 !! n] -- La comprobación es -- λ> quickCheck prop_especiales -- +++ OK, passed 100 tests. -- Comparación de eficiencia -- ========================= -- La comparación es -- Comparación -- λ> especiales1 !! 2000 -- 2390 -- (3.38 secs, 4,724,497,192 bytes) -- λ> especiales2 !! 2000 -- 2390 -- (1.91 secs, 4,303,415,512 bytes) -- λ> especiales3 !! 2000 -- 2390 -- (0.01 secs, 4,209,664 bytes) |
El código se encuentra en GitHub.
5. Método de bisección para aproximar raíces de funciones
El método de bisección para calcular un cero de una función en el intervalo [a,b] se basa en el teorema de Bolzano:
“Si f(x) es una función continua en el intervalo [a, b], y si, además, en los extremos del intervalo la función f(x) toma valores de signo opuesto (f(a) * f(b) < 0), entonces existe al menos un valor c en (a, b) para el que f(c) = 0".
El método para calcular un cero de la función f en el intervalo [a,b] con un error menor que e consiste en tomar el punto medio del intervalo c = (a+b)/2 y considerar los siguientes casos:
- Si |f(c)| < e, hemos encontrado una aproximación del punto que anula f en el intervalo con un error aceptable.
- Si f(c) tiene signo distinto de f(a), repetir el proceso en el intervalo [a,c].
- Si no, repetir el proceso en el intervalo [c,b].
Definir la función
|
1 |
biseccion :: (Double -> Double) -> Double -> Double -> Double -> Double |
tal que (biseccion f a b e) es una aproximación del punto del intervalo [a,b] en el que se anula la función f, con un error menor que e, calculada mediante el método de la bisección. Por ejemplo,
|
1 2 3 4 |
biseccion (\x -> x^2 - 3) 0 5 0.01 == 1.7333984375 biseccion (\x -> x^3 - x - 2) 0 4 0.01 == 1.521484375 biseccion cos 0 2 0.01 == 1.5625 biseccion (\x -> log (50-x) - 4) (-10) 3 0.01 == -5.125 |
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 |
import Test.QuickCheck (Arbitrary, Gen, Property, (==>), arbitrary, choose, sized, quickCheck) -- 1ª solución -- =========== biseccion1 :: (Double -> Double) -> Double -> Double -> Double -> Double biseccion1 f a b e | abs (f c) < e = c | f a * f c < 0 = biseccion1 f a c e | otherwise = biseccion1 f c b e where c = (a+b)/2 -- 2ª solución -- =========== biseccion2 :: (Double -> Double) -> Double -> Double -> Double -> Double biseccion2 f a b e = aux a b where aux a' b' | abs (f c) < e = c | f a' * f c < 0 = aux a' c | otherwise = aux c b' where c = (a'+b')/2 -- Comprobación de equivalencia -- ============================ newtype Polinomio = P [Int] deriving Show valorPolinomio :: Polinomio -> Double -> Double valorPolinomio (P cs) x = sum [fromIntegral c * x^n | (c,n) <- zip cs [0..]] polinomioArbitrario :: Int -> Gen Polinomio polinomioArbitrario 0 = return (P []) polinomioArbitrario n = do c <- choose (-10,10) (P xs) <- polinomioArbitrario (n `div` 2) return (P (c:xs)) instance Arbitrary Polinomio where arbitrary = sized polinomioArbitrario -- La propiedad es prop_biseccion :: Polinomio -> Double -> Double -> Double -> Property prop_biseccion p a b e = f a * f b < 0 && e > 0 ==> biseccion1 f a b e =~ biseccion2 f a b e where f = valorPolinomio p x =~ y = abs (x - y) < 0.001 -- La comprobación es -- λ> quickCheckWith (stdArgs {maxDiscardRatio=30}) prop_biseccion -- +++ OK, passed 100 tests; 2156 discarded. |
El código se encuentra en GitHub.