Boris Novgorodov

You are right, but also note that condition"if Abs(c) <= eps" really does not refer to "Not meaningful" 

Only simultaneous zero sums like "if Abs(c)  + Abs(s) <= eps" should be considered as extra case

5 hours ago, dummzeuch said:

>And some more in that thread, unfortunately I don't read Russian, so I don't know what the result of the tests were.)

That page test says that function with  $BE6EB50C constant (from GR32_Math.pas) is faster 

 

There is List of Delphi language features and version in which they were introduced/deprecated

You show decimal representation while version components correspond to nibbles of hex representation

  11264dec = $2C00, so your version number is 2.12.0.0  (2.$C.0.0)

Quick demonstration how space partition methods could help to get significant gain. Code uses simple 64x64 grid with cells 32x32. Each cell contains indexes of points belonging to this cell. At first  brute-force method counts (doubled) number of point pairs with distance<= 4. The second code piece checks only cells neighbouring to the cell containing examined point. For 64K points gain is about 200  times. Complexity is O(N^2) vs  O(N^3/2) (for proper cell size and uniform point distribution).  More advanced structures like kd-tree might provide better performance for complex cases too (worse distribution etc).

var
  Wdt, Hgt, NPts: Integer;
  i, j, k, r, m: Integer;
  t: DWord;
  sqd, maxsqd, cnt, gx, gy: Integer;
  P: TArray<TPoint>;
  Grid: array[0..63, 0..63] of TArray<Integer>;
begin
  Wdt := 2048;
  Hgt := 2048;
  NPts := 64 * 1024;
  RandSeed := 42;

  SetLength(P, NPts);
  for i := 0 to NPts - 1 do begin
    P[i] := Point(Random(Wdt), Random(Hgt));
    Grid[P[i].Y shr 5, P[i].X shr 5] := Grid[P[i].Y shr 5, P[i].X shr 5] + [i];
  end;

  maxsqd := 16;
  cnt := 0;
  t := GetTickCount;
  for j := 0 to NPts - 1 do
     for i := 0 to NPts - 1 do begin
       sqd := Sqr(P[i].X - P[j].X) + Sqr(P[i].Y - P[j].Y);
       if sqd <= maxsqd then
         Inc(cnt);
     end;
  Memo1.Lines.Add(Format('%d %d', [cnt, GetTickCount - t]));

  cnt := 0;
  t := GetTickCount;
  for i := 0 to NPts - 1 do begin
    gx := P[i].X shr 5;
    gy := P[i].Y shr 5;
    for j := Max(0, gy - 1) to Min(gy + 1, 63) do
       for k := Max(0, gx - 1) to Min(gx + 1, 63) do
          for m := 0 to High(Grid[j, k]) do begin
             r := Grid[j, k, m];
             sqd := Sqr(P[i].X - P[r].X) + Sqr(P[i].Y - P[r].Y);
             if sqd <= maxsqd then
                Inc(cnt);
          end;
  end;
  Memo1.Lines.Add(Format('%d %d', [cnt, GetTickCount - t]));

count      time (milliseconds)

115406   11466
115406    62