Ответ на предновогоднюю задачку

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

Задача действительно интересная.
На данный момент мне прислали два её решения, причем у каждого оказался свой индивидуальный подход.

Примерное время на решение, около 3-4 часов (один вечер).
Именно столько было затрачено каждым из решивших задачу программистов, включая меня.

Ну и ответ на введенное максимальное пороговое число (10 в 12 степени) будет: $259814D6C9AAF914221E (это вам для самопроверки).

Пора перейти к самой сути.

0. Математика длинных чисел:

Вспомним саму задачу:
Пользователь вводит некое положительное число N в диапазоне от единицы до 10 в 12 степени включительно. Нужно посчитать сумму остатков от деления по модулю в диапазоне от 1 до N (N mod 1 + N mod 2 + N mod 3 +...+ N mod N).

Даже в самом первом приближении результат такой суммы не получится сохранить ни в одном из известных Delphi типов, включая UInt64.

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

Проблема в том, что работать это будет очень долго (примерно 3 дня при вводе максимального числа). Да и в действительности такое суммирование можно сделать гораздо проще:

type
  UInt128 = record
    Lo, Hi: UInt64;
  end;

var
  VeriBigInteger: UInt128;

procedure Adder(X: Int64);
var
  Tmp: UInt128;
begin
  Tmp := VeriBigInteger;
  Inc(VeriBigInteger.Lo, X);
  if (VeriBigInteger.Lo < Tmp.Lo) or (VeriBigInteger.Lo < X) then
    Inc(VeriBigInteger.Hi);
end;

Идея, описанная данном алгоритме простая, попробую ее описать для начала на переменных типа Byte.

Давайте напишем вот такой код:

program ByteOverflow;

{$APPTYPE CONSOLE}

var
  A: Byte;
begin
  A := 255;
  Inc(A);
  Writeln(A);
end.

Нам выведется число ноль, хотя по факту реальное значение суммы равно 256.
Так получилось из-за того что переменная типа Byte просто не может в себе разместить столь большое число (255 максимум).
Но, мы все-же сможем выполнить такое сложение и учесть возможность переполнения.
Так как мы знаем что проводим только суммирование, можно сделать вывод: как только значение переменной А стало меньше значения, хранившегося в ней ранее, или стало меньше числа, с которым производилось суммирование - произошла операция переполнения.

Что происходит при переполнении?
В более старший разряд (который мы не можем хранить) добавляется один бит, и вот этот бит нам и нужно учесть (256 = 1 0000 0000 в двоичной системе счисления)

Пишем такой код для проверки:

program ByteOverflow2;

{$APPTYPE CONSOLE}

uses
  SysUtils;

type
  MyUInt16 = record
    Lo, Hi: Byte;
  end;

var
  Z: MyUInt16;

procedure PrintMyUInt16;
begin
  Write('Result: 0x');
  Write(IntToHex(Z.Hi, 2));
  Write(IntToHex(Z.Lo, 2));
  Writeln;
end;

procedure Adder(X: Byte);
var
  Tmp: MyUInt16;
begin
  Tmp := Z;
  Inc(Z.Lo, X);
  // проверяем, не произошло ли переполнения
  if (Z.Lo < Tmp.Lo) or (Z.Lo < X) then
    // если произошло, увеличиваем старший байт
    Inc(Z.Hi);
end;

var
  A: Byte;
begin
  Adder(255);
  Adder(1);
  PrintMyUInt16;
  Readln;
end.

Что выведет нам в результате число $100 (256 в десятичной, что и ожидалось).

Почему к старшему байту прибавляется только единица.
Ну тут все просто: нет таких двух чисел одинаковой размерности, сумма которых вызвала бы увеличение старшего байта на 2 и больше.
К примеру 255 + 255, как максимальные значения байта, дадут в итоге всего лишь 510, что в двоичном представлении будет: 1 1111 1110.

Теперь, если этот момент понятен, обратите внимание на самый первый вариант кода - практически полная идентичность, за исключением того, что в этом случае мы работаем уже не с Byte, а с UInt64, все остальное по такому же принципу.

1. Определение порога арифметических прогрессий.

Теперь надо разобраться, что из себя представляют все эти числа, которые нам придется суммировать.

Строим график от 1 до N, где N = 1000.

Опираясь на график видим, что, по сути, мы имеем дело с набором арифметических прогрессий, причем самая правая (с конца) идет от нуля (нижний предел прогрессии pMinLimit) с шагом pIncrement (изначально равным единице) ровно до середины максимального значения N - 1 (верхний предел прогрессии pMaxLimit).

Х = 1000, 999, 998, 997...503, 502, 501
Y = 0, 1, 2, 3...497, 498, 499

Как только Х достигает 500, начинается новая арифметическая прогрессия, но с шагом pIncrement увеличенным на единицу, причем максимальный порог рассчитывается уже не до половины N - 1, а до третьей части от N минус какое-то число (назовем его пока UnknownValue).

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

Начнем с вычисления верхнего предела прогрессии для каждого шага и посмотрим на значения вершин:

499, 332, 247, 196, 165, 142, 118, 104, 91, 90, 76 и т.д.

Вот тут я их выделил:

Изначально я предположил что эти числа получаются по следующей формуле:

pMaxLimit := N div (pIncrement + 1) - pIncrement;

Т.е. на первом шаге pIncrement равен единице, делим N на увеличенное значение шага (на два в первом случае) и от результат вычитаем само значение шага.

Это работает для следующих чисел:
pIncrement = 1, результат (1000 / 2 - 1) = 499
pIncrement = 3, результат (1000 / 4 - 3) = 247
pIncrement = 4, результат (1000 / 5 - 4) = 196
pIncrement = 7, результат (1000 / 8 - 7) = 118
pIncrement = 9, результат (1000 / 10 - 9) = 91

Но не работает на следующих:
pIncrement = 2, результат (1000 / 3 - 2) = 331, а должно быть 332, разница в 1
pIncrement = 5, результат (1000 / 6 - 5) = 161, а должно быть 165, разница в 4
pIncrement = 6, результат (1000 / 7 - 6) = 136, а должно быть 142, разница в 6
pIncrement = 8, результат (1000 / 9 - 8) = 103, а должно быть 104, разница в 1
pIncrement = 10, результат (1000 / 11 - 10) = 80, а должно быть 90, разница в 10
pIncrement = 11, результат (1000 / 12 - 11) = 72, а должно быть 76, разница в 4

и т.д.

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

pIncrement = 2, Difference = (1000 mod 3) = 1
pIncrement = 5, Difference = (1000 mod 6) = 4
pIncrement = 6, Difference = (1000 mod 7) = 6
pIncrement = 8, Difference = (1000 mod 9) = 1
pIncrement = 10, Difference = (1000 mod 11) = 10
pIncrement = 11, Difference = (1000 mod 12) = 4

Ну а для тех значений pIncrement, для которых поправок не требуется (1, 3, 4, 7, 9) результат деления по модулю всегда равен нулю.

Таким образом финальная формула получения максимального значения арифметической прогрессии стала выглядеть вот так:

  pMaxLimit := Limit div (pIncrement + 1) -
    (pIncrement - (Limit mod (pIncrement + 1)));

Нижний лимит прогрессии определяется банально:

pMinLimit := Limit mod nCount;

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

function summ(a1, an: Int64; n: Integer): Int64;
begin
  // если количество элементов ряда равно одному, возвращаем значение как есть
  if n = 1 then
    Result := a1
  else
    // в противном случае считаем сумму арифметической прогрессии
    Result := (n * (a1 + an)) shr 1;
end;

2. Делим прогрессию на блоки

Теперь надо разобраться как правильно суммировать результаты сумм прогрессий.
Дело в том что для числа 10 в 12 степени сумма самой правой прогрессии (с шагом 1) не влезет в диапазон Int64 и будет равна следующему числу: 0x1A78 4379D963 7EF6BC00.

Для этого определимся сколько всего остатков деления по модулю влезет в Int64 без переполнения, и напишем такой тест:

procedure Tst;
var
  I, A, Z: Int64;
begin
  I := 10000000000000 div 2 - 1;
  A := 0;
  Z := 0;
  while I > 0 do
  begin
    Inc(Z);
    A := A + I;
    if A < 0 then
      Writeln(Z);
    Dec(I);
  end;
end;

Зная что работать мы будем с рядами прогрессий, где самая большая будет самой правой, с шагом pIncrement = 1, просто ручками просуммируем значение всего ряда, начиная от самых старших значений.
Как только значение А станет отрицательным, нам выведется число: 1844675.
Подстрахуемся и будем считать что количество таких элементов в два раза меньше и возьмем размер буфера в 900 тысяч элементов.

Зная что сумма ряда равна сумме частей ряда мы можем дробить длинные ряды на порции по 900 тысяч элементов и считать сумму от них, не выходя за переполнение.

К примеру ряд (1 + 2 + 3 + 4) = 10, если мы раздробим его на 2 части (1 + 2) и (3 + 4), то их суммы так-же совпадут (что логично :).

Осталось написать код.

3. Финальное решение задачи.

program mod_sum;

{$APPTYPE CONSOLE}

{$R *.res}

uses
  Windows,
  SysUtils;

type
  UInt128 = record
    Lo, Hi: UInt64;
  end;

var
  VeryBigInteger: UInt128;

procedure Adder(X: Int64);
var
  Tmp: UInt128;
begin
  Tmp := VeryBigInteger;
  Inc(VeryBigInteger.Lo, X);
  if (VeryBigInteger.Lo < Tmp.Lo) or (VeryBigInteger.Lo < X) then
    Inc(VeryBigInteger.Hi);
end;

function summ(a1, an: UInt64; n: Integer): UInt64;
begin
  // если количество элементов ряда равно одному, возвращаем значение как есть
  if n = 1 then
    Result := a1
  else
    // в противном случае считаем сумму арифметической прогрессии
    Result := (n * (a1 + an)) shr 1;
end;

procedure Calc(Limit: Int64);
const
  MaxProgressionSize = 900000;
var
  nCount,             // количество не обработанных элементов
  pMinLimit,          // минимальное значение прогрессии
  pMaxLimit,          // максимальное значение прогрессии
  pIncrement,         // дистанция между элементами ряда
  pCount,             // количество элементов ряда
  pSumm,              // сумма элементов ряда
  TmpMaxLimit: Int64; // временное максимальное значение прогрессии
  Start: DWORD;
begin
  // первичная инициализация
  Start := GetTickCount;
  ZeroMemory(@VeryBigInteger, SizeOf(UInt128));
  nCount := Limit;
  pIncrement := 1;
  pMinLimit := 0;
  pMaxLimit := Limit div (pIncrement + 1) -
    (pIncrement - (Limit mod (pIncrement + 1)));
  while nCount > 0 do
  begin
    // рассчитываем количество элементов ряда
    pCount := 1 + (pMaxLimit - pMinLimit) div pIncrement;
    // если это число больше чем максимальное кол-во элементов,
    // у которых сумма остатков от деления по модулю влезет в Int64
    // то считаем сумму ряда частями
    while pCount > MaxProgressionSize do
    begin
      // рассчитываем временный верхний предел ряда
      TmpMaxLimit := pMinLimit + (pIncrement * (MaxProgressionSize - 1));
      // считаем сумму части
      pSumm := summ(pMinLimit, TmpMaxLimit, MaxProgressionSize);
      // суммируем
      Adder(pSumm);
      // правим кол-во оставшихся элементов ряда
      Dec(pCount, MaxProgressionSize);
      // правим кол-во всех оставшихся элементов
      Dec(nCount, MaxProgressionSize);
      // сдвигаем нижний предел ряда
      pMinLimit := TmpMaxLimit + pIncrement;
    end;

    // считаем сумму ряда
    pSumm := summ(pMinLimit, pMaxLimit, pCount);
    // суммируем
    Adder(pSumm);
    // правим кол-во всех оставшихся элементов
    Dec(nCount, pCount);
    // если все посчитали - выходим
    if nCount = 0 then
      Break;
    // рассчитываем новый нижний предел ряда
    pMinLimit := Limit mod nCount;
    // рассчитываем новый максимальный предел ряда
    Inc(pIncrement);
    pMaxLimit := Limit div (pIncrement + 1) -
      (pIncrement - (Limit mod (pIncrement + 1)));
    // который, к слову сказать, может стать отрицательным
    if pMaxLimit <= 0 then
      // в таком случае выставляем его в ноль,
      // а pCount на след итерации станет равным единице (что есть правильно)
      pMaxLimit := 0;
  end;
  // ну и выводим результат
  Writeln('Result: 0x', IntToHex(VeryBigInteger.Hi, 16), IntToHex(VeryBigInteger.Lo, 16));
  Writeln('Time elapsed: ', GetTickCount - Start);
end;

var
  N: int64;
begin
  try
    repeat
      Writeln('Input N or zero (0) for exit: ');
      Readln(N);
      if N > 1000000000000 then
      begin
        Writeln('Too much value, repeat input.');
        Continue;
      end;
      if N < 0 then
      begin
        Writeln('Wrong value, repeat input.');
        Continue;
      end;
      if N <> 0 then
      begin
        Writeln('calculate...');
        Calc(N);
      end;
    until N = 0;
  except
    on E: Exception do
      Writeln(E.ClassName, ': ', E.Message);
  end;
end.

Выглядит в итоге вот так:

В отличие от самого первого варианта решения, это работает всего за 2 секунды, но...

4. Можно ли оптимизировать?

Конечно можно.
Это решение от Бориса Новгородова (MBo):

program mod_sum2;

{$APPTYPE CONSOLE}

{$R *.res}

uses
  Windows,
  SysUtils;

type
  TUInt128 = record
    case Integer of
      0: (Lo, Hi: UInt64);
      1: (LoLo, LoHi, HiLo, HiHi: Cardinal);
      2: (Bytes: array [0 .. 15] of Byte);
  end;

  procedure Add128(var Sum, Addend: TUInt128);
  begin
    if Sum.Lo >= UInt64($FFFFFFFFFFFFFFFF) - Addend.Lo then
      Inc(Sum.Hi);
    Sum.Lo := Sum.Lo + Addend.Lo;
    Sum.Hi := Sum.Hi + Addend.Hi;
  end;

  procedure Sub128(var Sum, Addend: TUInt128);
  begin
    if Sum.Lo < Addend.Lo then
      Dec(Sum.Hi);
    Sum.Lo := Sum.Lo - Addend.Lo;
    Sum.Hi := Sum.Hi - Addend.Hi;
  end;

  //Sq = Sqr(Value.Lo)
  procedure Sqr128(var Sq, Value: TUInt128);
  var
    ad: TUInt128;
  begin
    Sq.Hi := UInt64(Value.LoHi) * Value.LoHi;
    Sq.Lo := UInt64(Value.LoLo) * Value.LoLo;
    ad.Lo := UInt64(Value.LoHi) * Value.LoLo;
    ad.Hi := ad.LoHi;
    ad.Lo := ad.Lo shl 32;
    Add128(Sq, ad);
    Add128(Sq, ad);
  end;

  function SumOfMods128(N: Int64): TUInt128;
  var
    Minus, Divider: UInt64;
    Even: Boolean;
    V, Sq: TUInt128;
  begin
    Result.Lo := 0;
    Result.Hi := 0;
    Minus := 1;
    Divider := 2;
    Even := True;
    V.Hi := 0;
    while Minus < N do begin
      V.Lo := (N - Minus) div Divider;
      Sqr128(Sq, V);
      if Even then
        Add128(Result, Sq)
      else
        Sub128(Result, Sq);
      Minus := Minus + Divider;
      Divider := Divider + 1;
      Even := not Even;
    end;
  end;

var
  N: Int64;
  Res128: TUInt128;
  Start: DWORD;
begin
  Start := GetTickCount;
  N := 1000000000000;
  Res128 := SumOfMods128(N);
  Writeln('Result: 0x', IntToHex(Res128.Hi, 16), IntToHex(Res128.Lo, 16));
  Writeln('Time elapsed: ', GetTickCount - Start);
  Readln;
end.

Работает почти в 3 раза быстрее моего варианта.
Как работает и почему - я не объясню, я не математик :)

Единственно дам его пояснение:

Вот начало ряда сумм S(N) и модули, из которых эти суммы складываются.

N   S     mods
1   0      0
2   0      0  0
3   1      0  1  0
4   1      0  0  1  0
5   4      0  1  2  1  0
6   3      0  0  0  2  1  0
7   8      0  1  1  3  2  1  0
8   8      0  0  2  0  3  2  1  0
9  12      0  1  0  1  4  3  2  1  0
10 13      0  0  1  2  0  4  3  2  1  0
11 22      0  1  2  3  1  5  4  3  2  1  0
12 17      0  0  0  0  2  0  5  4  3  2  1  0
13 28      0  1  1  1  3  1  6  5  4  3  2  1  0

Видно, что правая часть представляет собой арифм. последовательность, сумму которой я уже приводил, а вот левая похитрее. Компактного математического выражения не нашел, но эмпирическое исследование привело к следующему:
главный член суммы резко растет на каждом втором (нечетном) числе, рост квадратичный, и это квадрат выражения (N-1) div 2
Q1 = Sqr((N-1) div 2)
Теперь, если вывести разницу Q1 - S(N), то можно увидеть, что она чаще всего скачет на каждом третьем числе, а сдвиг от начала будет уже 3
Q2 = -Sqr((N-3) div 3)
Далее аналогично
Q3 = Sqr((N-6) div 4)
Q4 = -Sqr((N-10) div 5)
и т.д. Числа сдвига - "треугольные", вида k(k-1)/2, где k - знаменатель.

Итого с использованием длинной арифметики получается функция, которая отличается от такой же для короткой арифметики только типом результата и заменой Sqr на TUint128.Sqr128 (если Multiply реализовать, то в обоих случаях можно одинаково X*X)

function SumOfMods128(N: Int64): TUInt128;
var
  Minus, Divider: UInt64;
  Even: Boolean;
begin
  Result := 0;
  Minus := 1;
  Divider := 2;
  Even := True;
  while Minus < N do begin
    if Even then
      Result := Result + TUint128.Sqr128((N - Minus) div Divider)
    else
      Result := Result - TUint128.Sqr128((N - Minus) div Divider);
    Minus := Minus + Divider;
    Divider := Divider + 1;
    Even := not Even;
  end;
end;

Прониклись? Тогда идем дальше...

5. А можно ли еще ускорить?

Конечно можно, и вот вариант решения от Александра Шарахова, более известного как Sha.

Если не знаете кто это такой, то откройте system.pas и найдите там вот такой блок текста:

(* ***** BEGIN LICENSE BLOCK *****
 *
 * The function PosEx is licensed under the CodeGear license terms.
 *
 * The initial developer of the original code is Fastcode
 *
 * Portions created by the initial developer are Copyright (C) 2002-2004
 * the initial developer. All Rights Reserved.
 *
 * Contributor(s): Aleksandr Sharahov
 *
 * ***** END LICENSE BLOCK ***** *)

Выглядит его решение вот так:

program mod_sum3;

{$APPTYPE CONSOLE}

{$R *.res}

uses
  Windows,
  SysUtils;

type
  TSum= array[0..3] of int64;

procedure Zero128(var sum: TSum);
begin;
  sum[0]:=0;
  sum[1]:=0;
  sum[2]:=0;
  sum[3]:=0;
end;

procedure Add128Progression(var sum: TSum; val, count: int64); //sum:=sum + val * count shr 1;
var
  t: int64;
begin;
  if val and 1=0 then val:=val shr 1 else count:=count shr 1;
  t:=int64(Int64Rec(val).Lo) * Int64Rec(count).Lo;
  sum[0]:=sum[0] + Int64Rec(t).Lo;
  sum[1]:=sum[1] + Int64Rec(t).Hi;
  t:=int64(Int64Rec(val).Lo) * Int64Rec(count).Hi;
  sum[1]:=sum[1] + Int64Rec(t).Lo;
  sum[2]:=sum[2] + Int64Rec(t).Hi;
  t:=int64(Int64Rec(val).Hi) * Int64Rec(count).Lo;
  sum[1]:=sum[1] + Int64Rec(t).Lo;
  sum[2]:=sum[2] + Int64Rec(t).Hi;
  t:=int64(Int64Rec(val).Hi) * Int64Rec(count).Hi;
  sum[2]:=sum[2] + Int64Rec(t).Lo;
  sum[3]:=sum[3] + Int64Rec(t).Hi;
end;

procedure Norm128(var sum: TSum);
begin;
  sum[1]:=sum[1] + Int64Rec(sum[0]).Hi; Int64Rec(sum[0]).Hi:=0;
  sum[2]:=sum[2] + Int64Rec(sum[1]).Hi; Int64Rec(sum[1]).Hi:=0;
  sum[3]:=sum[3] + Int64Rec(sum[2]).Hi; Int64Rec(sum[2]).Hi:=0;
                                        Int64Rec(sum[3]).Hi:=0;
end;

procedure Add128(var sum: TSum; val: int64); //sum:=sum + val;
begin;
  sum[0]:=sum[0] + Int64Rec(val).Lo;
  sum[1]:=sum[1] + Int64Rec(val).Hi;
end;

procedure Add128Norm(var sum: TSum; val: int64); //sum:=norm(sum + val);
begin;
  Add128(sum, val);
  Norm128(sum);
end;

procedure Not128Norm(var sum: TSum); //sum:=norm(not sum);
begin;
  Norm128(sum);
  Int64Rec(sum[0]).Lo:=not Int64Rec(sum[0]).Lo;
  Int64Rec(sum[1]).Lo:=not Int64Rec(sum[1]).Lo;
  Int64Rec(sum[2]).Lo:=not Int64Rec(sum[2]).Lo;
  Int64Rec(sum[3]).Lo:=not Int64Rec(sum[3]).Lo;
end;

procedure Neg128(var sum: TSum); //sum:=-sum;
begin;
  Not128Norm(sum);
  Add128(sum, 1);
end;

procedure Sum128(var sum: TSum; n: int64);
var
  m, c: cardinal;
  i, t: int64;
begin;
  Zero128(sum); t:=0;
  c:=trunc(sqrt(n+0.0));
  m:=c+1;
  dec(c, ord(int64(c)*c=n));
  while c>1 do begin;
    i:=n div c;
    t:=t + (cardinal(n)-cardinal(i)*c);
    Add128Progression(sum, i+m, i-m+1);
    dec(c);
    end;
  Neg128(sum);
  Add128Progression(sum, n-m, n-m+1);
  Add128Norm(sum,t);
  end;

var
  m, n: int64;
  sum: TSum;
  t: DWORD;
begin
  t := GetTickCount;
  N := 1000000000000;
  Sum128(sum, n);
  t:=GetTickCount-t;
  Writeln(Format('%dms   %.8x %.8x %.8x %.8x',[t, sum[3], sum[2], sum[1], sum[0]]));
  Readln;
end.

Обоснование:

Делим слагаемые на 2 группы: [1 .. trunc(sqrt(N))] и [trunc(sqrt[N))+1 .. N].

В первой группе слагаемых мало, складываем их все непосредственно.
Заметим, что и их количество, и каждое слагаемое не превышает sqrt(N).
Значит, их сумма не превысит N.

Вторая группа разбивается на sqrt(N) непересекающихся подгрупп. В первую войдут все номера больше N div 2. Во вторую – больше N div 3, не вошедшие в первую подгруппу. В третью – большие N div 4, не вошедшие в первые две, и т.д.

Это те самые пики, но считать их неудобно.
Удобнее считать сумму ПЕРЕСЕКАЮШИХСЯ подгрупп. Первая – это вся правая часть от trunc(sqrt(N))+1 до N, вторая – до N div 2 включительно, третья – до N div 3 включительно.

В этом случае мы считаем общую сумму как прогрессию 0+1+2+…, а потом просто отнимаем превышения результата для каждой пересекающейся подгруппы. Так, например, для всех чисел от trunc(sqrt(N))+1 до N div 2 (это первая итерация) результат будет завышен по крайней мере (k+m) * (k-m+1) shr 1, где(k+m) - сумма первого и последнего элементов прогрессии, (k-m+1) - количество элементов в этой подгруппе. И т.д.

Ну а работает этот код примерно в районе 73 (!!!) миллисекунд :)
Единственное, что могу сказать по поводу этого варианта решения задачи:

На сем откланиваюсь, а исходники всех трех примеров можно забрать тут.
Всем удачного Нового Года

---

© Александр (Rouse_) Багель

Январь, 2015