//CREATED BY: mystical, 2024
/*
This is a search program named multiperiodphotonsearcher v0.1, which can be used to find
photons with arbitrarily large periods. It can be set up to search for a specific period 
(using the feature of forbidden periods) or it can be set up to search all periods less 
than max period. Interestingly, the time taken to search should only be linear in max period. 

The main idea used here is similar to that of photonsrc by LaundryPizza03 (https://conwaylife.com/forums/viewtopic.php?f=9&t=4649),
but I came up with most of the ideas used in this program independently.

The code is slightly user unfriendly for now, and you'll have to set some variables like max period 
& width in the source code itself. some parameters like forbidden periods and the rulespace can 
be specified as command line arguments.

You can compile this program on linux it using the command 

    >>> g++ multi_period_photon_searcherv0.1.cpp -O3 -march=native -Wall -Wextra -o mpps

and run it using (example search in B27(8)/S24(8) with forbidden periods 4 and 6)

    >>> ./mpps B27/S24 B278/S248 f 4 6

Certain parameters need to be set directly in the part of the code where it says 'SET PARAMETERS HERE'
*/

// see more of my ca related programs at https://github.com/Aaron-d-k
// This program is not yet on my github, but I will put it there when I get the time

#include <bits/stdc++.h>
#include <immintrin.h>

using namespace std;
using ll = long long; 
using pll = pair<ll,ll>;
#define e0 first
#define e1 second
constexpr ll INFTY = 2e18;

// SET PARAMETERS HERE
// vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
    constexpr int MAX_PERIOD = 50;
    constexpr size_t W = 10;
    constexpr ll MAX_ITER = INFTY;
    constexpr ll VERBOSITY = 2;// can be -1, 0, 1 or 2
    #undef EVEN // for even symmetric searches, use #define EVEN
    enum {BFS,DFS}; constexpr auto search_style = BFS; // dfs usually doesn't do very well
    constexpr int CACHE_TILL_PERIOD = 100; // keeping this very high doesn't impact performance much
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

static_assert(W<=16, "Sorry, widths above 16 not supported yet.");

//TODO: prune dead ends, period reduction checking



struct node;
using otca_rule = uint32_t;

size_t create_tree_node(array<uint64_t,MAX_PERIOD+1>& rows, int P, size_t parent);
void set_rule(otca_rule rule);
constexpr uint64_t rulelookup(uint64_t s, uint64_t v);
bool evolve(int gens);
string node_to_rle(size_t nodeidx);
int find_depth(size_t currnode);


void tryextend(size_t nodeidx, deque<size_t>& extensions);
void rulesearch(ll verbosity, ll forbidhighmultiples=INFTY, ll maxiter=INFTY, ll minprintperiod=2);

otca_rule parse_rule(const string& rule);
string int_to_rule(otca_rule rule);
vector<otca_rule> getrulespace(otca_rule minrule, otca_rule maxrule);


constexpr uint64_t CELL_POSN = ((1ull<<(4ull*W))-1ull)/0xFull;

// All the partials found are represented as a large tree
// This does mean that the memory usage keeps steadily increasing over time
// till the computer ooms or the search ends.
// Pruning dead ends, when implemented, should reduce memory usage.
vector<node> tree;

vector<int> forbidden_periods;
set<int> seen_periods;

struct fingerprint
{
    vector<uint16_t> data;
    bool operator<(const fingerprint& other) const
    {
        return data < other.data;
    }

    bool operator==(const fingerprint& other) const
    {
        return data == other.data;
    }

    void or_at(uint16_t d, size_t pos)
    {
        data[pos/16] |= d << (pos%16);
        if ((pos/16)+1<data.size()) data[(pos/16)+1] |= d >> (16-(pos%16));
    }
};

// credits: http://graphics.stanford.edu/~seander/bithacks.html
uint16_t reverselowbits(uint16_t v, size_t num_bits)
{
    // swap odd and even bits
    v = ((v >> 1) & 0x5555) | ((v & 0x5555) << 1);
    // swap consecutive pairs
    v = ((v >> 2) & 0x3333) | ((v & 0x3333) << 2);
    // swap nibbles ... 
    v = ((v >> 4) & 0x0F0F) | ((v & 0x0F0F) << 4);
    // swap bytes
    v = ((v >> 8) & 0x00FF) | ((v & 0x00FF) << 8);
    return v>>(16-num_bits);
}

struct node
{
    size_t parent;
    vector<uint16_t> row;

    int period() const
    {
        return row.size();
    }

    // converts the compressed stored representation (1 bit per cell) into an expanded 
    // form with 4 bits per cell.
    void expand(array<uint64_t,MAX_PERIOD>& r1, array<uint64_t,MAX_PERIOD>& r2) const
    {
        const vector<uint16_t>& prevrow=tree[parent].row;
        for (int i = 0; i < MAX_PERIOD; i++)
        {
            r1[i]=_pdep_u64(uint64_t(prevrow[i%(prevrow.size())]),CELL_POSN);
        }

        for (int i = 0; i < MAX_PERIOD; i++)
        {
            r2[i]=_pdep_u64(uint64_t(row[i%(row.size())]),CELL_POSN);
        }
    }

    bool is_terminated() const
    {
        for (uint16_t s : row)
        {
            if (s!=0) return false;
        }
        for (uint16_t s : tree[parent].row)
        {
            if (s!=0) return false;
        }
        return true;
    }

    // used to generate a unique number for each set of 2 rows to help avoid repeating partials
    fingerprint getfingerprint() const
    {
        const size_t P1 = period();
        const size_t P2 = tree[parent].period();
        const vector<uint16_t>& prevrow=tree[parent].row;

        assert(P1>=P2 && P1%P2==0);
        assert(P1<=CACHE_TILL_PERIOD);

        size_t minstart = 0;
        for (size_t i = 1; i < P1; i++)
        {
            if (pair<uint16_t,uint16_t>{ row[i], prevrow[i%P2] }
                < pair<uint16_t,uint16_t>{ row[minstart], prevrow[minstart%P2] }
            ) {
                minstart=i;
            }
        }

        bool saveflipped=false;
#ifndef EVEN
        if (reverselowbits(row[minstart],W)<row[minstart])
        {
           saveflipped = true;
        }
#endif

        fingerprint outp = {vector<uint16_t>(1+(P1*W*2)/16,uint16_t(0)), };
        for (size_t i = 0; i < P1; i++)
        {
            if (saveflipped)
            {
                outp.or_at(reverselowbits(row[(i+minstart)%P1],W), i*W*2);
                outp.or_at(reverselowbits(prevrow[(i+minstart)%P2],W), i*W*2+W);
            }
            else
            {
                outp.or_at(row[(i+minstart)%P1], i*W*2);
                outp.or_at(prevrow[(i+minstart)%P2], i*W*2+W);
            }
        }
        return outp;
    }
};

size_t create_tree_node(array<uint64_t,MAX_PERIOD+1>& rows, int P, size_t parent)
{
    tree.push_back(node{parent,vector<uint16_t>(P)});
    for (int i = 0; i < P; i++)
    {
        tree.back().row[i] = _pext_u64(rows[i], CELL_POSN);
    }
    
    return tree.size()-1;
}

array<bool,18> birth_survival;
string rule_name;

void set_rule(otca_rule rule)
{
    birth_survival.fill(0);
    for (ll i = 0; i < 18; i++)
    {
        if ((rule>>i)&1) {
            birth_survival[i]=1;
        }
    }
    rule_name=int_to_rule(rule);
}



constexpr uint64_t rulelookup(uint64_t s, uint64_t v)
{
    return birth_survival[s+v*9];
}

array<uint64_t,MAX_PERIOD> s1;
array<uint64_t,MAX_PERIOD> s2; 
array<uint64_t,MAX_PERIOD+1> s3;

// evolves a candidate extension (s[3][0]) of the photon for specified number of gens.
// returns false if it fails to work in some way (eg it escapes the bounding box)
bool evolve(int gens)
{
    constexpr uint64_t cellmask = 0xF;

// Profiling says that 95% of the time is being spent here, but Im not too sure
// how I can make it faster
    for (int i = 0; i < gens; i++)
    {
        uint64_t sum = s1[i]+(s1[i]<<4)+(s1[i]>>4)
        +(s2[i]<<4)+(s2[i]>>4)
        +s3[i]+(s3[i]<<4)+(s3[i]>>4)
        #ifdef EVEN
        +((s1[i]+s2[i]+s3[i])&cellmask)
        #endif
        ;

        s3[i+1]=0;
        for (int j = 0; j < int(W); j++)
        {
            s3[i+1] |= rulelookup(cellmask&(sum>>(j*4)),cellmask&(s2[i]>>(j*4)))<<(j*4);
        }

        // check if pattern has breached the bounding box
        #ifndef EVEN
        if (0!=rulelookup((s1[i]+s2[i]+s3[i])&cellmask,0)) return false;
        #endif

        if (0!=rulelookup(((s1[i]+s2[i]+s3[i])>>(W*4-4))&cellmask,0)) return false;
    }

    return true;
}

string node_to_rle(size_t nodeidx)
{
    ostringstream ss;
    if (tree[nodeidx].is_terminated()) 
    {
        ss << "#C +++++++ p"<< tree[nodeidx].period()<<" ship +++++++\n";
    }
    else 
    {
        ss << "#C ---- p"<<tree[nodeidx].period()<<" partial ----\n";
    }

    ss << "x = 0, y = 0, rule = " << rule_name << "\n";
    while (nodeidx!=0)
    {
        uint64_t currrow = tree[nodeidx].row[0];
        for (size_t b = 0; b < W; b++)
        {
            ss << (((currrow>>b)&1)?'o':'b');
        }
        nodeidx=tree[nodeidx].parent;
        ss << ((nodeidx==0)?"!\n":"$\n");
    }
    return ss.str();
}

int find_depth(size_t currnode)
{
    int depth=0;
    while (currnode!=0)
    {
        currnode=tree[currnode].parent;
        depth++;
    }
    return depth;
}

array<set<fingerprint>,CACHE_TILL_PERIOD+1> lowperiodhashes;

void tryextend(size_t nodeidx, deque<size_t>& extensions)
{
    node& currnode = tree[nodeidx];
    int P = currnode.period();
    for (int p_f : forbidden_periods) 
    {
        if (P%p_f==0) return;
    }
    currnode.expand(s1,s2);
    for (uint64_t i = 0; i < 1ull<<W; i++)
    {
        s3[0] =  _pdep_u64(i,CELL_POSN);
        if (!evolve(MAX_PERIOD-(MAX_PERIOD%P))) continue;

        for (int tryperiod = P; tryperiod <= MAX_PERIOD; tryperiod+=P)
        {
            if (s3[tryperiod]==s3[0])
            {
                for (int p_f : forbidden_periods) 
                {
                    if (tryperiod%p_f==0) goto next_continuation_trial;
                } 
                
                size_t newnode = create_tree_node(s3, tryperiod, nodeidx);

                //Note that the duplicate checking means all but first ship of a certain period are rejected...
                //shouldn't be too hard to change if necessary.
                if (tryperiod<=CACHE_TILL_PERIOD)
                {
                    fingerprint f = tree.back().getfingerprint();
                    if (!lowperiodhashes[tryperiod].insert(f).e1)
                    {
                        tree.pop_back();
                        goto next_continuation_trial;
                    }
                }
                
                extensions.push_back(newnode);
                goto next_continuation_trial;
            }
        }
        next_continuation_trial:; 
    }
}

void rulesearch(ll verbosity, ll forbidhighmultiples, ll maxiter, ll minprintperiod)
{
    seen_periods.clear();
    tree.clear();
    lowperiodhashes={};
    tree.push_back({0,vector<uint16_t>{0,}});

    deque<size_t> extensions;
    extensions.push_back(0);

    ll iter=0;
    ll recorddepth=3;
    auto lastprinttime=chrono::steady_clock::now();
    while (!extensions.empty() && iter++ < maxiter)
    {
        size_t curr;
        if (search_style==BFS)
        {
            curr = extensions.front();
            extensions.pop_front();
        }
        else if (search_style==DFS)
        {
            curr = extensions.back();
            extensions.pop_back();
        }
        else throw runtime_error("not a valid search style yet...");

        if (verbosity>0 && iter%100==0) 
        {
            using namespace std::literals::chrono_literals;
            auto currtime= chrono::steady_clock::now();
            if ((currtime-lastprinttime)>1s) 
            {
                int depth = find_depth(curr);
                clog << "\rdepth = " << depth << ", queue = " << extensions.size() << ", tree size = " << tree.size() << ", iter = " << iter;
                if (verbosity>1 && recorddepth<depth) 
                {
                    clog << "\n" << node_to_rle(curr) << "\n"; 
                    recorddepth=depth; 
                }
                clog.flush();
                lastprinttime=currtime;
            }
        }

        if (curr!=0 && tree[curr].is_terminated())
        {
            // don't want to show ships from periods shown already
            if (seen_periods.insert(tree[curr].period()).e1) 
            {
                if (minprintperiod<=tree[curr].period()) 
                {
                    cout << node_to_rle(curr);
                    cout.flush();
                }
                if (forbidhighmultiples <= tree[curr].period())
                {
                    forbidden_periods.push_back(tree[curr].period());
                }
            }
        }
        else
        {
            tryextend(curr,extensions);
        }
    }
}

otca_rule parse_rule(const string& rule)
{
    otca_rule outp=0;
    if (0==rule.size() || rule[0]!='B') throw runtime_error("rule format error:1");

    ll idx=1;
    while (idx<ll(rule.size()) && '0'<=rule[idx]&&rule[idx]<='8')
    {
        outp |= 1<<(rule[idx]-'0');
        idx++;
    }

    if (idx+1>=ll(rule.size()) || rule[idx]!='/' || rule[idx+1]!='S') throw runtime_error("rule format error:2");
    idx+=2;

    while (idx<ll(rule.size()) && '0'<=rule[idx]&&rule[idx]<='8')
    {
        outp |= 1<<(rule[idx]-'0'+9);
        idx++;
    }

    if (idx!=ll(rule.size())) throw runtime_error("rule format error:3");

    return outp;
}

string int_to_rule(otca_rule rule)
{
    string o="B";
    for (ll i = 0; i < 9; i++)
    {
        if ((rule>>i)&1) {
            o+= i+'0';
        }
    }
    o+="/S";
    for (ll i = 0; i < 9; i++)
    {
        if ((rule>>(i+9))&1) {
            o+= i+'0';
        }
    }
    return o;
}

vector<otca_rule> getrulespace(otca_rule minrule, otca_rule maxrule)
{
    if ((minrule&~maxrule)!=0) throw runtime_error("invalid rulespace");

    otca_rule variabletrans = minrule^maxrule;
    ll numtrans = __builtin_popcount(variabletrans);
    vector<otca_rule> outp(1<<numtrans,0);
    for (otca_rule i = 0; i < (1ull<<numtrans); i++)
    {
        outp[i]= minrule|_pdep_u32(i, variabletrans);
    }
    return outp;
}

int main(int argc, char *argv[])
{
    cin.tie(nullptr);
    ios_base::sync_with_stdio(false);

   if (argc<3) throw runtime_error("too few args");

    otca_rule minrule = parse_rule(argv[1]);
    otca_rule maxrule = parse_rule(argv[2]);

    if (argc >= 4 && argv[3][0]=='f')
    {
        for (int i = 4; i < argc; i++)
        {
            forbidden_periods.push_back(stoi(argv[i]));
        }
    }

    auto lastprinttime=chrono::steady_clock::now()-2s;
    for (otca_rule r : getrulespace(minrule,maxrule))
    {
        set_rule(r);
        using namespace std::literals::chrono_literals;
        auto currtime= chrono::steady_clock::now();
        if (VERBOSITY>=1 || (VERBOSITY==0&&(currtime-lastprinttime)>1s))
        {
            clog<<"\nSearching rule " << rule_name << "..." << endl;
            lastprinttime=currtime;
        }
        rulesearch(VERBOSITY,INFTY,MAX_ITER,2); 
    }
    cout << endl;
}