Le théorème des quatre carrés de Lagrange nous dit que tout nombre naturel peut être représenté comme la somme de quatre nombres carrés. Votre tâche consiste à écrire un programme qui fait cela.

Entrée: un nombre naturel (inférieur à 1 milliard)

Sortie: quatre nombres dont les carrés totalisent ce nombre (l'ordre n'a pas d'importance)

Remarque: vous n'avez pas à effectuer de recherche par force brute! Détails ici et ici . S'il existe une fonction qui banalise ce problème (je vais le déterminer), elle n'est pas autorisée. Les fonctions principales automatisées et la racine carrée sont autorisées. S'il y a plus d'une représentation, tout va bien. Si vous avez choisi de faire de la force brute, elle doit s'exécuter dans un délai raisonnable (3 minutes)

échantillon d'entrée

123456789

échantillon de sortie (soit est très bien)

10601 3328 2 0
10601 3328 2

CJam, 50 octets

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

Ma troisième (et dernière, je le promets) réponse. Cette approche est fortement basée sur la réponse de Primo .

Essayez-le en ligne dans l' interpréteur CJam .

Usage

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

Contexte

1. Après avoir vu l'algorithme mis à jour de primo, j'ai dû voir comment une implémentation CJam marquerait:

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

   Seulement 58 octets! Cet algorithme fonctionne en temps presque constant et ne présente pas beaucoup de variation pour différentes valeurs de N. Changeons cela ...

2. Au lieu de commencer à floor(sqrt(N))et décrémenter, nous pouvons commencer à 1et incrémenter. Cela économise 4 octets.

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

3. Au lieu d'exprimer en Ntant que 4**a * b, nous pouvons l'exprimer comme p**(2a) * b- où pest le plus petit facteur premier de N- pour économiser 1 octet de plus.

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

4. La modification précédente nous permet de modifier légèrement l'implémentation (sans toucher à l'algorithme lui-même): Au lieu de diviser Npar p**(2a)et de multiplier la solution par p**a, nous pouvons directement restreindre les solutions possibles à des multiples de p**a. Cela économise 2 octets supplémentaires.

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

5. Ne pas restreindre le premier entier à des multiples de p**asauvegarde un octet supplémentaire.

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

Algorithme final

1. Trouvez aet btel que N = p**(2a) * b, où bn'est pas un multiple de p**2et pest le plus petit facteur premier 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. Si N - j**2 - k**2 - l**2 - m**2 > 0, réglez j = j + 1et revenez à l'étape 3.

Cela peut être implémenté comme suit:

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

Repères

J'ai exécuté les 5 versions sur tous les entiers positifs jusqu'à 999 999 999 sur mon Intel Core i7-3770, mesuré le temps d'exécution et compté les itérations nécessaires pour trouver une solution.

Le tableau suivant montre le nombre moyen d'itérations et le temps d'exécution pour un seul entier:

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. À seulement 4 itérations et 6,6 micro secondes par entier, l'algorithme de primo est incroyablement rapide.

2. Commencer par floor(sqrt(N))est plus logique, car cela nous laisse des valeurs plus petites pour la somme des trois carrés restants. Comme prévu, commencer à 1 est beaucoup plus lent.

3. Ceci est un exemple classique d'une bonne idée mal mise en œuvre. Pour réellement réduire la taille du code, nous nous appuyons sur mF, qui factorise l'entier N. Bien que la version 3 nécessite moins d'itérations que la version 2, elle est beaucoup plus lente en pratique.

4. Bien que l'algorithme ne change pas, la version 4 est beaucoup plus lente. En effet, il effectue une division en virgule flottante supplémentaire et une multiplication entière à chaque itération.

5. Pour l'entrée N = p**(2a) ** b, l'algorithme 5 nécessitera des (k - 1) * p**a + 1itérations, où kest le nombre d'itérations requis par l'algorithme 4. Si k = 1ou a = 0, cela ne fait aucune différence.

   Cependant, toute entrée du formulaire 4**a * (4**c * (8 * d + 7) + 1)peut fonctionner assez mal. Pour la valeur de départ j = p**a, N - 4**a = 4**(a + c) * (8 * d + 7)elle ne peut donc pas être exprimée comme une somme de trois carrés. Ainsi, k > 1et au moins des p**aitérations sont nécessaires.

   Heureusement, l'algorithme original de primo est incroyablement rapide et N < 1,000,000,000. Le pire cas que j'ai pu trouver à la main est celui 265,289,728 = 4**10 * (4**1 * (7 * 8 + 7) + 1)qui nécessite 6 145 itérations. Le temps d'exécution est inférieur à 300 ms sur ma machine. En moyenne, cette version est 13,5 fois plus lente que l'implémentation de l'algorithme de primo.

FRACTRAN: 156 98 fractions

Puisqu'il s'agit d'un problème classique de théorie des nombres, quelle meilleure façon de résoudre ce problème que d'utiliser des nombres!

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

Prend en entrée la forme 2 ⁿ × 193 et ​​sort 3 ^a × 5 ^b × 7 ^c × 11 ^d . Peut fonctionner en 3 minutes si vous avez un très bon interprète. Peut être.

... d'accord, pas vraiment. Cela semblait être un problème tellement amusant à faire dans FRACTRAN que j'ai dû l' essayer. Évidemment, ce n'est pas une bonne solution à la question car cela ne fait pas les exigences de temps (ça force brutalement) et c'est à peine même joué au golf, mais je pensais que je posterais ceci ici parce que ce n'est pas tous les jours qu'une question Codegolf peut se faire en FRACTRAN;)

Allusion

Le code est équivalent au pseudo-Python suivant:

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

Trois méthodes sont présentées. Seule la première approche répond aux exigences de temps.

1-FindInstance (51 caractères)

Cela renvoie une seule solution l'équation.

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

Exemples et délais

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

Cela fonctionne également, mais est trop lent pour répondre aux exigences de vitesse.

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

Range[0, Floor@Sqrt@n]^2est l'ensemble de tous les carrés inférieur à la racine carrée de n(le plus grand carré possible dans la partition).

{4}nécessite que les partitions entières se ncomposent de 4 éléments de l'ensemble de carrés susmentionné.

1, dans la fonction IntegerPartitionsrenvoie la première solution.

[[1]]supprime les accolades extérieures; la solution a été renvoyée sous la forme d'un ensemble d'un élément.

f[123456]

{348, 44, 20, 4}

3-PowerReprésentations

PowerRepresentations renvoie toutes les solutions au problème des 4 carrés. Il peut également résoudre des sommes d'autres pouvoirs.

PowersRepresentations renvoie, en moins de 5 secondes, les 181 façons d'exprimer 123456789 comme la somme de 4 carrés:

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

Cependant, il est beaucoup trop lent pour d'autres sommes.

Perl - 116 octets 87 octets (voir la mise à jour ci-dessous)

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

En comptant le shebang comme un octet, des nouvelles lignes ont été ajoutées pour une raison horizontale.

Quelque chose d'une combinaison de code-golf - soumission de code le plus rapide .

La complexité moyenne (pire?) Des cas semble être O (log n) O (n ^0,07 ) . Rien de ce que j'ai trouvé ne tourne plus lentement que 0,001 s, et j'ai vérifié toute la plage de 900000000 à 999999999 . Si vous trouvez quelque chose qui prend beaucoup plus de temps que cela, ~ 0,1 s ou plus, faites-le moi savoir.

Exemple d'utilisation

$ 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

Les deux derniers semblent être les scénarios les plus défavorables pour les autres soumissions. Dans les deux cas, la solution présentée est littéralement la toute première chose vérifiée. Car 123456789c'est le second.

Si vous souhaitez tester une plage de valeurs, vous pouvez utiliser le script suivant:

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

Meilleur lorsqu'il est redirigé vers un fichier. La plage 1..1000000prend environ 14 secondes sur mon ordinateur (71000 valeurs par seconde) et la plage 999000000..1000000000prend environ 20 secondes (50000 valeurs par seconde), conformément à la complexité moyenne O (log n) .

Mise à jour

Edit : Il s'avère que cet algorithme est très similaire à celui qui a été utilisé par les calculatrices mentales pendant au moins un siècle .

Depuis la publication initiale, j'ai vérifié chaque valeur sur la plage de 1 à 1000000000 . Le comportement du «pire des cas» a été démontré par la valeur 699731569 , qui a testé un grand total de 190 combinaisons avant d'arriver à une solution. Si vous considérez 190 comme une petite constante - et je le fais certainement - le pire des cas sur la plage requise peut être considéré comme O (1) . Autrement dit, aussi rapide que la recherche de la solution à partir d'une table géante, et en moyenne, peut-être plus rapidement.

Une autre chose cependant. Après 190 itérations, rien de plus grand que 144400 n'a même dépassé le premier passage. La logique de la traversée en largeur est sans valeur - elle n'est même pas utilisée. Le code ci-dessus peut être raccourci un peu:

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

Qui n'effectue que la première passe de la recherche. Nous devons cependant confirmer qu'aucune valeur inférieure à 144400 n'a nécessité la deuxième passe:

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

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

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

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

En bref, pour la gamme 1..1000000000 , une solution à temps quasi constant existe, et vous la regardez.

Mise à jour mise à jour

@Dennis et moi avons apporté plusieurs améliorations à cet algorithme. Vous pouvez suivre les progrès dans les commentaires ci-dessous, et la discussion ultérieure, si cela vous intéresse. Le nombre moyen d'itérations pour la plage requise est passé d'un peu plus de 4 à 1,229 , et le temps nécessaire pour tester toutes les valeurs pour 1..1000000000 a été amélioré de 18 m 54 s à 2 m 41 s. Le pire des cas nécessitait auparavant 190 itérations; le pire des cas maintenant, 854382778 , n'a besoin que de 21 .

Le code Python final est le suivant:

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

Celui-ci utilise deux tables de correction pré-calculées, l'une de 10 Ko, l'autre de 253 Ko. Le code ci-dessus inclut les fonctions de générateur pour ces tables, bien qu'elles devraient probablement être calculées au moment de la compilation.

Une version avec des tables de correction de taille plus modeste peut être trouvée ici: http://codepad.org/1ebJC2OV Cette version nécessite une moyenne de 1.620 itérations par terme, avec le pire des cas de 38 , et toute la gamme fonctionne en environ 3m 21s. Un peu de temps est compensé, en utilisant au niveau du bit andpour la correction b , plutôt que modulo.

Améliorations

Les valeurs paires sont plus susceptibles de produire une solution que les valeurs impaires.
L'article de calcul mental lié à précédemment note que si, après avoir supprimé tous les facteurs de quatre, la valeur à décomposer est paire, cette valeur peut être divisée par deux, et la solution reconstruite:

Bien que cela puisse être logique pour le calcul mental (les valeurs plus petites ont tendance à être plus faciles à calculer), cela n'a pas beaucoup de sens algorithmiquement. Si vous prenez 256 4 -tuples aléatoires et examinez la somme des carrés modulo 8 , vous constaterez que les valeurs 1 , 3 , 5 et 7 sont chacune atteintes en moyenne 32 fois. Cependant, les valeurs 2 et 6 sont chacune atteintes 48 fois. La multiplication des valeurs impaires par 2 trouvera une solution, en moyenne, dans 33% d' itérations en moins. La reconstruction est la suivante:

Il faut veiller à ce que a et b aient la même parité, ainsi que c et d , mais si une solution a été trouvée, un bon ordre est garanti d'exister.

Les chemins impossibles n'ont pas besoin d'être vérifiés.
Après avoir sélectionné la deuxième valeur, b , il peut déjà être impossible pour une solution d'exister, étant donné les résidus quadratiques possibles pour un modulo donné. Au lieu de vérifier de toute façon ou de passer à l'itération suivante, la valeur de b peut être «corrigée» en la décrémentant de la plus petite quantité qui pourrait éventuellement conduire à une solution. Les deux tables de correction stockent ces valeurs, l'une pour b et l'autre pour c . L'utilisation d'un modulo plus élevé (plus précisément, l'utilisation d'un modulo avec relativement moins de résidus quadratiques) se traduira par une meilleure amélioration. La valeur a n'a besoin d'aucune correction; en modifiant n pour qu'il soit pair, toutes les valeurs dea sont valides.

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

Après avoir réduit l'entrée Npour qu'elle ne soit pas divisible par 4, elle doit être exprimable sous la forme d'une somme de quatre carrés où l'un d'eux est soit la plus grande valeur possible, a=int(N**0.5)soit un de moins que cela, ne laissant qu'un petit reste pour la somme des trois autres carrés prendre soin de. Cela réduit considérablement l'espace de recherche.

Voici une preuve plus tard ce code trouve toujours une solution. Nous souhaitons trouver une asorte qui n-a^2soit la somme de trois carrés. D'après le théorème des trois carrés de Legendre , un nombre est la somme de trois carrés, sauf s'il s'agit de la forme 4^j(8*k+7). En particulier, ces nombres sont soit 0 soit 3 (modulo 4).

Nous montrons que deux consécutifs apeuvent faire la quantité de restes N-a^2ont une telle forme pour les deux valeurs consécutives .. Nous pouvons le faire en faisant simplement une table de aet Nmodulo 4, notant que N%4!=0parce que nous avons extrait tous les pouvoirs de 4 sur N.

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

Parce que pas deux années consécutives adonner (N-a*a)%4 in [0,3], l' un d'entre eux est sûr à utiliser. Donc, nous utilisons avidement le plus grand possible navec n^2<=N, et n-1. Depuis N<(n+1)^2, le reste N-a^2à représenter comme une somme de trois carrés est au plus (n+1)^2 -(n-1)^2, ce qui équivaut 4*n. Il suffit donc de vérifier uniquement les valeurs jusqu'à 2*sqrt(n), ce qui correspond exactement à la plage R.

On pourrait encore optimiser le temps d'exécution en s'arrêtant après une seule solution, en calculant plutôt qu'en itérant pour la dernière valeur det en recherchant uniquement parmi les valeurs b<=c<=d. Mais, même sans ces optimisations, le pire cas que j'ai pu trouver s'est terminé en 45 secondes sur ma machine.

La chaîne "for x in R" est regrettable. Il peut probablement être raccourci par substitution ou remplacement de chaîne en itérant sur un seul index qui code (a, b, c, d). L'importation d'itertools s'est avérée inutile.

Edit: Changé à int(2*n**.5)+1partir 2*int(n**.5)+2de faire plus propre argument, même nombre de caractères.

CJam , 91 90 74 71 octets

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

Compact, mais plus lent que mon autre approche.

Essayez-le en ligne! Collez le code , saisissez l'entier souhaité dans Entrée et cliquez sur Exécuter .

Contexte

Ce message a commencé comme une réponse GolfScript de 99 octets . Bien qu'il y ait encore place à amélioration, GolfScript manque de fonction sqrt intégrée. J'ai gardé la version GolfScript jusqu'à la révision 5 , car elle était très similaire à la version CJam.

Cependant, les optimisations depuis la révision 6 nécessitent des opérateurs qui ne sont pas disponibles dans GolfScript, donc au lieu de publier des explications distinctes pour les deux langues, j'ai décidé d'abandonner la version moins compétitive (et beaucoup plus lente).

L'algorithme implémenté calcule les nombres par force brute:

1. Pour l'entrée m, trouver Net Wtel que m = 4**W * N.

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

3. Set i = i + 1.

4. Trouvez des entiers j, k, ltels que i = 257**2 * j + 257 * k + l, où k, l < 257.

5. Vérifiez si d = N - j**2 - k**2 - l**2c'est un carré parfait.

6. Si ce n'est pas le cas, et revenez à l'étape 3.

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

Exemples

$ 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

Les horaires correspondent à un Intel Core i7-4700MQ.

Comment ça fonctionne

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

Ceci est basé sur l'algorithme de la page Wikipedia, qui est une force brute 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 octets

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

Rapide, mais long. La nouvelle ligne peut être supprimée.

Essayez-le en ligne. Notez que l'interprète en ligne a un délai de 5 secondes, il peut donc ne pas fonctionner pour tous les numéros.

Contexte

L'algorithme tire parti du théorème de Legendre à trois carrés , qui stipule que tout nombre naturel n qui n'est pas de la forme

peut être exprimée comme la somme de trois carrés.

L'algorithme fait ce qui suit:

1. Exprimez le nombre comme 4**i * j.

2. Trouvez le plus grand entier ktel que k**2 <= jet j - k**2vérifie l'hypothèse du théorème de Legendre à trois carrés.

3. Set i = 0.

4. Vérifiez si j - k**2 - (i / 252)**2 - (i % 252)**2c'est un carré parfait.

5. Si ce n'est pas le cas, incrémentez iet revenez à l'étape 4.

Exemples

$ 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

Les horaires correspondent à un Intel Core i7-4700MQ.

Comment ça fonctionne

~              # 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.

Rév.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();}

C'est encore plus de mémoire que rev 0. Même principe: construisez une table avec une valeur positive pour toutes les sommes possibles de 2 carrés (et zéro pour les nombres qui ne sont pas des sommes de deux carrés), puis recherchez-la.

Dans cette version, utilisez un tableau de shortau lieu de charpour stocker les résultats, afin que je puisse stocker la racine de l'une des paires de carrés dans la table au lieu d'un simple indicateur. Cela simplifie considérablement la fonction p(pour décoder la somme de 2 carrés) car il n'y a pas besoin de boucle.

Windows a une limite de 2 Go sur les tableaux. Je peux contourner ce short s[15<<26]qui est un tableau de 1006632960 éléments, assez pour se conformer aux spécifications. Malheureusement, la taille totale d'exécution du programme est toujours supérieure à 2 Go et (malgré des ajustements des paramètres du système d'exploitation), je n'ai pas pu exécuter cette taille (bien que cela soit théoriquement possible.) Le mieux que je puisse faire est short s[14<<26](939524096 éléments.) m*mDoit être strictement inférieur à cela (30651 ^ 2 = 939483801.) Néanmoins, le programme fonctionne parfaitement et devrait fonctionner sur tout système d'exploitation qui n'a pas cette restriction.

Code non golfé

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

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

Ceci est une bête affamée de mémoire. Il prend un tableau de 1 Go, calcule toutes les sommes possibles de 2 carrés et stocke un indicateur pour chacun dans le tableau. Ensuite, pour l'entrée utilisateur z, il recherche dans le tableau deux sommes de 2 carrés a et za.

la fonction preconstitue ensuite les carrés originaux qui ont été utilisés pour faire les sommes de 2 carrés aet les z-aimprime, le premier de chaque paire sous forme d'entier, le second sous forme de double (s'il doit s'agir de tous les entiers, deux caractères supplémentaires sont nécessaires, t> m=t.)

Le programme prend quelques minutes pour construire le tableau des sommes des carrés (je pense que cela est dû à des problèmes de gestion de la mémoire, je vois l'allocation de mémoire augmenter lentement au lieu de sauter comme on pourrait s'y attendre.) Cependant, une fois cela fait, produit des réponses très rapidement (si plusieurs nombres doivent être calculés, le programme à scanfpartir de maintenant peut être mis en boucle.

code non golfé

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

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

Exemple de sortie

Le premier est par la question. Le second a été choisi comme difficile à rechercher. Dans ce cas, le programme doit rechercher jusqu'à 8192 ^ 2 + 8192 ^ 2 = 134217728, mais ne prend que quelques secondes une fois la table créée.

123456789
0 2.000000 3328 10601.000000

805306368
8192 8192.000000 8192 24576.000000

Mathematica, 138 caractères

Il s'avère donc que cela produit des résultats négatifs et imaginaires pour certaines entrées comme indiqué par edc65 (par exemple, 805306368), donc ce n'est pas une solution valide. Je vais le laisser pour l'instant, et peut-être, si je déteste vraiment mon temps, je vais y retourner et essayer de le réparer.

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, sans tremblement:

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

Je n'ai pas trop regardé les algorithmes, mais je pense que c'est la même idée. Je viens juste de trouver la solution évidente et je l'ai modifiée jusqu'à ce qu'elle fonctionne. Je l'ai testé pour tous les nombres entre 1 et un milliard et ... ça marche. Le test ne prend qu'environ 100 secondes sur ma machine.

Le bon côté à ce sujet est que, puisque b, c et d sont définis avec des affectations retardées :=, ils n'ont pas besoin d'être redéfinis lorsque a est décrémenté. Cela a sauvé quelques lignes supplémentaires que j'avais auparavant. Je pourrais le jouer plus loin et imbriquer les parties redondantes, mais voici le premier projet.

Oh, et vous l'exécutez en tant que S@123456789et vous pouvez le tester avec {S@#, Total[(S@#)^2]} & @ 123456789ou # == Total[(S@#)^2]&[123456789]. Le test exhaustif est

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

J'ai utilisé une instruction Print [] auparavant, mais cela l'a beaucoup ralentie, même si elle n'est jamais appelée. Allez comprendre.

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]

Force brute simple sur des carrés pré-calculés.

Il a besoin de l' -Ooption de compilation (j'ai ajouté 3 caractères pour cela). Cela prend moins d'une minute pour le pire des cas 999950883.

Testé uniquement sur GHC.

C: 198 caractères

Je peux probablement le réduire à un peu plus de 100 caractères. Ce que j'aime dans cette solution, c'est la quantité minimale d'ordure, juste une simple boucle for, faisant ce qu'une boucle for devrait faire (ce qui est fou).

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

Et fortement embelli:

#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: Ce n'est pas assez rapide pour toutes les entrées, mais je reviendrai avec une autre solution. Je vais laisser ce désordre d'opération ternaire rester à partir de maintenant.

Rév 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);}

Merci à @Dennis pour les améliorations. Le reste de la réponse ci-dessous n'est pas mis à jour à partir de la rév. A.

Rév 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);}

Beaucoup plus rapide que mon autre réponse et avec beaucoup moins de mémoire!

Cela utilise http://en.wikipedia.org/wiki/Legendre%27s_three-square_theorem . Tout nombre non de la forme suivante peut être exprimé comme la somme de 3 carrés (j'appelle cela la forme interdite):

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

Notez que tous les nombres carrés impairs sont de la forme (8b+1)et tous les nombres carrés pairs sont superficiellement de la forme 4b. Cependant, cela masque le fait que tous les nombres pairs sont de la forme 4^a*(odd square)==4^a*(8b+1). Par conséquent, 2^x-(any square number < 2^(x-1))pour impair xsera toujours de la forme interdite. Par conséquent, ces nombres et leurs multiples sont des cas difficiles, c'est pourquoi tant de programmes ici divisent les puissances de 4 dans un premier temps.

Comme indiqué dans la réponse de @ xnor, N-a*ane peut pas être de la forme interdite pour 2 valeurs consécutives de a. Ci-dessous, je présente une forme simplifiée de son tableau. En plus du fait qu'après division par 4 N%4ne peut pas être égal à 0, notez qu'il n'y a que 2 valeurs possibles pour (a*a)%4.

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

Donc, nous voulons éviter que les valeurs de (N-a*a)cela puissent être de la forme interdite, à savoir celles où (N-a*a)%4est 3 ou 0. Comme on peut le voir, cela ne peut pas se produire pour la même chose Nà la fois impaire et pair(a*a) .

Donc, mon algorithme fonctionne comme ceci:

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

J'aime particulièrement la façon dont je fais l'étape 3. J'ajoute 3 à N, de sorte que le décrément est requis si (3+N-a*a)%4 =3 ou 2. (mais pas 1 ou 0.) Divisez cela par 2 et le travail peut être fait par une expression assez simple .

Code non golfé

Notez la seule forboucle qet l' utilisation de la division / modulo pour dériver les valeurs de bet cde lui. J'ai essayé d'utiliser acomme diviseur au lieu de 256 pour économiser des octets, mais parfois la valeur de an'était pas correcte et le programme se bloquait, peut-être indéfiniment. 256 était le meilleur compromis que je peux utiliser >>8au lieu de /256pour la division.

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

Production

Une bizarrerie intéressante est que si vous entrez un nombre carré, N-(a*a)= 0. Mais le programme détecte que 0%4= 0 et diminue au carré suivant. En conséquence, les entrées de nombres carrés sont toujours décomposées en un groupe de carrés plus petits, sauf si elles sont de la forme 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 - 175191176 173 caractères

Force brute, mais rapide.

Édition rapide mais pas suffisante pour une entrée désagréable. J'ai dû ajouter une première étape de réduction par multiplications de 4.

Edit 2 Débarrassez-vous de la fonction, sortez dans la boucle puis forcez la sortie

Modifier 3 0 n'est pas une entrée valide

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

Non golfé:

v = prompt();

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

Exemple de sortie

123456789
11111,48,10,8

805306368
16384,16384,16384,0