Menu Close

Etiqueta: sqrt

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

Primo suma de dos cuadrados

Definir la sucesión

   primosSumaDe2Cuadrados :: [Integer]

cuyos elementos son los números primos que se pueden escribir como sumas de dos cuadrados. Por ejemplo,

   λ> take 20 primosSumaDe2Cuadrados
   [2,5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181]
   λ> primosSumaDe2Cuadrados !! (2*10^5)
   5803241

En el ejemplo anterior,

  • 13 está en la sucesión porque es primo y 13 = 2²+3².
  • 11 no está en la sucesión porque no se puede escribir como suma de dos cuadrados (en efecto, 11-1=10, 11-2²=7 y 11-3²=2 no son cuadrados).
  • 20 no está en la sucesión porque, aunque es suma de dos cuadrados (20=4²+2²), no es primo.

Soluciones

import Data.Numbers.Primes (primes)
 
-- 1ª solución
-- ===========
 
primosSumaDe2Cuadrados1 :: [Integer]
primosSumaDe2Cuadrados1 =
    [n | n <- primes,
         esSumaDe2Cuadrados n]
 
esSumaDe2Cuadrados :: Integer -> Bool
esSumaDe2Cuadrados n = not (null (sumas2cuadrados n))
 
sumas2cuadrados :: Integer -> [(Integer,Integer)]
sumas2cuadrados n = 
    [(x,y) | x <- [1..ceiling (sqrt (fromIntegral n))],
             let z = n - x^2,
             esCuadrado z, 
             let y = raiz z,
             x <= y]
 
-- (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 x = floor (sqrt (fromIntegral x))
 
-- 2ª solución
-- ===========
 
primosSumaDe2Cuadrados2 :: [Integer]
primosSumaDe2Cuadrados2 =
    2 : [n | n <- primes,
             n `mod` 4 == 1]
 
-- Comparación de eficiencia
-- =========================
 
--    λ> primosSumaDe2Cuadrados1 !! (2*10^3)
--    38153
--    (3.03 secs, 568,239,744 bytes)
--    λ> primosSumaDe2Cuadrados2 !! (2*10^3)
--    38153
--    (0.05 secs, 20,017,912 bytes)

Referencias

Sumas de potencias de 3 primos

Los primeros números de la forma p²+q³+r⁴, con p, q y r primos son

   28 = 2² + 2³ + 2⁴
   33 = 3² + 2³ + 2⁴
   47 = 2² + 3³ + 2⁴
   49 = 5² + 2³ + 2⁴

Definir la sucesión

   sumas3potencias :: [Integer]

cuyos elementos son los números que se pueden escribir de la forma p²+q³+r⁴, con p, q y r primos. Por ejemplo,

   λ> take 15 sumas3potencias
   [28,33,47,49,52,68,73,92,93,98,112,114,117,133,138]
   λ> sumas3potencias !! 234567
   8953761

Soluciones

import Data.Numbers.Primes (primes)
 
-- 1ª solución
-- ===========
 
sumas3potencias1 :: [Integer]
sumas3potencias1 = filter esTripleSuma [1..]
 
-- (esTripleSuma n) se verifica si n se puede escribir de la forma
-- p²+q³+r⁴, con p, q y r primos. Por ejemplo, 
--    esTripleSuma 33  ==  True
--    esTripleSuma 45  ==  False
esTripleSuma :: Integer -> Bool
esTripleSuma = not . null . triplesSumas  
 
-- (triplesSumas n) es la lista de las ternas de primos (p,q,r) tales
-- que p²+q³+r⁴ = n. Por ejemplo,
--    triplesSumas 33   ==  [(3,2,2)]
--    triplesSumas 145  ==  [(11,2,2),(2,5,2)]
--    triplesSumas 45   ==  []
triplesSumas :: Integer -> [(Integer,Integer,Integer)]
triplesSumas n =  
    [(p,q,r) 
    | r <- xs, 
      q <- xs, 
      let p2 = n - q^3 - r^4,
      let p = (floor . sqrt . fromIntegral) p2,
      p*p == p2,
      isPrime p]
    where xs = takeWhile (< (ceiling $ sqrt $ fromIntegral n)) primes
 
-- 2ª solución
-- ===========
 
sumas3potencias2 :: [Integer]
sumas3potencias2 = sumaOrdenadaInf as (sumaOrdenadaInf bs cs)
 
sumaOrdenadaInf:: [Integer] -> [Integer] -> [Integer]
sumaOrdenadaInf xs ys = mezclaTodas [map (+x) ys | x <- xs]
 
mezclaTodas :: Ord a => [[a]] -> [a]
mezclaTodas = foldr1 xmezcla
    where xmezcla (x:xs) ys = x : mezcla xs ys
 
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
 
as = map (^2) primes
bs = map (^3) primes
cs = map (^4) primes

Número de representaciones de n como suma de dos cuadrados

Sea n un número natural cuya factorización prima es

   n = 2^a * p(1)^b(1) *...* p(n)^b(n) * q(1)^c(1) *...* q(m)^c(m)

donde los p(i) son los divisores primos de n congruentes con 3 módulo 4 y los q(j) son los divisores primos de n congruentes con 1 módulo 4 Entonces, el número de forma de descomponer n como suma de dos cuadrados es 0, si algún b(i) es impar y es el techo (es decir, el número entero más próximo por exceso) de

   ((1+c(1)) *...* (1+c(m))) / 2

en caso contrario. Por ejemplo, el número 2^3*(3^9*7^8)*(5^3*13^6) no se puede descomponer como sumas de dos cuadrados (porque el exponente de 3 es impar) y el número 2^3*(3^2*7^8)*(5^3*13^6) tiene 14 descomposiciones como suma de dos cuadrados (porque los exponentes de 3 y 7 son pares y el techo de ((1+3)*(1+6))/2 es 14).

Definir la función

   nRepresentaciones :: Integer -> Integer

tal que (nRepresentaciones n) es el número de formas de representar n como suma de dos cuadrados. Por ejemplo,

   nRepresentaciones (2^3*3^9*5^3*7^8*13^6)        ==  0
   nRepresentaciones (2^3*3^2*5^3*7^8*13^6)        ==  14
   head [n | n <- [1..], nRepresentaciones n > 8]  ==  71825

Usando la función representaciones del ejercicio anterior, comprobar con QuickCheck la siguiente propiedad

   prop_nRepresentaciones :: Integer -> Property
   prop_nRepresentaciones n =
       n > 0 ==>
         nRepresentaciones2 n == genericLength (representaciones n)

Soluciones

import Test.QuickCheck
import Data.List (genericLength, group)
import Data.Numbers.Primes (primeFactors)
 
-- 1ª definición
-- =============
 
nRepresentaciones1 :: Integer -> Integer
nRepresentaciones1 n =
    ceiling (fromIntegral (product (map aux (factorizacion n))) / 2)
    where aux (p,e) | p == 2         = 1
                    | p `mod` 4 == 3 = if even e then 1 else 0
                    | otherwise      = e+1
 
-- (factorizacion n) es la factorización prima de n. Por ejemplo,
--    factorizacion 600  ==  [(2,3),(3,1),(5,2)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
    map (\xs -> (head xs, genericLength xs)) (group (primeFactors n))
 
-- 2ª definición
-- =============
 
nRepresentaciones2 :: Integer -> Integer
nRepresentaciones2 n =
    (1 + product (map aux (factorizacion n))) `div`  2
    where aux (p,e) | p == 2         = 1
                    | p `mod` 4 == 3 = if even e then 1 else 0
                    | otherwise      = e+1
 
 
-- Comprobación de la propiedad
-- ============================
 
-- La propiedad es
prop_nRepresentaciones :: Integer -> Property
prop_nRepresentaciones n =
    n > 0 ==>
      nRepresentaciones2 n == genericLength (representaciones n)
 
representaciones :: Integer -> [(Integer,Integer)]
representaciones n =
    [(x,raiz z) | x <- [0..raiz (n `div` 2)], 
                  let z = n - x*x,
                  esCuadrado z]
 
esCuadrado :: Integer -> Bool
esCuadrado x = x == y * y
    where y = raiz x
 
raiz :: Integer -> Integer
raiz x = floor (sqrt (fromIntegral x))
 
-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=7}) prop_representacion
--    +++ OK, passed 100 tests.

Referencias

La sucesión [nRepresentaciones n | n <- [0..]] es la A000161 de la OEIS.

Representaciones de un número como suma de dos cuadrados

Definir la función

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

tal que (representaciones n) es la lista de pares de números naturales (x,y) tales que n = x² + y² 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)]

Comprobar con QuickCheck que un número natural n se puede representar como suma de dos cuadrados si, y sólo si, en la factorización prima de n todos los exponentes de sus factores primos congruentes con 3 módulo 4 son pares.

Soluciones

import Test.QuickCheck
import Data.List (genericLength, group)
import Data.Numbers.Primes (primeFactors)
 
-- 1ª solución
-- ===========
 
representaciones1 :: Integer -> [(Integer,Integer)]
representaciones1 n =
    [(x,y) | x <- [0..n], y <- [x..n], n == x*x + y*y]
 
-- 2ª solución
-- ===========
 
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 x = floor (sqrt (fromIntegral x))
 
-- 3ª solución
-- ===========
 
representaciones3 :: Integer -> [(Integer, Integer)]
representaciones3 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)
 
-- Comparación de eficiencia
-- =========================
 
--    λ> representaciones1 2000
--    [(8,44),(20,40)]
--    (4.08 secs, 613,941,832 bytes)
--    λ> representaciones2 2000
--    [(8,44),(20,40)]
--    (0.01 secs, 0 bytes)
--    λ> representaciones3 2000
--    [(8,44),(20,40)]
--    (0.01 secs, 0 bytes)
--    
--    λ> length (representaciones2 1000000000000)
--    7
--    (3.95 secs, 760,612,080 bytes)
--    λ> length (representaciones3 1000000000000)
--    7
--    (3.41 secs, 470,403,832 bytes)
 
-- Comprobación de la propiedad
-- ============================
 
-- La propiedad es
prop_representacion :: Integer -> Property
prop_representacion n =
    n > 0 ==>
      not (null (representaciones2 n)) == 
          all (\(p,e) -> p `mod` 4 /= 3 || even e) (factorizacion n)
 
-- (factorizacion n) es la factorización prima de n. Por ejemplo,
--    factorizacion 600  ==  [(2,3),(3,1),(5,2)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
    map (\xs -> (head xs, genericLength xs)) (group (primeFactors n))
 
-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=7}) prop_representacion
--    +++ OK, passed 100 tests.

Factorizaciones de números de Hilbert

Un número de Hilbert es un entero positivo de la forma 4n+1. Los primeros números de Hilbert son 1, 5, 9, 13, 17, 21, 25, 29, 33, 37, 41, 45, 49, 53, 57, 61, 65, 69, …

Un primo de Hilbert es un número de Hilbert n que no es divisible por ningún número de Hilbert menor que n (salvo el 1). Los primeros primos de Hilbert son 5, 9, 13, 17, 21, 29, 33, 37, 41, 49, 53, 57, 61, 69, 73, 77, 89, 93, 97, 101, 109, 113, 121, 129, 133, 137, …

Definir la función

   factorizacionesH :: Integer -> [[Integer]]

tal que (factorizacionesH n) es la listas de primos de Hilbert cuyo producto es el número de Hilbert n. Por ejemplo,

  factorizacionesH  25  ==  [[5,5]]
  factorizacionesH  45  ==  [[5,9]]
  factorizacionesH 441  ==  [[9,49],[21,21]]

Comprobar con QuickCheck que todos los números de Hilbert son factorizables como producto de primos de Hilbert (aunque laa factorización, como para el 441, puede no ser única).

Soluciones

import Test.QuickCheck
 
-- numerosH es la sucesión de los números de Hilbert. Por ejemplo,
--    take 15 numerosH  ==  [1,5,9,13,17,21,25,29,33,37,41,45,49,53,57]
numerosH :: [Integer]
numerosH = [1,5..]
 
-- primosH es la sucesión de primos de Hilbert. Por ejemplo,
--    take 15 primosH  ==  [5,9,13,17,21,29,33,37,41,49,53,57,61,69,73]
primosH :: [Integer]
primosH = filter esPrimoH (tail numerosH) 
    where esPrimoH n = all noDivideAn [5,9..m]
              where noDivideAn x = n `mod` x /= 0
                    m            = ceiling (sqrt (fromIntegral n))
 
factorizacionesH :: Integer -> [[Integer]]
factorizacionesH n = aux n primosH
    where aux n (x:xs)
              | x == n         = [[n]]
              | x > n          = []
              | n `mod` x == 0 = map (x:) (aux (n `div` x) (x:xs))
                                 ++ aux n xs
              | otherwise      = aux n xs
 
-- La propiedad es
prop_factorizable :: Integer -> Property
prop_factorizable n =
    n > 0 ==> not (null (factorizacionesH (1 + 4 * n)))
 
-- La comprobación es
--    λ> quickCheck prop_factorizable
--    +++ OK, passed 100 tests.

Referencias

Basado en el artículo Failure of unique factorization (A simple example of the failure of the fundamental theorem of arithmetic) de R.J. Lipton en el blog Gödel’s Lost Letter and P=NP.

Otras referencias

Números primos de Hilbert

Un número de Hilbert es un entero positivo de la forma 4n+1. Los primeros números de Hilbert son 1, 5, 9, 13, 17, 21, 25, 29, 33, 37, 41, 45, 49, 53, 57, 61, 65, 69, …

Un primo de Hilbert es un número de Hilbert n que no es divisible por ningún número de Hilbert menor que n (salvo el 1). Los primeros primos de Hilbert son 5, 9, 13, 17, 21, 29, 33, 37, 41, 49, 53, 57, 61, 69, 73, 77, 89, 93, 97, 101, 109, 113, 121, 129, 133, 137, …

Definir la sucesión

   primosH :: [Integer]

tal que sus elementos son los primos de Hilbert. Por ejemplo,

   take 15 primosH  == [5,9,13,17,21,29,33,37,41,49,53,57,61,69,73]
   primosH !! 20000 == 203221

Soluciones

-- 1ª definición
-- =============
 
-- numerosH es la sucesión de los números de Hilbert. Por ejemplo,
--    take 15 numerosH  ==  [1,5,9,13,17,21,25,29,33,37,41,45,49,53,57]
numerosH :: [Integer]
numerosH = [1,5..]
 
-- (divisoresH n) es la lista de los números de Hilbert que dividen a
-- n. Por ejemplo,
--   divisoresH 117  ==  [1,9,13,117]
--   divisoresH  21  ==  [1,21]
divisoresH :: Integer -> [Integer]
divisoresH n = [x | x <- takeWhile (<=n) numerosH,
                    n `mod` x == 0]
 
primosH1 :: [Integer]
primosH1 = [n | n <- tail numerosH,
                divisoresH n == [1,n]]
 
-- 2ª definición
-- =============
 
primosH2 :: [Integer]
primosH2 = filter esPrimoH (tail numerosH) 
    where esPrimoH n = all noDivideAn [5,9..m]
              where noDivideAn x = n `mod` x /= 0
                    m            = ceiling (sqrt (fromIntegral n))
 
-- Comparación de eficiencia
-- =========================
 
--    λ> primosH1 !! 2000
--    16957
--    (6.93 secs, 750,291,352 bytes)
--    λ> primosH2 !! 2000
--    16957
--    (0.13 secs, 18,066,288 bytes)

Raíces enteras de los números primos

Definir la sucesión

   raicesEnterasDePrimos :: [Integer]

cuyos elementos son las partes enteras de las raíces cuadradas de los números primos. Por ejemplo,

   λ> take 30 raicesEnterasDePrimos
   [1,1,2,2,3,3,4,4,4,5,5,6,6,6,6,7,7,7,8,8,8,8,9,9,9,10,10,10,10,10]
   λ> raicesEnterasDePrimos !!  9963
   322
   λ> raicesEnterasDePrimos !!  9964
   323

Comprobar con QuickCheck que la diferencia entre dos términos consecutivos de la sucesión es como máximo igual a 1.

Soluciones

import Data.Numbers.Primes (primes)
import Test.QuickCheck
 
-- 1ª solución
-- ===========
 
raicesEnterasDePrimos1 :: [Integer]
raicesEnterasDePrimos1 = map raizEntera primes
 
-- (raizEntera x) es la parte entera de la raíz cuadrada de x. Por
-- ejemplo,
--    raizEntera  8  ==  2
--    raizEntera  9  ==  3
--    raizEntera 10  ==  3
raizEntera :: Integer -> Integer
raizEntera = floor . sqrt . fromIntegral 
 
-- 2ª solución
-- ===========
 
raicesEnterasDePrimos2 :: [Integer]
raicesEnterasDePrimos2 = map raizEntera2 primes
 
raizEntera2 :: Integer -> Integer
raizEntera2 n = aux 1
    where aux k | k*k > n   = k-1
                | otherwise = aux (k+1)
 
-- 3º solución
-- ===========
 
raicesEnterasDePrimos3 :: [Integer]
raicesEnterasDePrimos3 = aux primes [1..]
    where aux (p:ps) (x:xs) | p > x*x   = aux (p:ps) xs
                            | otherwise = (x-1) : aux ps (x:xs)
 
 
-- Comparación de eficiencia
--    ghci> raicesEnterasDePrimos1 !! 400000
--    2408
--    (2.86 secs, 1177922500 bytes)
--    ghci> raicesEnterasDePrimos2 !! 400000
--    2408
--    (3.08 secs, 1177432260 bytes)
--    ghci> raicesEnterasDePrimos3 !! 400000
--    2408
--    (3.88 secs, 1260772112 bytes)
 
-- En lo sucesivo usaremos la 1ª definición
raicesEnterasDePrimos :: [Integer]
raicesEnterasDePrimos = raicesEnterasDePrimos3
 
-- La propiedad es
prop_raicesEnterasDePrimos :: Int -> Property
prop_raicesEnterasDePrimos n =
    n >= 0 ==> 
    raicesEnterasDePrimos !! (n+1) - raicesEnterasDePrimos !! n <= 1
 
-- La comprobación es
--    λ> quickCheck prop_raicesEnterasDePrimos
--    +++ OK, passed 100 tests.

Números libres de cuadrados

Un número es libre de cuadrados si no es divisible el cuadrado de ningún entero mayor que 1. Por ejemplo, 70 es libre de cuadrado porque sólo es divisible por 1, 2, 5, 7 y 70; en cambio, 40 no es libre de cuadrados porque es divisible por 2².

Definir la función

   libreDeCuadrados :: Integer -> Bool

tal que (libreDeCuadrados x) se verifica si x es libre de cuadrados. Por ejemplo,

   libreDeCuadrados 70                 ==  True
   libreDeCuadrados 40                 ==  False
   libreDeCuadrados 510510             ==  True
   libreDeCuadrados (((10^10)^10)^10)  ==  False

Soluciones

-- 1ª definición:
libreDeCuadrados :: Integer -> Bool
libreDeCuadrados x = x == product (divisoresPrimos x)
 
-- (divisoresPrimos x) es la lista de los divisores primos de x. Por
-- ejemplo,  
--    divisoresPrimos 40  ==  [2,5]
--    divisoresPrimos 70  ==  [2,5,7]
divisoresPrimos :: Integer -> [Integer]
divisoresPrimos x = [n | n <- divisores x, primo n]
 
-- (divisores n) es la lista de los divisores del número n. Por ejemplo,
--    divisores 30  ==  [1,2,3,5,6,10,15,30]  
divisores :: Integer -> [Integer]
divisores n = [x | x <- [1..n], n `mod` x == 0]
 
-- (primo n) se verifica si n es primo. Por ejemplo,
--    primo 30  == False
--    primo 31  == True  
primo :: Integer -> Bool
primo n = divisores n == [1, n]
 
-- 2ª definición
libreDeCuadrados2 :: Integer -> Bool
libreDeCuadrados2 n = 
    null [x | x <- [2..n], rem n (x^2) == 0]
 
-- 3ª definición
libreDeCuadrados3 :: Integer -> Bool
libreDeCuadrados3 n = 
    null [x | x <- [2..floor (sqrt (fromIntegral n))], 
              rem n (x^2) == 0]
 
-- Comparación de eficiencia
-- =========================
 
--    λ> libreDeCuadrados 510510
--    True
--    (0.76 secs, 89,522,360 bytes)
--    λ> libreDeCuadrados2 510510
--    True
--    (1.78 secs, 371,826,320 bytes)
--    λ> libreDeCuadrados3 510510
--    True
--    (0.01 secs, 0 bytes)

Perímetro más frecuente de triángulos rectángulos

El grado perimetral de un número p es la cantidad de tres triángulos rectángulos de lados enteros cuyo perímetro es p. Por ejemplo, el grado perimetral de 120 es 3 ya que sólo hay 3 triángulos rectángulos de lados enteros cuyo perímetro es 120: {20,48,52}, {24,45,51} y {30,40,50}.

Definir la función

   maxGradoPerimetral :: Int -> (Int,[Int])

tal que (maxGradoPerimetral n) es el par (m,ps) tal que m es el máximo grado perimetral de los números menores o iguales que n y ps son los perímetros, menores o iguales que n, cuyo grado perimetral es m. Por ejemplo,

   maxGradoPerimetral   50  ==  (1,[12,24,30,36,40,48])
   maxGradoPerimetral  100  ==  (2,[60,84,90])
   maxGradoPerimetral  200  ==  (3,[120,168,180])
   maxGradoPerimetral  400  ==  (4,[240,360])
   maxGradoPerimetral  500  ==  (5,[420])
   maxGradoPerimetral  750  ==  (6,[720])
   maxGradoPerimetral  839  ==  (6,[720])
   maxGradoPerimetral  840  ==  (8,[840])
   maxGradoPerimetral 1500  ==  (8,[840,1260])
   maxGradoPerimetral 2000  ==  (10,[1680])
   maxGradoPerimetral 3000  ==  (12,[2520])

Soluciones

import Data.List
import Data.Array (accumArray, assocs)
 
-- 1ª solución                                                      --
-- ===========
 
maxGradoPerimetral1 :: Int -> (Int,[Int])
maxGradoPerimetral1 p = (m,[x | (n,x) <- ts, n == m])
    where ts    = [(length (triangulos x),x) | x <- [1..p]] 
          (m,_) = maximum ts 
 
-- (triangulos p) es el conjunto de triángulos rectángulos de perímetro
-- p. Por ejemplo,
--    triangulos 120  ==  [(20,48,52),(24,45,51),(30,40,50)]
triangulos :: Int -> [(Int,Int,Int)]
triangulos p = 
    [(a,b,c) | a <- [1..q],
               b <- [a..q],
               let c = p-a-b,
               a*a+b*b == c*c]
    where q = p `div` 2
 
-- 2ª solución                                                      --
-- ===========
 
maxGradoPerimetral2 :: Int -> (Int,[Int])
maxGradoPerimetral2 p = (m,[x | (n,x) <- ts, n == m])
    where ts    = [(n,x) | (x,n) <- numeroPerimetrosTriangulos p, n > 0]
          (m,_) = maximum ts 
 
-- (numeroPerimetrosTriangulos p) es la lista formado por los números de
-- 1 a p y la cantidad de triángulos rectángulos enteros cuyo perímetro
-- es dicho número. Por ejemplo,
--    ghci>  [(p,n) | (p,n) <- numeroPerimetrosTriangulos 70, n > 0]
--    [(12,1),(24,1),(30,1),(36,1),(40,1),(48,1),(56,1),(60,2),(70,1)]
numeroPerimetrosTriangulos :: Int -> [(Int,Int)] 
numeroPerimetrosTriangulos p = 
    assocs (accumArray (\x _ -> 1+x) 0 (1,p) (perimetrosTriangulos p))
 
-- (perimetrosTriangulos p) es la lista formada por los perímetros y los
-- lados de los triángulos rectángulos enteros cuyo perímetro es menor o
-- igual que p. Por ejemplo,
--    ghci> perimetrosTriangulos 70
--    [(12,(3,4,5)),   (30,(5,12,13)),(24,(6,8,10)),  (56,(7,24,25)),
--     (40,(8,15,17)), (36,(9,12,15)),(60,(10,24,26)),(48,(12,16,20)),
--     (60,(15,20,25)),(70,(20,21,29))]
perimetrosTriangulos :: Int -> [(Int,(Int,Int,Int))]
perimetrosTriangulos p =
    [(q,(a,b,c)) | a <- [1..p1],
                   b <- [a..p1],
                   esCuadrado (a*a+b*b), 
                   let c = raizCuadradaE (a*a+b*b), 
                   let q = a+b+c,
                   q <= p]
    where p1 = p `div` 2
 
-- (esCuadrado n) se verifica si n es un cuadrado. Por ejemplo,
--    esCuadrado 25  ==  True
--    esCuadrado 27  ==  False
esCuadrado :: Int -> Bool
esCuadrado n = a*a == n
    where a = raizCuadradaE n
 
-- (raizCuadradaE n) es la raíz cuadrada entera de n. Por ejemplo,
--    raizCuadradaE 25  ==  5
--    raizCuadradaE 27  ==  5
--    raizCuadradaE 35  ==  5
--    raizCuadradaE 36  ==  6
raizCuadradaE :: Int -> Int
raizCuadradaE = floor . sqrt . fromIntegral
 
-- 3ª solución                                                      --
-- ===========
 
maxGradoPerimetral3 :: Int -> (Int,[Int])
maxGradoPerimetral3 p = (m,[x | (n,x) <- ts, n == m])
    where ts    = [(n,x) | (x,n) <- numeroPerimetrosTriangulos2 p, n > 0]
          (m,_) = maximum ts 
 
-- (numeroPerimetrosTriangulos2 p) es la lista formado por los números de
-- 1 a p y la cantidad de triángulos rectángulos enteros cuyo perímetro
-- es dicho número. Por ejemplo,
--    ghci>  [(p,n) | (p,n) <- numeroPerimetrosTriangulos2 70, n > 0]
--    [(12,1),(24,1),(30,1),(36,1),(40,1),(48,1),(56,1),(60,2),(70,1)]
numeroPerimetrosTriangulos2 :: Int -> [(Int,Int)] 
numeroPerimetrosTriangulos2 p = 
    [(head xs, length xs) | xs <- group (sort (perimetrosTriangulos2 p))]
 
-- (perimetrosTriangulos2 p) es la lista formada por los perímetros de
-- los triángulos rectángulos enteros cuyo perímetro es menor o igual
-- que p. Por ejemplo, 
--    perimetrosTriangulos2 70  ==  [12,30,24,56,40,36,60,48,60,70]
perimetrosTriangulos2 :: Int -> [Int]
perimetrosTriangulos2 p =
    [q | a <- [1..p1],
         b <- [a..p1],
         esCuadrado (a*a+b*b), 
         let c = raizCuadradaE (a*a+b*b), 
         let q = a+b+c,
         q <= p]
    where p1 = p `div` 2
 
-- Comparación de eficiencia
-- =========================
 
-- La comparación es
--    ghci> maxGradoPerimetral1 1000
--    (8,[840])
--    (120.08 secs, 21116625136 bytes)
--    ghci> maxGradoPerimetral2 1000
--    (8,[840])
--    (0.66 secs, 132959056 bytes)
--    ghci> maxGradoPerimetral3 1000
--    (1000,[1])
--    (0.66 secs, 133443816 bytes)