Menu Close

PFH: La semana en Exercitium (25 de junio de 2022)

Esta semana he publicado en Exercitium las soluciones de los siguientes problemas:

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

   menor :: Integer -> Integer

tal que (menor n) es el menor número con 2^n divisores. Por ejemplo,

   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

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

   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

   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

   integral 0 1 (^3) 0.01  ==  0.24998750000000042

Otros ejemplos son

   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

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

   suma :: Integer -> Integer

tal que (suma n) es la suma 1·1! + 2·2! + 3·3! + ... + n·n!. Por ejemplo,

   suma 1  ==  1
   suma 2  ==  5
   suma 3  ==  23
   suma 4  ==  119
   suma 5  ==  719
   take 9 (show (suma 70000))  ==  "823780458"

Soluciones

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

   mcm(1,2,3,4,5) = 60 = mcm(1,2,3,4,5,6)

Definir la sucesión

   especiales :: [Integer]

cuyos términos son los números especiales. Por ejemplo,

   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

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

   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,

   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

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.

PFH