Lagranges Vierquadrat-Theorem besagt, dass jede natürliche Zahl als die Summe von vier Quadratzahlen dargestellt werden kann. Ihre Aufgabe ist es, ein Programm zu schreiben, das dies tut.

Input: Eine natürliche Zahl (unter 1 Milliarde)

Ausgabe: Vier Zahlen, deren Quadrate sich zu dieser Zahl addieren (Reihenfolge spielt keine Rolle)

Hinweis: Sie müssen keine Brute-Force-Suche durchführen! Details hier und hier . Wenn es eine Funktion gibt, die dieses Problem trivialisiert (ich werde feststellen), ist dies nicht zulässig. Automatisierte Primfunktionen und Quadratwurzel sind erlaubt. Wenn es mehr als eine Darstellung gibt, ist jede in Ordnung. Wenn Sie sich für brachiale Gewalt entscheiden, muss sie innerhalb einer angemessenen Zeit (3 Minuten) ausgeführt werden.

Probeneingabe

123456789

Beispielausgabe (entweder ist in Ordnung)

10601 3328 2 0
10601 3328 2

CJam, 50 Bytes

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

Meine dritte (und letzte, ich verspreche es) Antwort. Dieser Ansatz basiert stark auf der Antwort von primo .

Probieren Sie es online im CJam-Interpreter aus .

Verwendung

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

Hintergrund

1. Primo aktualisierten Algorithmus Nachdem ich, ich hatte , um zu sehen , wie ein CJam Implementierung würde punkten:

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

   Nur 58 Bytes! Dieser Algorithmus wird in nahezu konstanter Zeit ausgeführt und weist für verschiedene Werte von keine großen Schwankungen auf N. Lassen Sie uns das ändern ...

2. Anstatt bei zu beginnen floor(sqrt(N))und zu dekrementieren, können wir bei zu beginnen 1und zu inkrementieren. Das spart 4 Bytes.

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

3. Anstatt Nals 4**a * bauszudrücken, können wir es als p**(2a) * b- wo pist der kleinste Primfaktor von N- ausdrücken , um 1 weiteres Byte zu sparen.

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

4. Die vorherige Änderung ermöglicht es uns , die Umsetzung etwas zu ändern (ohne den Algorithmus selbst zu berühren): Anstelle der Unterteilung Ndurch p**(2a)und multiplizieren Sie die Lösung durch p**a, können wir direkt die möglichen Lösungen auf ein Vielfaches von beschränken p**a. Dies spart 2 weitere Bytes.

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

5. Wenn Sie die erste Ganzzahl nicht auf ein Vielfaches von beschränken, wird p**aein zusätzliches Byte gespeichert.

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

Endgültiger Algorithmus

1. Finde aund bso, dass N = p**(2a) * b, wo bnicht ein Vielfaches von p**2und pder kleinste Primfaktor von ist N.

2. Einstellen j = p**a.

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

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

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

6. Wenn ja N - j**2 - k**2 - l**2 - m**2 > 0, stellen Sie ein j = j + 1und fahren Sie mit Schritt 3 fort.

Dies kann wie folgt implementiert werden:

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

Ich habe auf meinem Intel Core i7-3770 alle 5 Versionen über alle positiven ganzen Zahlen bis zu 999.999.999 ausgeführt, die Ausführungszeit gemessen und die für die Lösungsfindung erforderlichen Iterationen gezählt.

Die folgende Tabelle zeigt die durchschnittliche Anzahl der Iterationen und die Ausführungszeit für eine einzelne Ganzzahl:

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. Mit nur 4 Iterationen und 6,6 Mikrosekunden pro Ganzzahl ist der Algorithmus von primo unglaublich schnell.

2. Ab floor(sqrt(N))ist sinnvoller, da wir dann für die Summe der verbleibenden drei Quadrate kleinere Werte haben. Wie erwartet ist der Start bei 1 viel langsamer.

3. Dies ist ein klassisches Beispiel für eine schlecht umgesetzte gute Idee. Um die Codegröße tatsächlich zu reduzieren , verlassen wir uns darauf mF, dass die ganze Zahl faktorisiert wird N. Obwohl Version 3 weniger Iterationen erfordert als Version 2, ist es in der Praxis viel langsamer.

4. Obwohl sich der Algorithmus nicht ändert, ist Version 4 viel langsamer. Dies liegt daran, dass in jeder Iteration eine zusätzliche Gleitkommadivision und eine ganzzahlige Multiplikation durchgeführt wird.

5. Für die Eingabe benötigt der N = p**(2a) ** bAlgorithmus 5 (k - 1) * p**a + 1Iterationen, wobei kdie Anzahl der Iterationen ist, die der Algorithmus 4 benötigt. Wenn k = 1oder a = 0, macht dies keinen Unterschied.

   Jede Eingabe des Formulars 4**a * (4**c * (8 * d + 7) + 1)kann jedoch sehr schlecht ausgeführt werden. Für den Startwert j = p**a, N - 4**a = 4**(a + c) * (8 * d + 7)so kann es nicht als eine Summe von drei Quadraten ausgedrückt werden. Somit sind k > 1zumindest p**aIterationen erforderlich.

   Zum Glück ist der ursprüngliche Algorithmus von primo unglaublich schnell und N < 1,000,000,000. Der schlechteste Fall, den ich von Hand finden konnte 265,289,728 = 4**10 * (4**1 * (7 * 8 + 7) + 1), erfordert 6.145 Iterationen. Die Ausführungszeit auf meinem Computer beträgt weniger als 300 ms. Im Durchschnitt ist diese Version 13,5-mal langsamer als die Implementierung des primo-Algorithmus.

FRACTRAN: 156 98 Fraktionen

Da dies ein klassisches Problem der Zahlentheorie ist, gibt es keinen besseren Weg, dies zu lösen, als Zahlen zu verwenden!

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

Nimmt die Eingabe der Form 2 ⁿ × 193 auf und gibt 3 ^a × 5 ^b × 7 ^c × 11 ^{d aus} . Könnte in 3 Minuten laufen, wenn Sie einen wirklich guten Dolmetscher haben. Vielleicht.

... okay, nicht wirklich. Dies schien in FRACTRAN ein so lustiges Problem zu sein, dass ich es versuchen musste. Offensichtlich ist dies keine richtige Lösung für die Frage, da es nicht den Zeitbedarf (es ist brachial) und es ist kaum Golf, aber ich dachte, ich würde dies hier posten, weil es nicht jeden Tag eine Codegolf-Frage ist kann in FRACTRAN gemacht werden;)

Hinweis

Der Code entspricht dem folgenden 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

Es werden drei Methoden gezeigt. Nur der erste Ansatz erfüllt den Zeitbedarf.

1-FindInstance (51 Zeichen)

Dies gibt eine einzelne Lösung der Gleichung zurück.

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

Beispiele und Timings

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

Dies funktioniert auch, ist jedoch zu langsam, um die Geschwindigkeitsanforderungen zu erfüllen.

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

Range[0, Floor@Sqrt@n]^2ist die Menge aller Quadrate, die kleiner ist als die Quadratwurzel von n(das größtmögliche Quadrat in der Partition).

{4}erfordert, dass die ganzzahligen Partitionen naus 4 Elementen aus der oben genannten Menge von Quadraten bestehen.

1, innerhalb der Funktion IntegerPartitionsgibt die erste Lösung zurück.

[[1]]entfernt die äußeren Klammern; Die Lösung wurde als Satz von einem Element zurückgegeben.

f[123456]

{348, 44, 20, 4}

3-PowerRepresentations

PowerRepresentations gibt alle Lösungen für das 4-Quadrate-Problem zurück. Es kann auch nach Summen anderer Mächte aufgelöst werden.

PowersRepresentations gibt in weniger als 5 Sekunden die 181 Möglichkeiten zurück, um 123456789 als die Summe von 4 Quadraten auszudrücken:

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

Für andere Beträge ist es jedoch viel zu langsam.

Perl - 116 Bytes 87 Bytes (siehe Update unten)

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

Zählt man den Shebang als ein Byte, werden Zeilenumbrüche hinzugefügt, um die horizontale Vernunft zu gewährleisten.

So etwas wie eine Kombination aus Code-Golf -Code- Einreichung mit dem schnellsten Code .

Die durchschnittliche (schlechteste?) Fallkomplexität scheint O (log n) O (n ^0,07 ) zu sein . Nichts, was ich gefunden habe, läuft langsamer als 0,001s und ich habe den gesamten Bereich von 900000000 - 999999999 überprüft . Wenn Sie etwas finden, das wesentlich länger dauert (~ 0,1s oder länger), lassen Sie es mich bitte wissen.

Beispielnutzung

$ 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

Die letzten beiden scheinen Worst-Case-Szenarien für andere Einreichungen zu sein. In beiden Fällen wird die gezeigte Lösung buchstäblich als erstes überprüft. Für 123456789ist es das zweite.

Wenn Sie einen Wertebereich testen möchten, können Sie das folgende Skript verwenden:

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

Am besten in eine Datei umleiten. 1..1000000Auf meinem Computer dauert der Bereich ungefähr 14 Sekunden (71000 Werte pro Sekunde), und der Bereich 999000000..1000000000dauert ungefähr 20 Sekunden (50000 Werte pro Sekunde), was der durchschnittlichen Komplexität von O (log n) entspricht .

Aktualisieren

Bearbeiten : Es stellt sich heraus, dass dieser Algorithmus einem sehr ähnlich ist, der von mentalen Taschenrechnern seit mindestens einem Jahrhundert verwendet wird .

Seit dem ursprünglichen Posten habe ich jeden Wert im Bereich von 1..1000000000 überprüft . Das "Worst-Case" -Verhalten wurde durch den Wert 699731569 angezeigt , bei dem insgesamt 190 Kombinationen getestet wurden, bevor eine Lösung gefunden wurde. Wenn Sie 190 als kleine Konstante betrachten - und das tue ich sicherlich -, kann das Worst-Case-Verhalten im erforderlichen Bereich als O (1) betrachtet werden . Das ist so schnell wie das Nachschlagen der Lösung von einem riesigen Tisch und im Durchschnitt möglicherweise sogar schneller.

Eine andere Sache. Nach 190 Iterationen hat alles, was größer als 144400 ist, noch nicht einmal den ersten Durchgang geschafft. Die Logik für die erste Durchquerung ist wertlos - sie wird nicht einmal verwendet. Der obige Code kann einiges gekürzt werden:

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

Welches führt nur den ersten Durchgang der Suche. Wir müssen jedoch bestätigen, dass es keine Werte unter 144400 gibt , für die der zweite Durchgang erforderlich ist:

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

Kurz gesagt, für den Bereich 1..1000000000 gibt es eine nahezu konstante Zeitlösung , die Sie sich ansehen .

Aktualisiertes Update

@Dennis und ich haben diesen Algorithmus mehrfach verbessert. Sie können den Fortschritt in den Kommentaren und der anschließenden Diskussion verfolgen, wenn dies Sie interessiert. Die durchschnittliche Anzahl der Iterationen für den erforderlichen Bereich ist von etwas mehr als 4 auf 1,229 gesunken , und die Zeit, die zum Testen aller Werte für 1..1000000000 erforderlich ist, wurde von 18 auf 2 Millionen Sekunden auf 41 Sekunden gesenkt . Der schlimmste Fall erforderte zuvor 190 Iterationen. Der schlimmste Fall, 854382778 , benötigt nur noch 21 .

Der endgültige Python-Code lautet wie folgt:

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

Dabei werden zwei vorberechnete Korrekturtabellen verwendet, eine mit einer Größe von 10 KB und eine mit einer Größe von 253 KB. Der obige Code enthält die Generatorfunktionen für diese Tabellen, obwohl diese wahrscheinlich zur Kompilierungszeit berechnet werden sollten.

Eine Version mit kleineren Korrekturtabellen finden Sie hier: http://codepad.org/1ebJC2OV Diese Version erfordert durchschnittlich 1.620 Iterationen pro Term, wobei der schlechteste Fall 38 ist. Die gesamte Reichweite beträgt ca. 3 m 21 s. Ein wenig Zeit wird wettgemacht, indem andfür die b- Korrektur bitweise anstelle von Modulo verwendet wird.

Verbesserungen

Gerade Werte ergeben eher eine Lösung als ungerade Werte.
Der Artikel zur mentalen Berechnung, der mit dem vorherigen Artikel verknüpft ist, stellt fest, dass, wenn nach dem Entfernen aller vier Faktoren der zu zerlegende Wert gerade ist, dieser Wert durch zwei geteilt und die Lösung rekonstruiert werden kann:

Obwohl dies für die mentale Berechnung sinnvoll sein kann (kleinere Werte sind in der Regel einfacher zu berechnen), ist es algorithmisch wenig sinnvoll. Wenn Sie 256 zufällige 4- Tupel nehmen und die Summe der Quadrate Modulo 8 untersuchen , werden Sie feststellen, dass die Werte 1 , 3 , 5 und 7 jeweils durchschnittlich 32- mal erreicht werden. Die Werte 2 und 6 werden jedoch jeweils 48 Mal erreicht. Wenn Sie ungerade Werte mit 2 multiplizieren , erhalten Sie im Durchschnitt eine Lösung mit 33% weniger Iterationen. Die Rekonstruktion ist die folgende:

Es muss darauf geachtet werden, dass a und b sowie c und d die gleiche Parität haben , aber wenn überhaupt eine Lösung gefunden wurde, ist eine ordnungsgemäße Reihenfolge garantiert.

Unmögliche Pfade müssen nicht überprüft werden.
Nach der Auswahl des zweiten Wertes, b , kann es bereits unmöglich sein, dass eine Lösung existiert, wenn die möglichen quadratischen Reste für ein gegebenes Modulo gegeben sind. Anstatt es trotzdem zu überprüfen oder mit der nächsten Iteration fortzufahren, kann der Wert von b "korrigiert" werden, indem er um den kleinsten Betrag verringert wird, der möglicherweise zu einer Lösung führen könnte. In den beiden Korrekturtabellen werden diese Werte gespeichert, einer für b und der andere für c . Die Verwendung eines höheren Modulos (genauer gesagt eines Moduls mit relativ weniger quadratischen Resten) führt zu einer besseren Verbesserung. Der Wert a muss nicht korrigiert werden. durch Ändern von n auf gerade werden alle Werte vona sind gültig.

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

Nachdem wir die Eingabe so reduziert haben N, dass sie nicht durch 4 teilbar ist, muss sie als Summe von vier Quadraten ausgedrückt werden, wobei eines der Quadrate entweder der größtmögliche Wert a=int(N**0.5)oder ein kleinerer Wert ist , sodass nur ein kleiner Rest für die Summe der drei anderen Quadrate übrig bleibt sich kümmern um. Dies reduziert den Suchraum erheblich.

Hier ist ein Beweis später findet dieser Code immer eine Lösung. Wir möchten ein aso finden, n-a^2das die Summe von drei Quadraten ist. Nach dem Drei-Quadrate-Satz von Legendre ist eine Zahl die Summe von drei Quadraten, sofern es sich nicht um die Form handelt 4^j(8*k+7). Insbesondere sind solche Zahlen entweder 0 oder 3 (Modulo 4).

Wir zeigen, dass keine zwei aufeinanderfolgenden a den Restbetrag N-a^2für beide aufeinanderfolgenden Werte eine solche Form haben können. Wir können dies tun, indem wir einfach eine Tabelle von aund Nmodulo 4 erstellen, wobei N%4!=0wir beachten, dass wir alle Potenzen von 4 aus extrahiert haben N.

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

Da keine zwei aufeinanderfolgenden ageben (N-a*a)%4 in [0,3], ist einer von ihnen sicher zu bedienen. Also nutzen wir gierig das größtmögliche nmitn^2<=N und n-1. Da N<(n+1)^2ist der Rest N-a^2, der als Summe von drei Quadraten dargestellt werden soll, höchstens (n+1)^2 -(n-1)^2gleich 4*n. Es genügt also, nur Werte bis zu überprüfen 2*sqrt(n), die genau dem Bereich entsprechen R.

Man könnte die Laufzeit weiter optimieren, indem man nach einer einzelnen Lösung stoppt, anstatt nach dem letzten Wert zu iterieren dund nur zwischen den Werten zu suchen b<=c<=d. Aber selbst ohne diese Optimierungen war die schlimmste Instanz, die ich finden konnte, in 45 Sekunden auf meinem Computer beendet.

Die Kette von "für x in R" ist unglücklich. Sie kann möglicherweise durch Zeichenfolgensubstitution oder durch Iteration über einen einzelnen Index, der (a, b, c, d) codiert, verkürzt werden. Der Import von itertools hat sich nicht gelohnt.

Bearbeiten: Geändert in int(2*n**.5)+1von 2*int(n**.5)+2, um das Argument sauberer zu machen, bei gleicher Zeichenanzahl.

CJam , 91 90 74 71 Bytes

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

Kompakt, aber langsamer als mein anderer Ansatz.

Probieren Sie es online! Fügen Sie den Code ein , geben Sie die gewünschte Ganzzahl in die Eingabe ein und klicken Sie auf Ausführen .

Hintergrund

Dieser Beitrag begann als 99-Byte-GolfScript-Antwort . Während es noch Verbesserungspotential gab, fehlt GolfScript die eingebaute sqrt-Funktion. Ich habe die GolfScript-Version bis beibehalten 5 , da sie der CJam-Version sehr ähnlich war.

Die Optimierungen seit Revision 6 erfordern jedoch Operatoren, die in GolfScript nicht verfügbar sind. Daher habe ich mich entschlossen, die weniger wettbewerbsfähige (und viel langsamere) Version zu streichen, anstatt separate Erklärungen für beide Sprachen zu veröffentlichen.

Der implementierte Algorithmus berechnet die Zahlen mit Brute Force:

1. Für die Eingabe mfinden Nund Wso m = 4**W * N.

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

3. Einstellen i = i + 1.

4. Finden Sie ganze Zahlen, j, k, lso dass i = 257**2 * j + 257 * k + l, wo k, l < 257.

5. Überprüfen Sie, ob d = N - j**2 - k**2 - l**2es sich um ein perfektes Quadrat handelt.

6. Wenn nicht, fahren Sie mit Schritt 3 fort.

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

Beispiele

$ 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

Die Timings entsprechen einem Intel Core i7-4700MQ.

Wie es funktioniert

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

Dies basiert auf dem Algorithmus auf der Wikipedia-Seite, bei dem es sich um eine O (n) Brute Force handelt.

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

Schnell, aber langwierig. Der Zeilenumbruch kann entfernt werden.

Probieren Sie es online aus. Beachten Sie, dass der Online-Dolmetscher ein Zeitlimit von 5 Sekunden hat und daher möglicherweise nicht für alle Nummern funktioniert.

Hintergrund

Der Algorithmus nutzt den Drei-Quadrate-Satz von Legendre , der besagt, dass jede natürliche Zahl n nicht von der Form ist

kann als die Summe von drei Quadraten ausgedrückt werden .

Der Algorithmus führt Folgendes aus:

1. Drücken Sie die Nummer als aus 4**i * j.

2. Finden Sie die größte Ganzzahl kso, dass k**2 <= jund j - k**2die Hypothese von Legendres Drei-Quadrat-Theorem erfüllt.

3. Einstellen i = 0.

4. Überprüfen Sie, ob j - k**2 - (i / 252)**2 - (i % 252)**2es sich um ein perfektes Quadrat handelt.

5. Ist dies nicht der Fall, erhöhen Sie den Wert iund fahren Sie mit Schritt 4 fort.

Beispiele

$ 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

Die Timings entsprechen einem Intel Core i7-4700MQ.

Wie es funktioniert

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

Dies ist noch speicherhungriger als Rev. 0. Gleiches Prinzip: Erstellen Sie eine Tabelle mit einem positiven Wert für alle möglichen Summen von 2 Quadraten (und Null für die Zahlen, die keine Summen von zwei Quadraten sind), und durchsuchen Sie sie dann.

Verwenden Sie in dieser Version ein Array von shortanstelle von char, um die Treffer zu speichern, damit ich die Wurzel eines der beiden Quadrate in der Tabelle anstelle eines Flags speichern kann. Dies vereinfacht die Funktion p(zum Decodieren der Summe von 2 Quadraten) erheblich, da keine Schleife erforderlich ist.

Windows hat eine Beschränkung von 2 GB für Arrays. Ich kann short s[15<<26]das umgehen, womit sich ein Array von 1006632960 Elementen ergibt, das ausreicht, um der Spezifikation zu entsprechen. Leider beträgt die Gesamtgröße der Programmlaufzeit immer noch mehr als 2 GB, und ich konnte diese Größe (obwohl dies theoretisch möglich ist) short s[14<<26](939524096 Elemente) trotz Optimierung der Betriebssystemeinstellungen nicht überschreiten.m*m Muß Streng weniger (30651 ^ 2 = 939483801.) Trotzdem läuft das Programm einwandfrei und sollte auf jedem Betriebssystem ohne diese Einschränkung funktionieren.

Ungolfed Code

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

Dies ist eine erinnerungshungrige Bestie. Es benötigt ein 1-GB-Array, berechnet alle möglichen Summen von 2 Quadraten und speichert für jedes ein Flag im Array. Dann durchsucht es für die Benutzereingabe z das Array nach zwei Summen von 2 Quadraten a und za.

Die Funktion setzt pdann die ursprünglichen Quadrate, aus denen die Summe von 2 Quadraten gebildet wurde, erneut zusammen aund gibt z-asie aus, wobei das erste jedes Paares als Ganzzahl und das zweite als Doppelzahl (wenn es sich um ganze Zahlen handeln muss , werden zwei weitere Zeichen benötigt). t> m=t.)

Das Programm benötigt ein paar Minuten, um die Tabelle der Quadratsummen zu erstellen (ich glaube, dies liegt an Problemen mit der Speicherverwaltung. Ich sehe, dass die Speicherzuweisung langsam ansteigt, anstatt wie erwartet zu springen.) Sobald dies erledigt ist Es liefert sehr schnell Antworten (wenn mehrere Zahlen berechnet werden sollen, kann das Programm von nun scanfan in eine Schleife gestellt werden.

ungolfed code

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)  

Beispielausgabe

Das erste ist pro die Frage. Der zweite war schwer zu finden. In diesem Fall muss das Programm bis zu 8192 ^ 2 + 8192 ^ 2 = 134217728 suchen, dauert jedoch nach dem Erstellen der Tabelle nur einige Sekunden.

123456789
0 2.000000 3328 10601.000000

805306368
8192 8192.000000 8192 24576.000000

Mathematica, 138 Zeichen

Es stellt sich also heraus, dass dies für bestimmte Eingaben negative und imaginäre Ergebnisse liefert, wie von edc65 (z. B. 805306368) hervorgehoben. Dies ist also keine gültige Lösung. Ich lasse es fürs Erste und wenn ich meine Zeit wirklich hasse, gehe ich zurück und versuche es zu reparieren.

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

Oder unsquished:

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

Ich habe mir die Algorithmen nicht so genau angesehen, aber ich gehe davon aus, dass dies die gleiche Idee ist. Ich habe gerade die offensichtliche Lösung gefunden und sie optimiert, bis sie funktioniert hat. Ich habe es für alle Zahlen zwischen 1 und 1 Milliarde getestet und ... es funktioniert. Der Test dauert auf meinem Rechner nur ca. 100 Sekunden.

Das Schöne daran ist, dass, da b, c und d mit verzögerten Zuweisungen definiert sind :=, sie nicht neu definiert werden müssen, wenn a dekrementiert wird. Dies ersparte ein paar zusätzliche Zeilen, die ich zuvor hatte. Ich könnte weiter Golf spielen und die überflüssigen Teile verschachteln, aber hier ist der erste Entwurf.

Oh, und du läufst es so S@123456789und kannst es mit {S@#, Total[(S@#)^2]} & @ 123456789oder testen # == Total[(S@#)^2]&[123456789]. Der ausführliche Test ist

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

Ich habe vorher eine Print [] -Anweisung verwendet, aber das hat sie sehr verlangsamt, obwohl sie nie aufgerufen wird. Stelle dir das vor.

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]

Einfache rohe Gewalt über vorberechneten Quadraten.

Es benötigt die -OKompilierungsoption (dazu habe ich 3 Zeichen hinzugefügt). Im schlimmsten Fall dauert es weniger als 1 Minute 999950883.

Nur auf GHC getestet.

C: 198 Zeichen

Ich kann es wahrscheinlich auf etwas mehr als 100 Zeichen reduzieren. Was ich an dieser Lösung mag, ist die minimale Menge an Junk, nur eine einfache for-Schleife, die das tut, was eine for-Schleife tun sollte (was verrückt sein soll).

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

Und stark verschönert:

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

Bearbeiten: Es ist nicht schnell genug für alle Eingaben, aber ich werde mit einer anderen Lösung zurück sein. Ich lasse dieses ternäre Operations-Chaos ab sofort bestehen.

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

Vielen Dank an @Dennis für die Verbesserungen. Der Rest der Antwort unten ist nicht von Rev. A aktualisiert.

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

Viel schneller als meine andere Antwort und mit viel weniger Gedächtnis!

Dies verwendet http://en.wikipedia.org/wiki/Legendre%27s_three-square_theorem . Jede Zahl, die nicht die folgende Form hat, kann als Summe von 3 Quadraten ausgedrückt werden (ich nenne dies die verbotene Form):

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

Beachten Sie, dass alle ungeraden Quadratzahlen von der Form sind (8b+1)und alle geraden Quadratzahlen oberflächlich von der Form sind 4b. Dies verbirgt jedoch die Tatsache, dass alle geraden quadratischen Zahlen von der Form sind 4^a*(odd square)==4^a*(8b+1). Als Ergebnis 2^x-(any square number < 2^(x-1))für ungerade xwird immer die verbotene Form sein. Daher sind diese Zahlen und ihre Vielfachen schwierige Fälle, weshalb so viele der Programme hier als ersten Schritt Potenzen von 4 aufteilen.

Wie in der Antwort von @ xnor angegeben, N-a*akann nicht für 2 aufeinanderfolgende Werte von die verbotene Form sein a. Unten stelle ich eine vereinfachte Form seiner Tabelle vor. Zusätzlich zu der Tatsache, dass nach der Division durch 4 N%4nicht gleich 0 sein kann, ist zu beachten, dass es nur 2 mögliche Werte für gibt (a*a)%4.

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

Wir wollen also vermeiden, (N-a*a)dass Werte , die möglicherweise von der verbotenen Form sind, nämlich solche mit (N-a*a)%43 oder 0. Wie zu sehen ist, kann dies Nbei ungeraden und geraden Werten nicht für dasselbe auftreten (a*a).

Mein Algorithmus funktioniert also so:

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

Ich mag besonders die Art und Weise, wie ich Schritt 3 mache. Ich addiere 3 zu N, so dass die Dekrementierung erforderlich ist, wenn (3+N-a*a)%4 =3 oder 2. (aber nicht 1 oder 0.) Teilen Sie dies durch 2 und die ganze Arbeit kann durch einen ziemlich einfachen Ausdruck erledigt werden .

Ungolfed Code

Beachten Sie die einzelne forSchleife qund die Verwendung von Division / Modulo die Werte abzuleiten bund cdaraus. Ich habe versucht, aanstelle von 256 Bytes einen Divisor zu verwenden, aber manchmal stimmte der Wert von anicht, und das Programm hing möglicherweise auf unbestimmte Zeit. 256 war der beste Kompromiss, den ich >>8anstelle /256für die Teilung verwenden kann.

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

Ausgabe

Eine interessante Eigenheit ist, dass bei Eingabe einer Quadratzahl N-(a*a)= 0. Aber das Programm erkennt, dass 0%4= 0 und dekrementiert zum nächsten Quadrat nach unten. Infolgedessen werden Quadratzahleneingaben immer in eine Gruppe kleinerer Quadrate zerlegt, es sei denn, sie haben die 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 Zeichen

Brachiale Kraft, aber schnell.

Schnelle Bearbeitung, aber nicht genug für schlechte Eingaben. Ich musste einen ersten Schritt der Reduzierung um ein Vielfaches von 4 hinzufügen.

Edit 2 Funktion entfernen, Ausgabe innerhalb der Schleife, dann Exit-Contition erzwingen

Edit 3 0 ist keine gültige Eingabe

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

Beispielausgabe

123456789
11111,48,10,8

805306368
16384,16384,16384,0