Теорема Лагранжа о четырех квадратах говорит нам, что любое натуральное число может быть представлено как сумма четырех квадратных чисел. Ваша задача - написать программу, которая делает это.

Ввод: натуральное число (ниже 1 миллиарда)

Вывод: четыре числа, чьи квадраты суммируются с этим числом (порядок не имеет значения)

Примечание: вам не нужно искать грубой силой! Подробности здесь и здесь . Если есть функция, которая тривиализирует эту проблему (я буду определять), это не разрешено. Автоматические простые функции и квадратный корень разрешены. Если есть более одного представления, то все в порядке. Если вы решили применить грубую силу, она должна быть запущена в течение разумного времени (3 минуты)

образец ввода

123456789

образец вывода (либо в порядке)

10601 3328 2 0
10601 3328 2

CJam, 50 байтов

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

Мой третий (и последний, я обещаю) ответ. Этот подход в значительной степени основан на ответе primo .

Попробуйте онлайн в интерпретаторе CJam .

использование

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

Фон

1. После просмотра обновленного алгоритма primo я должен был увидеть, как будет работать реализация CJam:

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

   Всего 58 байт! Этот алгоритм работает в почти постоянное время и не демонстрирует большой разброс при различных значениях N. Давайте изменим это ...

2. Вместо того, чтобы начинать с floor(sqrt(N))и уменьшать, мы можем начинать с 1и увеличивать. Это экономит 4 байта.

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

3. Вместо того, чтобы выражать Nкак 4**a * b, мы можем выразить это как p**(2a) * b- где pнаименьший простой фактор N- чтобы сохранить еще 1 байт.

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

4. Предыдущая модификация позволяет слегка изменить реализацию (не касаясь сам алгоритм): Вместо деления Nна p**(2a)и умножить на решении, p**aмы можем непосредственно ограничить возможные решения кратных p**a. Это экономит еще 2 байта.

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

5. Не ограничивая первое целое число кратным, p**aсохраняет дополнительный байт.

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

Окончательный алгоритм

1. Найти aи bтакое, что N = p**(2a) * b, где bне кратно p**2и pявляется наименьшим простым множителем 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. Если N - j**2 - k**2 - l**2 - m**2 > 0, установите j = j + 1и вернитесь к шагу 3.

Это может быть реализовано следующим образом:

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

Ориентиры

Я выполнил все 5 версий для всех натуральных чисел до 999 999 999 на моем Intel Core i7-3770, измерил время выполнения и посчитал количество итераций, необходимых для поиска решения.

В следующей таблице показано среднее количество итераций и время выполнения для одного целого числа:

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. Всего за 4 итерации и 6,6 микросекунды на целое число алгоритм primo невероятно быстр.

2. Начиная с floor(sqrt(N))более логично, так как мы получаем меньшие значения для суммы оставшихся трех квадратов. Как и ожидалось, начиная с 1 намного медленнее.

3. Это классический пример плохо реализованной хорошей идеи. Чтобы уменьшить размер кода, мы полагаемся mF, что факторизирует целое число N. Хотя версия 3 требует меньше итераций, чем версия 2, на практике она намного медленнее.

4. Хотя алгоритм не меняется, версия 4 намного медленнее. Это потому, что он выполняет дополнительное деление с плавающей запятой и целочисленное умножение в каждой итерации.

5. Для ввода N = p**(2a) ** bалгоритму 5 потребуются (k - 1) * p**a + 1итерации, где kчисло итераций, необходимое алгоритму 4. Если k = 1или a = 0, это не имеет значения.

   Тем не менее, любой ввод формы 4**a * (4**c * (8 * d + 7) + 1)может работать довольно плохо. Для начального значения j = p**a, N - 4**a = 4**(a + c) * (8 * d + 7)поэтому его нельзя выразить как сумму трех квадратов. Таким образом, k > 1и хотя бы p**aитерации не требуется.

   К счастью, оригинальный алгоритм Primo невероятно быстр и N < 1,000,000,000. Худший случай, который я мог найти вручную, - 265,289,728 = 4**10 * (4**1 * (7 * 8 + 7) + 1)это 6,145 итераций. Время выполнения составляет менее 300 мс на моей машине. В среднем эта версия в 13,5 раз медленнее, чем реализация алгоритма primo.

ФРАКТРАН: 156 98 фракций

Поскольку это классическая проблема теории чисел, что может быть лучше, чем использовать числа!

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

Принимает на входе форму 2 ⁿ × 193 и выводит 3 ^a × 5 ^b × 7 ^c × 11 ^d . Бегите за 3 минуты, если у вас действительно хороший переводчик. Может быть.

... хорошо, не совсем. Это казалось такой забавной проблемой во ФРАКТРАНЕ, что мне пришлось ее попробовать. Очевидно, что это не правильное решение вопроса, поскольку оно не предъявляет требований к времени (это грубые силы), и это даже не игра в гольф, но я подумал, что я опубликую это здесь, потому что не каждый день вопрос Codegolf можно сделать во ФРАКТРАНЕ;)

намек

Код эквивалентен следующему псевдопифону:

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

Три метода показаны. Только первый подход соответствует требованию времени.

1-FindInstance (51 символ)

Это возвращает единственное решение уравнения.

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

Примеры и сроки

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

Это также работает, но слишком медленно, чтобы удовлетворить требования к скорости.

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

Range[0, Floor@Sqrt@n]^2является набором всех квадратов, меньших квадратного корня из n(максимально возможный квадрат в разбиении).

{4}требует целочисленных разбиений, nсостоящих из 4 элементов из вышеупомянутого набора квадратов.

1, внутри функции IntegerPartitionsвозвращает первое решение.

[[1]]удаляет внешние брекеты; решение было возвращено как набор из одного элемента.

f[123456]

{348, 44, 20, 4}

3-PowerRepresentations

PowerRepresentations возвращает все решения проблемы 4 квадратов. Это может также решить для сумм других полномочий.

PowersRepresentations возвращает менее чем за 5 секунд 181 способ выражения 123456789 в виде суммы четырех квадратов:

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

Тем не менее, это слишком медленно для других сумм.

Perl - 116 байт 87 байт (см. Обновление ниже)

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

Считая Шебанг одним байтом, добавлены новые строки для горизонтального рассудка.

Что-то вроде комбинации код-гольф, подача самого быстрого кода .

Кажется, что средняя (наихудшая?) Сложность случая составляет O (log n) O (n ^0,07 ) . Ничто из того, что я нашел, не работает медленнее, чем 0,001 с, и я проверил весь диапазон от 900000000 до 999999999 . Если вы обнаружите что-либо, что занимает значительно больше времени, ~ 0,1 с или более, пожалуйста, дайте мне знать.

Образец использования

$ 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

Последние два из них, как представляется, являются наихудшими сценариями для других представлений. В обоих случаях показанное решение является буквально самой первой проверенной вещью. Ведь 123456789это второй.

Если вы хотите проверить диапазон значений, вы можете использовать следующий скрипт:

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

Лучше всего, когда передано в файл. Диапазон 1..1000000занимает около 14 секунд на моем компьютере (71000 значений в секунду), а диапазон - 999000000..1000000000около 20 секунд (50000 значений в секунду), что соответствует средней сложности O (log n) .

Обновить

Редактировать : Оказывается, этот алгоритм очень похож на тот, который использовался ментальными калькуляторами в течение по крайней мере века .

С момента первоначальной публикации я проверил все значения в диапазоне от 1..1000000000 . Поведение «наихудшего случая» было продемонстрировано значением 699731569 , которое проверяло в общей сложности 190 комбинаций, прежде чем прийти к решению. Если вы считаете 190 малой константой - и я, конечно, знаю - поведение наихудшего случая в требуемом диапазоне можно считать O (1) . То есть, так же быстро, как поиск решения с гигантского стола, и в среднем, вполне возможно, быстрее.

Другое дело, хотя. После 190 итераций все, что больше 144400 , даже после первого прохода не вышло . Логика обхода в ширину ничего не стоит - она ​​даже не используется. Приведенный выше код может быть несколько сокращен:

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

Который выполняет только первый проход поиска. Нам нужно подтвердить, что нет значений ниже 144400, которые требовали бы второго прохода, хотя:

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

Короче говоря, для диапазона 1..1000000000 существует решение с почти постоянным временем, и вы смотрите на него.

Обновленное обновление

@Dennis и я сделали несколько улучшений этого алгоритма. Вы можете следить за прогрессом в комментариях ниже и последующим обсуждением, если это вас интересует. Среднее число итераций для требуемого диапазона уменьшилось с чуть более 4 до 1,229 , а время, необходимое для проверки всех значений для 1..1000000000 , было увеличено с 18 м 54 с до 2 м 41 с. В худшем случае ранее требовалось 190 итераций; в худшем случае, 854382778 , нужен только 21 .

Окончательный код Python выглядит следующим образом:

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

При этом используются две предварительно вычисленные таблицы поправок: одна размером 10 КБ, другая - 253 КБ. Приведенный выше код включает функции генератора для этих таблиц, хотя, вероятно, их следует вычислять во время компиляции.

Версию с таблицами коррекции более скромного размера можно найти здесь: http://codepad.org/1ebJC2OV. Эта версия требует в среднем 1,620 итераций в семестр, наихудший случай - 38 , а весь диапазон работает примерно за 3 м 21 с. Немного времени компенсируется, используя побитовое andдля коррекции b , а не по модулю.

улучшения

Четные значения с большей вероятностью дают решение, чем нечетные.
В статье о ментальных вычислениях, связанной с ранее сделанными заметками, отмечается, что если после удаления всех четырех факторов значение, которое нужно разложить, является четным, это значение можно разделить на два, и решение можно восстановить:

Хотя это может иметь смысл для умственного вычисления (меньшие значения, как правило, легче вычислить), это не имеет большого смысла алгоритмически. Если вы возьмете 256 случайных 4- кортежей и изучите сумму квадратов по модулю 8 , вы обнаружите, что значения 1 , 3 , 5 и 7 достигаются в среднем по 32 раза. Однако значения 2 и 6 достигаются по 48 раз. Умножив нечетные значения на 2, вы найдете решение, в среднем, на 33% меньше итераций. Реконструкция заключается в следующем:

Необходимо позаботиться о том, чтобы a и b имели одинаковое соотношение, равно как и c и d , но если решение вообще было найдено, правильное упорядочение гарантировано существует.

Невозможные пути не нужно проверять.
После выбора второго значения b уже может оказаться невозможным существование решения, учитывая возможные квадратичные вычеты для любого заданного модуля. Вместо того, чтобы все равно проверять или переходить к следующей итерации, значение b можно «исправить», уменьшив его на наименьшее значение, которое может привести к решению. В двух таблицах поправок хранятся эти значения: одно для b , а другое для c . Использование более высокого модуля (точнее, использования модуля с относительно меньшим количеством квадратичных остатков) приведет к лучшему улучшению. Значение a не нуждается в коррекции; изменив n на четное, все значенияА действительны.

Питон 3 (177)

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

После того, как мы уменьшим входное значение, Nчтобы оно не делилось на 4, оно должно быть выражено в виде суммы четырех квадратов, где один из них является либо наибольшим возможным значением, a=int(N**0.5)либо одним меньшим, чем, оставляя лишь небольшой остаток для суммы трех других квадратов. заботиться о. Это значительно уменьшает пространство поиска.

Вот доказательство позже, этот код всегда находит решение. Мы хотим найти aтак, что n-a^2это сумма трех квадратов. Из теоремы Лежандра о трех квадратах число представляет собой сумму трех квадратов, если это не форма 4^j(8*k+7). В частности, такими числами являются либо 0, либо 3 (по модулю 4).

Мы показываем, что никакие два последовательных элемента не aмогут заставить оставшуюся сумму N-a^2иметь такую ​​форму для обоих последовательных значений. Мы можем сделать это, просто составив таблицу aи Nпо модулю 4, отметив это, N%4!=0потому что мы извлекли все степени 4 из N.

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

Поскольку нет двух последовательных aуступок (N-a*a)%4 in [0,3], один из них является безопасным для использования. Итак, мы жадно используем максимально возможное nс n^2<=N, и n-1. Поскольку N<(n+1)^2остаток, который N-a^2должен быть представлен как сумма трех квадратов, равен не более (n+1)^2 -(n-1)^2, что равно 4*n. Таким образом, достаточно проверить только значения до 2*sqrt(n), что в точности соответствует диапазону R.

Можно было бы дополнительно оптимизировать время выполнения, останавливаясь после одного решения, вычисляя, а не повторяя последнее значение d, и выполняя поиск только среди значений b<=c<=d. Но даже без этих оптимизаций худший случай, который я смог найти, закончился за 45 секунд на моей машине.

Цепочка "для х в R" вызывает сожаление. Вероятно, его можно сократить с помощью подстановки или замены строки путем итерации по одному индексу, который кодирует (a, b, c, d). Импортировать itertools оказалось не стоит.

Изменить: изменено int(2*n**.5)+1с, 2*int(n**.5)+2чтобы сделать аргумент чище, то же количество символов.

CJam , 91 90 74 71 байт

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

Компактный, но медленнее, чем мой другой подход.

Попробуйте онлайн! Вставьте код , введите нужное целое число в поле ввода и нажмите « Выполнить» .

Фон

Этот пост начинался как 99-байтовый ответ GolfScript . Хотя возможности для улучшения еще есть, в GolfScript отсутствует встроенная функция sqrt. Я сохранил версию GolfScript до 5- й версии , так как она была очень похожа на версию CJam.

Однако для оптимизации, начиная с версии 6, требуются операторы, которых нет в GolfScript, поэтому вместо публикации отдельных объяснений для обоих языков я решил отказаться от менее конкурентоспособной (и гораздо более медленной) версии.

Реализованный алгоритм вычисляет числа грубой силой:

1. Для ввода mнайди Nи Wтакое что m = 4**W * N.

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

3. Set i = i + 1.

4. Найти целые числа j, k, lтакие, что i = 257**2 * j + 257 * k + l, где k, l < 257.

5. Проверьте, если d = N - j**2 - k**2 - l**2это идеальный квадрат.

6. Если это не так, вернитесь к шагу 3.

7. Версия для печати 2**W * j, 2**W * k, 2**W * l, 2**W * sqrt(m).

Примеры

$ 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

Сроки соответствуют Intel Core i7-4700MQ.

Как это устроено

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

С 228

Это основано на алгоритме на странице Википедии, который представляет собой грубую силу 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 байтов

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

Быстро, но долго. Новая строка может быть удалена.

Попробуйте онлайн. Обратите внимание, что онлайн-переводчик имеет ограничение в 5 секунд, поэтому он может работать не для всех номеров.

Фон

Алгоритм использует теорему Лежандра о трех квадратах , которая утверждает, что каждое натуральное число n , которое не имеет форму

может быть выражена как сумма трех квадратов.

Алгоритм выполняет следующие действия:

1. Выразите число как 4**i * j.

2. Найти наибольшее целое число kтакое, что k**2 <= jи j - k**2удовлетворяет гипотезе теоремы Лежандра о трех квадратах.

3. Set i = 0.

4. Проверьте, если j - k**2 - (i / 252)**2 - (i % 252)**2это идеальный квадрат.

5. Если это не так, увеличьте значение iи вернитесь к шагу 4.

Примеры

$ 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

Сроки соответствуют Intel Core i7-4700MQ.

Как это устроено

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

Ред. 1: С, 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();}

Это еще больше потребляет память, чем версия 0. Тот же принцип: создайте таблицу с положительным значением для всех возможных сумм из 2 квадратов (и ноль для тех чисел, которые не являются суммами из двух квадратов), а затем найдите ее.

В этой версии используйте массив shortвместо вместо того, charчтобы хранить попадания, чтобы я мог сохранить корень одного из пары квадратов в таблице, а не просто флаг. Это значительно упрощает функцию p(для декодирования суммы 2 квадратов), поскольку нет необходимости в цикле.

Windows имеет ограничение в 2 ГБ для массивов. Я могу обойти то, с short s[15<<26]чем является массив из 1006632960 элементов, достаточно, чтобы соответствовать спецификации. К сожалению, общий размер программы во время выполнения еще более 2 ГБ и (несмотря на настройки параметров ОС) Я не был в состоянии работать над этим размер (хотя теоретически это возможно.) Лучшее , что я могу сделать , это short s[14<<26](939524096 элементы.) m*mДолжны быть строго меньше этого (30651 ^ 2 = 939483801.) Тем не менее, программа работает отлично и должна работать на любой ОС, которая не имеет этого ограничения.

Код без правил

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

Это голодный зверь памяти. Он занимает массив 1 ГБ, вычисляет все возможные суммы из 2 квадратов и сохраняет флаг для каждого в массиве. Затем для пользовательского ввода z он ищет в массиве две суммы из 2 квадратов a и za.

функция pзатем reconsitutes первоначальные квадраты , которые были использованы , чтобы сделать суммы 2 квадратов aи z-aи выводит их, первую из каждой пары в виде целого числа, второй , как двойные (если она есть , чтобы быть все целыми числами , которые необходимы еще два символа, t> m=t.)

Программа занимает пару минут, чтобы построить таблицу сумм квадратов (я думаю, что это связано с проблемами управления памятью, я вижу, что выделение памяти увеличивается медленно, а не подпрыгивает, как можно было ожидать). Однако, как только это будет сделано выдает ответы очень быстро (если нужно рассчитать несколько чисел, программа scanfможет быть зациклена).

Негольфированный код

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)  

Пример вывода

Первый вопрос. Второй был выбран как трудный для поиска. В этом случае программа должна выполнить поиск до 8192 ^ 2 + 8192 ^ 2 = 134217728, но это займет всего несколько секунд после создания таблицы.

123456789
0 2.000000 3328 10601.000000

805306368
8192 8192.000000 8192 24576.000000

Mathematica, 138 символов

Таким образом, получается, что это приводит к отрицательным и мнимым результатам для определенных входных данных, как указано в edc65 (например, 805306368), так что это не является правильным решением. Я оставлю это пока, и, возможно, если я действительно ненавижу свое время, я вернусь и попытаюсь это исправить.

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

Или непокоренный

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

Я не слишком внимательно смотрел на алгоритмы, но я ожидаю, что это та же идея. Я просто придумал очевидное решение и подправил его, пока оно не сработало. Я проверил это для всех чисел от 1 до одного миллиарда, и ... это работает. Тест занимает около 100 секунд на моей машине.

Приятный момент в том, что, поскольку b, c и d определены с отложенными присваиваниями :=, их не нужно переопределять при уменьшении a. Это сохранило несколько дополнительных строк, которые у меня были раньше. Я мог бы поиграть в гольф дальше и вложить лишние части, но вот первый проект.

О, и вы запускаете его как, S@123456789и вы можете проверить это с помощью {S@#, Total[(S@#)^2]} & @ 123456789или # == Total[(S@#)^2]&[123456789]. Исчерпывающий тест

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

Раньше я использовал оператор Print [], но это сильно замедляло его, хотя он никогда не вызывался. Пойди разберись.

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]

Простая грубая сила по заранее рассчитанным квадратам.

Он нуждается в -Oопции компиляции (я добавил 3 символов для этого). Для худшего случая 999950883 требуется менее 1 минуты.

Только проверено на GHC.

C: 198 символов

Я, вероятно, могу сжать его до чуть более 100 символов. Что мне нравится в этом решении, так это минимальное количество мусора, просто цикл for, делающий то, что должен делать цикл for (что должно быть сумасшедшим).

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

И сильно претенциозно:

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

Изменить: Это не достаточно быстро для всех ввода, но я вернусь с другим решением. Я оставлю этот троичный беспорядок на месте.

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

Спасибо @Dennis за улучшения. Остальная часть ответа ниже не обновляется с версии А.

Ред. А: С, 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);}

Гораздо быстрее, чем мой другой ответ и с гораздо меньшим количеством памяти!

Для этого используется http://en.wikipedia.org/wiki/Legendre%27s_three-square_theorem . Любое число не следующей формы может быть выражено как сумма 3 квадратов (я называю это запрещенной формой):

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

Обратите внимание, что все нечетные квадратные числа имеют форму, (8b+1)а все четные квадратные числа имеют внешнюю форму 4b. Однако это скрывает тот факт, что все четные квадратные числа имеют форму 4^a*(odd square)==4^a*(8b+1). В результате 2^x-(any square number < 2^(x-1))для нечетного xвсегда будет иметь запрещенную форму. Следовательно, эти числа и их кратные являются сложными случаями, поэтому многие программы здесь делят степени 4 в качестве первого шага.

Как указано в ответе @ xnor, N-a*aнельзя использовать запрещенную форму для двух последовательных значений a. Ниже я представлю упрощенную форму своей таблицы. В дополнение к тому факту, что после деления на 4 N%4не может быть равен 0, обратите внимание, что есть только 2 возможных значения для (a*a)%4.

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

Итак, мы хотим избежать значений, (N-a*a)которые могут иметь запрещенную форму, а именно те, где (N-a*a)%43 или 0. Как можно видеть, это не может происходить для одного и того же Nс нечетным и четным (a*a).

Итак, мой алгоритм работает так:

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

Мне особенно нравится способ, которым я делаю шаг 3. Я добавляю 3 к N, так что декремент необходим, если (3+N-a*a)%4 =3 или 2. (но не 1 или 0.) Разделите это на 2, и вся работа может быть выполнена довольно простым выражением ,

Код без правил

Обратите внимание на один forцикл qи использование деления / по модулю для получения значений bи cиз него. Я пытался использовать aв качестве делителя вместо 256 для сохранения байтов, но иногда значение aбыло неправильным, и программа зависала, возможно, на неопределенный срок. 256 был лучшим компромиссом, который я могу использовать >>8вместо /256деления.

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

Выход

Интересно, что если вы введете квадратное число, N-(a*a)= 0. Но программа обнаруживает, что 0%4= 0 и уменьшается до следующего квадрата вниз. В результате входные данные квадрата всегда разбиваются на группу меньших квадратов, если они не имеют форму 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 символов

Грубая сила, но быстро.

Редактировать быстро, но недостаточно для какого-то неприятного ввода. Я должен был добавить первый шаг сокращения умножением на 4.

Редактировать 2 Избавиться от функции, вывод внутри цикла, затем принудительно завершить работу

Редактировать 3 0 неверный ввод

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

Пример вывода

123456789
11111,48,10,8

805306368
16384,16384,16384,0