// PatFromSegments_VS2015_TotalCorrect.h/.cpp (single-file) - C++11/C++14 friendly (VS2015 v140)
// - U,V from bounding box (rectangular tile)
// - Normalizes to (0,0)
// - Deduplicates input (unordered_map)
// - Groups by near-parallel directions (epsAngleDeg), but writes measured angle from canonical u (no snapping)
// - Period + spacing computed with validated approxGCD
// - Merges overlaps on same generator line (unifies segments)
// - Handles segments longer than period => continuous line (empty dash list in .pat)
// - Basepoint anchored on generator line and wrapped ONLY along u (keeps it on the line)
// - Stable sort via quantized keys
//
// NOTE: assumes no segments cross tile borders (as you stated).

#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <fstream>
#include <iomanip>
#include <limits>
#include <map>
#include <numeric>
#include <string>
#include <tuple>
#include <utility>
#include <vector>
#include <unordered_map>

namespace patgen {

// VS2015-safe constant (avoid inline constexpr)
static const double PI = 3.14159265358979323846;

struct Vec2d { double x=0, y=0; };
struct Seg2d { Vec2d p1, p2; };

static inline Vec2d sub(const Vec2d& a, const Vec2d& b){ return Vec2d{a.x-b.x, a.y-b.y}; }
static inline Vec2d add(const Vec2d& a, const Vec2d& b){ return Vec2d{a.x+b.x, a.y+b.y}; }
static inline Vec2d mul(const Vec2d& a, double s){ return Vec2d{a.x*s, a.y*s}; }
static inline double dot(const Vec2d& a, const Vec2d& b){ return a.x*b.x + a.y*b.y; }
static inline double norm(const Vec2d& a){ return std::sqrt(dot(a,a)); }

static inline double rad2deg(double r){ return r * 180.0 / PI; }
static inline double wrap360(double a){
  while (a < 0) a += 360.0;
  while (a >= 360.0) a -= 360.0;
  return a;
}

// VS2015-safe quantize (avoid llround portability quirks)
static inline long long quantize(double v, double eps){
  double q = v / eps;
  return (long long)(q + (q >= 0 ? 0.5 : -0.5));
}

static inline double wrapRange(double v, double L){
  if (L <= 0) return v;
  v = std::fmod(v, L);
  if (v < 0) v += L;
  return v;
}

// Euclide approssimato + validazione "divide davvero"
static inline double approxGCD(double a, double b, double eps){
  a = std::abs(a); b = std::abs(b);
  if (a < eps) return b;
  if (b < eps) return a;

  const double a0 = a, b0 = b;

  while (b > eps) {
    double r = std::fmod(a, b);
    a = b;
    b = std::abs(r);
  }

  double g = a;
  if (g < eps) return std::min(a0, b0);

  // validate divisibility within eps
  auto divides = [&](double x)->bool{
    double r = std::fmod(x, g);
    r = std::min(r, g - r);
    return r <= eps;
  };

  if (!divides(a0) || !divides(b0)) {
    return std::min(a0, b0);
  }
  return g;
}

static inline std::vector<double> compactDash(const std::vector<double>& in, double eps){
  std::vector<double> out;
  if (in.empty()) return out;
  out.reserve(in.size());

  double cur = in[0];
  for (size_t i=1; i<in.size(); ++i) {
    const double v = in[i];
    if ((cur >= 0 && v >= 0) || (cur <= 0 && v <= 0)) cur += v;
    else {
      if (std::abs(cur) > eps) out.push_back(cur);
      cur = v;
    }
  }
  if (std::abs(cur) > eps) out.push_back(cur);

  // merge cyclic ends if same sign
  if (out.size() >= 2) {
    double a = out.front(), b = out.back();
    if ((a >= 0 && b >= 0) || (a <= 0 && b <= 0)) {
      out.front() = a + b;
      out.pop_back();
    }
  }
  return out;
}

struct PatLine {
  double angleDeg = 0.0;
  double x0=0, y0=0;
  double dx=0, dy=0;
  std::vector<double> dashes; // empty => continuous
};

struct BuildResult {
  double tileW = 0.0;
  double tileH = 0.0;
  std::vector<PatLine> lines;
};

class PatFromSegments {
public:
  struct Settings {
    double epsCoord    = 1e-6;
    double epsProj     = 1e-6;
    double epsAngleDeg = 0.15;
    double minSegLen   = 1e-9;
    int    decimals    = 8;
    bool   deduplicateInput = true;
  };

  explicit PatFromSegments(const Settings& s) : s_(s) {}

  template <typename T>
  void setSegments(const std::vector<std::array<T,4>>& segsP1P2xy) {
    segs_.clear();
    segs_.reserve(segsP1P2xy.size());
    for (size_t i=0; i<segsP1P2xy.size(); ++i) {
      Seg2d s;
      s.p1 = Vec2d{ (double)segsP1P2xy[i][0], (double)segsP1P2xy[i][1] };
      s.p2 = Vec2d{ (double)segsP1P2xy[i][2], (double)segsP1P2xy[i][3] };
      segs_.push_back(s);
    }
  }

  void setSegments(const std::vector<Seg2d>& segs) { segs_ = segs; }

  const std::string& lastWarning() const { return lastWarning_; }

  BuildResult build() const {
    lastWarning_.clear();
    BuildResult res;
    if (segs_.empty()) return res;

    // ---- 1) Bounding box
    double minx =  std::numeric_limits<double>::infinity();
    double miny =  std::numeric_limits<double>::infinity();
    double maxx = -std::numeric_limits<double>::infinity();
    double maxy = -std::numeric_limits<double>::infinity();

    for (size_t i=0; i<segs_.size(); ++i) {
      const Seg2d& s = segs_[i];
      minx = std::min(minx, std::min(s.p1.x, s.p2.x));
      miny = std::min(miny, std::min(s.p1.y, s.p2.y));
      maxx = std::max(maxx, std::max(s.p1.x, s.p2.x));
      maxy = std::max(maxy, std::max(s.p1.y, s.p2.y));
    }

    const double rawW = maxx - minx;
    const double rawH = maxy - miny;
    const double W = std::max(s_.epsCoord, rawW);
    const double H = std::max(s_.epsCoord, rawH);

    res.tileW = W;
    res.tileH = H;

    if (rawW < s_.epsCoord || rawH < s_.epsCoord) {
      lastWarning_ = "Degenerate tile bounding box (W or H ~ 0). Pattern may be empty.";
    }

    // ---- 2) Normalize to origin
    Vec2d shift{minx, miny};
    std::vector<Seg2d> segs;
    segs.reserve(segs_.size());
    for (size_t i=0; i<segs_.size(); ++i) {
      Seg2d s = segs_[i];
      s.p1 = sub(s.p1, shift);
      s.p2 = sub(s.p2, shift);
      if (norm(sub(s.p2, s.p1)) >= s_.minSegLen) segs.push_back(s);
    }
    if (segs.empty()) return res;

    // ---- 3) Optional dedup input (unordered_map)
    if (s_.deduplicateInput) {
      struct SegKey {
        long long x1,y1,x2,y2;
        bool operator==(const SegKey& o) const {
          return x1==o.x1 && y1==o.y1 && x2==o.x2 && y2==o.y2;
        }
      };
      struct SegKeyHash {
        std::size_t operator()(const SegKey& k) const {
          // 64-bit-ish mix
          std::size_t h = 1469598103934665603ULL;
          auto mix = [&](long long v){
            h ^= (std::size_t)v + 0x9e3779b97f4a7c15ULL + (h<<6) + (h>>2);
          };
          mix(k.x1); mix(k.y1); mix(k.x2); mix(k.y2);
          return h;
        }
      };

      std::unordered_map<SegKey, Seg2d, SegKeyHash> uniq;
      uniq.reserve(segs.size() * 2);

      for (size_t i=0; i<segs.size(); ++i) {
        Seg2d s = segs[i];
        // stabilize endpoints
        if (s.p2.x < s.p1.x || (s.p2.x == s.p1.x && s.p2.y < s.p1.y)) std::swap(s.p1, s.p2);

        SegKey k{
          quantize(s.p1.x, s_.epsCoord),
          quantize(s.p1.y, s_.epsCoord),
          quantize(s.p2.x, s_.epsCoord),
          quantize(s.p2.y, s_.epsCoord)
        };
        uniq.emplace(k, s);
      }

      segs.clear();
      segs.reserve(uniq.size());
      for (typename std::unordered_map<SegKey, Seg2d, SegKeyHash>::const_iterator it = uniq.begin();
           it != uniq.end(); ++it) {
        segs.push_back(it->second);
      }
      if (segs.empty()) return res;
    }

    // ---- helpers
    auto unitDir = [&](const Seg2d& s)->Vec2d {
      Vec2d v = sub(s.p2, s.p1);
      double L = norm(v);
      if (L < s_.minSegLen) return Vec2d{1,0};
      v = mul(v, 1.0/L);
      // canonical sign
      if (v.x < 0 || (std::abs(v.x) < s_.epsCoord && v.y < 0)) v = mul(v, -1.0);
      return v;
    };
    auto angleFromU = [&](const Vec2d& u)->double {
      return wrap360(rad2deg(std::atan2(u.y, u.x)));
    };

    // ---- 4) Group by near-parallel direction
    struct DirGroup {
      Vec2d u;
      double angleDeg;
      std::vector<Seg2d> segs;
    };
    std::vector<DirGroup> groups;
    groups.reserve(16);

    for (size_t i=0; i<segs.size(); ++i) {
      const Seg2d& sg = segs[i];
      Vec2d u = unitDir(sg);
      double a = angleFromU(u);

      bool placed = false;
      for (size_t j=0; j<groups.size(); ++j) {
        DirGroup& g = groups[j];
        double da = std::abs(g.angleDeg - a);
        da = std::min(da, 360.0 - da);
        double da180 = std::abs(wrap360(g.angleDeg + 180.0) - a);
        da180 = std::min(da180, 360.0 - da180);
        if (std::min(da, da180) <= s_.epsAngleDeg) {
          g.segs.push_back(sg);
          placed = true;
          break;
        }
      }
      if (!placed) {
        DirGroup ng;
        ng.u = u;
        ng.angleDeg = a;
        ng.segs.push_back(sg);
        groups.push_back(ng);
      }
    }

    Vec2d Uvec{W, 0.0};
    Vec2d Vvec{0.0, H};

    std::vector<PatLine> out;

    // ---- 5) For each direction group
    for (size_t gi=0; gi<groups.size(); ++gi) {
      const DirGroup& g = groups[gi];
      const Vec2d u = g.u;
      const Vec2d n = Vec2d{-u.y, u.x};

      // spacing between parallels = gcd(|dot(n,U)|,|dot(n,V)|)
      double spU = std::abs(dot(n, Uvec));
      double spV = std::abs(dot(n, Vvec));
      if (spU < s_.epsCoord && spV < s_.epsCoord) continue;

      double spacing = 0.0;
      if (spU < s_.epsCoord) spacing = spV;
      else if (spV < s_.epsCoord) spacing = spU;
      else {
        spacing = approxGCD(spU, spV, s_.epsCoord);
        if (spacing < s_.epsCoord) spacing = std::min(spU, spV);
      }
      Vec2d delta = mul(n, spacing);

      // period along u = gcd(|dot(u,U)|,|dot(u,V)|)
      double pu = std::abs(dot(u, Uvec));
      double pv = std::abs(dot(u, Vvec));
      double period = 0.0;
      if (pu < s_.epsProj && pv < s_.epsProj) period = std::max(W, H);
      else if (pu < s_.epsProj) period = pv;
      else if (pv < s_.epsProj) period = pu;
      else {
        period = approxGCD(pu, pv, s_.epsProj);
        if (period < s_.epsProj) period = std::min(pu, pv);
      }

      // group by generator line offset c=dot(n,p)
      struct LineBucket { double cSum=0; int cCnt=0; std::vector<Seg2d> segs; };
      std::map<long long, LineBucket> byLine;

      for (size_t si=0; si<g.segs.size(); ++si) {
        const Seg2d& sg = g.segs[si];
        double c = dot(n, sg.p1);
        long long qc = quantize(c, s_.epsCoord);
        LineBucket& b = byLine[qc];
        b.cSum += c;
        b.cCnt += 1;
        b.segs.push_back(sg);
      }

      // each generator line
      for (std::map<long long, LineBucket>::iterator it = byLine.begin(); it != byLine.end(); ++it) {
        LineBucket& bucket = it->second;
        if (bucket.segs.empty()) continue;

        const double cMean = (bucket.cCnt > 0) ? (bucket.cSum / bucket.cCnt) : 0.0;

        // basepoint anchored on line, wrapped ONLY along u
        Vec2d base0 = mul(n, cMean);
        double tBase = dot(u, base0);
        double tWrapped = wrapRange(tBase, period);
        Vec2d base = add(mul(n, cMean), mul(u, tWrapped));

        struct Interval { double a,b; };
        std::vector<Interval> intervals;
        intervals.reserve(bucket.segs.size() * 2);

        auto modP = [&](double t){
          t = std::fmod(t, period);
          if (t < 0) t += period;
          return t;
        };

        for (size_t k=0; k<bucket.segs.size(); ++k) {
          const Seg2d& sg = bucket.segs[k];
          double t1 = dot(u, sg.p1);
          double t2 = dot(u, sg.p2);
          if (t2 < t1) std::swap(t1, t2);
          double L = t2 - t1;
          if (L < s_.epsProj) continue;

          // covers whole period -> continuous
          if (L >= period - s_.epsProj) {
            intervals.clear();
            intervals.push_back(Interval{0.0, period});
            break;
          }

          double a = modP(t1);
          double b = modP(t2);
          if (b >= a) intervals.push_back(Interval{a,b});
          else {
            intervals.push_back(Interval{0.0,b});
            intervals.push_back(Interval{a,period});
          }
        }

        if (intervals.empty()) continue;

        std::sort(intervals.begin(), intervals.end(),
                  [](const Interval& i1, const Interval& i2){ return i1.a < i2.a; });

        // merge overlaps
        std::vector<Interval> merged;
        merged.reserve(intervals.size());
        Interval cur = intervals.front();
        for (size_t mi=1; mi<intervals.size(); ++mi) {
          const Interval& in = intervals[mi];
          if (in.a <= cur.b + s_.epsProj) cur.b = std::max(cur.b, in.b);
          else { merged.push_back(cur); cur = in; }
        }
        merged.push_back(cur);

        // build cyclic dash/gap
        std::vector<double> dashes;
        double t = 0.0;
        for (size_t mi=0; mi<merged.size(); ++mi) {
          const Interval& m = merged[mi];
          double gap = m.a - t;
          if (gap > s_.epsProj) dashes.push_back(-gap);
          double dash = m.b - m.a;
          if (dash > s_.epsProj) dashes.push_back(dash);
          t = m.b;
        }
        double closingGap = period - t;
        if (closingGap > s_.epsProj) dashes.push_back(-closingGap);

        // rotate to start with dash if possible
        if (!dashes.empty() && dashes.front() < 0) {
          for (size_t r=0; r<dashes.size(); ++r) {
            std::rotate(dashes.begin(), dashes.begin()+1, dashes.end());
            if (dashes.front() > 0) break;
          }
        }

        dashes = compactDash(dashes, s_.epsProj);

        // continuous line in .pat => empty dash list
        if (dashes.size() == 1 && dashes[0] >= period - s_.epsProj) dashes.clear();

        PatLine pl;
        pl.angleDeg = g.angleDeg;
        pl.x0 = base.x;
        pl.y0 = base.y;
        pl.dx = delta.x;
        pl.dy = delta.y;
        pl.dashes.swap(dashes);
        out.push_back(pl);
      }
    }

    // stable sort via quantized keys
    struct Key { long long qa,qy,qx; };
    auto makeKey = [&](const PatLine& p)->Key{
      Key k;
      k.qa = quantize(p.angleDeg, s_.epsAngleDeg);
      k.qy = quantize(p.y0, s_.epsCoord);
      k.qx = quantize(p.x0, s_.epsCoord);
      return k;
    };

    std::sort(out.begin(), out.end(),
      [&](const PatLine& a, const PatLine& b){
        Key ka = makeKey(a);
        Key kb = makeKey(b);
        if (ka.qa != kb.qa) return ka.qa < kb.qa;
        if (ka.qy != kb.qy) return ka.qy < kb.qy;
        return ka.qx < kb.qx;
      });

    if (out.empty() && (rawW < s_.epsCoord || rawH < s_.epsCoord)) {
      if (lastWarning_.empty()) lastWarning_ = "Degenerate tile: could not build any .pat lines.";
    }

    res.lines.swap(out);
    return res;
  }

  bool writePat(const std::string& filename,
                const std::string& patName,
                const std::string& description) const
  {
    BuildResult br = build();
    if (br.lines.empty()) return false;

    std::ofstream out(filename.c_str(), std::ios::out | std::ios::trunc);
    if (!out) return false;

    out << "*" << patName << "," << description << "\n";
    out << std::fixed << std::setprecision(s_.decimals);

    for (size_t i=0; i<br.lines.size(); ++i) {
      const PatLine& l = br.lines[i];
      out << l.angleDeg << ", "
          << l.x0 << ", " << l.y0 << ", "
          << l.dx << ", " << l.dy;

      for (size_t j=0; j<l.dashes.size(); ++j) {
        double d = l.dashes[j];
        if (std::abs(d) < s_.epsProj) d = 0.0;
        out << ", " << d;
      }
      out << "\n";
    }
    return true;
  }

private:
  Settings s_;
  std::vector<Seg2d> segs_;
  mutable std::string lastWarning_;
};

} // namespace patgen

// ---------------------
// Example usage
// ---------------------
/*
#include "PatFromSegments_VS2015_TotalCorrect.h"

int main() {
  using namespace patgen;

  PatFromSegments::Settings s;
  s.epsCoord = 1e-5;
  s.epsProj  = 1e-5;
  s.epsAngleDeg = 0.2;
  s.decimals = 8;

  PatFromSegments conv(s);

  std::vector<std::array<float,4>> segs = {
    {0.f,0.f,  1.f,0.f},
    {1.f,0.f,  2.f,0.f},
    {0.f,1.f,  2.f,1.f},
  };

  conv.setSegments(segs);

  if (!conv.writePat("tile.pat", "TILEPAT", "Converted from ARC+ 2D segments")) {
    // optional:
    // std::cout << conv.lastWarning() << std::endl;
  }
  return 0;
}
*/