Персональна сторінка вчителя математики Ірини Кирдей

Як відомо, будь-яке просте числу виду 4k+1 можна представити у вигляді суми двох квадратів, до того одним єдиним способом.

N = P^2 + Q^2

Нижче наведено конструкцію Лежандра для пошуку P і Q для заданого N.

Конструкція Лежандра передбачає розклад у ланцюговий дріб квадратний корінь з числа N.
Цей розклад маэ певний цикл m, після яого члени послідовності почнуть повторюватись в оберненому порядку.

Отже нехай
P(0) = 0; Q(0) = 1
A(i) = [(sqrt(N)+P(i))/Q(i)];
P(i) = -P(i-1)+A(i-1)*Q(i-1);
Q(i) = (N - P(i)^2) / Q(i-1);
де [] - округлення вниз.

Обчислення продовжуються, доки не отримаємо Q(j)=1, тоді елементи P и Q з індексами (j+1)/2 - шукані.
Наприклад для числа 29 період дорівнює 5: 2,1,1,2,10
А от для числа 1297540301, наприклад, період виявився рівним 14023.
При цьому
P(7012) = 35870
Q(7012) = 3299
35870^2 + 3299^2 = 1297540301

Програма на Паскалі для розкладу довільного просто N виду 4k+1 на суму двох квадратів:

program legendre;

{$APPTYPE CONSOLE}

uses
SysUtils,
IniFiles,
Classes,
Math,
StrUtils;

procedure WriteElapsed(td: TDateTime);
var s: string;
begin
s := DateTimeToStr(td);
writeln(StringOfChar(chr(8), Length(s)), s);
end;

procedure DecomposePrime(Prime: Int64);
var
Found : boolean;
P, Q, A, i: int64;
sqrtN: Extended;
Elapsed: TDateTime;
begin
Elapsed := Now;
WriteElapsed(Elapsed);
Found := False;
i := 0;
P := 0;
Q := 1;
sqrtN := Sqrt(Prime + 0.0);
while not Found do begin
A := Trunc((sqrtN + P)/Q);
P := -P + A * Q;
Q := Round((Prime - P*P) / Q);
i := i + 1;
if P*P + Q*Q = Prime then begin
Writeln(IntToStr(Prime) + ' = ' + IntToStr(P) + '^2 + ' + IntToStr(Q) + '^2');
Elapsed := Now;
WriteElapsed(Elapsed);
Writeln(IntToStr(i) + ' loops');
Found := true;
end;
end;
end;

begin
DecomposePrime(StrToInt64(ParamStr(1)));
end.

legendre.exe 1844674407370954349
06.02.2013 1:43:39
1844674407370954349 = 401327318^2 + 1297540285^2
06.02.2013 1:43:41
45187956 loops

Час виконання 2 секунди

Інтрига полягає у тому, що 2 секунди - це не межа можливостей.
Існує конструкція, яка дозволяє отримати цей розклад майже миттєво...
Якщо цікаво - пишіть мені на @ukr.net