Twierdzenie Lagrange'a o czterech kwadratach mówi nam, że dowolna liczba naturalna może być reprezentowana jako suma czterech liczb kwadratowych. Twoim zadaniem jest napisanie programu, który to robi.

Dane wejściowe: liczba naturalna (poniżej 1 miliarda)

Wynik: cztery liczby, których kwadraty sumują się do tej liczby (kolejność nie ma znaczenia)

Uwaga: nie musisz szukać brutalnej siły! Szczegóły tutaj i tutaj . Jeśli istnieje funkcja, która trywializuje ten problem (ustalę), nie jest dozwolona. Automatyczne funkcje pierwsze i pierwiastek kwadratowy są dozwolone. Jeśli jest więcej niż jedna reprezentacja, wszystko jest w porządku. Jeśli zdecydujesz się na brutalną siłę, musi ona działać w rozsądnym czasie (3 minuty)

przykładowe dane wejściowe

123456789

próbka wyjściowa (każda z nich jest w porządku)

10601 3328 2 0
10601 3328 2

CJam, 50 bajtów

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

Moja trzecia (i ostatnia, obiecuję) odpowiedź. To podejście opiera się w dużej mierze na odpowiedzi primo .

Wypróbuj online w interpretatorze CJam .

Stosowanie

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

tło

1. Po obejrzeniu zaktualizowanego algorytmu primo musiałem zobaczyć, jak uzyskałaby implementacja CJam:

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

   Tylko 58 bajtów! Ten algorytm działa w prawie stałym czasie i nie wykazuje dużej zmienności dla różnych wartości N. Zmieńmy to ...

2. Zamiast zaczynać floor(sqrt(N))i zmniejszać, możemy zacząć od 1przyrostu. To oszczędza 4 bajty.

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

3. Zamiast wyrażać Njako 4**a * b, możemy wyrazić to jako p**(2a) * b- gdzie pjest najmniejszy czynnik główny N- aby zaoszczędzić 1 bajt.

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

4. Poprzedniego modyfikacja pozwala nam nieznacznie zmienić realizację (bez dotykania samego algorytmu): Zamiast dzielenia Nprzez p**(2a)i pomnożyć przez rozwiązanie p**a, możemy bezpośrednio ograniczają możliwych rozwiązań do wielokrotności p**a. To oszczędza 2 kolejne bajty.

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

5. Nie ograniczanie pierwszej liczby całkowitej do wielokrotności p**azapisuje dodatkowy bajt.

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

Ostateczny algorytm

1. Znajdź ai btakie, że N = p**(2a) * bgdzie bnie jest wielokrotnością p**2i pjest najmniejszym czynnikiem podstawowym 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. Jeśli N - j**2 - k**2 - l**2 - m**2 > 0, ustaw j = j + 1i wróć do kroku 3.

Można to zaimplementować w następujący sposób:

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

Benchmarki

Uruchomiłem wszystkie 5 wersji dla wszystkich liczb całkowitych dodatnich do 999,999,999 na moim rdzeniu Intel Core i7-3770, zmierzyłem czas wykonania i policzyłem iteracje wymagane do znalezienia rozwiązania.

Poniższa tabela pokazuje średnią liczbę iteracji i czas wykonania dla pojedynczej liczby całkowitej:

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. Przy zaledwie 4 iteracjach i 6,6 mikro sekundy na liczbę całkowitą algorytm primo jest niewiarygodnie szybki.

2. Rozpoczęcie od floor(sqrt(N))ma większy sens, ponieważ pozostawia nam to mniejsze wartości dla sumy pozostałych trzech kwadratów. Zgodnie z oczekiwaniami, rozpoczęcie od 1 jest znacznie wolniejsze.

3. To klasyczny przykład źle zaimplementowanego dobrego pomysłu. Aby faktycznie zmniejszyć rozmiar kodu, polegamy na tym mF, co rozkłada liczbę całkowitą N. Chociaż wersja 3 wymaga mniej iteracji niż wersja 2, w praktyce jest znacznie wolniejsza.

4. Chociaż algorytm się nie zmienia, wersja 4 jest znacznie wolniejsza. Wynika to z tego, że wykonuje ona dodatkowy podział zmiennoprzecinkowy i mnożenie liczb całkowitych w każdej iteracji.

5. Dla danych wejściowych N = p**(2a) ** balgorytm 5 będzie wymagał (k - 1) * p**a + 1iteracji, gdzie kjest liczba wymaganych algorytmów 4. Jeśli k = 1lub a = 0, to nie ma znaczenia.

   Jednak każde wprowadzanie formularza 4**a * (4**c * (8 * d + 7) + 1)może działać dość źle. W odniesieniu do wartości wyjściowej j = p**a, N - 4**a = 4**(a + c) * (8 * d + 7)tak że nie może zostać wyrażona jako suma trzech pól. Zatem wymagane są k > 1przynajmniej p**aiteracje.

   Na szczęście oryginalny algorytm primo jest niesamowicie szybki i N < 1,000,000,000. Najgorszy przypadek, jaki udało mi się znaleźć ręcznie, to 265,289,728 = 4**10 * (4**1 * (7 * 8 + 7) + 1)6,145 iteracji. Czas realizacji jest mniejszy niż 300 ms na moim komputerze. Ta wersja jest średnio 13,5 razy wolniejsza niż implementacja algorytmu primo.

FRACTRAN: 156 98 frakcji

Ponieważ jest to klasyczny problem teorii liczb, czy jest lepszy sposób na rozwiązanie tego niż używanie liczb!

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

Pobiera dane wejściowe formularza 2 ⁿ × 193 i dane wyjściowe 3 ^a × 5 ^b × 7 ^c × 11 ^d . Może biec za 3 minuty, jeśli masz naprawdę dobrego tłumacza. Może.

... ok, niezupełnie. Wydawało się, że jest to tak fajny problem w FRACTRAN, że musiałem go wypróbować. Oczywiście nie jest to właściwe rozwiązanie tego pytania, ponieważ nie wymaga czasu (to brutalne siły) i prawie nie gra w golfa, ale pomyślałem, że opublikuję to tutaj, ponieważ nie codziennie pytanie Codegolfa można to zrobić we FRACTRAN;)

Wskazówka

Kod jest równoważny następującemu pseudo-Pythonowi:

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

Pokazane są trzy metody. Tylko pierwsze podejście spełnia wymagania czasowe.

1-FindInstance (51 znaków)

Zwraca to jedno równanie równania.

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

Przykłady i terminy

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

Działa to również, ale jest zbyt wolne, aby spełnić wymagania dotyczące prędkości.

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

Range[0, Floor@Sqrt@n]^2to zbiór wszystkich kwadratów mniejszy niż pierwiastek kwadratowy z n(największy możliwy kwadrat w partycji).

{4}wymaga partycji całkowitych nskładających się z 4 elementów z wyżej wspomnianego zestawu kwadratów.

1, w ramach funkcji IntegerPartitionszwraca pierwsze rozwiązanie.

[[1]]usuwa zewnętrzne szelki; rozwiązanie zostało zwrócone jako zestaw jednego elementu.

f[123456]

{348, 44, 20, 4}

3-PowerRepresentations

PowerRepresentations zwraca wszystkie rozwiązania problemu 4 kwadratów. Może również rozwiązać sumy innych mocy.

PowersRepresentations zwraca, w mniej niż 5 sekund, 181 sposobów wyrażenia 123456789 jako sumę 4 kwadratów:

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

Jest jednak o wiele za wolny na inne kwoty.

Perl - 116 bajtów 87 bajtów (patrz aktualizacja poniżej)

#!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"

Licząc shebang jako jeden bajt, dodano nowe znaki dla poziomego zdrowia psychicznego.

Coś z kombinacji kodu-kodu najszybszego kodu .

Średnia (najgorsza?) Złożoność przypadku wydaje się być O (log n) O (n ^0,07 ) . Nic, co znalazłem, działa wolniej niż 0,001s i sprawdziłem cały zakres od 900000000 do 999999999 . Jeśli znajdziesz coś, co trwa znacznie dłużej, około 0,1 s lub więcej, daj mi znać.

Przykładowe użycie

$ 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

Ostatnie dwa z nich wydają się najgorszym scenariuszem dla innych wniosków. W obu przypadkach pokazane rozwiązanie jest dosłownie pierwszą sprawdzoną rzeczą. Bo 123456789to drugi.

Jeśli chcesz przetestować zakres wartości, możesz użyć następującego skryptu:

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

Najlepiej gdy jest przesyłany do pliku. Zakres 1..1000000zajmuje około 14 sekund na moim komputerze (71000 wartości na sekundę), a zakres 999000000..1000000000zajmuje około 20 sekund (50000 wartości na sekundę), zgodnie ze średnią złożonością O (log n) .

Aktualizacja

Edycja : Okazuje się, że ten algorytm jest bardzo podobny do tego, który jest używany przez kalkulatory umysłowe od co najmniej stulecia .

Od czasu pierwotnego opublikowania sprawdziłem każdą wartość z zakresu od 1..1000000000 . Zachowanie „najgorszego przypadku” wykazało wartość 699731569 , która przetestowała w sumie 190 kombinacji przed znalezieniem rozwiązania. Jeśli uważasz 190 za małą stałą - a ja na pewno tak - najgorsze zachowanie w wymaganym zakresie można uznać za O (1) . Oznacza to, że tak szybko, jak wyszukiwanie rozwiązania z gigantycznego stołu, i średnio całkiem możliwe, że szybciej.

Ale inna sprawa. Po 190 iteracjach cokolwiek większego niż 144400 nawet nie przekroczyło pierwszego przejścia. Logika pierwszego przejścia jest bezwartościowa - nawet nie jest używana. Powyższy kod można nieco skrócić:

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

Który wykonuje tylko pierwszy przebieg wyszukiwania. Musimy jednak potwierdzić, że nie ma żadnych wartości poniżej 144400, które wymagałyby drugiego przejścia:

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

Krótko mówiąc, dla zakresu 1..1000000000 istnieje rozwiązanie o prawie stałym czasie i na to patrzysz.

Zaktualizowana aktualizacja

@Dennis i ja wprowadziliśmy kilka ulepszeń tego algorytmu. Możesz śledzić postępy w komentarzach poniżej i późniejszą dyskusję, jeśli Cię to interesuje. Średnia liczba iteracji dla wymaganego zakresu spadła z nieco ponad 4 do 1,229 , a czas potrzebny do przetestowania wszystkich wartości dla 1..1000000000 został skrócony z 18m 54s do 2m 41s. Najgorszy przypadek wymagał wcześniej 190 iteracji; najgorszy teraz, 854382778 , potrzebuje tylko 21 .

Ostateczny kod Pythona jest następujący:

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

Wykorzystuje dwie wstępnie obliczone tabele korekcji, jedna o wielkości 10 kb, druga 253 kb. Powyższy kod zawiera funkcje generatora dla tych tabel, chociaż prawdopodobnie powinny one zostać obliczone w czasie kompilacji.

Wersję ze skromniejszymi tabelami korekcji można znaleźć tutaj: http://codepad.org/1ebJC2OV Ta wersja wymaga średnio 1.620 iteracji na semestr, w najgorszym przypadku 38 , a cały zakres trwa około 3m 21s. Trochę czasu składa się na, za pomocą bitowego anddla b korekty, zamiast modulo.

Ulepszenia

Nawet wartości parzyste częściej dają rozwiązanie niż nieparzyste.
Artykuł dotyczący obliczeń mentalnych powiązany z wcześniejszymi uwagami stwierdza, że ​​jeśli po usunięciu wszystkich czterech czynników wartość do rozkładu jest parzysta, wartość tę można podzielić na dwa, a rozwiązanie odtworzyć:

Chociaż może to mieć sens w obliczeniach umysłowych (mniejsze wartości wydają się być łatwiejsze do obliczenia), nie ma to większego sensu algorytmicznego. Jeśli weźmiesz 256 losowych 4- kropek i zbadasz sumę kwadratów modulo 8 , przekonasz się, że wartości 1 , 3 , 5 i 7 są osiągane średnio 32 razy. Jednak wartości 2 i 6 są osiągane 48 razy. Pomnożenie nieparzystych wartości przez 2 znajdzie rozwiązanie średnio w 33% mniej iteracji. Rekonstrukcja jest następująca:

Należy zwrócić szczególną uwagę na to, i b mają taką samą parzystości oraz C i D , ale jeśli roztwór został znaleziony w wszystkim właściwej kolejności gwarantuje istnieć.

Niemożliwe ścieżki nie muszą być sprawdzane.
Po wybraniu drugiej wartości, b , rozwiązanie może już być niemożliwe, biorąc pod uwagę możliwe reszty kwadratowe dla dowolnego danego modułu. Zamiast sprawdzania i tak dalej lub przechodzenia do następnej iteracji, wartość b można „skorygować”, zmniejszając ją o najmniejszą wartość, która może doprowadzić do rozwiązania. Dwie tabele korekcji przechowują te wartości, jedną dla b , a drugą dla c . Użycie wyższego modułu (dokładniej, użycie modułu z relatywnie mniejszą liczbą reszt kwadratowych) spowoduje lepszą poprawę. Wartość a nie wymaga korekty; przez zmianę n na parzystą, wszystkie wartościa są ważne.

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

Po zmniejszeniu nakładu, Naby nie był podzielny przez 4, musi on być wyrażony jako suma czterech kwadratów, przy czym jeden z nich ma albo największą możliwą wartość, a=int(N**0.5)albo jeden mniejszą, pozostawiając jedynie niewielką pozostałą sumę dla trzech pozostałych kwadratów zaopiekować się. To znacznie zmniejsza przestrzeń wyszukiwania.

Oto dowód później, że ten kod zawsze znajduje rozwiązanie. Chcemy znaleźć atak, że n-a^2jest to suma trzech kwadratów. Z Twierdzenia Trzech Kwadratów Legendre'a liczba jest sumą trzech kwadratów, chyba że jest to forma 4^j(8*k+7). W szczególności takimi liczbami są 0 lub 3 (moduł 4).

Pokazujemy, że żadne dwa kolejne nie amogą sprawić, że pozostała ilość N-a^2ma taki kształt dla obu kolejnych wartości. Możemy to zrobić, po prostu tworząc tabelę ai Nmodulo 4, zauważając, że N%4!=0ponieważ wydobyliśmy wszystkie moce 4 z N.

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

Ponieważ nie ma dwóch kolejnych apodań (N-a*a)%4 in [0,3], jeden z nich jest bezpieczny w użyciu. Tak więc łapczywie używamy największej możliwej nz n^2<=N, i n-1. Ponieważ N<(n+1)^2reszta N-a^2reprezentowana jako suma trzech kwadratów jest co najwyżej (n+1)^2 -(n-1)^2równa 4*n. Wystarczy więc sprawdzić tylko wartości do 2*sqrt(n), czyli dokładnie tego zakresu R.

Można jeszcze bardziej zoptymalizować czas działania, zatrzymując się po pojedynczym rozwiązaniu, obliczając zamiast iterować ostatnią wartość di wyszukując tylko wartości b<=c<=d. Ale nawet bez tych optymalizacji najgorszy przypadek, jaki udało mi się skończyć w 45 sekund na moim komputerze.

Łańcuch „for x in R” jest niefortunny. Prawdopodobnie można go skrócić przez podstawienie ciągu lub zastąpienie przez iterację pojedynczego indeksu, który koduje (a, b, c, d). Importowanie narzędzi itertools nie było tego warte.

Edycja: Zmieniono na int(2*n**.5)+1z, 2*int(n**.5)+2aby argument był czystszy, ta sama liczba znaków.

CJam , 91 90 74 71 bajtów

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

Kompaktowy, ale wolniejszy niż moje inne podejście.

Wypróbuj online! Wklej kod , wpisz żądaną liczbę całkowitą w polu Input i kliknij Uruchom .

tło

Ten post rozpoczął się jako 99-bajtowa odpowiedź GolfScript . Podczas gdy wciąż było miejsce na ulepszenia, GolfScript nie ma wbudowanej funkcji sqrt. Zachowałem wersję GolfScript do wersji 5 , ponieważ była bardzo podobna do wersji CJam.

Jednak optymalizacje od wersji 6 wymagają operatorów, które nie są dostępne w GolfScript, więc zamiast publikować osobne objaśnienia dla obu języków, postanowiłem zrezygnować z mniej konkurencyjnej (i znacznie wolniejszej) wersji.

Zaimplementowany algorytm oblicza liczby za pomocą brutalnej siły:

1. W przypadku danych wejściowych mznajdź Ni Wtakie m = 4**W * N.

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

3. Set i = i + 1.

4. Znajdź liczby całkowite, j, k, ltakie jak i = 257**2 * j + 257 * k + l, gdzie k, l < 257.

5. Sprawdź, czy d = N - j**2 - k**2 - l**2jest to idealny kwadrat.

6. Jeśli nie, i wróć do kroku 3.

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

Przykłady

$ 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

Czasy odpowiadają procesorowi Intel Core i7-4700MQ.

Jak to działa

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

Jest to oparte na algorytmie na stronie Wikipedii, który jest brutalną siłą 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 bajtów

~{.[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>`

Szybki, ale długi. Nowa linia może zostać usunięta.

Wypróbuj online. Pamiętaj, że tłumacz online ma limit 5 sekund, więc może nie działać dla wszystkich numerów.

tło

Algorytm korzysta z twierdzenia trzech kwadratów Legendre'a , które stwierdza, że ​​każda liczba naturalna n, która nie jest formą

można wyrazić jako sumę trzech kwadratów.

Algorytm wykonuje następujące czynności:

1. Wyraź liczbę jako 4**i * j.

2. Znajdź największą liczbę całkowitą ktakie, że k**2 <= ji j - k**2spełnia hipotezę o trzech kwadratowy twierdzenie Legendre za.

3. Set i = 0.

4. Sprawdź, czy j - k**2 - (i / 252)**2 - (i % 252)**2jest to idealny kwadrat.

5. Jeśli nie, zwiększ wartość ii wróć do kroku 4.

Przykłady

$ 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

Czasy odpowiadają procesorowi Intel Core i7-4700MQ.

Jak to działa

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

Jest to nawet bardziej wymagające pod względem pamięci niż rev 0. Ta sama zasada: zbuduj tabelę z wartością dodatnią dla wszystkich możliwych sum 2 kwadratów (i zero dla liczb, które nie są sumami dwóch kwadratów), a następnie przeszukaj ją.

W tej wersji użyj tablicy shortzamiast chardo przechowywania trafień, dzięki czemu mogę przechowywać pierwiastek jednego z dwóch kwadratów w tabeli zamiast samej flagi. Upraszcza to pznacznie funkcję (do dekodowania sumy 2 kwadratów), ponieważ nie ma potrzeby stosowania pętli.

System Windows ma limit 2 GB na tablice. Mogę short s[15<<26]obejść ten, z którym jest tablica 1006632960 elementów, wystarczająca do spełnienia specyfikacji. Niestety, całkowita wielkość Runtime program jest nadal ponad 2 GB i (mimo szczypanie ustawienia OS) nie były w stanie przejechać ten rozmiar (choć jest to teoretycznie możliwe). Najlepsze, co mogę zrobić, to short s[14<<26](939524096 elementy.) m*mMusi być ściśle mniej niż to (30651 ^ 2 = 939483801.) Niemniej jednak program działa idealnie i powinien działać na każdym systemie operacyjnym, który nie ma tego ograniczenia.

Nieskluczony kod

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

To bestia głodna pamięci. Zajmuje tablicę 1 GB, oblicza wszystkie możliwe sumy 2 kwadratów i przechowuje flagę dla każdego w tablicy. Następnie dla wpisu z użytkownika, przeszukuje tablicę w poszukiwaniu dwóch sum 2 kwadratów a i za.

pnastępnie funkcja ponownie przelicza oryginalne kwadraty, które zostały użyte do utworzenia sumy 2 kwadratów ai z-adrukuje je, pierwsza z każdej pary jako liczba całkowita, druga jako liczba podwójna (jeśli muszą to być wszystkie liczby całkowite, potrzebne są jeszcze dwa znaki, t> m=t.)

Program zajmuje kilka minut, aby zbudować tabelę sum kwadratów (myślę, że jest to spowodowane problemami z zarządzaniem pamięcią, widzę, że przydzielanie pamięci rośnie powoli zamiast podskakiwać, jak można by się spodziewać.) Jednak po zakończeniu tej operacji generuje odpowiedzi bardzo szybko (jeśli kilka liczb ma być obliczonych, program od początku scanfmożna umieścić w pętli.

nieprzypisany kod

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)  

Przykładowe dane wyjściowe

Pierwszy dotyczy pytania. Drugi został wybrany jako trudny do wyszukania. W takim przypadku program musi szukać aż do 8192 ^ 2 + 8192 ^ 2 = 134217728, ale zajmuje to tylko kilka sekund po zbudowaniu tabeli.

123456789
0 2.000000 3328 10601.000000

805306368
8192 8192.000000 8192 24576.000000

Mathematica, 138 znaków

Okazuje się więc, że daje to negatywne i urojone wyniki dla niektórych danych wejściowych, jak wskazał edc65 (np. 805306368), więc nie jest to prawidłowe rozwiązanie. Na razie zostawię to i może, jeśli naprawdę nie znoszę swojego czasu, wrócę i spróbuję to naprawić.

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

Lub niewykonane:

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

Nie patrzyłem zbyt mocno na algorytmy, ale spodziewam się, że to ten sam pomysł. Właśnie wymyśliłem oczywiste rozwiązanie i modyfikowałem je, aż zadziałało. Przetestowałem go dla wszystkich liczb od 1 do 1 miliarda i ... działa. Test trwa tylko około 100 sekund na moim komputerze.

Zaletą tego jest to, że ponieważ b, c i d są zdefiniowane z opóźnieniami, :=nie trzeba ich przedefiniowywać, gdy a jest zmniejszane. Dzięki temu zaoszczędziłem kilka dodatkowych linii, które miałem wcześniej. Mogę pograć w golfa dalej i zagnieździć zbędne części, ale oto pierwszy szkic.

Aha, uruchamiasz go jak S@123456789i możesz go przetestować za pomocą {S@#, Total[(S@#)^2]} & @ 123456789lub # == Total[(S@#)^2]&[123456789]. Wyczerpujący test to

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

Wcześniej korzystałem z instrukcji Print [], ale bardzo ją to spowolniło, mimo że nigdy się nie wywołuje. Domyśl.

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]

Prosta brutalna siła ponad wstępnie obliczonymi kwadratami.

Potrzebuje -Oopcji kompilacji (do tego dodałem 3 znaki). W najgorszym przypadku 999950883 zajmuje mniej niż 1 minutę.

Testowane tylko na GHC.

C: 198 znaków

Prawdopodobnie uda mi się to zmniejszyć do nieco ponad 100 znaków. To, co lubię w tym rozwiązaniu, to minimalna ilość śmieci, tylko zwykła pętla for-for, robiąca to, co powinna zrobić pętla for (co ma być szalone).

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

I mocno upiększony:

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

Edycja: Nie jest wystarczająco szybki dla wszystkich danych wejściowych, ale wrócę z innym rozwiązaniem. Pozwolę, by ten potrójny bałagan operacji został na razie.

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

Dzięki @Dennis za ulepszenia. Pozostała część odpowiedzi poniżej nie jest aktualizowana od wersji 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);}

Znacznie szybciej niż moja inna odpowiedź i przy znacznie mniejszej pamięci!

Wykorzystuje to http://en.wikipedia.org/wiki/Legendre%27s_three-square_theorem . Dowolną liczbę nie w poniższej formie można wyrazić jako sumę 3 kwadratów (nazywam to formą zabronioną):

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

Zauważ, że wszystkie nieparzyste liczby kwadratowe mają postać, (8b+1)a wszystkie parzyste liczby kwadratowe są powierzchownie formą 4b. Ukrywa to jednak fakt, że wszystkie parzyste liczby kwadratowe mają formę 4^a*(odd square)==4^a*(8b+1). W rezultacie 2^x-(any square number < 2^(x-1))nieparzyste xzawsze będą miały zabronioną formę. Dlatego te liczby i ich wielokrotności są trudnymi przypadkami, dlatego tak wiele programów tutaj dzieli siły 4 jako pierwszy krok.

Jak stwierdzono w odpowiedzi @ xnor, N-a*anie może mieć formy zabronionej dla 2 kolejnych wartości a. Poniżej przedstawiam uproszczoną formę jego stołu. Oprócz faktu, że po podzieleniu przez 4 N%4nie może być równe 0, należy pamiętać, że istnieją tylko 2 możliwe wartości dla (a*a)%4.

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

Chcemy więc unikać wartości, (N-a*a)które mogą mieć niedozwoloną formę, a mianowicie te, w których (N-a*a)%4jest 3 lub 0. Jak widać, nie może się tak zdarzyć Nzarówno z nieparzystymi, jak i parzystymi (a*a).

Mój algorytm działa w ten sposób:

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

Szczególnie podoba mi się sposób, w jaki robię krok 3. Dodam 3 do N, więc zmniejszenie jest wymagane, jeśli (3+N-a*a)%4 =3 lub 2. (ale nie 1 lub 0.) Podziel to przez 2, a całą pracę można wykonać za pomocą dość prostego wyrażenia .

Nieskluczony kod

Zwróć uwagę na pojedynczą forpętlę qi użycie dzielenia / modulo, aby uzyskać z niej wartości bi c. Próbowałem użyć ajako dzielnika zamiast 256, aby zapisać bajty, ale czasami wartość anie była właściwa i program zawiesił się, być może na czas nieokreślony. 256 był najlepszym kompromisem, jaki mogę zastosować >>8zamiast /256do podziału.

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

Wynik

Ciekawym dziwactwem jest to, że jeśli wprowadzisz liczbę kwadratową, N-(a*a)= 0. Ale program wykrywa, że 0%4= 0 i zmniejsza do następnego kwadratu w dół. W rezultacie liczby kwadratowe są zawsze rozkładane na grupę mniejszych kwadratów, chyba że mają formę 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 znaków

Brutalna siła, ale szybka.

Edycja Szybka, ale niewystarczająca do niektórych nieprzyjemnych danych wejściowych. Musiałem dodać pierwszy krok redukcji przez wielokrotność 4.

Edycja 2 Pozbądź się funkcji, wyjdź wewnątrz pętli, a następnie wymuś rywalizację o wyjście

Edytuj 3 0 niepoprawne dane wejściowe

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='')

Nie golfowany:

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

Przykładowe dane wyjściowe

123456789
11111,48,10,8

805306368
16384,16384,16384,0