Menu Close

Órbita con raíz entera (OME1997 P4)

El enunciado del problema 4 de la OME (Olimpiada Matemática Española) del 1997 es

Sea p un número primo. Determinar todos los enteros k tales que sqrt(k² – k*p) es natural.

Definir las funciones

   orbita        :: Integer -> [Integer]
   orbitaDePrimo :: Integer -> [Integer]

tales que

  • (orbita n) es la lista de todos los enteros k tales que sqrt(k² – k*n) es natural. Por ejemplo,
     take 4  (orbita 6)   == [0,-2,6,8]
     take 5  (orbita 36)  == [0,-12,36,48,-64]
     take 6  (orbita 9)   == [0,-3,9,12,-16,25]
     take 8  (orbita 27)  == [0,-9,27,36,-48,75,-169,196]
     take 10 (orbita 111) == [0,-37,111,148,-289,400,-972,1083,-3025,3136]
  • (orbitaDePrimo p) es la lista de todos los enteros k tales que sqrt(k² – k*p) es natural, suponiendo que p es un número primo. Por ejemplo,
     orbitaDePrimo 5                  == [0,-4,5,9]
     orbitaDePrimo (primes !! (10^6)) == [0,15485867,-59953011442489,59953026928356]

Soluciones

import Data.Numbers.Primes (primes)
 
orbita :: Integer -> [Integer]
orbita n =
  [k | k <- enteros,
       k^2 - k*n >= 0,
       esCuadrado (k^2 - k*n)]
 
-- entero es la lista de los números enteros. Por ejemplo,
--    λ> take 20 enteros
--    [0,-1,1,-2,2,-3,3,-4,4,-5,5,-6,6,-7,7,-8,8,-9,9,-10]
enteros :: [Integer]
enteros = 0 : concat [[-n,n] | n <- [1..]]
 
-- (esCuadrado x) se verifica si x es un cuadrado perfecto. Por
-- ejemplo,
--    esCuadrado 16  ==  True
--    esCuadrado 27  ==  False
esCuadrado :: Integer -> Bool
esCuadrado x =
  (raizEntera x)^2 == x
 
-- (raizEntera x) es el mayor entero cuyo cuadrado es menor o igual que
-- x. Por ejemplo,
--    raizEntera 16  ==  4
--    raizEntera 27  ==  5
raizEntera :: Integer -> Integer
raizEntera x = aux (1,x)
    where aux (a,b) | d == x    = c
                    | c == a    = c
                    | d < x     = aux (c,b)
                    | otherwise = aux (a,c)
              where c = (a+b) `div` 2
                    d = c^2
 
-- 1ª definición de orbitaDePrimo
-- ==============================
 
orbitaDePrimo1 :: Integer -> [Integer]
orbitaDePrimo1 2 = take 2 (orbita 2)
orbitaDePrimo1 p = take 4 (orbita p)
 
-- 2ª definición de orbitaDePrimo
-- ==============================
 
-- Basada en los siguientes cálculos
--    orbitaDePrimo1 2  == [0,2]
--    orbitaDePrimo1 3  == [0,-1,3,4]
--    orbitaDePrimo1 5  == [0,-4,5,9]
--    orbitaDePrimo1 7  == [0, 7,  -9, 16]
--    orbitaDePrimo1 11 == [0,11, -25, 36]
--    orbitaDePrimo1 13 == [0,13, -36, 49]
--    orbitaDePrimo1 17 == [0,17, -64, 81]
--    orbitaDePrimo1 19 == [0,19, -81,100]
--    orbitaDePrimo1 23 == [0,23,-121,144]
 
orbitaDePrimo2 :: Integer -> [Integer]
orbitaDePrimo2 p
  | p == 2    = [0,2]
  | p <= 5    = [0, -((p-1) `div` 2)^2, p, ((p+1) `div` 2)^2]
  | otherwise = [0, p, -((p-1) `div` 2)^2, ((p+1) `div` 2)^2]
 
-- Comprobación de equivalencia
-- ============================
 
-- La comprobación es
--    λ> and [orbitaDePrimo1 n == orbitaDePrimo2 n | n <- take 30 primes]
--    True
 
-- Comparación de eficiencia
-- =========================
 
-- La comparación es
--    λ> orbitaDePrimo1 (primes !! 100)
--    [0,547,-74529,75076]
--    (4.94 secs, 4,471,368,256 bytes)
--    λ> orbitaDePrimo2 (primes !! 100)
--    [0,547,-74529,75076]
--    (0.01 secs, 302,096 bytes)

Nuevas soluciones

  • En los comentarios se pueden escribir nuevas soluciones.
  • El código se debe escribir entre una línea con <pre lang="haskell"> y otra con </pre>

Una solución de “Órbita con raíz entera (OME1997 P4)

  1. j0sejuan

    (Por comodidad se usa aritmética real exacta, buscar directamente las raíces enteras sería mucho más raṕido)

    type Re = CReal 1
     
    {-
      Quitando las soluciones obvias
     
        0^2 - 0 p = 0^2
        p^2 - p p = 0^2
     
      Ha de ser
     
        k^2 - k p = z^2
     
      es decir
     
        k = (p +/- sqrt(p^2 + 4z^2)) / 2
     
      debe ser entonces
     
        p^2 + (2z)^2 = w^2
     
      que es una terna pitagórica primitiva por `p`.
     
      Ha de ser entonces
     
        p = m^2 - n^2
     
      una cota para m es que
     
        m^2 - (m - 1)^2 <= p
     
      de donde
     
        m <= (p + 1) / 2
     
      sacamos
     
        n = sqrt(m^2 - p)
     
      y es
     
        z = m n
     
      y de ahí k.
    -}
     
    órbitaPrima :: Integer -> [Integer]
    órbitaPrima p = [k | m <- [minM..maxM]
                       , let m2 = m^2
                       , m2 >= p
                       , let n = round $ sqrt (fromIntegral m2 - p')
                       , p == m^2 - n^2
                       , let z = m * n
                       , let r = round $ sqrt (p'**2 + 4 * (fromIntegral (z^2)))
                       , s <- [-1, 1]
                       , let k = (p + s * r) `div` 2
                       , k^2 - k * p == z^2]
      where p' = fromIntegral p :: Re
            minM = floor   $ sqrt p'
            maxM = ceiling $ (p' + 1) / 2
     
    {-
     
    *Main> órbitaPrima 15485867
    [-59953011442489,59953026928356]
     
    -}

Leave a Reply

Este sitio usa Akismet para reducir el spam. Aprende cómo se procesan los datos de tus comentarios.