El teorema de cuatro cuadrados de Lagrange nos dice que cualquier número natural puede representarse como la suma de cuatro números cuadrados. Su tarea es escribir un programa que haga esto.

Entrada: un número natural (menos de mil millones)

Salida: cuatro números cuyos cuadrados suman ese número (el orden no importa)

Nota: ¡No tiene que hacer una búsqueda de fuerza bruta! Detalles aquí y aquí . Si hay una función que trivializa este problema (lo determinaré), no está permitido. Se permiten funciones primarias automatizadas y raíz cuadrada. Si hay más de una representación, cualquiera está bien. Si elige hacer fuerza bruta, debe ejecutarse dentro de un tiempo razonable (3 minutos)

entrada de muestra

123456789

salida de muestra (cualquiera está bien)

10601 3328 2 0
10601 3328 2

CJam, 50 bytes

li:NmF0=~2/#:J(_{;)__*N\-[{_mqJ/iJ*__*@\-}3*])}g+p

Mi tercera (y última, lo prometo) respuesta. Este enfoque se basa en gran medida en la respuesta de primo .

Pruébelo en línea en el intérprete de CJam .

Uso

$ cjam 4squares.cjam <<< 999999999
[189 31617 567 90]

Antecedentes

1. Después de ver algoritmo de Primo actualizado, que tenía que ver cómo tendría una puntuación de aplicación Cjam:

   li{W):W;:N4md!}g;Nmqi)_{;(__*N\-[{_mqi__*@\-}3*])}g+2W#f*p
   

   ¡Solo 58 bytes! Este algoritmo se realiza en un tiempo casi constante y no presenta mucha variación para diferentes valores de N. Cambiemos eso ...

2. En lugar de comenzar en floor(sqrt(N))y disminuir, podemos comenzar en 1e incrementar. Esto ahorra 4 bytes.

   li{W):W;:N4md!}g;0_{;)__*N\-[{_mqi__*@\-}3*])}g+2W#f*p
   

3. En lugar de expresar Ncomo 4**a * b, podemos expresarlo como p**(2a) * b- donde pes el factor primo más pequeño de N- para guardar 1 byte más.

   li_mF0=~2/#:J_*/:N!_{;)__*N\-[{_mqi__*@\-}3*])}g+Jf*p
   

4. La modificación anterior nos permite cambiar ligeramente la implementación (sin tocar el algoritmo mismo): en lugar de dividir Npor p**(2a)y multiplicar la solución por p**a, podemos restringir directamente las posibles soluciones a múltiplos de p**a. Esto ahorra 2 bytes más.

   li:NmF0=~2/#:J!_{;J+__*N\-[{_mqJ/iJ*__*@\-}3*])}g+`
   

5. Sin restringir el primer entero a múltiplos de p**aguarda un byte adicional.

   li:NmF0=~2/#:J(_{;)__*N\-[{_mqJ/iJ*__*@\-}3*])}g+`
   

Algoritmo final

1. Encuentre ay btal que N = p**(2a) * b, donde bno sea múltiplo de p**2y psea ​​el factor primo más pequeño de N.

2. Conjunto j = p**a.

3. Conjunto k = floor(sqrt(N - j**2) / A) * A.

4. Conjunto l = floor(sqrt(N - j**2 - k**2) / A) * A.

5. Conjunto m = floor(sqrt(N - j**2 - k**2 - l**2) / A) * A.

6. Si N - j**2 - k**2 - l**2 - m**2 > 0, configurar j = j + 1y volver al paso 3.

Esto se puede implementar de la siguiente manera:

li:N          " Read an integer from STDIN and save it in “N”.                        ";
mF            " Push the factorization of “N”. Result: [ [ p1 a1 ] ... [ pn an ] ]    ";
0=~           " Push “p1” and “a1”. “p1” is the smallest prime divisor of “N”.        ";
2/#:J         " Compute p1**(a1/2) and save the result “J”.                           ";
(_            " Undo the first two instructions of the loop.                          ";
{             "                                                                       ";
  ;)_         " Pop and discard. Increment “J” and duplicate.                         ";
  _*N\-       " Compute N - J**2.                                                     ";
  [{          "                                                                       ";
    _mqJ/iJ*  " Compute K = floor(sqrt(N - J**2)/J)*J.                                ";
    __*@      " Duplicate, square and rotate. Result: K   K**2   N - J**2             ";
    \-        " Swap and subtract. Result: K   N - J**2 - K**2                        ";
  }3*]        " Do the above three times and collect in an array.                     ";
  )           " Pop the array. Result: N - J**2 - K**2 - L**2 - M**2                  ";
}g            " If the result is zero, break the loop.                                ";
+p            " Unshift “J” in [ K L M ] and print a string representation.           ";

Puntos de referencia

Ejecuté las 5 versiones en todos los enteros positivos hasta 999,999,999 en mi Intel Core i7-3770, medí el tiempo de ejecución y conté las iteraciones necesarias para encontrar una solución.

La siguiente tabla muestra el número promedio de iteraciones y el tiempo de ejecución para un solo entero:

Version               |    1    |    2    |    3    |    4    |    5
----------------------+---------+---------+---------+---------+---------
Number of iterations  |  4.005  |  28.31  |  27.25  |  27.25  |  41.80
Execution time [µs]   |  6.586  |  39.69  |  55.10  |  63.99  |  88.81

1. Con solo 4 iteraciones y 6.6 microsegundos por entero, el algoritmo de primo es increíblemente rápido.

2. Comenzar en floor(sqrt(N))tiene más sentido, ya que esto nos deja con valores más pequeños para la suma de los tres cuadrados restantes. Como se esperaba, comenzar en 1 es mucho más lento.

3. Este es un ejemplo clásico de una buena idea mal implementada. Para reducir realmente el tamaño del código, confiamos en mFque factoriza el número entero N. Aunque la versión 3 requiere menos iteraciones que la versión 2, es mucho más lenta en la práctica.

4. Aunque el algoritmo no cambia, la versión 4 es mucho más lenta. Esto se debe a que realiza una división de coma flotante adicional y una multiplicación entera en cada iteración.

5. Para la entrada N = p**(2a) ** b, el algoritmo 5 requerirá (k - 1) * p**a + 1iteraciones, donde kes la cantidad de iteraciones que requiere el algoritmo 4. Si k = 1o a = 0, esto no hace ninguna diferencia.

   Sin embargo, cualquier entrada del formulario 4**a * (4**c * (8 * d + 7) + 1)puede funcionar bastante mal. Para el valor de inicio j = p**a, N - 4**a = 4**(a + c) * (8 * d + 7), por lo que no se puede expresar como una suma de tres cuadrados. Por lo tanto, k > 1y al menos p**ase requieren iteraciones.

   Afortunadamente, el algoritmo original de primo es increíblemente rápido y N < 1,000,000,000. El peor caso que pude encontrar a mano es 265,289,728 = 4**10 * (4**1 * (7 * 8 + 7) + 1), que requiere 6,145 iteraciones. El tiempo de ejecución es inferior a 300 ms en mi máquina. En promedio, esta versión es 13.5 veces más lenta que la implementación del algoritmo de primo.

FRACTRAN: 156 98 fracciones

Como este es un problema clásico de la teoría de números, ¡qué mejor manera de resolverlo que usar números!

37789/221 905293/11063 1961/533 2279/481 57293/16211 2279/611 53/559 1961/403 53/299 13/53 1/13 6557/262727 6059/284321 67/4307 67/4661 6059/3599 59/83 1/59 14279/871933 131/9701 102037079/8633 14017/673819 7729/10057 128886839/8989 13493/757301 7729/11303 89/131 1/89 31133/2603 542249/19043 2483/22879 561731/20413 2483/23701 581213/20687 2483/24523 587707/21509 2483/24797 137/191 1/137 6215941/579 6730777/965 7232447/1351 7947497/2123 193/227 31373/193 23533/37327 5401639/458 229/233 21449/229 55973/24823 55973/25787 6705901/52961 7145447/55973 251/269 24119/251 72217/27913 283/73903 281/283 293/281 293/28997 293/271 9320827/58307 9831643/75301 293/313 28213/293 103459/32651 347/104807 347/88631 337/347 349/337 349/33919 349/317 12566447/68753 13307053/107143 349/367 33197/349 135199/38419 389/137497 389/119113 389/100729 383/389 397/383 397/39911 397/373 1203/140141 2005/142523 2807/123467 4411/104411 802/94883 397/401 193/397 1227/47477 2045/47959 2863/50851 4499/53743 241/409 1/241 1/239

Toma entrada de la forma 2 ⁿ × 193 y emite 3 ^a × 5 ^b × 7 ^c × 11 ^d . Podría correr en 3 minutos si tiene un intérprete realmente bueno. Tal vez.

... bueno, en realidad no. Esto parece ser un problema tan divertido que hacer en FRACTRAN que tenía que probarlo. Obviamente, esta no es una solución adecuada a la pregunta, ya que no requiere tiempo (fuerza bruta) y apenas se juega al golf, pero pensé en publicar esto aquí porque no todos los días es una pregunta de Codegolf se puede hacer en FRACTRAN;)

Insinuación

El código es equivalente al siguiente pseudo-Python:

a, b, c, d = 0, 0, 0, 0

def square(n):
    # Returns n**2

def compare(a, b):
    # Returns (0, 0) if a==b, (1, 0) if a<b, (0, 1) if a>b

def foursquare(a, b, c, d):
    # Returns square(a) + square(b) + square(c) + square(d)

while compare(foursquare(a, b, c, d), n) != (0, 0):
    d += 1

    if compare(c, d) == (1, 0):
        c += 1
        d = 0

    if compare(b, c) == (1, 0):
        b += 1
        c = 0
        d = 0

    if compare(a, b) == (1, 0):
        a += 1
        b = 0
        c = 0
        d = 0

Mathematica 61 66 51

Se muestran tres métodos. Solo el primer enfoque cumple el requisito de tiempo.

1-FindInstance (51 caracteres)

Esto devuelve una única solución a la ecuación.

FindInstance[a^2 + b^2 + c^2 + d^2 == #, {a, b, c, d}, Integers] &

Ejemplos y horarios

FindInstance[a^2 + b^2 + c^2 + d^2 == 123456789, {a, b, c, d}, Integers] // AbsoluteTiming

{0.003584, {{a -> 2600, b -> 378, c -> 10468, d -> 2641}}}

FindInstance[a^2 + b^2 + c^2 + d^2 == #, {a, b, c, d}, Integers] &[805306368]

{0.004437, {{a -> 16384, b -> 16384, c -> 16384, d -> 0}}}

Particiones de 2 enteros

Esto también funciona, pero es demasiado lento para cumplir con el requisito de velocidad.

f@n_ := Sqrt@IntegerPartitions[n, {4}, Range[0, Floor@Sqrt@n]^2, 1][[1]]

Range[0, Floor@Sqrt@n]^2es el conjunto de todos los cuadrados menores que la raíz cuadrada de n(el cuadrado más grande posible en la partición).

{4}requiere que las particiones enteras nconstan de 4 elementos del conjunto de cuadrados mencionado anteriormente.

1, dentro de la función IntegerPartitionsdevuelve la primera solución.

[[1]]elimina los tirantes externos; La solución se devolvió como un conjunto de un elemento.

f[123456]

{348, 44, 20, 4}

Representaciones de 3 poderes

PowerRepresentations devuelve todas las soluciones al problema de los 4 cuadrados. También puede resolver sumas de otros poderes.

PowersRepresentations devuelve, en menos de 5 segundos, las 181 formas de expresar 123456789 como la suma de 4 cuadrados:

n= 123456;
PowersRepresentations[n, 4, 2] //AbsoluteTiming

Sin embargo, es demasiado lento para otras sumas.

Perl - 116 bytes 87 bytes (ver actualización a continuación)

#!perl -p
$.<<=1,$_>>=2until$_&3;
{$n=$_;@a=map{$n-=$a*($a-=$_%($b=1|($a=0|sqrt$n)>>1));$_/=$b;$a*$.}($j++)x4;$n&&redo}
$_="@a"

Contando el shebang como un byte, se agregaron nuevas líneas para la cordura horizontal.

Algo de una combinación de código de golf más rápido de envío de código .

La complejidad del caso promedio (¿peor?) Parece ser O (log n) O (n ^0.07 ) . Nada de lo que he encontrado funciona más lento que 0.001s, y he verificado todo el rango de 900000000 a 999999999 . Si encuentra algo que lleva mucho más tiempo que eso, ~ 0.1s o más, hágamelo saber.

Uso de muestra

$ echo 123456789 | timeit perl four-squares.pl
11110 157 6 2

Elapsed Time:     0:00:00.000

$ echo 1879048192 | timeit perl four-squares.pl
32768 16384 16384 16384

Elapsed Time:     0:00:00.000

$ echo 999950883 | timeit perl four-squares.pl
31621 251 15 4

Elapsed Time:     0:00:00.000

Los dos últimos de estos parecen ser los escenarios del peor de los casos para otras presentaciones. En ambos casos, la solución que se muestra es, literalmente, lo primero que se verifica. por123456789 , es el segundo.

Si desea probar un rango de valores, puede usar el siguiente script:

use Time::HiRes qw(time);

$t0 = time();

# enter a range, or comma separated list here
for (1..1000000) {
  $t1 = time();
  $initial = $_;
  $j = 0; $i = 1;
  $i<<=1,$_>>=2until$_&3;
  {$n=$_;@a=map{$n-=$a*($a-=$_%($b=1|($a=0|sqrt$n)>>1));$_/=$b;$a*$i}($j++)x4;$n&&redo}
  printf("%d: @a, %f\n", $initial, time()-$t1)
}
printf('total time: %f', time()-$t0);

Mejor cuando se canaliza a un archivo. El rango 1..1000000toma alrededor de 14 segundos en mi computadora (71000 valores por segundo), y el rango 999000000..1000000000toma alrededor de 20 segundos (50000 valores por segundo), de acuerdo con la complejidad promedio de O (log n) .

Actualizar

Editar : Resulta que este algoritmo es muy similar a uno que ha sido utilizado por calculadoras mentales durante al menos un siglo .

Desde la publicación original, he verificado todos los valores en el rango de 1..1000000000 . El comportamiento del "peor de los casos" fue exhibido por el valor 699731569 , que probó un total de 190 combinaciones antes de llegar a una solución. Si considera que 190 es una constante pequeña, y ciertamente lo hago, el peor de los casos en el rango requerido puede considerarse O (1) . Es decir, tan rápido como buscar la solución desde una mesa gigante, y en promedio, posiblemente más rápido.

Sin embargo, otra cosa. Después de 190 iteraciones, cualquier cosa mayor que 144400 ni siquiera ha superado el primer paso. La lógica del primer recorrido transversal no tiene valor, ni siquiera se usa. El código anterior se puede acortar bastante:

#!perl -p
$.*=2,$_/=4until$_&3;
@a=map{$=-=$%*($%=$=**.5-$_);$%*$.}$j++,(0)x3while$=&&=$_;
$_="@a"

Que solo realiza el primer paso de la búsqueda. Sin embargo, debemos confirmar que no hay valores por debajo de 144400 que necesiten la segunda pasada:

for (1..144400) {
  $initial = $_;

  # reset defaults
  $.=1;$j=undef;$==60;

  $.*=2,$_/=4until$_&3;
  @a=map{$=-=$%*($%=$=**.5-$_);$%*$.}$j++,(0)x3while$=&&=$_;

  # make sure the answer is correct
  $t=0; $t+=$_*$_ for @a;
  $t == $initial or die("answer for $initial invalid: @a");
}

En resumen, para el rango 1..1000000000 , existe una solución de tiempo casi constante, y la está viendo.

Actualización actualizada

@Dennis y yo hemos realizado varias mejoras en este algoritmo. Puede seguir el progreso en los comentarios a continuación y la discusión posterior, si eso le interesa. El número promedio de iteraciones para el rango requerido se ha reducido de poco más de 4 a 1.229 , y el tiempo necesario para probar todos los valores para 1..1000000000 se ha mejorado de 18m 54s, a 2m 41s. El peor de los casos requería anteriormente 190 iteraciones; El peor de los casos ahora, 854382778 , solo necesita 21 .

El código final de Python es el siguiente:

from math import sqrt

# the following two tables can, and should be pre-computed

qqr_144 = set([  0,   1,   2,   4,   5,   8,   9,  10,  13,
                16,  17,  18,  20,  25,  26,  29,  32,  34,
                36,  37,  40,  41,  45,  49,  50,  52,  53,
                56,  58,  61,  64,  65,  68,  72,  73,  74,
                77,  80,  81,  82,  85,  88,  89,  90,  97,
                98, 100, 101, 104, 106, 109, 112, 113, 116,
               117, 121, 122, 125, 128, 130, 133, 136, 137])

# 10kb, should fit entirely in L1 cache
Db = []
for r in range(72):
  S = bytearray(144)
  for n in range(144):
    c = r

    while True:
      v = n - c * c
      if v%144 in qqr_144: break
      if r - c >= 12: c = r; break
      c -= 1

    S[n] = r - c
  Db.append(S)

qr_720 = set([  0,   1,   4,   9,  16,  25,  36,  49,  64,  81, 100, 121,
              144, 145, 160, 169, 180, 196, 225, 241, 244, 256, 265, 289,
              304, 324, 340, 361, 369, 385, 400, 409, 436, 441, 481, 484,
              496, 505, 529, 544, 576, 580, 585, 601, 625, 640, 649, 676])

# 253kb, just barely fits in L2 of most modern processors
Dc = []
for r in range(360):
  S = bytearray(720)
  for n in range(720):
    c = r

    while True:
      v = n - c * c
      if v%720 in qr_720: break
      if r - c >= 48: c = r; break
      c -= 1

    S[n] = r - c
  Dc.append(S)

def four_squares(n):
  k = 1
  while not n&3:
    n >>= 2; k <<= 1

  odd = n&1
  n <<= odd

  a = int(sqrt(n))
  n -= a * a
  while True:
    b = int(sqrt(n))
    b -= Db[b%72][n%144]
    v = n - b * b
    c = int(sqrt(v))
    c -= Dc[c%360][v%720]
    if c >= 0:
      v -= c * c
      d = int(sqrt(v))

      if v == d * d: break

    n += (a<<1) - 1
    a -= 1

  if odd:
    if (a^b)&1:
      if (a^c)&1:
        b, c, d = d, b, c
      else:
        b, c = c, b

    a, b, c, d = (a+b)>>1, (a-b)>>1, (c+d)>>1, (c-d)>>1

  a *= k; b *= k; c *= k; d *= k

  return a, b, c, d

Utiliza dos tablas de corrección calculadas previamente, una de 10 kb de tamaño y la otra de 253 kb. El código anterior incluye las funciones generadoras para estas tablas, aunque probablemente deberían calcularse en tiempo de compilación.

Puede encontrar una versión con tablas de corrección de tamaño más modesto aquí: http://codepad.org/1ebJC2OV Esta versión requiere un promedio de 1.620 iteraciones por término, con un peor caso de 38 , y el rango completo se ejecuta en aproximadamente 3 m 21 s. Se compensa un poco de tiempo, utilizando bitwise andpara la corrección b , en lugar del módulo.

Mejoras

Los valores pares tienen más probabilidades de producir una solución que los valores impares.
El artículo de cálculo mental vinculado a notas anteriores señala que si, después de eliminar todos los factores de cuatro, el valor a descomponer es par, este valor se puede dividir por dos, y la solución reconstruida:

Aunque esto podría tener sentido para el cálculo mental (los valores más pequeños tienden a ser más fáciles de calcular), no tiene mucho sentido algorítmicamente. Si toma 256 tuplas aleatorias de 4 y examina la suma de los cuadrados del módulo 8 , encontrará que los valores 1 , 3 , 5 y 7 se alcanzan en promedio 32 veces. Sin embargo, los valores 2 y 6 se alcanzan cada uno 48 veces. Multiplicando valores impares por 2 encontrará una solución, en promedio, en 33% menos iteraciones. La reconstrucción es la siguiente:

El cuidado necesita ser tomado que una y B tienen la misma paridad, así como c y d , pero si se encontró una solución en absoluto, una ordenación adecuada está garantizado para existir.

No es necesario verificar las rutas imposibles.
Después de seleccionar el segundo valor, b , puede que ya sea imposible que exista una solución, dados los posibles residuos cuadráticos para cualquier módulo dado. En lugar de verificar de todos modos, o pasar a la siguiente iteración, el valor de b puede 'corregirse' disminuyéndolo en la cantidad más pequeña que pueda conducir a una solución. Las dos tablas de corrección almacenan estos valores, uno para b y el otro para c . Usar un módulo más alto (más exactamente, usar un módulo con relativamente menos residuos cuadráticos) dará como resultado una mejor mejora. El valor a no necesita ninguna corrección; modificando n para que sea par, todos los valores deuna son válidos.

Pitón 3 (177)

N=int(input())
k=1
while N%4<1:N//=4;k*=2
n=int(N**.5)
R=range(int(2*n**.5)+1)
print([(a*k,b*k,c*k,d*k)for d in R for c in R for b in R for a in[n,n-1]if a*a+b*b+c*c+d*d==N][0])

Después de reducir la entrada Npara que no sea divisible por 4, debe ser expresable como una suma de cuatro cuadrados donde uno de ellos sea el valor más grande posible a=int(N**0.5)o uno menor que eso, dejando solo un pequeño resto para la suma de los otros tres cuadrados. para cuidar de. Esto reduce en gran medida el espacio de búsqueda.

Aquí hay una prueba de que este código siempre encuentra una solución. Deseamos encontrar un valor aque n-a^2sea ​​la suma de tres cuadrados. Del teorema de tres cuadrados de Legendre , un número es la suma de tres cuadrados a menos que sea la forma 4^j(8*k+7). En particular, tales números son 0 o 3 (módulo 4).

Mostramos que no hay dos consecutivos a puedan hacer que la cantidad sobrante N-a^2tenga esa forma para ambos valores consecutivos. Podemos hacerlo simplemente haciendo una tabla de ay Nmódulo 4, notando eso N%4!=0porque hemos extraído todas las potencias de 4 N.

  a%4= 0123
      +----
     1|1010
N%4= 2|2121  <- (N-a*a)%4
     3|3232

Debido a que no hay dos consecutiva adar (N-a*a)%4 in [0,3], uno de ellos es seguro de usar. Entonces, utilizamos con avidez la mayor cantidad posible nconn^2<=N , y n-1. Como N<(n+1)^2, el resto N-a^2que se representará como una suma de tres cuadrados es como máximo (n+1)^2 -(n-1)^2, lo que es igual 4*n. Por lo tanto, es suficiente verificar solo los valores hasta 2*sqrt(n), que es exactamente el rango R.

Uno podría optimizar aún más el tiempo de ejecución deteniéndose después de una única solución, calculando en lugar de iterar por el último valor dy buscando solo entre los valores b<=c<=d. Pero, incluso sin estas optimizaciones, la peor instancia que pude encontrar terminó en 45 segundos en mi máquina.

La cadena de "para x en R" es lamentable. Probablemente se puede acortar mediante la sustitución o sustitución de cadenas iterando sobre un índice único que codifica (a, b, c, d). Importar itertools resultó que no valía la pena.

Editar: se cambió a int(2*n**.5)+1desde 2*int(n**.5)+2para hacer un argumento más limpio, el mismo carácter cuenta.

CJam , 91 90 74 71 bytes

q~{W):W;:N4md!}gmqi257:B_**_{;)_[Bmd\Bmd]_N\{_*-}/mq_i@+\1%}g{2W#*}%`\;

Compacto, pero más lento que mi otro enfoque.

Pruébalo en línea! Pegue el Código , escriba el número entero deseado en Entrada y haga clic en Ejecutar .

Antecedentes

Esta publicación comenzó como una respuesta de 99 bytes de GolfScript . Si bien todavía había margen de mejora, GolfScript carece de la función sqrt incorporada. Mantuve la versión de GolfScript hasta la revisión 5 , ya que era muy similar a la versión de CJam.

Sin embargo, las optimizaciones desde la revisión 6 requieren operadores que no están disponibles en GolfScript, por lo que en lugar de publicar explicaciones separadas para ambos idiomas, decidí abandonar la versión menos competitiva (y mucho más lenta).

El algoritmo implementado calcula los números por fuerza bruta:

1. Para entrada m, busque Ny Wtal que m = 4**W * N.

2. Conjunto i = 257**2 * floor(sqrt(N/4)).

3. Conjunto i = i + 1.

4. Encuentra enteros j, k, ltales que i = 257**2 * j + 257 * k + l, donde k, l < 257.

5. Comprueba si d = N - j**2 - k**2 - l**2es un cuadrado perfecto.

6. Si no es así, y regrese al paso 3.

7. Imprimir 2**W * j, 2**W * k, 2**W * l, 2**W * sqrt(m).

Ejemplos

$ TIME='\n%e s' time cjam lagrange.cjam <<< 999999999
[27385 103 15813 14]
0.46 s
$ TIME='\n%e s' time cjam lagrange.cjam <<< 805306368
[16384 16384 0 16384]
0.23 s

Los tiempos corresponden a un Intel Core i7-4700MQ.

Cómo funciona

q~              " Read and interpret the input. ";
{
  W):W;         " Increment “W” (initially -1). ";
  :N            " Save the integer on the stack in “N”. ';
  4md!          " Push N / 4 and !(N % 4). ";
}g              " If N / 4 is an integer, repeat the loop.
mqi             " Compute floor(sqrt(N/4)). ";
257:B_**        " Compute i = floor(sqrt(N)) * 257**2. ";
_               " Duplicate “i” (dummy value). ";
{               " ";
  ;)_           " Pop and discard. Increment “i”. ";
  [Bmd\Bmd]     " Push [ j k l ], where i = 257**2 * j + 257 * k + l and k, l < 257. ";
  _N\           " Push “N” and swap it with a copy of [ j k l ]. ";
  {_*-}/        " Compute m = N - j**2 - k**2 - l**2. ";
  mq            " Compute sqrt(m). ";
  _i            " Duplicate sqrt(m) and compute floor(sqrt(m)). ";
  @+\           " Form [ j k l floor(sqrt(m)) ] and swap it with sqrt(m). ";
  1%            " Check if sqrt(m) is an integer. ";
}g              " If it is, we have a solution; break the loop. ";
{2W#*}%         " Push 2**W * [ j k l sqrt(m) ]. ";
`\;             " Convert the array into a string and discard “i”. ";

C, 228

Esto se basa en el algoritmo en la página de Wikipedia, que es una fuerza bruta O (n).

n,*l,x,y,i;main(){scanf("%d",&n);l=calloc(n+1,8);
for(x=0;2*x*x<=n;x++)for(y=x;(i=x*x+y*y)<=n;y++)l[2*i]=x,l[2*i+1]=y;
for(x=0,y=n;;x++,y--)if(!x|l[2*x+1]&&l[2*y+1]){
printf("%d %d %d %d\n",l[2*x],l[2*x+1],l[2*y],l[2*y+1]);break;}}

GolfScript, 133 130 129 bytes

~{.[4*{4/..4%1$!|!}do])\}:r~,(2\?:f;{{..*}:^~4-1??n*,}:v~)..
{;;(.^3$\-r;)8%!}do-1...{;;;)..252/@252%^@^@+4$\-v^@-}do 5$]{f*}%-4>`

Rápido, pero largo. La nueva línea se puede eliminar.

Pruébalo en línea.Tenga en cuenta que el intérprete en línea tiene un límite de tiempo de 5 segundos, por lo que podría no funcionar para todos los números.

Antecedentes

El algoritmo aprovecha el teorema de tres cuadrados de Legendre , que establece que cada número natural n que no tiene la forma

se puede expresar como la suma de tres cuadrados.

El algoritmo hace lo siguiente:

1. Exprese el número como 4**i * j.

2. Encontrar el mayor entero ktal que k**2 <= jy j - k**2satisface las hipótesis del Teorema de tres cuadrados de Legendre.

3. Conjunto i = 0.

4. Comprueba si j - k**2 - (i / 252)**2 - (i % 252)**2es un cuadrado perfecto.

5. Si no es así, incremente iy regrese al paso 4.

Ejemplos

$ TIME='%e s' time golfscript legendre.gs <<< 0
[0 0 0 0]
0.02 s
$ TIME='%e s' time golfscript legendre.gs <<< 123456789
[32 0 38 11111]
0.02 s
$ TIME='%e s' time golfscript legendre.gs <<< 999999999
[45 1 217 31622]
0.03 s
$ TIME='%e s' time golfscript legendre.gs <<< 805306368
[16384 0 16384 16384]
0.02 s

Los tiempos corresponden a un Intel Core i7-4700MQ.

Cómo funciona

~              # Interpret the input string. Result: “n”
{              #
  .            # Duplicate the topmost stack item.
  [            #
    4*         # Multiply it by four.
    {          #
      4/       # Divide by four.
      ..       # Duplicate twice.
      4%1$     # Compute the modulus and duplicate the number.
      !|!      # Push 1 if both are truthy.
    }do        # Repeat if the number is divisible by four and non-zero.
  ]            # Collect the pushed values (one per iteration) into an array.
  )\           # Pop the last element from the array and swap it with the array.
}:r~           # Save this code block as “r” and execute it.
,(2\?          # Get the length of the array, decrement it and exponentiate.
:f;            # Save the result in “f”.
               # The topmost item on the stack is now “j”, which is not divisible by 
               # four and satisfies that n = f**2 * j.
{              #
  {..*}:^~     # Save a code block to square a number in “^” and execute it.
  4-1??        # Raise the previous number to the power of 1/4.
               # The two previous lines compute (x**2)**(1/4), which is sqrt(abs(x)).
  n*,          # Repeat the string "\n" that many times and compute its length.
               # This casts to integer. (GolfScript doesn't officially support Rationals.)
}:v~           # Save the above code block in “v” and execute it.
)..            # Undo the first three instructions of the loop.
{              #
   ;;(         # Discard two items from the stack and decrement.
   .^3$\-      # Square and subtract from “n”.
   r;)8%!      # Check if the result satisfies the hypothesis of the three-square theorem.
}do            # If it doesn't, repeat the loop.
-1...          # Push 0 (“i”) and undo the first four instructions of the loop.
{              #
  ;;;)         # Discard two items from the stack and increment “i”.
  ..252/@252%  # Push the digits of “i” in base 252.
  ^@^@+4$\-    # Square both, add and subtract the result 
  v^@-         # Take square root, square and compare.
}do            # If the difference is a perfect square, break the loop.
5$]            # Duplicate the difference an collect the entire stack into an array.
{f*}%          # Multiply very element of the array by “f”.
-4>            # Reduce the array to its four last elements (the four numbers).
`              # Convert the result into a string.

Rev 1: C, 190

a,z,m;short s[15<<26];p(){m=s[a=z-a];printf("%d %f ",m,sqrt(a-m*m));}
main(){m=31727;for(a=m*m;--a;s[z<m*m?z:m*m]=a%m)z=a/m*(a/m)+a%m*(a%m);scanf("%d",&z);for(;a*!s[a]||!s[z-a];a++);p();p();}

Esto requiere aún más memoria que rev 0. El mismo principio: construir una tabla con un valor positivo para todas las sumas posibles de 2 cuadrados (y cero para aquellos números que no son sumas de dos cuadrados), luego búscalo.

En esta revisión, uso una matriz de en shortlugar de charalmacenar los golpes, para que pueda almacenar la raíz de uno de los cuadrados en la tabla en lugar de solo una bandera. Esto simplifica la función p(para decodificar la suma de 2 cuadrados) considerablemente ya que no hay necesidad de un bucle.

Windows tiene un límite de 2 GB en las matrices. Puedo sortear aquello con lo short s[15<<26]que hay una matriz de elementos 1006632960, suficiente para cumplir con la especificación. Por desgracia, el tamaño total de tiempo de ejecución del programa es de más de 2 GB y (a pesar de ajustar las configuraciones del sistema operativo) no he sido capaz de ejecutar más de este tamaño (aunque en teoría es posible.) Lo mejor que puedo hacer es short s[14<<26](939524096 elementos.) m*mDebe ser estrictamente menor que esto (30651 ^ 2 = 939483801.) Sin embargo, el programa se ejecuta perfectamente y debería funcionar en cualquier sistema operativo que no tenga esta restricción.

Código sin golf

a,z,m;
short s[15<<26];     
p(){m=s[a=z-a];printf("%d %f ",m,sqrt(a-m*m));}      
main(){       
 m=31727;             
 for(a=m*m;--a;s[z<m*m?z:m*m]=a%m)   //assignment to s[] moved inside for() is executed after the following statement. In this rev excessively large values are thrown away to s[m*m].
   z=a/m*(a/m)+a%m*(a%m);            //split a into high and low half, calculate h^2+l^2.                                  
 scanf("%d",&z); 
 for(;a*!s[a]||!s[z-a];a++);         //loop until s[a] and s[z-a] both contain an entry. s[0] requires special handling as s[0]==0, therefore a* is included to break out of the loop when a=0 and s[z] contains the sum of 2 squares.
 p();                                //print the squares for the first sum of 2 squares 
 p();}                               //print the squares for the 2nd sum of 2 squares (every time p() is called it does a=z-a so the two sums are exchanged.) 

Rev 0 C, 219

a,z,i,m;double t;char s[1<<30];p(){for(i=t=.1;(m=t)-t;i++)t=sqrt(a-i*i);printf("%d %f ",i-1,t);}
main(){m=1<<15;for(a=m*m;--a;){z=a/m*(a/m)+a%m*(a%m);s[z<m*m?z:0]=1;}scanf("%d",&z);for(;1-s[a]*s[z-a];a++);p();a=z-a;p();}

Esta es una bestia hambrienta de memoria. Se necesita una matriz de 1 GB, calcula todas las sumas posibles de 2 cuadrados y almacena una bandera para cada una en la matriz. Luego, para la entrada de usuario z, busca en la matriz dos sumas de 2 cuadrados ay za.

la función pluego reconstituye los cuadrados originales que se usaron para hacer las sumas de 2 cuadrados ay los z-aimprime, el primero de cada par como un entero, el segundo como un doble (si tiene que ser todos los enteros, se necesitan dos caracteres más, t> m=t.)

El programa tarda un par de minutos en construir la tabla de sumas de cuadrados (creo que esto se debe a problemas de administración de memoria, veo que la asignación de memoria aumenta lentamente en lugar de saltar como cabría esperar). Sin embargo, una vez hecho esto produce respuestas muy rápidamente (si se calculan varios números, el programa en scanfadelante se puede poner en un bucle.

código sin golf

a,z,i,m;
double t;
char s[1<<30];                              //handle numbers 0 up to 1073741823
p(){
 for(i=t=.1;(m=t)-t;i++)t=sqrt(a-i*i);      //where a contains the sum of 2 squares, search until the roots are found
 printf("%d %f ",i-1,t);}                   //and print them. m=t is used to evaluate the integer part of t. 

main(){       
 m=1<<15;                                   //max root we need is sqrt(1<<30);
 for(a=m*m;--a;)                            //loop m*m-1 down to 1, leave 0 in a
   {z=a/m*(a/m)+a%m*(a%m);s[z<m*m?z:0]=1;}  //split a into high and low half, calculate h^2+l^2. If under m*m, store flag, otherwise throw away flag to s[0]
 scanf("%d",&z);
 for(;1-s[a]*s[z-a];a++);                   //starting at a=0 (see above) loop until flags are found for sum of 2 squares of both (a) and (z-a)
 p();                                       //reconsitute and print the squares composing (a)
 a=z-a;                                     //assign (z-a) to a in order to...
 p();}                                      //reconsitute and print the squares composing (z-a)  

Salida de ejemplo

El primero es por la pregunta. El segundo fue elegido como uno difícil de buscar. En este caso, el programa tiene que buscar hasta 8192 ^ 2 + 8192 ^ 2 = 134217728, pero solo toma unos segundos una vez que se construye la tabla.

123456789
0 2.000000 3328 10601.000000

805306368
8192 8192.000000 8192 24576.000000

Mathematica, 138 caracteres

Resulta que esto produce resultados negativos e imaginarios para ciertas entradas como lo señala edc65 (por ejemplo, 805306368), por lo que esta no es una solución válida. Lo dejaré por ahora, y tal vez, si realmente odio mi tiempo, volveré e intentaré arreglarlo.

S[n_]:=Module[{a,b,c,d},G=Floor@Sqrt@#&;a=G@n;b:=G[n-a^2];c:=G[n-a^2-b^2];d:=G[n-a^2-b^2-c^2];While[Total[{a,b,c,d}^2]!=n,a-=1];{a,b,c,d}]

O sin desquitar:

S[n_] := Module[{a, b, c, d}, G = Floor@Sqrt@# &;
 a = G@n;
 b := G[n - a^2];
 c := G[n - a^2 - b^2];
 d := G[n - a^2 - b^2 - c^2];
 While[Total[{a, b, c, d}^2] != n, a -= 1];
 {a, b, c, d}
]

No analicé demasiado los algoritmos, pero espero que esta sea la misma idea. Se me ocurrió la solución obvia y la ajusté hasta que funcionó. Lo probé para todos los números entre 1 y mil millones y ... funciona. La prueba solo toma alrededor de 100 segundos en mi máquina.

Lo bueno de esto es que, dado que b, cyd están definidos con asignaciones retrasadas :=, no es necesario redefinirlos cuando se disminuye a. Esto salvó algunas líneas adicionales que tenía antes. Podría jugar más golf y anidar las partes redundantes, pero aquí está el primer borrador.

Ah, y lo ejecutas como S@123456789y puedes probarlo con {S@#, Total[(S@#)^2]} & @ 123456789o # == Total[(S@#)^2]&[123456789]. La prueba exhaustiva es

n=0;
AbsoluteTiming@ParallelDo[If[e != Total[(S@e)^2], n=e; Abort[]] &, {e, 1, 1000000000}]
n

Utilicé una declaración Print [] antes, pero eso la ralentizó mucho, aunque nunca se llama. Imagínate.

Haskell 123 + 3 = 126

main=getLine>>=print.f.read
f n=head[map(floor.sqrt)[a,b,c,d]|a<-r,b<-r,c<-r,d<-r,a+b+c+d==n]where r=[x^2|x<-[0..n],n>=x^2]

Fuerza bruta simple sobre cuadrados precalculados.

Necesita la -Oopción de compilación (agregué 3 caracteres para esto). Se tarda menos de 1 minuto en el peor de los casos 999950883.

Solo probado en GHC.

C: 198 caracteres

Probablemente pueda reducirlo a poco más de 100 caracteres. Lo que me gusta de esta solución es la cantidad mínima de basura, solo un bucle for, haciendo lo que debe hacer un bucle for (que es una locura).

i,a,b,c,d;main(n){for(scanf("%d",&n);a*a+b*b-n?a|!b?a*a>n|a<b?(--a,b=1):b?++b:++a:(a=b=0,--n,++i):c*c+d*d-i?c|!d?c*c>i|c<d?(--c,d=1):d?++d:++c:(a=b=c=d=0,--n,++i):0;);printf("%d %d %d %d",a,b,c,d);}

Y muy prettified:

#include <stdio.h>

int n, i, a, b, c, d;

int main() {
    for (
        scanf("%d", &n);
        a*a + b*b - n
            ? a | !b
                ? a*a > n | a < b
                    ? (--a, b = 1)
                    : b
                        ? ++b
                        : ++a
                : (a = b = 0, --n, ++i)
            : c*c + d*d - i
                ? c | !d
                    ? c*c > i | c < d
                        ? (--c, d = 1)
                        : d
                            ? ++d
                            : ++c
                    : (a = b = c = d = 0, --n, ++i)
                : 0;
    );
    printf("%d %d %d %d\n", a, b, c, d);
    return 0;
}

Editar: no es lo suficientemente rápido para todas las entradas, pero volveré con otra solución. Dejaré que este lío de operaciones ternarias permanezca hasta ahora.

Rev B: C, 179

a,b,c,d,m=1,n,q,r;main(){for(scanf("%d",&n);n%4<1;n/=4)m*=2;
for(a=sqrt(n),a-=(3+n-a*a)%4/2;r=n-a*a-b*b-c*c,d=sqrt(r),d*d-r;c=q%256)b=++q>>8;
printf("%d %d %d %d",a*m,b*m,c*m,d*m);}

Gracias a @Dennis por las mejoras. El resto de la respuesta a continuación no se actualiza desde la rev. A.

Rev A: C, 195

a,b,c,d,n,m,q;double r=.1;main(){scanf("%d",&n);for(m=1;!(n%4);n/=4)m*=2;a=sqrt(n);a-=(3+n-a*a)%4/2;
for(;(d=r)-r;q++){b=q>>8;c=q%256;r=sqrt(n-a*a-b*b-c*c);}printf("%d %d %d %d ",a*m,b*m,c*m,d*m);}

¡Mucho más rápido que mi otra respuesta y con mucho menos memoria!

Esto utiliza http://en.wikipedia.org/wiki/Legendre%27s_three-square_theorem . Cualquier número que no sea de la siguiente forma puede expresarse como la suma de 3 cuadrados (a esto le llamo la forma prohibida):

4^a*(8b+7), or equivalently 4^a*(8b-1)

Tenga en cuenta que todos los números cuadrados impares son de la forma (8b+1)y todos los números cuadrados pares son superficialmente de la forma 4b. Sin embargo, esto oculta el hecho de que todos los números cuadrados pares son de la forma 4^a*(odd square)==4^a*(8b+1). Como resultado 2^x-(any square number < 2^(x-1))para impar xsiempre será de la forma prohibida. Por lo tanto, estos números y sus múltiplos son casos difíciles, razón por la cual muchos de los programas aquí dividen las potencias de 4 como primer paso.

Como se indica en la respuesta de @ xnor, N-a*a no puede tener la forma prohibida durante 2 valores consecutivos de a. A continuación presento una forma simplificada de su tabla. Además del hecho de que después de la división por 4 N%4no puede ser igual a 0, tenga en cuenta que solo hay 2 valores posibles para (a*a)%4.

(a*a)%4= 01
        +--
       1|10
  N%4= 2|21  <- (N-a*a)%4
       3|32

Por lo tanto, queremos evitar que los valores (N-a*a)puedan ser de la forma prohibida, es decir, aquellos en los que (N-a*a)%4es 3 o 0. Como se puede ver, esto no puede suceder Ncon los pares e impares (a*a).

Entonces, mi algoritmo funciona así:

1. Divide out powers of 4
2. Set a=int(sqrt(N)), the largest possible square
3. If (N-a*a)%4= 0 or 3, decrement a (only once)
4. Search for b and c such that N-a*a-b*b-c*c is a perfect square

Particularmente me gusta la forma en que hago el paso 3. Agrego 3 a N , de modo que la disminución se requiere si (3+N-a*a)%4 =3 o 2. (pero no 1 o 0.) Divida esto entre 2 y todo el trabajo se puede hacer con una expresión bastante simple .

Código sin golf

Tenga en cuenta el forbucle único qy el uso de la división / módulo para derivar los valores de by cde él. Intenté usarlo acomo un divisor en lugar de 256 para guardar bytes, pero a veces el valor de ano era correcto y el programa se bloqueaba, posiblemente indefinidamente. 256 fue el mejor compromiso que puedo usar en >>8lugar de /256para la división.

a,b,c,d,n,m,q;double r=.1;
main(){
  scanf("%d",&n);
  for(m=1;!(n%4);n/=4)m*=2;
  a=sqrt(n);
  a-=(3+n-a*a)%4/2;
  for(;(d=r)-r;q++){b=q>>8;c=q%256;r=sqrt(n-a*a-b*b-c*c);}
  printf("%d %d %d %d ",a*m,b*m,c*m,d*m);}

Salida

Una peculiaridad interesante es que si ingresa un número cuadrado, N-(a*a)= 0. Pero el programa detecta que 0%4= 0 y disminuye al siguiente cuadrado hacia abajo. Como resultado, las entradas de números cuadrados siempre se descomponen en un grupo de cuadrados más pequeños a menos que sean de la forma 4^x.

999999999
31621 1 161 294

805306368
16384 0 16384 16384

999950883
31621 1 120 221

1
0 0 0 1

2
1 0 0 1

5
2 0 0 1

9
2 0 1 2

25
4 0 0 3

36
4 0 2 4

49
6 0 2 3

81
8 0 1 4

121
10 1 2 4

JavaScript - 175 191 176 173 caracteres

Fuerza bruta, pero rápida.

Edite rápido pero no lo suficiente para una entrada desagradable. Tuve que agregar un primer paso de reducción por multiplicaciones de 4.

Editar 2 Deshacerse de la función, salir dentro del bucle y luego forzar la contición de salida

Editar 3 0 no es una entrada válida

v=(p=prompt)();for(m=1;!(v%4);m+=m)v/=4;for(a=-~(q=Math.sqrt)(v);a--;)for(w=v-a*a,b=-~q(w);b--;)for(x=w-b*b,c=-~q(x);c--;)(d=q(x-c*c))==~~d&&p([m*a, m*b, m*c, m*d],a=b=c='')

Sin golf:

v = prompt();

for (m = 1; ! (v % 4); m += m) 
{
  v /= 4;
}
for (a = - ~Math.sqrt(v); a--;) /* ~ force to negative integer, changing sign lead to original value + 1 */
{
  for ( w = v - a*a, b = - ~Math.sqrt(w); b--;)
  {
    for ( x = w - b*b, c = - ~Math.sqrt(x); c--;)
    {
      (d = Math.sqrt(x-c*c)) == ~~d && prompt([m*a, m*b, m*c, m*d], a=b=c='') /* 0s a,b,c to exit loop */
    }
  }
}

Salida de ejemplo

123456789
11111,48,10,8

805306368
16384,16384,16384,0