module ray;
import std.math, std.string, std.stdio;

alias double FP;

FP delta; static this() { delta = sqrt(double.epsilon); }

FP max(FP a, FP b) { return (a>b)?a:b; }

align(16) struct Vec {
  FP x, y, z;
  Vec opAdd(ref Vec other) { return Vec(x + other.x, y + other.y, z + other.z); }
  Vec opSub(ref Vec other) { return Vec(x - other.x, y - other.y, z - other.z); }
  Vec opMul(FP other) { return Vec(other * x, other * y, other * z); }
  FP dot(ref Vec other) { return x * other.x + y * other.y + z * other.z; }
  FP length() { return sqrt(dot(*this)); }
  Vec unit() { return *this * (1.0 / length); }
  string toString() { return format("{", x, ", ", y, ", ", z, "}"); }
}

struct Hit {
  FP dist;
  Vec normal;
  static Hit opCall() { Hit res; res.dist = FP.infinity; return res; }
  bool none() { return dist == FP.infinity; }
}

class Scene {
  abstract void intersect(ref Vec, ref Hit);
  abstract bool intersect(ref Vec, ref Vec);
  abstract void bound(ref Sphere b);
}

final class Sphere : Scene {
  Vec center; FP radius;
  this(Vec c, FP r) { center = c; radius = r; }
  private FP ray_sphere(ref Vec dir) {
    auto b = center.dot(dir), disc = b*b - center.dot(center) + radius * radius;
    if (disc > 0) {
      auto d = sqrt(disc), t2 = b + d;
      if (t2 > 0) {
        auto t1 = b - d;
        return (t1 > 0) ? t1 : t2;
      }
    }
    return FP.infinity;
  }
  override void intersect(ref Vec dir, ref Hit hit) {
    auto dist = ray_sphere(dir);
    if (dist < hit.dist) {
      hit.dist = dist;
      auto n = dist * dir - center, il = 1.0 / sqrt(n.dot(n));
      hit.normal = il * n;
    }
  }
  override bool intersect(ref Vec orig, ref Vec dir) {
    auto v = center - orig, b = v.dot(dir), disc = b*b - v.dot(v) + radius*radius;
    return (disc < 0) ? false : (b + sqrt(disc) >= 0);
  }
  override void bound(ref Sphere b) { b.radius = max(b.radius, (center - b.center).length + radius); }
}

final class Group : Scene {
  Sphere b; Scene[] child;
  this(Sphere _b, Scene[] c) {
    b = _b; child = c;
  }
  override void intersect(ref Vec dir, ref Hit hit) {
    if (b.ray_sphere(dir) < hit.dist) {
      // foreach (sc; child) sc.intersect(dir, hit);
      for (int i = 0; i < child.length; ++i) child[i].intersect(dir, hit);
    }
  }
  override bool intersect(ref Vec orig, ref Vec dir) {
    if (!b.intersect(orig, dir)) return false;
    // foreach (sc; child) if (sc.intersect(orig, dir)) return true;
    for (int i = 0; i < child.length; ++i)
      if (child[i].intersect(orig, dir)) return true;
    return false;
  }
  override void bound(ref Sphere b) { foreach (sc; child) sc.bound(b); }
}

FP ray_trace(ref Vec neg_light, ref Vec dir, Scene sc) {
  auto hit = Hit();
  sc.intersect(dir, hit);
  if (hit.none) return 0;
  auto g = hit.normal.dot(neg_light);
  if (g < 0) return 0;
  auto p = hit.dist * dir + delta * hit.normal;
  return (sc.intersect(p, neg_light) ? 0 : g);
}

Scene create(int level, Vec c, FP r) {
  auto s = new Sphere(c, r);
  if (level == 1) return s;
  Scene[] child;
  child ~= s;
  auto rn = 3 * r / sqrt(12.0);
  for (int dz = -1; dz <= 1; dz += 2)
    for (int dx = -1; dx <= 1; dx += 2)
      child ~= create(level-1, c + rn * Vec(dx, 1,  dz), r/2);
  auto b2 = new Sphere(c + Vec(0, r, 0), r);
  foreach (ch; child) ch.bound(b2);
  return new Group(b2, child);
}

void main(string[] args) {
  int level = 9, n = 512, ss = 4;
  if (args.length == 3) {
    level = args[1].atoi(); n = args[2].atoi();
  }
  auto neg_light = Vec(1, 3, -2).unit, s = create(level, Vec(0, -1, 4), 1);
  writefln("P5\n", n, " ", n, "\n255");
  for (int y = n - 1; y >= 0; --y)
    for (int x = 0; x < n; ++x) {
      FP g = 0;
      for (int dx = 0; dx < ss; ++dx)
        for (int dy = 0; dy < ss; ++dy) {
          auto dir = Vec(x + dx * 1.0 / ss - n / 2.0, y + dy * 1.0 / ss - n / 2.0, n).unit;
          g += ray_trace(neg_light, dir, s);
        }
      printf("%c", cast(char) cast(ubyte) (0.5 + 255.0 * g / (ss * ss)));
    }
  return 0;
}