unit MathsLib;

{Copyright 2005, Gary Darby, Intellitech Systems Inc., www.DelphiForFun.org

 This program may be used or modified for any non-commercial purpose
 so long as this original notice remains in place.
 All other rights are reserved
 }

 {Assortment of math related functions and procedure in various states}


 {Revision Copyright 2006, Charles Doumar,  January 2006
 Added:
   Tprimes.BSPrime ... Binary search function to find index of prime in prime table.
   Tprimes.MaxPrimeInTable ... Returns Max Prime in prime table.
   Tprimes.GetNthPrime ... Returns Nth prime in table (returns -1 if not in table).
 Optimized:
    Tprimes.IsPrime ... Optimized for values greater than MaxPrimeInTable squared.
    Tprimes.NextPrime ... Speed up lookup of values within table range, added large value support
    Tprime.PrevPrime ... Speed up lookup of values within table range, added large value support
    Tprime.GetFactors ... Now returns factors for numbers greater than MaxPrimeInTable Squared
  }

interface

uses Classes, SysUtils, Windows, Dialogs, UBigIntsV2;

type
  intset = set of byte;

  TPoint64=record
    x,y:int64;
  end;

function GetNextPandigital(size: integer; var Digits: array of integer): boolean;
function IsPolygonal(T: int64; var rank: array of integer): boolean;
function GeneratePentagon(n: integer): integer;
function IsPentagon(p: integer): boolean;
function isSquare(const N: int64): boolean;
function isCube(const N: int64): boolean;


function isPalindrome(n: int64): boolean;
function GetEulerPhi(n: int64): int64;
//procedure SortStrDown(var s: string);
//procedure rotatestrleft(var s: string);


function IntPower(a, b: int64): int64;
function gcd2(a, b: int64): int64;
function GCDMany(A: array of integer): integer;
function LCMMany(A: array of integer): integer;
procedure ContinuedFraction(A: array of int64; const wholepart: integer;
  var numerator, denominator: int64);
function Factorial(n: int64): int64;

type
  TPrimes = class(TObject)
  protected
    function BSPrime(const n: int64; var index: integer): boolean;
  public
    Prime:    array of int64;  {Array of primes - 0th entry is not used}
    nbrprimes, nbrfactors, nbrcanonicalfactors, nbrdivisors: integer;
    Factors:  array of int64; {array of factors - 0th entry is not used}
    CanonicalFactors: array of TPoint64;
    Divisors: array of int64;
    function GetNextPrime(n: int64): int64;
    function GetPrevPrime(n: int64): int64;
    function IsPrime(n: int64): boolean;
    procedure GetFactors(const n: int64);  {get all prime factors}
    function MaxPrimeInTable: int64;
    function GetNthPrime(const n: integer): int64;
    procedure GetCanonicalFactors(const n: int64);  {get ccanonical prime factors}
    procedure GetDivisors(const n: int64);          {get all divisors}
    function Getnbrdivisors(n: int64): integer;

    constructor Create;
    destructor Destroy; override;
  end;

var
  Primes:        TPrimes;
  maxprimes:     integer = 50000; {initial size of primes array}
  maxfactors:    integer = 200;   {initial size of factors array}
  maxval:        int64 = 1000000000000;
  {10^12 - 10^6 is  max prime to be tested from tables}
  {primes up to the sqrt of maxval will be tabled}
  continuants:   array of int64;
  maxcontinuant: integer;

implementation

uses Math;

var
  nbrdiv: int64;

{************* Nbrfactors *************}
function nbrfactors(n: integer): integer;
var
  i: integer;
begin
  Result := 0;
  for i := 1 to trunc(0.0 + sqrt(n)) - 1 do
    if n mod i = 0 then
      Result := Result + 2;
  if n mod trunc(0.0 + sqrt(n)) = 0 then
    Inc(Result); {perfect square}
end;



{************ IsPrime **********}
function isprime(f: int64): boolean;
(*
var
  i:    int64;
  stop: int64;
*)
begin
  result:=primes.isprime(f);
end;

{************* IsSquare *************}
function isSquare(const N: int64): boolean;
var
  ex: extended;
  x:  int64;
begin
  ex     := sqrt(N + 0.0);
  x      := trunc(ex + 1e-10);
  Result := x * x = N;
end;

{************* IsCube *************}
function isCube(const N: int64): boolean;
var
  ex: extended;
  x:  int64;
begin
  ex     := (exp(ln(0.0 + N) / 3));
  x      := trunc(ex + 1e-6);
  Result := x * x * x = N;
end;

function gcd2(a, b: int64): int64;
  {return gcd of a and b}
  {Euclids method}
var
  g, z: int64;
begin
  g := b;
  if b <> 0 then
    while g <> 0 do
    begin
      z := a mod g;
      a := g;
      g := z;
    end;
  Result := a;
end;

function GCDMany(A: array of integer): integer;
  {Greatest common denominator}

var
  i: integer;
  {m:integer;}
  g: integer;
begin
  {   result:=0; }
   (*
   m:=minintvalue(a);
   g:=gcd2(m,a[0]);
   For i := 1 to length(a)-1 do g:=min(gcd2(m,a[i]),g);
   *)
  g := a[0];
  if length(a) >= 2 then
  begin
    g := gcd2(g, a[1]);
    if length(a) > 2 then
      for i := 2 to length(a) - 1 do
        g := gcd2(g, a[i]);
  end;
  Result := g;
end;

function LCMMany(A: array of integer): integer;
var
  i, x: integer;
begin
  x := a[0];
  for i := 1 to length(a) - 1 do
    x := (x * a[i]) div gcd2(x, a[i]);
  Result := x;
end;

function intpower(a, b: int64): int64;
var
  i: integer;
begin
  Result := 1;
  for i := 1 to b do
    Result := Result * a;
end;

function getEulerPhi(n: int64): int64;
var
  i: integer;
  p: int64;
  k: int64;
begin
  primes.getfactors(n);
  Result := 1;
  p      := primes.factors[1];
  k      := 0;
  for i := 1 to primes.nbrfactors do
  begin
    if p = primes.factors[i] then
      Inc(k)
    else
    begin
      Result := Result * (p - 1) * intpower(p, k - 1);
      k      := 1;
      p      := primes.factors[i];
    end;
  end;
  Result := Result * (p - 1) * intpower(p, k - 1);
end;


{**************** LoadCommaText **********}
procedure loadcommatext(list: TStringList; filename: string);
var
  f:    Textfile;
  line: string;
begin
  assignfile(f, filename);
  reset(f);
  readln(f, line);
  list.commatext := line;
  closefile(f);
end;


{**************** GetBigComboCount *************}
procedure GetBigCombocount(const r, n: integer; var ccount: TInteger);
{Return number of combinations -- n things taken r at a time
 without replacement}
{Return number of combinations -- n things taken r at a time
 without replacement}
var
  work: TInteger;
begin
  work := TInteger.Create;
  if (r > 0) and (r < n) then
  begin
    ccount.Assign(N);
    ccount.Factorial;
    work.Assign(r);
    work.Factorial;
    ccount.Divide(work);
    work.Assign(n - r);
    work.Factorial;
    ccount.Divide(work);
  end
  else if r = n then
    ccount.Assign(1)
  else
    ccount.Assign(0);
  work.Free;
end;



{*************** IsPandigital *******}
function ispandigital(const n, Base: int64;
  const includezero, exactlyOnce: boolean): boolean;
var
  i:      integer;
  Digits: array of integer;
begin
  SetLength(Digits, Base);
  for i := 0 to Base - 1 do
    Digits[i] := 0;
  Result := True;
  i := n;
  while i > 0 do
  begin
    Inc(Digits[i mod Base]);
    i := i div Base;
  end;
  {sure that all digit from 1 to max digit found occured exactly once and no higher digit occurred a all}
  for i := 0 to Base - 1 do
  begin
    if exactlyonce then
    begin
      if ((i = 0) and includezero and (Digits[i] <> 1)) or
        ((i > 0) and (Digits[i] <> 1)) then
      begin
        Result := False;