O teorema dos quatro quadrados de Lagrange nos diz que qualquer número natural pode ser representado como a soma de quatro números quadrados. Sua tarefa é escrever um programa que faça isso.

Entrada: um número natural (abaixo de 1 bilhão)

Saída: quatro números cujos quadrados somam esse número (a ordem não importa)

Nota: Você não precisa fazer uma pesquisa de força bruta! Detalhes aqui e aqui . Se houver uma função que trivialize esse problema (determinarei), isso não será permitido. Funções primárias automatizadas e raiz quadrada são permitidas. Se houver mais de uma representação, qualquer uma está correta. Se você optar por fazer força bruta, ela deverá ser executada dentro de um prazo razoável (3 minutos)

entrada de amostra

123456789

saída de amostra (qualquer uma é boa)

10601 3328 2 0
10601 3328 2

CJam, 50 bytes

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

Minha terceira (e final, prometo) resposta. Essa abordagem baseia-se fortemente na resposta do primo .

Experimente on-line no intérprete CJam .

Uso

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

fundo

1. Depois de ver o algoritmo atualizado do primo, tive que ver como a implementação do CJam seria pontuada:

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

   Apenas 58 bytes! Esse algoritmo é executado em tempo quase constante e não exibe muita variação para diferentes valores de N. Vamos mudar isso ...

2. Em vez de começar floor(sqrt(N))e diminuir, podemos começar em 1e incrementar. Isso economiza 4 bytes.

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

3. Em vez de expressar Ncomo 4**a * b, podemos expressá-lo como p**(2a) * b- onde pé o menor fator primordial de N- para economizar mais 1 byte.

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

4. A modificação anterior nos permite mudar ligeiramente a implementação (sem tocar o próprio algoritmo): Em vez de dividir Npor p**(2a)e multiplicar a solução por p**a, podemos restringir directamente as possíveis soluções para os múltiplos de p**a. Isso economiza mais 2 bytes.

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

5. Não restringir o primeiro número inteiro a múltiplos de p**asalva um byte adicional.

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

Algoritmo final

1. Encontre ae btal que N = p**(2a) * b, onde bnão é um múltiplo de p**2e pseja o menor fator primo de N.

2. Set j = p**a.

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

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

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

6. Se N - j**2 - k**2 - l**2 - m**2 > 0, defina j = j + 1e volte para a etapa 3.

Isso pode ser implementado da seguinte maneira:

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.           ";

Benchmarks

Executei todas as 5 versões em todos os números inteiros positivos até 999.999.999 no meu Intel Core i7-3770, medi o tempo de execução e contei as iterações necessárias para encontrar uma solução.

A tabela a seguir mostra o número médio de iterações e o tempo de execução para um único número inteiro:

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. Com apenas 4 iterações e 6,6 microssegundos por número inteiro, o algoritmo do primo é incrivelmente rápido.

2. Começando em floor(sqrt(N))faz mais sentido, pois isso nos deixa com valores menores para a soma dos três quadrados restantes. Como esperado, a partir de 1 é muito mais lento.

3. Este é um exemplo clássico de uma boa ideia mal implementada. Para realmente reduzir o tamanho do código, contamos com mFo fator de número inteiro N. Embora a versão 3 exija menos iterações que a versão 2, é muito mais lenta na prática.

4. Embora o algoritmo não mude, a versão 4 é muito mais lenta. Isso ocorre porque ele executa uma divisão de ponto flutuante adicional e uma multiplicação de números inteiros em cada iteração.

5. Para a entrada N = p**(2a) ** b, o algoritmo 5 exigirá (k - 1) * p**a + 1iterações, onde ké o número de iterações 4 que o algoritmo requer. Se k = 1ou a = 0, isso não faz diferença.

   No entanto, qualquer entrada do formulário 4**a * (4**c * (8 * d + 7) + 1)pode ter um desempenho muito ruim. Para o valor inicial j = p**a, N - 4**a = 4**(a + c) * (8 * d + 7), por isso não pode ser expresso como uma soma de três quadrados. Assim, k > 1e pelo menos p**aiterações são necessárias.

   Felizmente, o algoritmo original do primo é incrivelmente rápido e N < 1,000,000,000. O pior caso que pude encontrar à mão é o 265,289,728 = 4**10 * (4**1 * (7 * 8 + 7) + 1)que requer 6.145 iterações. O tempo de execução é inferior a 300 ms na minha máquina. Em média, esta versão é 13,5 vezes mais lenta que a implementação do algoritmo do primo.

FRACTRAN: 156 98 frações

Como esse é um problema clássico da teoria dos números, que melhor maneira de resolver isso do 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 na entrada o formato 2 ⁿ × 193 e produz 3 ^a × 5 ^b × 7 ^c × 11 ^d . Pode ser executado em 3 minutos, se você tem um realmente bom intérprete. Talvez.

... ok, não mesmo. Este parecia ser um problema tão divertido de se fazer no FRACTRAN que eu tive que tentar. Obviamente, essa não é uma solução adequada para a questão, pois não exige tempo (força bruta) e nem sequer joga golfe, mas pensei em publicá-la aqui, porque não é todo dia que uma pergunta do Codegolf pode ser feito no FRACTRAN;)

Sugestão

O código é equivalente ao seguinte 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

Três métodos são mostrados. Somente a primeira abordagem atende ao requisito de tempo.

1-FindInstance (51 caracteres)

Isso retorna uma única solução da equação.

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

Exemplos e horários

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}}}

2-IntegerPartitions

Isso funciona também, mas é muito lento para atender aos requisitos de velocidade.

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

Range[0, Floor@Sqrt@n]^2é o conjunto de todos os quadrados menores que a raiz quadrada de n(o maior quadrado possível na partição).

{4}requer que as partições inteiras nconsistam em 4 elementos do conjunto de quadrados mencionado acima.

1, dentro da função IntegerPartitionsretorna a primeira solução.

[[1]]remove os aparelhos externos; a solução foi retornada como um conjunto de um elemento.

f[123456]

{348, 44, 20, 4}

3-PowerRepresentations

PowerRepresentations retorna todas as soluções para o problema dos 4 quadrados. Também pode resolver somas de outros poderes.

PowersRepresentations retorna, em menos de 5 segundos, as 181 maneiras de expressar 123456789 como a soma de 4 quadrados:

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

No entanto, é muito lento para outras somas.

Perl - 116 bytes 87 bytes (veja a atualização abaixo)

#!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 o shebang como um byte, novas linhas foram adicionadas à sanidade horizontal.

Algo como uma combinação de envio de código mais rápido para código-golfe .

A complexidade média (pior?) Do caso parece ser O (log n) O (n ^0,07 ) . Nada que encontrei é mais lento que 0,001s e verifiquei todo o intervalo de 900000000 a 999999999 . Se você encontrar algo que demore significativamente mais do que isso, ~ 0,1s ou mais, entre em contato.

Uso da amostra

$ 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

Os dois últimos parecem ser os piores cenários para outros envios. Nos dois casos, a solução mostrada é literalmente a primeira coisa verificada. Pois 123456789é o segundo.

Se você deseja testar um intervalo de valores, pode usar o seguinte 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);

Melhor quando canalizado para um arquivo. O intervalo 1..1000000leva cerca de 14s no meu computador (71000 valores por segundo) e o intervalo 999000000..1000000000leva cerca de 20s (50000 valores por segundo), consistente com a complexidade média de O (log n) .

Atualizar

Edit : Acontece que este algoritmo é muito semelhante ao que tem sido usado por calculadoras mentais há pelo menos um século .

Desde a publicação original, verifiquei todos os valores no intervalo de 1..1000000000 . O comportamento do "pior caso" foi exibido pelo valor 699731569 , que testou um total geral de 190 combinações antes de chegar a uma solução. Se você considera 190 uma pequena constante - e eu certamente considero - o pior comportamento possível no intervalo necessário pode ser considerado O (1) . Ou seja, tão rápido quanto procurar a solução em uma mesa gigante e, em média, possivelmente mais rápido.

Outra coisa, porém. Após 190 iterações, qualquer coisa maior que 144400 nem chegou além da primeira passagem. A lógica do percurso pela primeira vez é inútil - nem sequer é usada. O código acima pode ser bastante reduzido:

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

Que realiza apenas a primeira passagem da pesquisa. Precisamos confirmar que não há valores abaixo de 144400 que precisem da segunda passagem, no entanto:

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");
}

Resumindo, para o intervalo de 1 a 1000000000 , existe uma solução de tempo quase constante e você está olhando para ela.

Atualização atualizada

@Dennis e eu fizemos várias melhorias neste algoritmo. Você pode acompanhar o progresso nos comentários abaixo e a discussão subsequente, se isso lhe interessar. O número médio de iterações para o intervalo necessário caiu de pouco mais de 4 para 1,229 , e o tempo necessário para testar todos os valores para 1..1000000000 foi aprimorado de 18m 54s para 2m 41s. O pior caso anteriormente exigia 190 iterações; o pior caso agora, 854382778 , precisa apenas de 21 .

O código final do Python é o seguinte:

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

Isso usa duas tabelas de correção pré-calculadas, uma com 10kb de tamanho e a outra com 253kb. O código acima inclui as funções de gerador para essas tabelas, embora provavelmente devam ser calculadas em tempo de compilação.

Uma versão com tabelas de correção de tamanho mais modesto pode ser encontrada aqui: http://codepad.org/1ebJC2OV Esta versão requer uma média de 1.620 iterações por termo, com o pior caso de 38 , e todo o intervalo é executado em cerca de 3m 21s. Um pouco de tempo é compensado, usando bit anda bit para correção b , em vez de módulo.

Melhorias

Valores pares têm mais probabilidade de produzir uma solução do que valores ímpares.
O artigo de cálculo mental, vinculado anteriormente, observa que se, após remover todos os fatores de quatro, o valor a ser decomposto for par, esse valor poderá ser dividido por dois e a solução reconstruída:

Embora isso possa fazer sentido para o cálculo mental (valores menores tendem a ser mais fáceis de calcular), não faz muito sentido algoritmicamente. Se você pegar 256 pares aleatórios de 4 e examinar a soma dos quadrados do módulo 8 , verá que os valores 1 , 3 , 5 e 7 são alcançados, em média, 32 vezes. No entanto, os valores 2 e 6 são atingidos 48 vezes. A multiplicação dos valores ímpares por 2 encontrará uma solução, em média, em 33% menos iterações. A reconstrução é a seguinte:

Cuidado deve ser tomado para que um e b têm a mesma paridade, bem como c e d , mas se uma solução foi encontrada em tudo, uma ordenação adequada é garantido que existe.

Caminhos impossíveis não precisam ser verificados.
Depois de selecionar o segundo valor, b , já pode ser impossível a existência de uma solução, dados os possíveis resíduos quadráticos para qualquer módulo. Em vez de verificar de qualquer maneira, ou passar para a próxima iteração, o valor de b pode ser 'corrigido', diminuindo-o pela menor quantidade possível de uma solução. As duas tabelas de correção armazenam esses valores, um para be o outro para c . Usar um módulo mais alto (mais precisamente, usar um módulo com relativamente menos resíduos quadráticos) resultará em uma melhoria melhor. O valor a não precisa de correção; modificando n para ser par, todos os valores dea são válidos.

Python 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])

Depois de reduzirmos a entrada Npara não ser divisível por 4, ela deve ser expressável como uma soma de quatro quadrados, onde um deles é o maior valor possível a=int(N**0.5)ou um menor que isso, deixando apenas um pequeno restante para a soma dos outros três quadrados para cuidar. Isso reduz muito o espaço de pesquisa.

Aqui está uma prova mais tarde, esse código sempre encontra uma solução. Queremos encontrar um apara que n-a^2seja a soma de três quadrados. No Teorema dos Três Quadrados de Legendre , um número é a soma de três quadrados, a menos que seja a forma 4^j(8*k+7). Em particular, esses números são 0 ou 3 (módulo 4).

Mostramos que não dois consecutivos apodem fazer com que a quantidade restante N-a^2tenha tal forma para os dois valores consecutivos. Podemos fazer isso simplesmente criando uma tabela ae o Nmódulo 4, observando isso N%4!=0porque extraímos todos os poderes de 4 de N.

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

Porque não há dois consecutivos adar (N-a*a)%4 in [0,3], um deles é seguro para uso. Portanto, avidamente usamos o maior possível ncom n^2<=N, e n-1. Desde então N<(n+1)^2, o restante N-a^2a ser representado como uma soma de três quadrados é no máximo (n+1)^2 -(n-1)^2, o que é igual 4*n. Portanto, basta verificar apenas valores até 2*sqrt(n), que é exatamente o intervalo R.

Pode-se otimizar ainda mais o tempo de execução, parando após uma única solução, computando em vez de repetir o último valor de pesquisando apenas entre valores b<=c<=d. Mas, mesmo sem essas otimizações, a pior instância que encontrei concluída em 45 segundos na minha máquina.

A cadeia de "para x em R" é lamentável. Provavelmente, pode ser reduzido pela substituição ou substituição de cadeias, iterando sobre um único índice que codifica (a, b, c, d). Importar ferramentas de controle acabou não valendo a pena.

Editar: Alterado int(2*n**.5)+1de 2*int(n**.5)+2para tornar o argumento mais limpo, com a mesma contagem de caracteres.

CJam , 91 90 74 71 bytes

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

Compacto, mas mais lento que minha outra abordagem.

Experimente online! Cole o código , digite o número inteiro desejado em Entrada e clique em Executar .

fundo

Esta postagem começou como uma resposta de 99 bytes do GolfScript . Embora ainda houvesse espaço para melhorias, o GolfScript não possui a função sqrt integrada. Eu mantive a versão GolfScript até a revisão 5 , pois era muito semelhante à versão CJam.

No entanto, as otimizações desde a revisão 6 exigem operadores que não estão disponíveis no GolfScript. Portanto, em vez de postar explicações separadas para os dois idiomas, decidi abandonar a versão menos competitiva (e muito mais lenta).

O algoritmo implementado calcula os números por força bruta:

1. Para entrada m, encontre Ne Wassim por diante m = 4**W * N.

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

3. Set i = i + 1.

4. Encontre números inteiros j, k, lcomo i = 257**2 * j + 257 * k + l, onde k, l < 257.

5. Verifique se d = N - j**2 - k**2 - l**2é um quadrado perfeito.

6. Caso contrário, volte ao passo 3.

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

Exemplos

$ 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

Os horários correspondem a um Intel Core i7-4700MQ.

Como 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

Isso é baseado no algoritmo da página da Wikipedia, que é uma força 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, mas demorado. A nova linha pode ser removida.

Experimente online. Observe que o intérprete on-line tem um limite de 5 segundos, portanto, pode não funcionar para todos os números.

fundo

O algoritmo tira proveito do teorema de três quadrados de Legendre , que afirma que todo número natural n que não é da forma

pode ser expresso como a soma de três quadrados.

O algoritmo faz o seguinte:

1. Expresse o número como 4**i * j.

2. Encontre o maior inteiro ktal que k**2 <= je j - k**2satisfaz a hipótese do teorema de três-quadrado de Legendre.

3. Set i = 0.

4. Verifique se j - k**2 - (i / 252)**2 - (i % 252)**2é um quadrado perfeito.

5. Caso contrário, aumente ie volte para a etapa 4.

Exemplos

$ 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

Os horários correspondem a um Intel Core i7-4700MQ.

Como 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();}

Isso tem ainda mais memória do que a rev 0. O mesmo princípio: construa uma tabela com um valor positivo para todas as somas possíveis de 2 quadrados (e zero para os números que não são somas de dois quadrados), depois procure-a.

Nesta revisão, use uma matriz de em shortvez de charpara armazenar os hits, para que eu possa armazenar a raiz de um dos quadrados da tabela em vez de apenas uma bandeira. Isso simplifica a função p(para decodificar a soma de 2 quadrados) consideravelmente, pois não há necessidade de um loop.

O Windows tem um limite de 2 GB em matrizes. Eu posso contornar o short s[15<<26]que é uma matriz de 1006632960 elementos, o suficiente para cumprir as especificações. Infelizmente, o tamanho total do tempo de execução do programa ainda é superior a 2 GB e (apesar de ajustar as configurações do SO), não fui capaz de executar esse tamanho (embora seja teoricamente possível). O melhor que posso fazer é short s[14<<26](939524096 elementos.) m*mDeve ser estritamente menor que isso (30651 ^ 2 = 939483801.) No entanto, o programa funciona perfeitamente e deve funcionar em qualquer sistema operacional que não possua essa restrição.

Código ungolfed

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();}

Este é um animal com fome de memória. É necessário um array de 1 GB, calcula todas as somas possíveis de 2 quadrados e armazena um sinalizador para cada um no array. Em seguida, para a entrada do usuário z, ele procura na matriz duas somas de 2 quadrados ae za.

a função pentão reconsitutes os quadrados originais que foram usadas para fazer as somas de 2 quadrados ae z-ae imprime-los, o primeiro de cada par como um número inteiro, o segundo como um duplo (se ele tem de ser todos os inteiros são necessários mais de dois caracteres, t> m=t.)

O programa leva alguns minutos para criar a tabela de somas de quadrados (acho que isso se deve a problemas de gerenciamento de memória, vejo a alocação de memória subindo lentamente, em vez de saltar como seria de esperar.) No entanto, uma vez feito isso, produz respostas muito rapidamente (se vários números devem ser calculados, o programa a partir de então scanfpode ser colocado em um loop.

código não destruído

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)  

Saída de exemplo

O primeiro é por pergunta. A segunda foi escolhida como difícil de pesquisar. Nesse caso, o programa deve procurar até 8192 ^ 2 + 8192 ^ 2 = 134217728, mas leva apenas alguns segundos depois que a tabela é criada.

123456789
0 2.000000 3328 10601.000000

805306368
8192 8192.000000 8192 24576.000000

Mathematica, 138 caracteres

Portanto, verifica-se que isso produz resultados negativos e imaginários para determinadas entradas, conforme apontado por edc65 (por exemplo, 805306368), portanto, essa não é uma solução válida. Vou deixar por enquanto, e talvez, se realmente odeio meu tempo, voltarei e tentarei consertar.

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}]

Ou não qualificado:

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}
]

Não olhei muito para os algoritmos, mas espero que seja a mesma ideia. Eu apenas vim com a solução óbvia e a aprimorei até que funcionasse. Eu testei para todos os números entre 1 e 1 bilhão e ... funciona. O teste leva apenas cerca de 100 segundos na minha máquina.

O bom disso é que, como b, c e d são definidos com atribuições atrasadas :=, eles não precisam ser redefinidos quando a é decrementado. Isso salvou algumas linhas extras que eu tinha antes. Eu poderia jogar mais e aninhar as partes redundantes, mas aqui está o primeiro rascunho.

Ah, e você o executa como S@123456789e pode testá-lo com {S@#, Total[(S@#)^2]} & @ 123456789ou # == Total[(S@#)^2]&[123456789]. O teste exaustivo é

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

Eu usei uma declaração Print [] antes, mas isso diminuiu bastante a velocidade, mesmo que nunca seja chamada. Vai saber.

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]

Força bruta simples sobre quadrados pré-calculados.

Ele precisa da -Oopção de compilação (adicionei 3 caracteres para isso). Demora menos de 1 minuto para o pior caso 999950883.

Testado apenas no GHC.

C: 198 caracteres

Provavelmente, posso apertar para pouco mais de 100 caracteres. O que eu mais gosto nessa solução é a quantidade mínima de lixo, apenas um loop for simples, fazendo o que um loop for deve fazer (o que é muito louco).

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);}

E fortemente prettificado:

#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;
}

Edit: Não é rápido o suficiente para todas as entradas, mas voltarei com outra solução. Vou deixar essa bagunça da operação ternária ficar a partir de agora.

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);}

Obrigado a @Dennis pelas melhorias. O restante da resposta abaixo não é atualizado da 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);}

Muito mais rápido que minha outra resposta e com muito menos memória!

Isso usa http://en.wikipedia.org/wiki/Legendre%27s_three-square_theorem . Qualquer número que não seja da seguinte forma pode ser expresso como a soma de 3 quadrados (eu chamo isso de forma proibida):

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

Observe que todos os números quadrados ímpares são da forma (8b+1)e todos os números quadrados pares são superficialmente da forma 4b. No entanto, isso oculta o fato de que todos os números pares são da forma 4^a*(odd square)==4^a*(8b+1). Como resultado, o 2^x-(any square number < 2^(x-1))ímpar xserá sempre da forma proibida. Portanto, esses números e seus múltiplos são casos difíceis, e é por isso que muitos dos programas aqui dividem as potências de 4 como primeiro passo.

Conforme declarado na resposta de @ xnor, N-a*anão pode ter a forma proibida por 2 valores consecutivos de a. Abaixo, apresento uma forma simplificada de sua tabela. Além do fato de que após a divisão por 4 N%4não pode ser igual a 0, observe que existem apenas 2 valores possíveis para (a*a)%4.

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

Portanto, queremos evitar valores (N-a*a)que podem ser da forma proibida, ou seja, aqueles onde (N-a*a)%4é 3 ou 0. Como pode ser visto, isso não pode ocorrer para o mesmo Ncom ímpares e pares (a*a).

Então, meu algoritmo funciona assim:

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

Eu particularmente gosto da maneira como passo o passo 3. Acrescento 3 a N, de modo que o decremento seja necessário se (3+N-a*a)%4 =3 ou 2. (mas não 1 ou 0.) Divida isso por 2 e todo o trabalho pode ser feito com uma expressão bastante simples .

Código ungolfed

Observe o forloop único qe use divisão / módulo para derivar os valores de be cpara ele. Tentei usar acomo divisor em vez de 256 para salvar bytes, mas às vezes o valor de anão estava certo e o programa travou, possivelmente por tempo indeterminado. 256 foi o melhor compromisso que posso usar em >>8vez da /256divisão.

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);}

Resultado

Uma peculiaridade interessante é que, se você digitar um número quadrado, N-(a*a)= 0. Mas o programa detecta isso 0%4= 0 e diminui para o próximo quadrado abaixo. Como resultado, as entradas de números quadrados são sempre decompostas em um grupo de quadrados menores, a menos que sejam da 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

Força bruta, mas rápida.

Edite rápido, mas não o suficiente para algumas entradas desagradáveis. Eu tive que adicionar uma primeira etapa de redução pela multiplicação de 4.

Edit 2 Livre-se da função, produza dentro do loop e force a saída da contição

Editar 3 0 não é uma 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='')

Ungolfed:

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 */
    }
  }
}

Saída de exemplo

123456789
11111,48,10,8

805306368
16384,16384,16384,0