// Standalone search harness that replicates SkEdgeClipper::clipCubic float
// arithmetic exactly, to find inputs that overflow fPoints[54]/fVerbs[18].
//
// Build:  clang++ -O2 -std=c++17 search.cpp -o search -lm
// Run:    ./search [seed] [iterations]

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstring>
#include <random>

struct Pt { float x, y; };

// ---------- valid_unit_divide (SkGeometry.cpp:54) ----------
static int valid_unit_divide(float numer, float denom, float* ratio) {
    if (numer < 0) { numer = -numer; denom = -denom; }
    if (denom == 0 || numer == 0 || numer >= denom) return 0;
    float r = numer / denom;
    if (std::isnan(r)) return 0;
    if (r == 0) return 0;
    *ratio = r;
    return 1;
}

// ---------- SkFindUnitQuadRoots (SkGeometry.cpp:95) ----------
static int find_unit_quad_roots(float A, float B, float C, float roots[2]) {
    if (A == 0) return valid_unit_divide(-C, B, roots);
    float* r = roots;
    double dr = (double)B * B - 4 * (double)A * C;
    if (dr < 0) return 0;
    dr = std::sqrt(dr);
    float R = (float)dr;
    if (!std::isfinite(R)) return 0;
    float Q = (B < 0) ? -(B - R) / 2 : -(B + R) / 2;
    r += valid_unit_divide(Q, A, r);
    r += valid_unit_divide(C, Q, r);
    if (r - roots == 2) {
        if (roots[0] > roots[1]) std::swap(roots[0], roots[1]);
        else if (roots[0] == roots[1]) r -= 1;
    }
    return (int)(r - roots);
}

// ---------- SkFindCubicExtrema (SkGeometry.cpp:454) ----------
static int find_cubic_extrema(float a, float b, float c, float d, float t[2]) {
    float A = d - a + 3 * (b - c);
    float B = 2 * (a - b - b + c);
    float C = b - a;
    return find_unit_quad_roots(A, B, C, t);
}

// ---------- SkChopCubicAt single t (SkGeometry.cpp:473) ----------
static inline float mix(float a, float b, float t) { return (b - a) * t + a; }

static void chop_cubic_at(const Pt src[4], Pt dst[7], float t) {
    if (t == 1) {
        memcpy(dst, src, sizeof(Pt) * 4);
        dst[4] = dst[5] = dst[6] = src[3];
        return;
    }
    float ax = src[0].x, ay = src[0].y;
    float bx = src[1].x, by = src[1].y;
    float cx = src[2].x, cy = src[2].y;
    float dx = src[3].x, dy = src[3].y;
    float abx = mix(ax, bx, t), aby = mix(ay, by, t);
    float bcx = mix(bx, cx, t), bcy = mix(by, cy, t);
    float cdx = mix(cx, dx, t), cdy = mix(cy, dy, t);
    float abcx = mix(abx, bcx, t), abcy = mix(aby, bcy, t);
    float bcdx = mix(bcx, cdx, t), bcdy = mix(bcy, cdy, t);
    float abcdx = mix(abcx, bcdx, t), abcdy = mix(abcy, bcdy, t);
    dst[0] = {ax, ay};
    dst[1] = {abx, aby};
    dst[2] = {abcx, abcy};
    dst[3] = {abcdx, abcdy};
    dst[4] = {bcdx, bcdy};
    dst[5] = {cdx, cdy};
    dst[6] = {dx, dy};
}

// ---------- SkChopCubicAt with two t values (SkGeometry.cpp:504) ----------
static void chop_cubic_at_2(const Pt src[4], Pt dst[10], float t0, float t1) {
    if (t1 == 1) {
        chop_cubic_at(src, dst, t0);
        dst[7] = dst[8] = dst[9] = src[3];
        return;
    }
    // emulate the SIMD version: lo lane = t0, hi lane = t1
    auto do_lane = [&](float T, Pt out[7]) { chop_cubic_at(src, out, T); };
    Pt lo[7], hi[7];
    do_lane(t0, lo);
    do_lane(t1, hi);
    // middle = unchecked_mix(abc, bcd, shuffle<2,3,0,1>(T))
    // i.e. lo middle uses t1, hi middle uses t0
    // abc/bcd come from each lane
    // lo lane: middle = mix(lo.abc, lo.bcd, t1)
    // hi lane: middle = mix(hi.abc, hi.bcd, t0)
    float lo_mx = mix(lo[2].x, lo[4].x, t1), lo_my = mix(lo[2].y, lo[4].y, t1);
    float hi_mx = mix(hi[2].x, hi[4].x, t0), hi_my = mix(hi[2].y, hi[4].y, t0);
    dst[0] = lo[0];
    dst[1] = lo[1];
    dst[2] = lo[2];
    dst[3] = lo[3];
    dst[4] = {lo_mx, lo_my};
    dst[5] = {hi_mx, hi_my};
    dst[6] = hi[3];
    dst[7] = hi[4];
    dst[8] = hi[5];
    dst[9] = hi[6];
}

// ---------- SkChopCubicAt with array of t (SkGeometry.cpp:541) ----------
static void chop_cubic_at_n(const Pt src[4], Pt dst[], const float t[], int n) {
    if (n == 0) { memcpy(dst, src, 4 * sizeof(Pt)); return; }
    Pt cur[4]; memcpy(cur, src, sizeof(cur));
    Pt* d = dst;
    int i = 0;
    for (; i < n - 1; i += 2) {
        float t0 = t[i], t1 = t[i + 1];
        if (i != 0) {
            float lastT = t[i - 1];
            t0 = std::clamp((t0 - lastT) / (1 - lastT), 0.f, 1.f);
            t1 = std::clamp((t1 - lastT) / (1 - lastT), 0.f, 1.f);
        }
        chop_cubic_at_2(cur, d, t0, t1);
        memcpy(cur, d + 6, 4 * sizeof(Pt));
        d += 6;
    }
    if (i < n) {
        float tt = t[i];
        if (i != 0) {
            float lastT = t[i - 1];
            float div = (tt - lastT) / (1 - lastT);
            if (std::isnan(div)) div = 0;
            tt = std::clamp(div, 0.f, 1.f);
        }
        chop_cubic_at(cur, d, tt);
    }
}

// ---------- flatten_double_cubic_extrema (SkGeometry.cpp:686) ----------
// coords are interleaved (stride 2). Sets coords[4]=coords[8]=coords[6].
static void flatten_x(Pt* dst, int idx) {
    // offset by idx*3 points; flatten X at indices 2,3,4 → dst[idx*3+2].x = dst[idx*3+4].x = dst[idx*3+3].x
    dst[2].x = dst[4].x = dst[3].x;
}
static void flatten_y(Pt* dst) {
    dst[2].y = dst[4].y = dst[3].y;
}

static int chop_at_y_extrema(const Pt src[4], Pt dst[10]) {
    float t[2];
    int roots = find_cubic_extrema(src[0].y, src[1].y, src[2].y, src[3].y, t);
    chop_cubic_at_n(src, dst, t, roots);
    if (roots > 0) {
        // flatten Y at first chop boundary (points 2,3,4)
        dst[2].y = dst[4].y = dst[3].y;
        if (roots == 2) {
            dst[5].y = dst[7].y = dst[6].y;
        }
    }
    return roots;
}

static int chop_at_x_extrema(const Pt src[4], Pt dst[10]) {
    float t[2];
    int roots = find_cubic_extrema(src[0].x, src[1].x, src[2].x, src[3].x, t);
    chop_cubic_at_n(src, dst, t, roots);
    if (roots > 0) {
        dst[2].x = dst[4].x = dst[3].x;
        if (roots == 2) {
            dst[5].x = dst[7].x = dst[6].x;
        }
    }
    return roots;
}

// ---------- Simulate clipMonoCubic point/verb count ----------
// We use a tall clip so chop_cubic_in_Y is a no-op, simplifying.
struct ClipRect { float l, t, r, b; };

static void count_mono_cubic(const Pt src[4], const ClipRect& clip,
                             bool canCullRight, int* nVerbs, int* nPts,
                             bool tallClip) {
    Pt pts[4];
    // sort_increasing_Y
    if (src[0].y > src[3].y) {
        for (int i = 0; i < 4; i++) pts[i] = src[3 - i];
    } else {
        memcpy(pts, src, sizeof(pts));
    }
    // completely above/below
    if (pts[3].y <= clip.t || pts[0].y >= clip.b) return;
    // chop_cubic_in_Y: assume tall clip so it's a no-op (we ensure this in search)
    // X-sort
    if (pts[0].x > pts[3].x) {
        std::swap(pts[0], pts[3]);
        std::swap(pts[1], pts[2]);
    }
    // wholly left
    if (pts[3].x <= clip.l) { *nVerbs += 1; *nPts += 2; return; }
    // wholly right
    if (pts[0].x >= clip.r) {
        if (!canCullRight) { *nVerbs += 1; *nPts += 2; }
        return;
    }
    // partially left
    if (pts[0].x < clip.l) {
        *nVerbs += 1; *nPts += 2;
        // we don't need to actually chop for counting; pts[3].x unchanged
    }
    // partially right
    if (pts[3].x > clip.r) {
        *nVerbs += 2; *nPts += 6;  // cubic + vline
    } else {
        *nVerbs += 1; *nPts += 4;  // cubic only
    }
}

// ---------- Simulate full clipCubic ----------
struct Result {
    int totalPieces;
    int nVerbs;
    int nPts;
    int yPieces;
    int xPieces[3];
};

static bool too_big(const Pt p[4]) {
    float minx = p[0].x, maxx = p[0].x, miny = p[0].y, maxy = p[0].y;
    for (int i = 1; i < 4; i++) {
        minx = std::min(minx, p[i].x); maxx = std::max(maxx, p[i].x);
        miny = std::min(miny, p[i].y); maxy = std::max(maxy, p[i].y);
    }
    const float lim = 1 << 22;
    return minx < -lim || miny < -lim || maxx > lim || maxy > lim;
}

static Result simulate(const Pt src[4], const ClipRect& clip, bool canCullRight) {
    Result r{};
    if (too_big(src)) return r;
    // bounds check vs clip.t/b
    float miny = std::min({src[0].y, src[1].y, src[2].y, src[3].y});
    float maxy = std::max({src[0].y, src[1].y, src[2].y, src[3].y});
    if (!(maxy > clip.t && miny < clip.b)) return r;

    Pt monoY[10];
    int countY = chop_at_y_extrema(src, monoY);
    r.yPieces = countY + 1;
    for (int y = 0; y <= countY; y++) {
        Pt monoX[10];
        int countX = chop_at_x_extrema(&monoY[y * 3], monoX);
        r.xPieces[y] = countX + 1;
        for (int x = 0; x <= countX; x++) {
            r.totalPieces++;
            count_mono_cubic(&monoX[x * 3], clip, canCullRight, &r.nVerbs, &r.nPts, true);
        }
    }
    return r;
}

// ---------- Search ----------
static uint32_t f2u(float f) { uint32_t u; memcpy(&u, &f, 4); return u; }
static float u2f(uint32_t u) { float f; memcpy(&f, &u, 4); return f; }

int main(int argc, char** argv) {
    uint64_t seed = argc > 1 ? std::strtoull(argv[1], nullptr, 0) : 12345;
    uint64_t iters = argc > 2 ? std::strtoull(argv[2], nullptr, 0) : 200000000ULL;
    std::mt19937_64 rng(seed);

    int bestPts = 0, bestVerbs = 0, bestPieces = 0;
    Pt bestP[4]; ClipRect bestC{};

    // Strategies for generating candidates
    auto rand_float = [&](float lo, float hi) {
        std::uniform_real_distribution<float> d(lo, hi);
        return d(rng);
    };
    auto rand_bits = [&]() -> float {
        // random float bits, but reject NaN/Inf and >4M magnitude
        for (;;) {
            uint32_t u = (uint32_t)rng();
            float f = u2f(u);
            if (std::isfinite(f) && std::fabs(f) < (1 << 22)) return f;
        }
    };

    // We use a tall clip so chop_cubic_in_Y is no-op; narrow in X.
    // canCullToTheRight = false (path not convex) gives more output, but
    // either works for the spanning case.
    for (uint64_t i = 0; i < iters; i++) {
        Pt p[4];
        ClipRect c;
        int strat = rng() % 6;
        if (strat == 0) {
            // Y: 2 extrema (e.g., 0,1,-1,0 pattern scaled), X: oscillating
            float ys = rand_float(1.f, 1000.f);
            float xs = rand_float(100.f, 1e6f);
            p[0] = {-xs * rand_float(0.5f, 1.5f), 0};
            p[1] = {xs * rand_float(1.5f, 3.f), ys};
            p[2] = {-xs * rand_float(1.5f, 3.f), -ys};
            p[3] = {xs * rand_float(0.5f, 1.5f), 0};
            c = {-0.5f, -ys * 2, 0.5f, ys * 2};
        } else if (strat == 1) {
            // Random bit patterns (like fuzzer)
            for (int k = 0; k < 4; k++) p[k] = {rand_bits(), rand_bits()};
            float cl = rand_float(-100, 100), cw = rand_float(0.001f, 10.f);
            c = {cl, -4e6f, cl + cw, 4e6f};
        } else if (strat == 2) {
            // Subnormal/tiny mix (like crbug/698714)
            auto tiny = [&]() {
                uint32_t u = (rng() & 0x007FFFFF) | ((rng() % 4) << 23);
                if (rng() & 1) u |= 0x80000000;
                return u2f(u);
            };
            float big = rand_float(100.f, 300.f);
            for (int k = 0; k < 4; k++) {
                p[k].x = (rng() & 1) ? tiny() : ((rng() & 1) ? big : -big) + tiny();
                p[k].y = (rng() & 1) ? tiny() : ((rng() & 1) ? big : -big) + tiny();
            }
            c = {0, -1e6f, rand_float(0.01f, 65.f), 1e6f};
        } else if (strat == 3) {
            // Near-cancellation: A = d-a+3(b-c) ~ 0 in both X and Y
            // For Y: pick a,b,c then d = a - 3(b-c) + epsilon
            float ya = rand_float(-1000, 1000);
            float yb = rand_float(-1000, 1000);
            float yc = rand_float(-1000, 1000);
            float yd = ya - 3 * (yb - yc) + u2f((rng() & 0xFF) | 0x33800000) * ((rng()&1)?1:-1);
            float xa = rand_float(-1e5f, 1e5f);
            float xb = rand_float(-1e5f, 1e5f);
            float xc = rand_float(-1e5f, 1e5f);
            float xd = xa - 3 * (xb - xc) + u2f((rng() & 0xFF) | 0x33800000) * ((rng()&1)?1:-1);
            p[0]={xa,ya}; p[1]={xb,yb}; p[2]={xc,yc}; p[3]={xd,yd};
            float cl = rand_float(-10, 10);
            c = {cl, -1e6f, cl + rand_float(0.001f, 1.f), 1e6f};
        } else if (strat == 4) {
            // Mutate best
            if (bestPts > 0) {
                memcpy(p, bestP, sizeof(p));
                c = bestC;
                int idx = rng() % 8;
                float* f = &p[0].x + idx;
                uint32_t u = f2u(*f);
                int shift = rng() % 24;
                u ^= (1u << shift);
                *f = u2f(u);
                if (!std::isfinite(*f) || std::fabs(*f) >= (1 << 22)) continue;
            } else continue;
        } else {
            // Y has 2 close extrema, X oscillates wildly
            float scale = rand_float(1e3f, 4e6f) * 0.9f;
            // Y: cubic with 2 extrema near each other
            float t1 = rand_float(0.2f, 0.5f);
            float t2 = rand_float(0.5f, 0.8f);
            // Hard to invert; just random
            for (int k=0;k<4;k++) p[k].y = rand_float(-100,100);
            p[0].x = -scale; p[1].x = scale*3; p[2].x = -scale*3; p[3].x = scale;
            // bit-twiddle X
            for (int k=0;k<4;k++) {
                uint32_t u = f2u(p[k].x);
                u += (int)(rng()%21)-10;
                p[k].x = u2f(u);
            }
            float cl = rand_float(-5,5);
            c = {cl, -1e6f, cl+rand_float(0.001f,1.f), 1e6f};
        }

        Result r = simulate(p, c, false);
        if (r.nPts > bestPts || r.totalPieces > bestPieces) {
            if (r.nPts > bestPts || (r.nPts == bestPts && r.totalPieces > bestPieces)) {
                bestPts = r.nPts; bestVerbs = r.nVerbs; bestPieces = r.totalPieces;
                memcpy(bestP, p, sizeof(p)); bestC = c;
                fprintf(stderr,
                    "[%llu] pieces=%d (Y=%d X=%d,%d,%d) verbs=%d pts=%d  "
                    "P0=(%a,%a) P1=(%a,%a) P2=(%a,%a) P3=(%a,%a) clip=[%a,%a,%a,%a]\n",
                    (unsigned long long)i, r.totalPieces, r.yPieces,
                    r.xPieces[0], r.xPieces[1], r.xPieces[2],
                    r.nVerbs, r.nPts,
                    p[0].x,p[0].y,p[1].x,p[1].y,p[2].x,p[2].y,p[3].x,p[3].y,
                    c.l,c.t,c.r,c.b);
            }
        }
        if (r.nPts > 54 || r.nVerbs > 18) {
            printf("OVERFLOW FOUND! pieces=%d verbs=%d pts=%d\n",
                   r.totalPieces, r.nVerbs, r.nPts);
            printf("P0 = {0x%08x, 0x%08x}  // %g, %g\n", f2u(p[0].x), f2u(p[0].y), p[0].x, p[0].y);
            printf("P1 = {0x%08x, 0x%08x}  // %g, %g\n", f2u(p[1].x), f2u(p[1].y), p[1].x, p[1].y);
            printf("P2 = {0x%08x, 0x%08x}  // %g, %g\n", f2u(p[2].x), f2u(p[2].y), p[2].x, p[2].y);
            printf("P3 = {0x%08x, 0x%08x}  // %g, %g\n", f2u(p[3].x), f2u(p[3].y), p[3].x, p[3].y);
            printf("clip = {%a, %a, %a, %a}\n", c.l, c.t, c.r, c.b);
            return 0;
        }
    }
    fprintf(stderr, "Done. Best: pieces=%d verbs=%d pts=%d\n", bestPieces, bestVerbs, bestPts);
    return 1;
}