Integer to polynomial

#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <cstring>
#include <iostream>
#include <cmath>
#include <cstdint>
#include <vector>
#include <limits.h>
#include <algorithm>
#include <tuple>

#include <thread>
#include <mutex>
#include <condition_variable>
#include <queue>
#include <functional>


//Math functions
int64_t intPow(int64_t base, int64_t expo){
    int64_t result = 1;
    while (expo){
        if (expo & 1) result *= base;
        base *= base;
        expo >>= 1;
    }
    return result;
}

int64_t digitLen(int64_t n){
    n = abs(n);
    if (n < 10) return 1;
    if (n < 100) return 2;
    if (n < 1000) return 3;
    if (n < 10000) return 4;
    if (n < 100000) return 5;
    if (n < 1000000) return 6;
    if (n < 10000000) return 7;
    if (n < 100000000) return 8;
    if (n < 1000000000) return 9;
    if (n < 10000000000) return 10;
    if (n < 100000000000) return 11;
    if (n < 1000000000000) return 12;
    if (n < 10000000000000) return 13;
    if (n < 100000000000000) return 14;
    if (n < 1000000000000000) return 15;
    if (n < 10000000000000000) return 16;
    if (n < 100000000000000000) return 17;
    if (n < 1000000000000000000) return 18;
    return 19;
}

int64_t int64Hash(const char* cstr){
    int64_t result = 0;
    while(cstr[0]){
        result = (result<<8) | (result>>56);
        result |= (uint8_t)cstr[0];
        cstr++;
    };
    return result;
}

void reverseCstr(char* cstr){
    size_t l = strlen(cstr)-1;
    for(size_t i=0; i<=l/2; i++){
        char buf=cstr[i];
        cstr[i]=cstr[l-i];
        cstr[l-i]=buf;
    }
}

void int64ToChar(int64_t value, char* cstr, size_t byteSize) {
    if(cstr==NULL){ return; }
    memset(cstr, 0, byteSize);
    size_t maxSize = byteSize-1;
    maxSize = maxSize < sizeof(int64_t) ? maxSize : sizeof(int64_t);
    memcpy(cstr, &value, maxSize);
    if(cstr[0]=='\0'){ return; }
    reverseCstr(cstr);
}

template <class T>
bool isPrime(T n){
    if(n<=1){return false;}
    if(n==2){return true;}
    if(n%2==0){return false;}
    T max = sqrt(n);
    for(T i=3; i <= max; i+=2)
        if (n % i == 0)
            return false;
    return true;
}

bool is_numeric(char const *string)
{
    return std::all_of(string, string+strlen(string),
                       [](unsigned char c) { return ::isdigit(c); });
}



class ThreadPool {
public:
    ThreadPool(size_t numThreads);
    ~ThreadPool();

    template<class F>
    void enqueue(F&& f);

    void waitUntilFinished();

private:
    std::vector<std::thread> workers;
    std::queue<std::function<void()>> tasks;

    std::mutex queueMutex;
    std::condition_variable condition;
    std::condition_variable finishedCondition;
    bool stop;
    int tasksRemaining = 0;
};

ThreadPool::ThreadPool(size_t numThreads) : stop(false) {
    for (size_t i = 0; i < numThreads; ++i) {
        workers.emplace_back([this] {
            while (true) {
                std::function<void()> task;
                {
                    std::unique_lock<std::mutex> lock(this->queueMutex);
                    this->condition.wait(lock, [this] { return this->stop || !this->tasks.empty(); });
                    if (this->stop && this->tasks.empty()) return;
                    task = std::move(this->tasks.front());
                    this->tasks.pop();
                }

                task();

                {
                    std::lock_guard<std::mutex> lock(this->queueMutex);
                    tasksRemaining--;
                    if (tasksRemaining == 0) {
                        finishedCondition.notify_all();
                    }
                }
            }
        });
    }
}

ThreadPool::~ThreadPool() {
    {
        std::unique_lock<std::mutex> lock(queueMutex);
        stop = true;
    }
    condition.notify_all();
    for (std::thread &worker : workers) {
        worker.join();
    }
}

template<class F>
void ThreadPool::enqueue(F&& f) {
    {
        std::unique_lock<std::mutex> lock(queueMutex);
        tasks.emplace(std::forward<F>(f));
        tasksRemaining++;
    }
    condition.notify_one();
}

void ThreadPool::waitUntilFinished() {
    std::unique_lock<std::mutex> lock(queueMutex);
    finishedCondition.wait(lock, [this] { return tasksRemaining == 0; });
}









//Polynomial
struct Polynom {

    //Term
    struct Term {
        int64_t mult;
        int64_t base;
        int64_t expo;

        Term(int64_t m = 1, int64_t b = 1, int64_t e = 1) : mult(m), base(b), expo(e){}

        bool operator==(const Term& other) const;
        std::string getString() const;
        int64_t value() const;
        int64_t symbolCount() const;
    };

    struct SearchCache {
        std::vector<Polynom> *storage; //Pointer to unique polynom vector in thread scope
        const Polynom *curPolynom; //Pointer of mainList[i]
        int maxSymbols;
        int maxBranches;
        int termCount;
        int64_t maxCoef;
        int64_t maxExpo;

        SearchCache(
            std::vector<Polynom> *localeStorage = nullptr,
            const Polynom *polynomCur = nullptr,
            int symbolsMax = 8,
            int branchesMax = 8,
            int countTerm = 3,
            int64_t coefMax = 10000,
            int64_t expoMax = 62
        ):
            storage(localeStorage),
            curPolynom(curPolynom),
            maxSymbols(symbolsMax),
            maxBranches(branchesMax),
            termCount(countTerm),
            maxCoef(coefMax),
            maxExpo(expoMax)
        {}
    };

    std::vector<Term> terms; //polynomial terms, form a*b^c
    int64_t target;
    int64_t remainder; //How much needs to be added to reach the target
    bool validEnd=1;

    //Not enough redundatn calls to make caching useful
    //  mutable int64_t cacheSymCount = -1;
    //  mutable std::string cacheString;

    void add(Term term){
        remainder+=term.value();
        terms.push_back(term);
    }

    void terminate(){
        if(remainder==0){return;}
        terms.push_back({-1*remainder,1,1});
        remainder=0;
    }

    std::string getString() const;
    int64_t value() const;
    int symbolCount() const;

    int64_t scoreTerm(const Term& term, int64_t remainder) const;
    void sortTerms(SearchCache cfg, std::vector<Term> *candiTerms) const;
    int searchTerms(SearchCache cfg) const;
    bool isPrimesOnly() const;
};

std::string Polynom::getString() const {
//  if(cacheString!=""){return cacheString;}
    std::string str;
    for(int i=0;i<terms.size();i++){
        if(i>0 && (terms[i].mult>=0) ){str+= "+";} //positive sign for next term
        str+= terms[i].getString();
    }
//  cacheString=str;
    return str;
}

int Polynom::symbolCount() const {
//  if(cacheSymCount!=-1){return cacheSymCount;}

    int count=0;
    for(int i=0; i<terms.size();i++){
        count+=terms[i].symbolCount();
    }
    count+=terms.size()-1; // sign and digit
//  cacheSymCount=count;
    return count;
}

int64_t Polynom::value() const {
    int64_t sum=0;
    for(int i=0;i<terms.size();i++){
        sum+=terms[i].value();
    }
    return sum;
}

bool Polynom::isPrimesOnly() const {
    for(int i=0;i<terms.size();i++){
        Term term = terms[i];
        if(
            ( !isPrime(term.mult) && term.mult!=1 ) || //if any term isn't prime or 1
            ( !isPrime(term.base) && term.base!=1 ) ||
            ( !isPrime(term.expo) && term.expo!=1 )
        ){
            return 0;
        };
    }
    return 1;
}


//Polynom::Term::
bool Polynom::Term::operator==(const Term& other) const {

    if (
        intPow(other.base,other.expo)==1 && //ignore base^expo==1
        intPow(base,expo)==1
    ) {
        return (mult == other.mult);
    }
    return (mult == other.mult) && (base == other.base) && (expo == other.expo);
}

std::string Polynom::Term::getString() const {
    std::string str="";
    if(
        mult!=1 ||
        ( mult==1 && (mult==1 && base ==1) )
    ){
        str+= std::to_string(mult);
    }
    if(base==1){return str;} //a

    if(mult!=1){
        str+= "*";
    }
    str+= std::to_string(base);
    if(expo==1){return str;} //a*b

    str+= "^" + std::to_string(expo);
    return str; //a*b^c
}

int64_t Polynom::Term::value() const{
    return mult*intPow(base,expo);
}

int64_t Polynom::Term::symbolCount() const { //remove trivial components and return term a*b^c as a string
    bool trivExpo = intPow(base,expo)==1;

    int symCount = 0;
    if(mult!=1){
        symCount+=digitLen(mult);
    }
    if( trivExpo ){
        if(symCount==0){return 1;} // 1*1^1 (1 digit)
        return symCount; //a
    }

    symCount += mult!=1; //multiply sign *
    symCount+=digitLen(base);
    if(expo!=1){
        symCount+=digitLen(expo);
    }

    return symCount; //b or b^c or a*b^c
}





//Polynomial::
//Fitness value for new term. Higher is better
int64_t Polynom::scoreTerm(const Term& term, int64_t remainder) const {
    int64_t termLength = term.symbolCount(); //how many chars the term has

    int64_t reduceWeight = 1+termLength; //Each new term adds plus or minus sign, so we prefer remainder reduction
    int64_t lengthWeight = termLength;

    int64_t newRemainder = remainder + term.value();
    int64_t charsReduced = digitLen(remainder) - digitLen(newRemainder); //how many chars the term reduces from remainder

    //weight
    int64_t score = 0;
    score += charsReduced * reduceWeight; //Higher = better
    score -= termLength * lengthWeight; //Higher = worse
    return score;
}

void Polynom::sortTerms(SearchCache cfg, std::vector<Term> *candiTerms) const {
    //Precompute scores
    std::vector<std::pair<Term, int>> scoredTerms;
    for (const auto& term : *candiTerms){
        scoredTerms.emplace_back(term, scoreTerm(term, cfg.curPolynom->remainder));
    }

    //sort highest score last
    std::sort(scoredTerms.begin(), scoredTerms.end(), [](const auto& a, const auto& b){
        return a.second > b.second;
    });

    Term prevTerm;
    candiTerms->clear();
    for (const auto& pair : scoredTerms){
        //Don't add terms with same value
        if( !(pair.first==prevTerm) ){
            candiTerms->push_back(pair.first);
            prevTerm=pair.first;
        }
    }
}

//Find candidates for polynomial term close to remainder
int Polynom::searchTerms(SearchCache cfg) const {
    if(cfg.curPolynom->remainder==0){return 0;}

    std::vector<Term> candiTerms; //candidate terms
    Term newTerm = {0};
    int64_t sign = cfg.curPolynom->remainder < 0 ? 1 : -1; //if remainder is negative, next coefficient should be positive

    for(int64_t base=2;base<=cfg.maxCoef;base++){
        for(int64_t expo=1;expo<=cfg.maxExpo;expo++){

            int64_t power = intPow(base,expo);
            if( power<=0 || ( power>cfg.maxCoef+abs(cfg.curPolynom->remainder) ) ){break;} //ignore overflown powers

            //Scale the power near reminder
            int64_t mult = (abs(cfg.curPolynom->remainder)+power/2)/power;
            if(mult<=0){mult=1;}
            newTerm={ sign*mult, base, expo };

            int64_t newTermDiff = abs(cfg.curPolynom->remainder+newTerm.value()); //new term difference from remainder
            if( newTermDiff > abs( cfg.curPolynom->remainder ) ){ break; } //break if the new term results in bigger abs(remainder)

            //If the new term meets the conditions, add it to candidate term list
            if(newTerm.symbolCount() <= cfg.maxSymbols){ //Limit term symbol count
                candiTerms.push_back(newTerm);
            }
        }
    }

    if(candiTerms.empty()){ return 1; }

    sortTerms(cfg,&candiTerms);

    //Add new branchest to main list. Indices 0 to maxBranches
    for (int i = 0; (i < cfg.maxBranches) && (i < candiTerms.size()); i++){
        Polynom bufCandi = *cfg.curPolynom;
        bufCandi.add(candiTerms[i]);
        cfg.storage->push_back(bufCandi);
    }

    return 0;
}






//Main
void sortPolynoms(std::vector<Polynom> *polys, uint maxCoef=10000) {

    //invalidate branches that didn't reach remainder=0
    for (Polynom& polynom : *polys) {
        if(abs(polynom.remainder) < maxCoef){
            polynom.terminate();
        }else
        if(polynom.remainder != 0){
            polynom.validEnd=0;
        }
    }

    //Remove incomplete polynomials
    polys->erase(std::remove_if(polys->begin(), polys->end(), [](const Polynom& p) {
        return !p.validEnd;
    }), polys->end());

    //Precompute strings and symbol counts
    struct PolynomWithString {
        Polynom poly;
        std::string cachedString;
        int symbolCount;

        PolynomWithString(const Polynom& p) : poly(p), cachedString(p.getString()), symbolCount(p.symbolCount()) {}
    };

    //Copy original polynomials to
    std::vector<PolynomWithString> cachedPolys;
    cachedPolys.reserve(polys->size());
    for (const auto& poly : *polys) {
        cachedPolys.emplace_back(poly);
    }

    // Sort using the cached values
    std::sort(cachedPolys.begin(), cachedPolys.end(), [](const PolynomWithString& a, const PolynomWithString& b) {
        if (a.symbolCount != b.symbolCount) {
            return a.symbolCount > b.symbolCount; // Sort by symbol count ascending
        }
        return a.cachedString > b.cachedString; // Sort by ASCII order
    });

    //Remove duplicates
    auto last = std::unique(cachedPolys.begin(), cachedPolys.end(), [](const PolynomWithString& a, const PolynomWithString& b) {
        return a.cachedString == b.cachedString; // Remove duplicates based on cached string
    });

    polys->clear(); //copy back to polys
    for (auto it = cachedPolys.begin(); it != last; ++it) {
        polys->push_back(it->poly);
    }
}



void printManual(char** argv){
        std::string path = argv[0];
        std::string programName = path.substr(path.find_last_of("/\\") + 1);
        std::cout << "This program searches polynomial for a given number \n";
        std::cout << "Usage    : ./" << programName << " [integer] [terms] [branches] [polynoms] [digits]\n";
        std::cout << "Argument : Default : Explanation\n";
        std::cout << "integer  : 2^31-1  : Number 0 to 2^63-10001. Non-numerics are hashed\n";
        std::cout << "terms    : 2       : Number of terms in polynomial\n";
        std::cout << "branches : 100     : Branch count per node in the tree search\n";
        std::cout << "polynoms : 20      : Number of displayed polynomials\n";
        std::cout << "digits   : Varies  : Number of digits in term constants\n";
}

int main(int argc, char** argv){
    if(argc>1){
        std::string firstArg = argv[1];
        if( firstArg=="-h" || firstArg=="-?" || firstArg=="--help" ){
            printManual(argv);
            return 0;
        }
    }

    const int64_t MAX_THREADS = 24;

    Polynom starter;
    int64_t defaultTarget = int64Hash("Therma");
    starter.target=defaultTarget;
    //starter.target = 3083;
    //starter.target = 13238717;
    //starter.target=4611686018427387904;

    if(argc>0){
        if(argv[1]){
            if(is_numeric(argv[1])){
                size_t num=std::strtol(argv[1], NULL, 10);
                starter.target=num;
            }
            else{
                starter.target=int64Hash(argv[1]);
            }
        }
    }

    starter.remainder = -starter.target;
    //starter.target = 128;

    std::vector<Polynom> mainList; //Main list of polynom candidates
    mainList.push_back(starter);

    int64_t MAX_COEF=pow(10,digitLen(starter.target)/2);
    MAX_COEF = MAX_COEF > 10000 ? 10000 : MAX_COEF;

    //config
    Polynom::SearchCache cache;
    cache.maxBranches=300;
    cache.termCount=2;
    cache.maxSymbols=digitLen(starter.target);
    cache.maxExpo=62;
    cache.maxCoef = MAX_COEF;
    size_t maxLeaves = 100000;
    size_t maxDisplay = 20;
    bool primalityCheck = 0;

    if(argc>1){
        if(argv[2]){
            size_t num=std::strtol(argv[2], NULL, 10);
            cache.termCount=num;
        }
    }
    if(argc>2){
        if(argv[3]){
            size_t num=std::strtol(argv[3], NULL, 10);
            if(num){
                cache.maxBranches=num;
            }
        }
    }
    if(argc>3){
        if(argv[4]){
            size_t num=std::strtol(argv[4], NULL, 10);
            if(num){
                maxDisplay=num;
            }
        }
    }
    if(argc>4){
        if(argv[5]){
            size_t num=std::strtol(argv[5], NULL, 10);
            if(num){
                MAX_COEF=num;
                cache.maxCoef = MAX_COEF;
            }
        }
    }

    size_t testedTerms = 0;

    ThreadPool pool(MAX_THREADS); // Adjust number of threads to hardware
    for(int j=0;j<cache.termCount;j++){

        const std::vector<Polynom> prevList=mainList; //Make a copy that the list won't change
        size_t prevSize=prevList.size();
        prevSize= std::min(maxLeaves,prevSize); //Limit maximum leaf node count

        std::mutex mainListMutex;
        std::vector<std::vector<Polynom>> threadStorage; //New branches storage. unique per thread
        threadStorage.resize(prevSize);


        //searchTerms function will be adding new elements each loop
        for(int i=0;i<prevSize;i++){
            if( prevList[i].validEnd==0 || prevList[i].remainder==0 ){ continue; } //Skip invalid or terminated branches

            Polynom::SearchCache localCache=cache; //Copy config for each thread. Usually better for small caches
            localCache.curPolynom=&prevList[i]; //Address of polynom copy that the thread will branch out
            localCache.storage=&threadStorage[i];

            // Enqueue work for the thread pool
            pool.enqueue([i, localCache, &mainList, &prevList, &mainListMutex] {
                bool err = localCache.curPolynom->searchTerms(localCache); //populater threadStorage aka localCache.storage
            });
        }
        pool.waitUntilFinished();

        //invalidate previous branches that weren't terminated and therefore become part of the trunk
        for (int i=0; i<prevSize; i++){
            if (!threadStorage[i].empty()) {//If valid
                size_t branchCount = std::min(threadStorage[i].size(), (size_t)cache.maxBranches); //Insert last N elements to maunList
                mainList.insert(mainList.end(), threadStorage[i].end() - branchCount, threadStorage[i].end());
            }else if(mainList[i].remainder!=0){
                mainList[i].validEnd = 0;
                if (i != 0) {
                    std::cout << "No branches found for " << mainList[i].getString() << "\n";
                }
            }

            if( abs(prevList[i].remainder) < cache.maxCoef ){
                mainList[i].terminate(); //if remainder is X digits or shorter, add it as constant +N*1^1
            }else
            if (mainList[i].remainder!=0){
                mainList[i].validEnd=0;
            }
        }

        std::cout << "Depth: " << j << " Branches: " << mainList.size() << "\n";
        testedTerms+=mainList.size();
    }

    sortPolynoms(&mainList,MAX_COEF);

    //Print
    int count = 0;

    size_t startPrint=0;
    if(maxDisplay!=-1 && !primalityCheck){
        startPrint=0;
        if(mainList.size()>maxDisplay){
            startPrint=mainList.size()-maxDisplay;
        }
    }
    for (size_t i = startPrint; i < mainList.size(); i++) {
        if( primalityCheck && !mainList[i].isPrimesOnly() ){ continue; }
        std::string bufStr=mainList[i].getString();
        std::cout << bufStr;
//      std::cout << " " << mainList[i].symbolCount();
        std::cout << "\n";
        count++;
    }
    if(mainList.size()>0){
        int64_t polyValue = mainList[0].value();
        char polyAscii[9];
        int64ToChar(polyValue,polyAscii,sizeof(polyAscii));
        std::cout << "Decimal: " << polyValue << "\n";
        std::cout << "Ascii: " << polyAscii << "\n";
    }
    std::cout << "Listed: " << count << "\n";
    std::cout << "Terms tested: " << testedTerms << "\n";

    if(argc<=1){
            printManual(argv);
    }

    return 0;
}