Menu Close

Etiqueta: sqrt

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