// Isolate the upper-13 optimization on identical supplied lower35 planes.
// The full seed pipeline lives in 08_seed_cpu.cpp / 08_seed_cuda.cu.
#include "common.hpp"
#include "../research/slime_solver_core.hpp"
#include <cuda_runtime.h>
#include <array>

namespace sr=slime_research;
void check(cudaError_t e) { if(e!=cudaSuccess) throw std::runtime_error(cudaGetErrorString(e)); }
struct Solution {
    uint64_t lower35,seed;
    bool operator<(const Solution& o) const { return std::tie(lower35,seed)<std::tie(o.lower35,o.seed); }
    bool operator==(const Solution& o) const { return lower35==o.lower35 && seed==o.seed; }
};
struct Emitter {
    uint64_t lower35,capacity;
    Solution* out;
    unsigned long long* count;
    __device__ void operator()(uint64_t seed) {
        const auto slot=atomicAdd(count,1ULL);
        if(slot<capacity) out[slot]={lower35,seed};
    }
};
template<bool analytic> __global__ void upperKernel(
    uint64_t lower35,uint64_t planes,const sr::CellTerm* cells,int n,
    Solution* out,unsigned long long* count,uint64_t capacity
) {
    const uint64_t stride=uint64_t(gridDim.x)*blockDim.x;
    const uint64_t total=analytic?planes:planes*8192;
    for(uint64_t i=uint64_t(blockIdx.x)*blockDim.x+threadIdx.x;i<total;i+=stride) {
        const uint64_t plane=analytic?i:i/8192;
        const uint64_t lower=(lower35+plane)&((1ULL<<35)-1);
        Emitter emit{lower,capacity,out,count};
        if constexpr(analytic) sr::solveUpper13(lower,cells,n,emit);
        else {
            const uint64_t seed=lower+((i%8192)<<35);
            // Same first-draw candidate model as analytic, then exact verification.
            bool candidate=true;
            for(int c=0;c<n;++c) {
                if((sr::nextStateFromTerm(seed,cells[c].term)>>17)%10!=0) { candidate=false;break; }
            }
            if(candidate && sr::verifyPattern(seed,cells,n)) emit(seed);
        }
    }
}

int main(int argc,char** argv) {
    try {
        uint64_t lower=16420645998ULL,planes=4096;
        int repeats=5;
        for(int i=1;i<argc;++i) {
            const std::string key=argv[i];
            if(++i>=argc) throw std::runtime_error("missing value");
            if(key=="--lower35") lower=show::positive(argv[i]);
            else if(key=="--planes") planes=show::positive(argv[i]);
            else if(key=="--repeats") repeats=std::stoi(argv[i]);
            else throw std::runtime_error("unknown option");
        }
        if(!planes || planes>1000000 || repeats<1 || lower>=(1ULL<<35))
            throw std::runtime_error("invalid planes/repeats/lower35");
        std::array<sr::CellTerm,25> terms;
        int cell=0;
        for(int dx=0;dx<5;++dx) for(int dz=0;dz<5;++dz)
            terms[cell++]={sr::coordinateTerm(51069+dx,419124+dz),sr::kSlime};
        sr::CellTerm* deviceTerms=nullptr;
        Solution* out=nullptr;
        unsigned long long* count=nullptr;
        constexpr uint64_t capacity=100000;
        check(cudaMalloc(&deviceTerms,sizeof(terms))); check(cudaMalloc(&out,capacity*sizeof(Solution)));
        check(cudaMalloc(&count,sizeof(unsigned long long)));
        check(cudaMemcpy(deviceTerms,terms.data(),sizeof(terms),cudaMemcpyHostToDevice));
        cudaEvent_t begin,end; check(cudaEventCreate(&begin)); check(cudaEventCreate(&end));
        std::array<std::vector<double>,2> times;
        std::vector<Solution> expected;
        bool first=true;
        int dev=0; cudaDeviceProp prop{};
        check(cudaGetDevice(&dev)); check(cudaGetDeviceProperties(&prop,dev));
        for(int r=-1;r<repeats;++r) for(int j=0;j<2;++j) {
            const int mode=(j+std::max(0,r))%2;
            check(cudaMemset(count,0,sizeof(unsigned long long)));
            check(cudaEventRecord(begin));
            const uint64_t work=mode?planes:planes*8192;
            const int blocks=int(std::min<uint64_t>((work+255)/256,prop.multiProcessorCount*32));
            if(mode) upperKernel<true><<<blocks,256>>>(lower,planes,deviceTerms,25,out,count,capacity);
            else upperKernel<false><<<blocks,256>>>(lower,planes,deviceTerms,25,out,count,capacity);
            check(cudaGetLastError()); check(cudaEventRecord(end)); check(cudaEventSynchronize(end));
            float elapsed=0; check(cudaEventElapsedTime(&elapsed,begin,end));
            unsigned long long size=0; check(cudaMemcpy(&size,count,sizeof(size),cudaMemcpyDeviceToHost));
            if(size>capacity) throw std::runtime_error("upper results truncated");
            std::vector<Solution> hits(size);
            if(size) check(cudaMemcpy(hits.data(),out,size*sizeof(Solution),cudaMemcpyDeviceToHost));
            std::sort(hits.begin(),hits.end());
            if(first) { expected=hits;first=false; }
            else if(hits!=expected) throw std::runtime_error("enumerated/analytic upper results differ");
            if(r>=0) times[mode].push_back(elapsed);
        }
        std::cout.precision(10);
        std::cout << "{\"lower35_start\":" << lower << ",\"planes\":" << planes
                  << ",\"candidate_model\":\"first_draw\",\"comparison_equal\":true,\"rows\":[";
        for(int m=0;m<2;++m) {
            if(m) std::cout << ',';
            std::cout << "{\"mode\":\"" << (m?"analytic_upper13":"enumerate_upper13")
                      << "\",\"kernel_samples_ms\":"; show::array(times[m]);
            std::cout << ",\"kernel_median_ms\":" << show::median(times[m]) << '}';
        }
        std::cout << "],\"solutions\":[";
        for(size_t i=0;i<expected.size();++i) {
            if(i) std::cout << ',';
            std::cout << '[' << expected[i].lower35 << ',' << expected[i].seed << ']';
        }
        std::cout << "]}\n";
        cudaEventDestroy(begin);cudaEventDestroy(end);cudaFree(count);cudaFree(out);cudaFree(deviceTerms);
    } catch(const std::exception& e) { std::cerr << e.what() << '\n';return 1; }
}