Skip to content

Instantly share code, notes, and snippets.

@jin-x
Last active June 6, 2021 09:57
Show Gist options
  • Save jin-x/a6b97330d027471d46556d6ead42b7bf to your computer and use it in GitHub Desktop.
Save jin-x/a6b97330d027471d46556d6ead42b7bf to your computer and use it in GitHub Desktop.
@jinxonik / UniLecs #100 (Delphi fast version)
program VegCutFast;
// Решение задачи "Овощная нарезка" (https://t.me/unilecs/359)
// Быстрая (умная) версия, использующая многопоточность на больших числах
{$APPTYPE CONSOLE}
// Рассуждения и описание алгоритма:
// http://telegra.ph/Reshenie-zadachi-ovoshchnaya-narezka-06-07
// Модуль Velthuis.BigIntegers: https://github.com/rvelthuis/BigNumbers
// Описание модуля: http://www.rvelthuis.de/programs/bigintegers.html
//{$DEFINE DEBUG}//{$IFDEF DEBUG}{$DEFINE DEBUG2}{$ENDIF}
uses
Winapi.Windows, System.SysUtils, System.Math, System.Threading,
Velthuis.BigIntegers;
// Функция Мёбиуса
function Mu(d: Integer): Integer;
var M, i, i2: Integer;
begin
M := 0;
i := 2;
repeat
i2 := i*i;
if i2 > d then Break;
if d mod i2 = 0 then Exit(0);
while d mod i = 0 do
begin
Inc(M);
d := d div i;
end;
Inc(i);
until False;
if d > 1 then Inc(M);
if Odd(M) then Result := -1
else Result := 1;
end;
// Получение числа неприводимых многочленов степени k
function Ik(k: Integer): BigInteger;
var MaxI, i: Integer;
begin
Result := 0;
MaxI := Trunc(Sqrt(k));
for i := 1 to MaxI do
if k mod i = 0 then
begin
case Mu(i) of
1: Result := Result + BigInteger.Pow(2, k div i);
-1: Result := Result - BigInteger.Pow(2, k div i);
end;
if (i <> MaxI) or (i*i <> k) then
case (Mu(k div i)) of
1: Result := Result + BigInteger.Pow(2, i);
-1: Result := Result - BigInteger.Pow(2, i);
end;
end;
Result := Result div k;
end;
// Определение количества вариантов нарезки для N элементов
// Однопоточная функция
function CountVegCut(N: Integer): BigInteger;
var MaxI, Count, i: Integer;
begin
Result := 0;
if N < 1 then Exit;
MaxI := Trunc(Sqrt(N));
Count := 0;
for i := 2 to MaxI do
if N mod i = 0 then
begin
Inc(Count);
Result := Result + Ik(i);
if (i <> MaxI) or (i*i <> N) then
Result := Result + Ik(N div i);
end;
if Count = 0 then Result := (BigInteger.Pow(2, N) - 2) div N + 2
else Result := Result + {Ik(1)=}2 + Ik(N);
end;
{$IFDEF DEBUG}
var CS: TRTLCriticalSection;
procedure SafeWriteLn(const S: String; Wait: Boolean = False);
begin
EnterCriticalSection(CS);
Write(S);
if Wait then ReadLn
else WriteLn;
LeaveCriticalSection(CS);
end;
{$ENDIF}
// Определение количества вариантов нарезки для N элементов
// Многопоточная функция
function CountVegCutMultithread(N: Integer): BigInteger;
var
DivList: array of Integer;
IkList: array of BigInteger;
Count, Factor{$IFDEF DEBUG}, Added{$ENDIF}: Integer;
MaxI, i: Integer;
begin
Result := 0;
if N < 1 then Exit;
{$IFDEF DEBUG}InitializeCriticalSection(CS);{$ENDIF}
// Получаем список делителей
{$IFDEF DEBUG}SafeWriteLn('> calc DivList...');{$ENDIF}
Count := 0;
SetLength(DivList, Count);
MaxI := Trunc(Sqrt(N));
for i := 1 to MaxI do
if N mod i = 0 then
begin
Inc(Count);
SetLength(DivList, Count);
DivList[Count-1] := i;
if (i <> MaxI) or (i*i <> N) then
begin
Inc(Count);
SetLength(DivList, Count);
DivList[Count-1] := N div i;
end;
end;
// Для простых чисел
if Count = 2 then
begin
{$IFDEF DEBUG}SafeWriteLn(Format('> Count = %d (prime!)'#13#10'> start to calc Result...', [Count]));{$ENDIF}
SetLength(DivList, 0);
Result := (BigInteger.Pow(2, N) - 2) div N + 2;
{$IFDEF DEBUG}SafeWriteLn(Format('>> result is ready!', [Added, Count]));{$ENDIF}
Exit;
end;
// Вычисляем Ik
{$IFDEF DEBUG}SafeWriteLn(Format('> Count = %d'#13#10'> start to calc Ik...', [Count]));{$ENDIF}
SetLength(IkList, Count);
TParallel.For(0, Count-1, procedure(j: Integer)
begin
{$IFDEF DEBUG}SafeWriteLn(Format('>> calc Ik(DivList[%d]=%d) starts', [j, DivList[j]]));{$ENDIF}
IkList[j] := Ik(DivList[j]);
{$IFDEF DEBUG2}SafeWriteLn(Format('>> calc Ik(DivList[%d]=%d) ends', [j, DivList[j]]));{$ENDIF}
end);
{$IFDEF DEBUG}SafeWriteLn('> press enter to start adding...', True); Added := 1;{$ENDIF}
// Складываем числа
Factor := 1;
repeat
Factor := Factor * 2;
MaxI := (Count+Factor div 2-1) div Factor-1;
if MaxI < 0 then Break;
{$IFDEF DEBUG}SafeWriteLn(Format('>> Factor=%d, i=0..%d', [Factor, MaxI]));{$ENDIF}
TParallel.For(0, MaxI, procedure(j: Integer)
var i1, i2: Integer;
begin
i1 := j*Factor;
i2 := i1+Factor div 2;
{$IFDEF DEBUG}SafeWriteLn(Format('>> adding DivList[%d]+[%d] starts', [i1, i2])); AtomicIncrement(Added);{$ENDIF}
IkList[i1] := IkList[i1] + IkList[i2];
{$IFDEF DEBUG2}SafeWriteLn(Format('>> adding DivList[%d]+[%d] ends', [i1, i2]));{$ENDIF}
end);
until False;
// Результат
{$IFDEF DEBUG}SafeWriteLn(Format('>> result is ready, added=%d of %d!', [Added, Count]));{$ENDIF}
Result := IkList[0];
SetLength(DivList, 0);
SetLength(IkList, 0);
{$IFDEF DEBUG}DeleteCriticalSection(CS);{$ENDIF}
end;
////////////////////////////////////////////////////////////////////////////////
// Преобразовать число в строку с разделением по 3 разряда
function IntToStrBy3(N: Integer; Sep: Char = ''''): String;
var Len, i: Integer;
begin
Result := IntToStr(N);
Len := Length(Result);
for i := 1 to (Len-1) div 3 do
Insert(Sep, Result, Len-i*3+1);
end;
const
MAX_DIGS = 2000; // Сколько цифр отображать
MAX_DIGS_ASK_HEX = 1*1000*1000; // Когда спрашивать о переводе в Hex
MAX_DIGS_ASK_CANCEL = 100*1000*1000; // Когда предлагать отмену вывода
MIN_MULTITHREAD = 5*1000*1000; // Минимальное значение N, при котором используется многопоточность
var
N: Integer;
Answer, Big, Reminder: BigInteger;
StrLen, StrLen2: Integer;
S: String;
Tick1, Tick2, TickFreq: Int64;
// Основной код
begin
try
// Ввод
WriteLn('[Vegetable Cuts] by Jin X (08.06.2018), fast version (smartthread) :)');
Write('Enter N (1..2147483647): ');
try
ReadLn(N);
except
WriteLn('Wrong number is specified!');
Exit;
end;
WriteLn('Calculating...');
try
QueryPerformanceFrequency(TickFreq);
QueryPerformanceCounter(Tick1);
if N < MIN_MULTITHREAD then Answer := CountVegCut(N)
else Answer := CountVegCutMultithread(N);
QueryPerformanceCounter(Tick2);
except
WriteLn('Oh no! Something goes wrong (exception is occured)... :(');
Exit;
end;
if Answer > 0 then StrLen := Ceil(BigInteger.Log10(Answer))
else StrLen := 1;
StrLen2 := StrLen;
if StrLen >= MAX_DIGS_ASK_HEX then
begin
Write('Number is too large! Output in hex, it''s faster (y/n)? [y]: ');
ReadLn(S);
if Copy(LowerCase(S), 1, 1) <> 'n' then
begin
BigInteger.Hex;
StrLen2 := Ceil(BigInteger.Log(Answer, BigInteger.Base));
end;
if StrLen >= MAX_DIGS_ASK_CANCEL then
begin
Write('OMG, it''s very very huge! Cancel output (y/n)? [y]: ');
ReadLn(S);
if Copy(LowerCase(S), 1, 1) <> 'n' then StrLen2 := -1
end
end;
if StrLen2 = -1 then S := ''
else
begin
WriteLn('Almost done...');
if StrLen2 >= MAX_DIGS then
begin
Big := BigInteger.Pow(BigInteger.Base, MAX_DIGS div 2-1);
BigInteger.DivMod(Answer, Big, Answer, Reminder);
S := (Answer div BigInteger.Pow(BigInteger.Base, StrLen2-(MAX_DIGS-3))).ToString;
S := S + '...' + Reminder.ToString;
end
else S := Answer.ToString;
end;
WriteLn;
if S <> '' then
begin
WriteLn(S);
WriteLn;
end;
WriteLn(IntToStrBy3(StrLen), ' digits, calculated in ', (Tick2-Tick1)/TickFreq:0:6, ' seconds');
finally
WriteLn('Press Enter to quit...');
ReadLn;
end;
end.
@jin-x
Copy link
Author

jin-x commented Jun 8, 2018

Download THIS FILE (save as link) and rename it to VegCut.zip.
You'll get more sources (including Python) and compiled EXE... :))

@jin-x
Copy link
Author

jin-x commented Jun 8, 2018

VerCut.pas

program VegCut;
// Решение задачи "Овощная нарезка" (https://t.me/unilecs/359)
// Однопоточная версия
{$APPTYPE CONSOLE}

// Рассуждения и описание алгоритма:
// http://telegra.ph/Reshenie-zadachi-ovoshchnaya-narezka-06-07

// Модуль Velthuis.BigIntegers: https://github.com/rvelthuis/BigNumbers
// Описание модуля: http://www.rvelthuis.de/programs/bigintegers.html

uses
  Winapi.Windows, System.SysUtils, System.Math, Velthuis.BigIntegers;

// Функция Мёбиуса
function Mu(d: Integer): Integer;
var M, i, i2: Integer;
begin
  M := 0;
  i := 2;
  repeat
    i2 := i*i;
    if i2 > d then Break;
    if d mod i2 = 0 then Exit(0);
    while d mod i = 0 do
    begin
      Inc(M);
      d := d div i;
    end;
    Inc(i);
  until False;
  if d > 1 then Inc(M);
  if Odd(M) then Result := -1
  else Result := 1;
end;

// Получение числа неприводимых многочленов степени k
function Ik(k: Integer): BigInteger;
var MaxI, i: Integer;
begin
  Result := 0;
  MaxI := Trunc(Sqrt(k));
  for i := 1 to MaxI do
    if k mod i = 0 then
    begin
      case Mu(i) of
        1: Result := Result + BigInteger.Pow(2, k div i);
        -1: Result := Result - BigInteger.Pow(2, k div i);
      end;
      if (i <> MaxI) or (i*i <> k) then
        case (Mu(k div i)) of
          1: Result := Result + BigInteger.Pow(2, i);
          -1: Result := Result - BigInteger.Pow(2, i);
        end;
    end;
  Result := Result div k;
end;

// Определение количества вариантов нарезки для N элементов
function CountVegCut(N: Integer): BigInteger;
var MaxI, Count, i: Integer;
begin
  Result := 0;
  if N < 1 then Exit;
  MaxI := Trunc(Sqrt(N));
  Count := 0;
  for i := 2 to MaxI do
    if N mod i = 0 then
    begin
      Inc(Count);
      Result := Result + Ik(i);
      if (i <> MaxI) or (i*i <> N) then
        Result := Result + Ik(N div i);
    end;
  if Count = 0 then Result := (BigInteger.Pow(2, N) - 2) div N + 2
  else Result := Result + {Ik(1)=}2 + Ik(N);
end;

////////////////////////////////////////////////////////////////////////////////

// Преобразовать число в строку с разделением по 3 разряда
function IntToStrBy3(N: Integer; Sep: Char = ''''): String;
var Len, i: Integer;
begin
  Result := IntToStr(N);
  Len := Length(Result);
  for i := 1 to (Len-1) div 3 do
    Insert(Sep, Result, Len-i*3+1);
end;

const
  MAX_DIGS = 2000;                      // Сколько цифр отображать
  MAX_DIGS_ASK_HEX = 1*1000*1000;       // Когда спрашивать о переводе в Hex
  MAX_DIGS_ASK_CANCEL = 100*1000*1000;  // Когда предлагать отмену вывода

var
  N: Integer;
  Answer, Big, Reminder: BigInteger;
  StrLen, StrLen2: Integer;
  S: String;
  Tick1, Tick2, TickFreq: Int64;

// Основной код
begin
  try
    // Ввод
    WriteLn('[Vegetable Cuts] by Jin X (08.06.2018)');
    Write('Enter N (1..2147483647): ');
    try
      ReadLn(N);
    except
      WriteLn('Wrong number is specified!');
      Exit;
    end;
  
    WriteLn('Calculating...');
    try
      QueryPerformanceFrequency(TickFreq);
      QueryPerformanceCounter(Tick1);
      Answer := CountVegCut(N);
      QueryPerformanceCounter(Tick2);
    except
      WriteLn('Oh no! Something goes wrong (exception is occured)... :(');
      Exit;
    end;
  
    if Answer > 0 then StrLen := Ceil(BigInteger.Log10(Answer))
    else StrLen := 1;
    StrLen2 := StrLen;
    if StrLen >= MAX_DIGS_ASK_HEX then
    begin
      Write('Number is too large! Output in hex, it''s faster (y/n)? [y]: ');
      ReadLn(S);
      if Copy(LowerCase(S), 1, 1) <> 'n' then
      begin
        BigInteger.Hex;
        StrLen2 := Ceil(BigInteger.Log(Answer, BigInteger.Base));
      end;
      if StrLen >= MAX_DIGS_ASK_CANCEL then
      begin
        Write('OMG, it''s very very huge! Cancel output (y/n)? [y]: ');
        ReadLn(S);
        if Copy(LowerCase(S), 1, 1) <> 'n' then StrLen2 := -1
      end
    end;
    if StrLen2 = -1 then S := ''
    else
    begin
      WriteLn('Almost done...');
      if StrLen2 >= MAX_DIGS then
      begin
        Big := BigInteger.Pow(BigInteger.Base, MAX_DIGS div 2-1);
        BigInteger.DivMod(Answer, Big, Answer, Reminder);
        S := (Answer div BigInteger.Pow(BigInteger.Base, StrLen2-(MAX_DIGS-3))).ToString;
        S := S + '...' + Reminder.ToString;
      end
      else S := Answer.ToString;
    end;
  
    WriteLn;
    if S <> '' then
    begin
      WriteLn(S);
      WriteLn;
    end;
    WriteLn(IntToStrBy3(StrLen), ' digits, calculated in ', (Tick2-Tick1)/TickFreq:0:6, ' seconds');
  finally
    WriteLn('Press Enter to quit...');
    ReadLn;
  end;
end.

@jin-x
Copy link
Author

jin-x commented Jun 12, 2018

Метод Пойа: VegCutPoyaFast.pas

program VegCutPoyaFast;
// Решение задачи "Овощная нарезка" (https://t.me/unilecs/359)
// Быстрая (умная) версия, использующая теорему Пойа
{$APPTYPE CONSOLE}

// Модуль Velthuis.BigIntegers: https://github.com/rvelthuis/BigNumbers
// Описание модуля: http://www.rvelthuis.de/programs/bigintegers.html

//{$DEFINE DEBUG}//{$IFDEF DEBUG}{$DEFINE DEBUG2}{$ENDIF}

uses
  Winapi.Windows, System.SysUtils, System.Math, System.Threading,
  Velthuis.BigIntegers;

// Функция Эйлера
function Euler(N: Integer): Integer;
var i: Integer;
begin
  Result := N;
  for i := 2 to Trunc(Sqrt(N)) do
    if N mod i = 0 then
    begin
      Result := Result - Result div i;
      while N mod i = 0 do
        N := N div i;
    end;
  if N > 1 then 
    Result := Result - Result div N;
end;

// Определение количества вариантов нарезки для N элементов
// Однопоточная функция
function CountVegCut(N: Integer): BigInteger;
var MaxI, i: Integer;
begin
  Result := 0;
  if N < 1 then Exit;
  MaxI := Trunc(Sqrt(N));
  for i := 1 to MaxI do
    if N mod i = 0 then
    begin
      Result := Result + BigInteger.Pow(2, i) * Euler(N div i);
      if (i <> MaxI) or (i*i <> N) then
        Result := Result + BigInteger.Pow(2, N div i) * Euler(i);
    end;
  Result := Result div N;
end;

{$IFDEF DEBUG}
var CS: TRTLCriticalSection;
procedure SafeWriteLn(const S: String; Wait: Boolean = False);
begin
  EnterCriticalSection(CS);
  Write(S);
  if Wait then ReadLn
  else WriteLn;
  LeaveCriticalSection(CS);
end;
{$ENDIF}

// Определение количества вариантов нарезки для N элементов
// Многопоточная функция
function CountVegCutMultithread(N: Integer): BigInteger;
var
  DivList: array of Integer;
  SumList: array of BigInteger;
  Count, Factor{$IFDEF DEBUG}, Added{$ENDIF}: Integer;  
  MaxI, i: Integer;
begin
  Result := 0;
  if N < 1 then Exit;
  {$IFDEF DEBUG}InitializeCriticalSection(CS);{$ENDIF}
  // Получаем список делителей
  {$IFDEF DEBUG}SafeWriteLn('> calc DivList...');{$ENDIF}
  Count := 0;
  SetLength(DivList, Count);
  MaxI := Trunc(Sqrt(N));
  for i := 1 to MaxI do
    if N mod i = 0 then
    begin
      Inc(Count);
      SetLength(DivList, Count);
      DivList[Count-1] := i;
      if (i <> MaxI) or (i*i <> N) then
      begin
        Inc(Count);
        SetLength(DivList, Count);
        DivList[Count-1] := N div i;
      end;
    end;
(*
  // Для простых чисел
  if Count = 2 then
  begin
    {$IFDEF DEBUG}SafeWriteLn(Format('> Count = %d (prime!)'#13#10'> start to calc Result...', [Count]));{$ENDIF}
    SetLength(DivList, 0);
    Result := (BigInteger.Pow(2, N) - 2) div N + 2;
  {$IFDEF DEBUG}SafeWriteLn(Format('>> result is ready!', [Added, Count]));{$ENDIF}
    Exit;
  end;
*)
  // Вычисляем значения с помощью функции Эйлера
  {$IFDEF DEBUG}SafeWriteLn(Format('> Count = %d'#13#10'> start to calc Ik...', [Count]));{$ENDIF}
  SetLength(SumList, Count);
  TParallel.For(0, Count-1, procedure(j: Integer)
    begin
      {$IFDEF DEBUG}SafeWriteLn(Format('>> calc Euler(DivList[%d]=%d) starts', [j, DivList[j]]));{$ENDIF}
      SumList[j] := BigInteger.Pow(2, DivList[j]) * Euler(N div DivList[j]);
      {$IFDEF DEBUG2}SafeWriteLn(Format('>> calc Euler(DivList[%d]=%d) ends', [j, DivList[j]]));{$ENDIF}
    end);
  {$IFDEF DEBUG}SafeWriteLn('> press enter to start adding...', True); Added := 1;{$ENDIF}
  // Складываем числа
  Factor := 1;
  repeat
    Factor := Factor * 2;
    MaxI := (Count+Factor div 2-1) div Factor-1;
    if MaxI < 0 then Break;
    {$IFDEF DEBUG}SafeWriteLn(Format('>> Factor=%d, i=0..%d', [Factor, MaxI]));{$ENDIF}
    TParallel.For(0, MaxI, procedure(j: Integer)
      var i1, i2: Integer;
      begin
        i1 := j*Factor;
        i2 := i1+Factor div 2;
        {$IFDEF DEBUG}SafeWriteLn(Format('>> adding DivList[%d]+[%d] starts', [i1, i2])); AtomicIncrement(Added);{$ENDIF}
        SumList[i1] := SumList[i1] + SumList[i2];
        {$IFDEF DEBUG2}SafeWriteLn(Format('>> adding DivList[%d]+[%d] ends', [i1, i2]));{$ENDIF}
      end);
  until False;
  // Результат
  {$IFDEF DEBUG}SafeWriteLn(Format('>> result is ready, added=%d of %d!', [Added, Count]));{$ENDIF}
  Result := SumList[0] div N;
  SetLength(DivList, 0);
  SetLength(SumList, 0);
  {$IFDEF DEBUG}DeleteCriticalSection(CS);{$ENDIF}
end;

////////////////////////////////////////////////////////////////////////////////

// Преобразовать число в строку с разделением по 3 разряда
function IntToStrBy3(N: Integer; Sep: Char = ''''): String;
var Len, i: Integer;
begin
  Result := IntToStr(N);
  Len := Length(Result);
  for i := 1 to (Len-1) div 3 do
    Insert(Sep, Result, Len-i*3+1);
end;

const
  MAX_DIGS = 2000;                      // Сколько цифр отображать
  MAX_DIGS_ASK_HEX = 1*1000*1000;       // Когда спрашивать о переводе в Hex
  MAX_DIGS_ASK_CANCEL = 100*1000*1000;  // Когда предлагать отмену вывода
  MIN_MULTITHREAD = 5*1000*1000;        // Минимальное значение N, при котором используется многопоточность

var
  N: Integer;
  Answer, Big, Reminder: BigInteger;
  StrLen, StrLen2: Integer;
  S: String;
  Tick1, Tick2, TickFreq: Int64;

// Основной код
begin
  try
    // Ввод
    WriteLn('[Vegetable Cuts] by Jin X (12.06.2018), fast Poya version (smartthread) :)');
    Write('Enter N (1..2147483647): ');
    try
      ReadLn(N);
    except
      WriteLn('Wrong number is specified!');
      Exit;
    end;
  
    WriteLn('Calculating...');
    try
      QueryPerformanceFrequency(TickFreq);
      QueryPerformanceCounter(Tick1);
      if N < MIN_MULTITHREAD then Answer := CountVegCut(N)
      else Answer := CountVegCutMultithread(N);
      QueryPerformanceCounter(Tick2);
    except
      WriteLn('Oh no! Something goes wrong (exception is occured)... :(');
      Exit;
    end;
  
    if Answer > 0 then StrLen := Ceil(BigInteger.Log10(Answer))
    else StrLen := 1;
    StrLen2 := StrLen;
    if StrLen >= MAX_DIGS_ASK_HEX then
    begin
      Write('Number is too large! Output in hex, it''s faster (y/n)? [y]: ');
      ReadLn(S);
      if Copy(LowerCase(S), 1, 1) <> 'n' then
      begin
        BigInteger.Hex;
        StrLen2 := Ceil(BigInteger.Log(Answer, BigInteger.Base));
      end;
      if StrLen >= MAX_DIGS_ASK_CANCEL then
      begin
        Write('OMG, it''s very very huge! Cancel output (y/n)? [y]: ');
        ReadLn(S);
        if Copy(LowerCase(S), 1, 1) <> 'n' then StrLen2 := -1
      end
    end;
    if StrLen2 = -1 then S := ''
    else
    begin
      WriteLn('Almost done...');
      if StrLen2 >= MAX_DIGS then
      begin
        Big := BigInteger.Pow(BigInteger.Base, MAX_DIGS div 2-1);
        BigInteger.DivMod(Answer, Big, Answer, Reminder);
        S := (Answer div BigInteger.Pow(BigInteger.Base, StrLen2-(MAX_DIGS-3))).ToString;
        S := S + '...' + Reminder.ToString;
      end
      else S := Answer.ToString;
    end;
  
    WriteLn;
    if S <> '' then
    begin
      WriteLn(S);
      WriteLn;
    end;
    WriteLn(IntToStrBy3(StrLen), ' digits, calculated in ', (Tick2-Tick1)/TickFreq:0:6, ' seconds');
  finally
    WriteLn('Press Enter to quit...');
    ReadLn;
  end;
end.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment