coprime.pas

program coprime;

{$mode objfpc}{$H+}

uses
  SysUtils;

function Gcd(a, b: Int64): Int64;
var
  x, y, r: Int64;
begin
  x := Abs(a);
  y := Abs(b);

  while y <> 0 do
  begin
    r := x mod y;
    x := y;
    y := r;
  end;

  Result := x;
end;

function EulerPhi(n: Int64): Int64;
var
  resultValue, remaining, factor: Int64;
begin
  { Euler's phi starts at n and is adjusted once for each distinct prime
    divisor p of n. Among the numbers 1..n, every p-th number shares that
    factor with n, so the current count is multiplied by (1 - 1/p):
    resultValue := resultValue - resultValue div p. }
  resultValue := n;
  remaining := n;
  factor := 2;

  while factor * factor <= remaining do
  begin
    if remaining mod factor = 0 then
    begin
      resultValue := resultValue - resultValue div factor;

      { Remove all copies of this prime factor. Phi needs each distinct prime
        divisor only once, and shrinking remaining also shortens the search. }
      while remaining mod factor = 0 do
        remaining := remaining div factor;
    end;

    if factor = 2 then
      Inc(factor)
    else
      Inc(factor, 2);
  end;

  if remaining > 1 then
    { Anything left after trial division is a prime factor larger than the
      square root of the reduced value, so it contributes one final phi term. }
    resultValue := resultValue - resultValue div remaining;

  Result := resultValue;
end;

function FormatProbability(value: Double): String;
begin
  Str(value:0:16, Result);

  while (Length(Result) > 0) and (Result[Length(Result)] = '0') do
    Delete(Result, Length(Result), 1);

  if (Length(Result) > 0) and (Result[Length(Result)] = '.') then
    Result := Result + '0';
end;

procedure PrintCoprimePairProbabilitiesUsingPhi(upToN: Int64);
var
  i, coprimePairs: Int64;
  probability: Double;
begin
  { In an i by i grid of ordered pairs (a, b), growing from (i - 1) to i adds
    the new row and column for value i. The new coprime pairs are exactly
    (i, k) and (k, i) for each 1 <= k < i with gcd(i, k) = 1. There are
    phi(i) such k values, so each step adds 2 * phi(i). The diagonal (i, i)
    is not coprime for i > 1, while (1, 1) is the initial single coprime pair. }
  coprimePairs := 1;
  {Writeln('1: ', FormatProbability(1.0));}

  if upToN >= 2 then
    for i := 2 to upToN do
    begin
      coprimePairs := coprimePairs + 2 * EulerPhi(i);
      probability := coprimePairs / (i * i);
      {Writeln(i, ': ', FormatProbability(probability));}
    end;
    Writeln(FormatProbability(probability));

end;

procedure Usage(const programName: String);
begin
  Writeln('Usage: ', programName, ' N');
  Writeln('N must be a natural number.');
end;

var
  n: Int64;
  code: Integer;

begin
  if ParamCount <> 1 then
  begin
    Usage(ParamStr(0));
    Halt(1);
  end;

  Val(ParamStr(1), n, code);
  if (code <> 0) or (n < 1) then
  begin
    Usage(ParamStr(0));
    Halt(1);
  end;

  PrintCoprimePairProbabilitiesUsingPhi(n);
end.