Devouring Arrow simulation

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <iostream>
#include <vector>
#include <algorithm>
#include <random>
using namespace std;

const double EPS = 1e-6;
const double INF = 1e+12;

// CHARACTER CONFIG
const double BonusPierceChance = 0.25; // Ninth Cirri Satchel
const int GuaranteedPierce = 2; // Buriza-do Kyanon
const double CritChance = 0.6; // CC for Spray of Teeth
const int Trials = 1000000; // number of simulations

// TARGET CONFIG
struct Config {
  double distance;
  double radius;
  double size;
  int count;
} config = {
  22, // distance from the target
  1, // target spread (max distance from the center)
  2.5, // target size (collision radius)
  1, // number of targets
};

// SKILL DATA (do not change)
const double BaseDamage = 1.55;
const double Speed = 1.3;
const double SeekInterval = 0.3;
const double SeekIntervalDelta = 0.3;
const double SearchInitial = 11;
const double SearchAfterPierce = 20;
const double PierceChance = 0.35;
const double MaxTravelDistance = 120;
const double MaxDistanceFromCast = 90;
const double B_ConeSize = 50;
const double C_BonusDamage = 0.7;
const double E_Damage = 0.6;
const double E_Range = 10;
const double D_PierceChanceBonus = 0.15;

// SKILL RUNES
struct {
  char sym;
  char const* name;
} Runes[] = {
  { 'x', "Unruned" },
  { 'd', "Puncturing Arrow" },
  { 'a', "Serrated Arrow" },
  { 'b', "Shatter Shot" },
  { 'c', "Devouring Arrow" },
  { 'e', "Spray of Teeth" },
};

// vector class
struct Vector {
  double x, y;
  Vector(double x = 0, double y = 0)
    : x(x), y(y)
  {}
  Vector(Vector const& v)
    : x(v.x), y(v.y)
  {}

  double len() const {
    return sqrt(x * x + y * y);
  }
  double norm() const {
    return sqrt(x * x + y * y);
  }
  double len2() const {
    return x * x + y * y;
  }
  double operator ~() const {
    return sqrt(x * x + y * y);
  }

  Vector& normalize() {
    double n = norm();
    x /= n;
    y /= n;
    return *this;
  }
  Vector normalized() const {
    double n = norm();
    return Vector(x / n, y / n);
  }

  Vector operator -() const {
    return Vector(-x, -y);
  }
  Vector operator +() const {
    return Vector(x, y);
  }
  Vector perp() const {
    return Vector(-y, x);
  }

  Vector& operator =(Vector const& rhs) {
    x = rhs.x;
    y = rhs.y;
    return *this;
  }
  Vector& operator +=(Vector const& rhs) {
    x += rhs.x;
    y += rhs.y;
    return *this;
  }
  Vector& operator -=(Vector const& rhs) {
    x -= rhs.x;
    y -= rhs.y;
    return *this;
  }
  Vector& operator *=(double rhs) {
    x *= rhs;
    y *= rhs;
    return *this;
  }
  Vector& operator /=(double rhs) {
    x /= rhs;
    y /= rhs;
    return *this;
  }

  friend Vector operator +(Vector const& lhs, Vector const& rhs) {
    return Vector(lhs.x + rhs.x, lhs.y + rhs.y);
  }
  friend Vector operator -(Vector const& lhs, Vector const& rhs) {
    return Vector(lhs.x - rhs.x, lhs.y - rhs.y);
  }
  friend Vector operator *(Vector const& lhs, double rhs) {
    return Vector(lhs.x * rhs, lhs.y * rhs);
  }
  friend Vector operator *(double lhs, Vector const& rhs) {
    return Vector(lhs * rhs.x, lhs * rhs.y);
  }
  friend Vector operator /(Vector const& lhs, double rhs) {
    return Vector(lhs.x / rhs, lhs.y / rhs);
  }
  friend double operator &(Vector const& lhs, Vector const& rhs) {
    return lhs.x * rhs.x + lhs.y * rhs.y;
  }
  friend double operator ^(Vector const& lhs, Vector const& rhs) {
    return lhs.x * rhs.y - lhs.y * rhs.x;
  }
};

mt19937_64 engine;
uniform_real_distribution<double> rng(0, 1);
uniform_real_distribution<double> seekRng(SeekInterval, SeekInterval + SeekIntervalDelta);

const double PI = 2.0 * acos(0);
// rotate a vector `deg' degrees CCW
Vector rotate(Vector& v, double deg) {
  deg *= PI / 180.0;
  double cs = cos(deg);
  double sn = sin(deg);
  return Vector(v.x * cs - v.y * sn, v.x * sn + v.y * cs);
}

// target setup
struct Setup {
  Vector pos, dir; // current projectile position/direction
  vector<Vector> targets; // target positions
  double size; // target sizes

  // generate random targets based on the config
  Setup(Config const& cfg) {
    uniform_real_distribution<double> distr(-cfg.radius, cfg.radius);
    for (int i = 0; i < cfg.count; ++i) {
      Vector cur;
      do {
        cur.x = distr(engine);
        cur.y = distr(engine);
      } while (cur.len2() > cfg.radius * cfg.radius);
      targets.push_back(cur);
    }
    pos = Vector(-cfg.distance, 0);
    dir = Vector(1, 0);
    size = cfg.size;
  }

  // find the closest target to the current position
  Vector const& closest() {
    size_t index = 0;
    double dist2 = (targets[0] - pos).len2();
    for (size_t i = 1; i < targets.size(); ++i) {
      double cur2 = (targets[i] - pos).len2();
      if (cur2 < dist2) {
        dist2 = cur2;
        index = i;
      }
    }
    return targets[index];
  }

  // get the distance until collision with the given target,
  // or INF if the ray does not intersect the target
  double impact(Vector const& target) {
    Vector to = target - pos;
    double dy = abs(dir ^ to);
    if (dy > size - EPS) return INF;
    double dx = dir & to;
    double dt = sqrt(size * size - dy * dy);
    if (dx - dt < EPS) return INF;
    return dx - dt;
  }
  // get the time until next collision, or INF if no collisions
  double impact() {
    double best = impact(targets[0]);
    for (size_t i = 1; i < targets.size(); ++i) {
      best = min(best, impact(targets[i]));
    }
    return (best < INF ? best / Speed : best);
  }
  // get the number of targets within `range' yards of current position
  int aoe(double range) const {
    int count = 0;
    range += size;
    for (auto& v : targets) {
      if ((v - pos).len2() < range * range) {
        ++count;
      }
    }
    return count;
  }
};

double run(Setup& cfg, char rune) {
  double impact = cfg.impact(); // time of next collision
  double finish = MaxTravelDistance / Speed; // projectile expiration time
  double seek = seekRng(engine) * 60; // time until next direction change
  double prev = 0; // previous event time
  double seekRange = SearchInitial; // current seek range (11 before first hit, 20 after)
  double total = 0; // total damage
  double damage = BaseDamage; // current damage
  double chance = PierceChance + (rune == 'd' ? D_PierceChanceBonus : 0) + BonusPierceChance; // pierce chance
  int index = 0; // number of pierces
  while (finish > impact || finish > seek) {
    if (impact < seek) {
      // next event = collision
      total += damage; // add current damage
      seekRange = SearchAfterPierce; // now seek range = 20
      cfg.pos += cfg.dir * Speed * (impact - prev); // update projectile position
      prev = impact; // update last event time
      if (rune == 'e') {
        // add Spray of Teeth damage (TODO: currently assumes it can crit as well)
        total += E_Damage * CritChance * cfg.aoe(E_Range);
      }
      // check if we should pierce (Buriza gives 1-2 guaranteed pierces; otherwise do rng check)
      if (!(rune && index < GuaranteedPierce) && rng(engine) > chance) break;
      ++index;
      if (rune == 'b' && index == 1) {
        // fire spray of teeth on the first pierce
        // (assuming the extra projectiles are unruned and do not benefit from Buriza)
        Setup left(cfg);
        left.dir = rotate(cfg.dir, B_ConeSize);
        total += run(left, 0);
        Setup right(cfg);
        right.dir = rotate(cfg.dir, -B_ConeSize);
        total += run(right, 0);
      }
      if (rune == 'c') {
        // apply Devouring Arrow damage
        damage += C_BonusDamage * BaseDamage;
      }
      impact = prev + cfg.impact(); // calculate next collision time
    } else {
      // next event = seek
      cfg.pos += cfg.dir * Speed * (seek - prev); // update projectile position
      prev = seek; // update last event time
      auto to = cfg.closest(); // find closest target
      if ((cfg.pos - to).len2() < seekRange * seekRange) {
        cfg.dir = (to - cfg.pos).normalized(); // change direction
        impact = prev + cfg.impact(); // calculate next collision time
      }
      seek += seekRng(engine) * 60; // get next redirect time
    }
  }
  return total;
}

// average result over a bunch of tests
double run_avg(Config const& cfg, char rune) {
  double avg = 0;
  for (int i = 0; i < Trials; ++i) {
    avg += run(Setup(cfg), rune);
  }
  return avg / Trials;
}

// find the best distance to shoot from
pair<double, double> run_best(Config& cfg, char rune, double from, double to, double delta) {
  cfg.distance = from;
  pair<double, double> best(cfg.distance, run_avg(cfg, rune));
  for (cfg.distance = from + delta; cfg.distance <= to; cfg.distance += delta) {
    pair<double, double> cur(cfg.distance, run_avg(cfg, rune));
    if (cur.second > best.second) {
      best = cur;
    }
  }
  return best;
}

template<class X, class Y>
ostream& operator <<(ostream& os, pair<X, Y> const& p) {
  return os << p.first << " " << p.second;
}

int main() {
  for (auto& rune : Runes) {
    cout << rune.name << ": " << run_avg(config, rune.sym) << endl;
  }
  return 0;
}