Menu Close

Etiqueta: sqrt

El teorema de Navidad de Fermat

El 25 de diciembre de 1640, en una carta a Mersenne, Fermat demostró la conjetura de Girard: todo primo de la forma 4n+1 puede expresarse de manera única como suma de dos cuadrados. Por eso es conocido como el Teorema de Navidad de Fermat.

Definir las funciones

   representaciones :: Integer -> [(Integer,Integer)]
   primosImparesConRepresentacionUnica :: [Integer]
   primos4nM1 :: [Integer]

tales que

  • (representaciones n) es la lista de pares de números naturales (x,y) tales que n = x^2 + y^2 con x <= y. Por ejemplo,
     representaciones  20           ==  [(2,4)]
     representaciones  25           ==  [(0,5),(3,4)]
     representaciones 325           ==  [(1,18),(6,17),(10,15)]
     representaciones 100000147984  ==  [(0,316228)]
     length (representaciones (10^10))    ==  6
     length (representaciones (4*10^12))  ==  7
  • primosImparesConRepresentacionUnica es la lista de los números primos impares que se pueden escribir exactamente de una manera como suma de cuadrados de pares de números naturales (x,y) con x <= y. Por ejemplo,
     λ> take 20 primosImparesConRepresentacionUnica
     [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]
  • primos4nM1 es la lista de los números primos que se pueden escribir como uno más un múltiplo de 4 (es decir, que son congruentes con 1 módulo 4). Por ejemplo,
     λ> take 20 primos4nM1
     [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]

El teorema de Navidad de Fermat afirma que un número primo impar p se puede escribir exactamente de una manera como suma de dos cuadrados de números naturales p = x² + y^2 (con x <= y) si, y sólo si, p se puede escribir como uno más un múltiplo de 4 (es decir, que es congruente con 1 módulo 4).

Comprobar con QuickCheck el teorema de Navidad de Fermat; es decir, que para todo número n, los n-ésimos elementos de primosImparesConRepresentacionUnica y de primos4nM1 son iguales.

Soluciones

import Data.Numbers.Primes (primes)
import Test.QuickCheck
 
-- 1ª definición de representaciones
-- =================================
 
representaciones :: Integer -> [(Integer,Integer)]
representaciones n =
  [(x,y) | x <- [0..n], y <- [x..n], n == x*x + y*y]
 
-- 2ª definición de representaciones
-- =================================
 
representaciones2 :: Integer -> [(Integer,Integer)]
representaciones2 n =
  [(x,raiz z) | x <- [0..raiz (n `div` 2)] 
              , let z = n - x*x
              , esCuadrado z]
 
-- (esCuadrado x) se verifica si x es un número al cuadrado. Por
-- ejemplo,
--    esCuadrado 25  ==  True
--    esCuadrado 26  ==  False
esCuadrado :: Integer -> Bool
esCuadrado x = x == y * y
  where y = raiz x
 
-- (raiz x) es la raíz cuadrada entera de x. Por ejemplo,
--    raiz 25  ==  5
--    raiz 24  ==  4
--    raiz 26  ==  5
raiz :: Integer -> Integer 
raiz 0 = 0
raiz 1 = 1
raiz x = aux (0,x)
    where aux (a,b) | d == x    = c
                    | c == a    = a
                    | d < x     = aux (c,b)
                    | otherwise = aux (a,c) 
              where c = (a+b) `div` 2
                    d = c^2
 
-- 3ª definición de representaciones
-- =================================
 
representaciones3 :: Integer -> [(Integer,Integer)]
representaciones3 n =
  [(x,raiz3 z) | x <- [0..raiz3 (n `div` 2)] 
               , let z = n - x*x
               , esCuadrado3 z]
 
-- (esCuadrado x) se verifica si x es un número al cuadrado. Por
-- ejemplo,
--    esCuadrado3 25  ==  True
--    esCuadrado3 26  ==  False
esCuadrado3 :: Integer -> Bool
esCuadrado3 x = x == y * y
  where y = raiz3 x
 
-- (raiz3 x) es la raíz cuadrada entera de x. Por ejemplo,
--    raiz3 25  ==  5
--    raiz3 24  ==  4
--    raiz3 26  ==  5
raiz3 :: Integer -> Integer
raiz3 x = floor (sqrt (fromIntegral x))
 
-- 4ª definición de representaciones
-- =================================
 
representaciones4 :: Integer -> [(Integer, Integer)]
representaciones4 n = aux 0 (floor (sqrt (fromIntegral n)))
  where aux x y
          | x > y     = [] 
          | otherwise = case compare (x*x + y*y) n of
                          LT -> aux (x + 1) y
                          EQ -> (x, y) : aux (x + 1) (y - 1)
                          GT -> aux x (y - 1)
 
-- Equivalencia de las definiciones de representaciones
-- ====================================================
 
-- La propiedad es
prop_representaciones_equiv :: (Positive Integer) -> Bool
prop_representaciones_equiv (Positive n) =
  representaciones  n == representaciones2 n &&
  representaciones2 n == representaciones3 n &&
  representaciones3 n == representaciones4 n
 
-- La comprobación es
--    λ> quickCheck prop_representaciones_equiv
--    +++ OK, passed 100 tests.
 
-- Comparación de eficiencia de las definiciones de representaciones
-- =================================================================
 
--    λ> representaciones 3025
--    [(0,55),(33,44)]
--    (2.86 secs, 1,393,133,528 bytes)
--    λ> representaciones2 3025
--    [(0,55),(33,44)]
--    (0.00 secs, 867,944 bytes)
--    λ> representaciones3 3025
--    [(0,55),(33,44)]
--    (0.00 secs, 173,512 bytes)
--    λ> representaciones4 3025
--    [(0,55),(33,44)]
--    (0.00 secs, 423,424 bytes)
--    
--    λ> length (representaciones2 (10^10))
--    6
--    (3.38 secs, 2,188,903,544 bytes)
--    λ> length (representaciones3 (10^10))
--    6
--    (0.10 secs, 62,349,048 bytes)
--    λ> length (representaciones4 (10^10))
--    6
--    (0.11 secs, 48,052,360 bytes)
--
--    λ> length (representaciones3 (4*10^12))
--    7
--    (1.85 secs, 1,222,007,176 bytes)
--    λ> length (representaciones4 (4*10^12))
--    7
--    (1.79 secs, 953,497,480 bytes)
 
-- Definición de primosImparesConRepresentacionUnica
-- =================================================
 
primosImparesConRepresentacionUnica :: [Integer]
primosImparesConRepresentacionUnica =
  [x | x <- tail primes
     , length (representaciones4 x) == 1]
 
-- Definición de primos4nM1
-- ========================
 
primos4nM1 :: [Integer]
primos4nM1 = [x | x <- primes
                , x `mod` 4 == 1]
 
-- Teorema de Navidad de Fermat
-- ============================
 
-- La propiedad es
prop_teoremaDeNavidadDeFermat :: Positive Int -> Bool
prop_teoremaDeNavidadDeFermat (Positive n) =
  primosImparesConRepresentacionUnica !! n == primos4nM1 !! n
 
-- La comprobación es
--    λ> quickCheck prop_teoremaDeNavidadDeFermat
--    +++ OK, passed 100 tests.

Pensamiento

Dijo Dios: brote la nada
Y alzó su mano derecha,
hasta ocultar su mirada.
Y quedó la nada hecha.

Antonio Machado

Cálculo de pi usando la fórmula de Vieta

La fórmula de Vieta para el cálculo de pi es la siguiente
Calculo_de_pi_usando_la_formula_de_Vieta

Definir las funciones

   aproximacionPi :: Int -> Double
   errorPi :: Double -> Int

tales que

  • (aproximacionPi n) es la aproximación de pi usando n factores de la fórmula de Vieta. Por ejemplo,
     aproximacionPi  5  ==  3.140331156954753
     aproximacionPi 10  ==  3.1415914215112
     aproximacionPi 15  ==  3.141592652386592
     aproximacionPi 20  ==  3.1415926535886207
     aproximacionPi 25  ==  3.141592653589795
  • (errorPi x) es el menor número de factores de la fórmula de Vieta necesarios para obtener pi con un error menor que x. Por ejemplo,
     errorPi 0.1        ==  2
     errorPi 0.01       ==  4
     errorPi 0.001      ==  6
     errorPi 0.0001     ==  7
     errorPi 1e-4       ==  7
     errorPi 1e-14      ==  24
     pi                 ==  3.141592653589793
     aproximacionPi 24  ==  3.1415926535897913

Soluciones

-- 1ª definición de aproximacionPi
aproximacionPi :: Int -> Double
aproximacionPi n = product [2 / aux x | x <- [0..n]]
  where
    aux 0 = 1
    aux 1 = sqrt 2
    aux n = sqrt (2 + aux (n-1))
 
-- 2ª definición de aproximacionPi
aproximacionPi2 :: Int -> Double
aproximacionPi2 n = product [2/x | x <- 1 : xs] 
  where xs = take n $ iterate (\x -> sqrt (2+x)) (sqrt 2)
 
-- 3ª definición de aproximaxionPi
aproximacionPi3 :: Int -> Double
aproximacionPi3 n =  product (2 : take n (map (2/) xs))
  where xs = sqrt 2 : [sqrt (2 + x) | x <- xs]
 
-- 1ª definición de errorPi
errorPi :: Double -> Int
errorPi x = head [n | n <- [1..]
                    , abs (pi - aproximacionPi n) < x]
 
-- 2ª definición de errorPi
errorPi2 :: Double -> Int
errorPi2 x = until aceptable (+1) 1
  where aceptable n = abs (pi - aproximacionPi n) < x

Pensamiento

El tiempo que la barba me platea,
cavó mis ojos y agrandó mi frente,
va siendo en mi recuerdo transparente,
y mientras más al fondo, más clarea.

Antonio Machado

Suma de primos menores

La suma de los primos menores que 10 es 2 + 3 + 5 + 7 = 17.

Definir la función

   sumaPrimosMenores :: Integer -> Integer

tal que (sumaPrimosMenores n) es la suma de los primos menores que n. Por ejemplo,

   sumaPrimosMenores 10        ==  17
   sumaPrimosMenores (5*10^5)  ==  9914236195

Nota: Este ejercicio está basado en el problema 10 del Proyecto Euler

Soluciones

import Data.Numbers.Primes (primes)
 
-- 1ª solución
-- ===========
 
sumaPrimosMenores :: Integer -> Integer
sumaPrimosMenores n = sum (takeWhile (<n) primos)
 
-- primos es la lista de los números primos. Por ejemplo,
--    λ> take 12 primos
--    [2,3,5,7,11,13,17,19,23,29,31,37]
primos :: [Integer]
primos = 2 : filter esPrimo [3,5..]
 
esPrimo :: Integer -> Bool
esPrimo n = null [x | x <- [2..(ceiling . sqrt . fromIntegral) n]
                    , n `mod` x == 0]
 
-- 2ª solución
-- ===========
 
sumaPrimosMenores2 :: Integer -> Integer
sumaPrimosMenores2 n = sum (takeWhile (<n) primos2)
 
primos2 :: [Integer]
primos2 = 2 : filter esPrimo2 [3,5..]
 
esPrimo2 :: Integer -> Bool
esPrimo2 x =
  all ((/= 0) . mod x)
  (takeWhile (<= floor (sqrt (fromIntegral x))) primos2)
 
-- 3ª solución
-- ===========
 
sumaPrimosMenores3 :: Integer -> Integer
sumaPrimosMenores3 n = sum (takeWhile (<n) primes)
 
-- Comparación de eficiencia
-- =========================
 
--    λ> sumaPrimosMenores (2*10^5)
--    1709600813
--    (2.56 secs, 1,522,015,240 bytes)
--    λ> sumaPrimosMenores2 (2*10^5)
--    1709600813
--    (0.56 secs, 376,951,456 bytes)
--    λ> sumaPrimosMenores3 (2*10^5)
--    1709600813
--    (0.07 secs, 62,321,888 bytes)

Pensamiento

El movimiento no es nada esencial. La fuerza puede ser inmóvil (lo es en su estado de pureza); mas no por ello deja de ser activa.

Antonio Machado

Mayor divisor primo

Los divisores primos de 13195 son 5, 7, 13 y 29. Por tanto, el mayor divisor primo de 13195 es 29.

Definir la función

   mayorDivisorPrimo :: Integer -> Integer

tal que (mayorDivisorPrimo n) es el mayor divisor primo de n. Por ejemplo,

   mayorDivisorPrimo 13195            ==  29
   mayorDivisorPrimo 152416333181401  ==  12345701

Nota: Este ejercicio está basado en el problema 3 del Proyecto Euler

Soluciones

import Data.Numbers.Primes (primeFactors)
 
-- 1ª solución  (sin librerías auxiliares)
-- =======================================
 
mayorDivisorPrimo :: Integer -> Integer
mayorDivisorPrimo = last . divisoresPrimos 
 
-- (divisoresPrimos n) es la lista de los divisores primos de n. Por
-- ejemplo, 
--    divisoresPrimos 13195  ==  [5,7,13,29]
divisoresPrimos :: Integer -> [Integer]
divisoresPrimos 0 = []
divisoresPrimos 1 = []
divisoresPrimos n = m : divisoresPrimos (n `div` m)
  where m = menorDivisorPrimo n 
 
-- (menorDivisorPrimo n) es el menor divisor primo de n. Por ejemplo, 
--    menorDivisorPrimo 24  ==  2
--    menorDivisorPrimo 25  ==  5
--    menorDivisorPrimo 29  ==  29
menorDivisorPrimo :: Integer -> Integer
menorDivisorPrimo x =
  head [y | y <- 2 : [3,5..(ceiling . sqrt . fromIntegral) x] ++ [x]
          , x `mod` y == 0]
 
-- 2ª solución (con la librería Data.Numbers.Primes)
-- =================================================
 
mayorDivisorPrimo2 :: Integer -> Integer
mayorDivisorPrimo2 = last . primeFactors
 
-- Comparación de eficiencia
-- =========================
 
--   λ> mayorDivisorPrimo 152416333181401
--   12345701
--   (1.96 secs, 1,630,201,856 bytes)
--   λ> mayorDivisorPrimo2 152416333181401
--   12345701
--   (2.01 secs, 5,445,284,432 bytes)

Pensamiento

“Un programa de ordenador es una demostración.” ~ Igor Rivin

Caminos minimales en un árbol numérico

En la librería Data.Tree se definen los tipos de árboles y bosques como sigue

   data Tree a   = Node a (Forest a)
   type Forest a = [Tree a]

Se pueden definir árboles. Por ejemplo,

   ej = Node 3 [Node 5 [Node 9 []], Node 7 []]

Y se pueden dibujar con la función drawTree. Por ejemplo,

   λ> putStrLn (drawTree (fmap show ej))
   3
   |
   +- 5
   |  |
   |  `- 9
   |
   `- 7

Los mayores divisores de un número x son los divisores u tales que u > 1 y existe un v tal que 1 < v < u y u.v = x. Por ejemplo, los mayores divisores de 24 son 12, 8 y 6.

El árbol de los predecesores y mayores divisores de un número x es el árbol cuya raíz es x y los sucesores de cada nodo y > 1 es el conjunto formado por y-1 junto con los mayores divisores de y. Los nodos con valor 1 no tienen sucesores. Por ejemplo, el árbol de los predecesores y mayores divisores del número 6 es

       6
      / \
     5   3 
     |   |
     4   2
    / \  |
   3   2 1 
   |   | 
   2   1
   |
   1

Definir las siguientes funciones

   mayoresDivisores :: Int -> [Int]
   arbol            :: Int -> Tree Int
   caminos          :: Int -> [[Int]]
   caminosMinimales :: Int -> [[Int]]

tales que
+ (mayoresDivisores x) es la lista de los mayores divisores de x. Por ejemplo,

     mayoresDivisores 24  ==  [12,8,6]
     mayoresDivisores 16  ==  [8,4]
     mayoresDivisores 10  ==  [5]
     mayoresDivisores 17  ==  []
  • (arbol x) es el árbol de los predecesores y mayores divisores del número x. Por ejemplo,
     λ> putStrLn (drawTree (fmap show (arbol 6)))
     6
     |
     +- 5
     |  |
     |  `- 4
     |     |
     |     +- 3
     |     |  |
     |     |  `- 2
     |     |     |
     |     |     `- 1
     |     |
     |     `- 2
     |        |
     |        `- 1
     |
     `- 3
        |
        `- 2
           |
           `- 1
  • (caminos x) es la lista de los caminos en el árbol de los predecesores y mayores divisores del número x. Por ejemplo,
     λ> caminos 6
     [[6,5,4,3,2,1],[6,5,4,2,1],[6,3,2,1]]
  • (caminosMinimales x) es la lista de los caminos en de menor longitud en el árbol de los predecesores y mayores divisores del número x. Por ejemplo,
     λ> caminosMinimales 6
     [[6,3,2,1]]
     λ> caminosMinimales 17
     [[17,16,4,2,1]]
     λ> caminosMinimales 50
     [[50,25,5,4,2,1],[50,10,9,3,2,1],[50,10,5,4,2,1]]

Soluciones

import Data.Tree
import Test.QuickCheck
 
mayoresDivisores :: Int -> [Int]
mayoresDivisores x =
  [max u v | u <- [2..floor (sqrt (fromIntegral x))]
           , x `mod` u == 0
           , let v = x `div` u]  
 
arbol :: Int -> Tree Int
arbol 1 = Node 1 []
arbol x = Node x (arbol (x-1) : [arbol y | y <- mayoresDivisores x])
 
caminos :: Int -> [[Int]]
caminos = caminosArbol . arbol
 
--    λ> caminosArbol (arbol 6)
--    [[6,5,4,3,2,1],[6,5,4,2,1],[6,3,2,1]]
caminosArbol :: Tree a -> [[a]]
caminosArbol (Node x []) = [[x]]
caminosArbol (Node x as) = [x:ys | ys <- caminosBosque as]
 
caminosBosque :: Forest a -> [[a]]
caminosBosque = concatMap caminosArbol
 
caminosMinimales :: Int -> [[Int]]
caminosMinimales x = [ys | ys <- yss, length ys == m]
  where yss = caminos x
        m   = minimum (map length yss)

Pensamiento

Tras el vivir y el soñar,
está lo que más importa:
despertar.

Antonio Machado