Poisson Disk Generator (3 стр)

import std.getopt:getopt;
import std.math:abs,pow,trunc;
import std.stdio:File,writef,writefln;
import heapsort:HeapSort;

struct pair
 {
  ulong N;
  double coord;
  bool opEquals(ref const pair p)
   {
    return this.coord == p.coord;
   }
  byte opCmp(ref const pair p)
   {
    return (this.coord < p.coord) ? -1 : (this.coord > p.coord) ? 1 : 0;
   }
 }

struct point
 {
  double x;
  double y;
  bool IsOn = true;
  bool opEquals(ref const point p)
   {
    return (this.x == p.x) && (this.y == p.y);
   }
 }

double[] random(ulong N)
 {
  const ubyte L=16;
  
  double sum(ref double[L] a)
   {
    double s = 0;
    for (ubyte k = 0; k < a.length; k++) s += a[k];
    return s;
   }
  double NextRand(double a, double m)
   {
    a *= m;
    return a - trunc(a);
   }
  double NextRandSum(ref double[L] r)
   {
    double s = sum(r)*1.000000001;
    return s - trunc(s);
   }
  double[L] r = 0;
  double a = 0.5;
  double m = pow(pow(2.0,5.0), pow(2.0,0.5));
  for (ulong k = 0; k < 100000; k++) a = NextRand(a, m);
  for (ulong k = 0; k < L; k++) r[k] = (a = NextRand(a, m));
  auto rnd = new double[N];
  ulong z = 0;
  while (z < N)
   {
    for (ulong k = 0; k < L; k++)
     {
      rnd[z++] = (r[k] = NextRandSum(r));
      if (z == N) break;
     }
   }
  return rnd;
 }

ulong[] nearests(ref pair[] pairs, ulong N, double coord, double r)
 {
  ulong[] nearest;
  double LB = coord - r;
  double RB = coord + r;
  ulong pos = pairs.length;
  ulong step = pos / 2;
  while (step > 0)
   {
    if (LB <= pairs[pos - step].coord)
     {
      pos -= step;
      while (pos < step)
       {
        step /= 2;
       }
     }
    else
     {
      step /= 2;
     }
   }
  while ((pos < pairs.length) && (LB < pairs[pos].coord) && (pairs[pos].coord < RB))
   {
    if (pairs[pos].N != N) nearest ~= [pairs[pos].N];
    pos++;
   }
  return nearest;
 }

void PointsClean(ref point[] points, ref pair[] X, ulong k, double r)
 {
  double r2 = r * r;
  double dx, dy;
  foreach (ulong m; nearests(X, k, points[k].x, r))
   {
    if (points[m].IsOn && ((dx = points[m].x - points[k].x) * dx +
                           (dy = points[m].y - points[k].y) * dy < r2)) points[m].IsOn = false;
   }
 }

void main(string[] args)
 {
  ulong N;
  double r;
  getopt(args, "N", &N, "r", &r);
  auto points = new point[N];
  auto X = new pair[N];
  auto data = random(2 * N);
  for (ulong k = 0; k < N; k++)
   {
    X[k].coord = (points[k].x = data[k]);
    points[k].y = data[k + N];
    X[k].N = k;
   }
  HeapSort(X);
  auto PointsFile = File("points.txt","w");
  for (ulong k = 0; k < N; k++)
   {
    if (points[k].IsOn)
     {
      PointsClean(points, X, k, r);
      PointsFile.writefln("%.16f %.16f", points[k].x, points[k].y);
     }
   }
  PointsFile.close();
 }