/*
   Copyright (C) 2022 Stephen Hewitt, cambridgeclarion.org

   This file is part of a programme to research fast Fermat factorisation algorithms.

   This is free software: you can redistribute it and/or modify
   it under the terms of version 3 of the GNU General Public License as
   published by the Free Software Foundation.

   It is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public Licence
   along with this programme.  If not, see <https://www.gnu.org/licenses/>.

   A copy of the licence is also available at
   https://www.cambridgeclarion.org/download/GPLv3.txt
   It has the following cryptographic hash
   SHA256 3972dc9744f6499f0f9b2dbf76696f2ae7ad8af9b23dde66d6af86c9dfb36986

*/

/*
    This programme is for research on Fermat factorisation algorithms.

    For each of 24 different categories of a number to factorise it prints a list of possible optimum filter values.

    For more information see the articles on www.cambridgeclarion.org
    "Empirical explorations of faster Fermat factorisation" by Stephen Hewitt
*/

#include <cstdio>
#include <list>
#include <cmath>

#define ARRAYSIZE(a) (sizeof a/sizeof a[0])

/*
    TR or table ratio is the number of elements that will be needed in a table to implement this filter.
    It is equal to the number of values (within its input range) that pass the filter.
    It is called a ratio because if two filters F1 and F2 with table ratios TR1 and TR2 are combined to make a filter Fc,
    then the table ratio of the combined filter is:

    TRc = TR1 * TR2

    It also follows from the definition of tr that rr = input_range/tr
*/

class Filter {
public:
    unsigned long long input_range;
    double rr;
    unsigned long long tr;
bool operator < (const Filter & f2) { return rr < f2.rr; }
};

#include "series.gh"

class Odd_prime_data
{
public:
unsigned base;
unsigned series_A_designator_mask;
unsigned series_B_designator_mask;
const unsigned long *tr_a_series;
unsigned max_digits;
unsigned requested_digits;
};

const Odd_prime_data db[] = {
2, 0,0,0,0,0,
3, 0xa0, 0xb0, tr_base3_series_A, ARRAYSIZE(tr_base3_series_A), 5,
5, 0xa00, 0xb00, tr_base5_series_A, ARRAYSIZE(tr_base5_series_A), 3,
7, 0xa000, 0xb000, tr_base7_series_A, ARRAYSIZE(tr_base7_series_A), 2,
// 11, 0xa0000, 0xb0000, tr_base11_series_A, ARRAYSIZE(tr_base11_series_A),
};

/* The designator mask is a sort of pun in hex.
   An integer is used to designate the category of N, with each nibble representing one prime base.
   When the integer is printed in hex it will be for example "ABA3",
   meaning series A in base 7, series B in base 5, series A in base 3 and series 3 in base 2.
   The masks can be simply ORed together to form this.
*/

class Base2_data
{
public:
unsigned series_designator_mask;
const unsigned long *tr_series;
unsigned max_digits;
unsigned requested_digits;
};

/* The following is slightly conceptually different to the db[]
   - because here each row is for a different series, rather than a different base
   Note some of the terms in the series could (clearly) be eliminated by hand, because they result in no increase in rr
   - but rely on pareto_purge() to remove these for now.
   Motivation is we can then simply use the auto-generated series in each of the three arrays, without hand intervention
*/

const Base2_data base2_db[] = {
1, tr_base2_series_001, ARRAYSIZE(tr_base2_series_001), 10,
5, tr_base2_series_101, ARRAYSIZE(tr_base2_series_101), 5,
3, tr_base2_series_11, ARRAYSIZE(tr_base2_series_11), 3,
};

typedef std::list<Filter> Filters;

void show_factors(unsigned long long number)
{
    static unsigned bases[] = {11,7,5,3,2};
    unsigned digits[ARRAYSIZE(bases)];

    if (!number)
        return;

    unsigned bases_printed = 0;
    for (unsigned i = 0; i < ARRAYSIZE(bases); i++)
    {
        unsigned factor = bases[i];
        digits[i] = 0;
        while (0 == (number % factor))
        {
            factor *= bases[i];
            digits[i]++;
        }
        if (digits[i])
        {
            if (bases_printed > 0)
                putchar(' ');
            printf("%u", bases[i]);
            if (digits[i] > 1)
                printf("^%u", digits[i]);
            bases_printed++;
        }
    }
}

int status_message(int status, const char* message)
{
    fprintf(stderr, "series: error %d %s\n", status, message);
    return status;
}

/* Delete every filter in the passed @filters list where there is another that has a better rr and no worse tr
   Do by sorting into ascending rr, so the worst rr is at the start
   - then for each entry, if there is one further up the list with no worse tr, the current one on the list can be purged
*/

void pareto_purge(Filters& filters)
{
    filters.sort();
    auto it = filters.begin();
    while (it != filters.end())
    {
        auto here = it++;
        for (auto it2 = it; it2 != filters.end(); it2++)
        {
            if (it2->tr <= here->tr)
            {
                filters.erase(here);
                break;
            }
        }
    }
}

void print_filter(const Filters& filters, unsigned long N_classification)
{
    printf("starting new N classification %lX\n", N_classification);
    for (auto filter: filters)
    {
        if (filter.input_range == 1)
            continue;
        printf("%lX rr=%.2f; size=%llu; ", N_classification, filter.rr, filter.tr);
        printf("q=%.3f; ", log(filter.rr)/log(filter.tr));
        printf(" M=%llu ", filter.input_range);
        show_factors(filter.input_range);
        putchar('\n');
    }
}

/*
    This recursive function models the tree of combinations by calling itself.

    It adds another layer (another prime base) of filtering to the passed list of filters
    and passes the resulting list to itself to add the next layer
    until all the allowed layers (prime bases) have been done, whereupon it instead prints the list.

    There are two potential conceptual trees here:

    1. the possibilities at each base for the category of N (series A or B)
    2. the choices at each base of how many digits to include in the filter.

    The list of filters that is passed to it is for a specific classification of N (up to the prime base considered)
    From the one list, it derives two lists for the next base: one for series A and one for series B.
    It calls itself for each of these two lists.

    For example it might be called say with base 3 filters with N classified as A1
    It will then generate filters based on base 5, covering the two cases:
    N classified as AA1
    N classified as BA1

    Note the conceptual tree of filter choices can be pruned by the pareto_purge() before reaching the leaves.
    This is because of the monotonic way in which filters combine.
*/

void odd_prime_filter(const Filters& filters_so_far, unsigned base_index, unsigned long N_classification)
{
    unsigned base = db[base_index].base;
    unsigned base_power = base;

    /* first A series */
    Filters derived_a = filters_so_far;

    /* for each allowed number (k) of LSDs in this base (as defined in database db[])
       for every filter on the list so far,
       create a new filter that adds filtering using k LSDs of this base to the existing filter
       The derived list consists of every filter on the original list plus all the newly created filters
    */
    for (unsigned k = 1; k <= db[base_index].requested_digits; k++)
    {
        for (auto it = filters_so_far.begin(); it != filters_so_far.end(); it++)
        {
            Filter f;
            f.input_range = it->input_range*base_power;
            f.tr = it->tr * db[base_index].tr_a_series[k - 1];
            f.rr = (double)f.input_range/f.tr;
            derived_a.push_front(f);
        }
        base_power *= base;
    }
    pareto_purge(derived_a);

    /* now do B series - only two cases to consider: with and without this base */
    Filters derived_b = filters_so_far;

    /* for every filter on the list so far,
       add another filter filtering also using the first LSD of this base
    */
    for (auto it = filters_so_far.begin(); it != filters_so_far.end(); it++)
    {
        Filter f;
        f.input_range = it->input_range*base;
        f.tr = it->tr * (base - 1)/2; // from empirical discovery August 2022, www.cambridgeclarion.org/60.html
        f.rr = (double)f.input_range/f.tr;
        derived_b.push_front(f);
    }
    pareto_purge(derived_b);

    if (base_index + 1 < ARRAYSIZE(db))
    {
        odd_prime_filter(derived_a, base_index + 1, N_classification | db[base_index].series_A_designator_mask);
        odd_prime_filter(derived_b, base_index + 1, N_classification | db[base_index].series_B_designator_mask);
    }
    else
    {
        print_filter(derived_a, N_classification | db[base_index].series_A_designator_mask);
        print_filter(derived_b, N_classification | db[base_index].series_B_designator_mask);
    }
}

void even_prime_filter(const Filters& filters_so_far, unsigned long N_classification)
{
    /* for each of the three series or N classifications */
    for (unsigned ss = 0; ss < ARRAYSIZE(base2_db); ss++)
    {
        Filters derived = filters_so_far;

        /* for each allowed number (k) of LSDs in this series (except k=1 which we know is pointless)
           for every filter on the list so far,
           add this filter filtering using k LSDs of base 2 in this series
        */
        unsigned base_power = 4;
        for (unsigned k = 2; k <= base2_db[ss].requested_digits; k++)
        {
            for (auto it = filters_so_far.begin(); it != filters_so_far.end(); it++)
            {
                Filter f;
                f.input_range = it->input_range*base_power;
                f.tr = it->tr * base2_db[ss].tr_series[k - 1];
                f.rr = (double)f.input_range/f.tr;
                derived.push_front(f);
            }
            base_power *= 2;
        }
        pareto_purge(derived);
        odd_prime_filter(derived, 1, N_classification | base2_db[ss].series_designator_mask);
    }
}

int main(int argc, char* argv[])
{
    Filter unity; // dummy filter which has no effect
    Filters filters;

    unity.input_range = 1;
    unity.tr = 1;
    unity.rr = 1;

    filters.push_back(unity);

    if (argc != 1)
        return status_message(9, "bad invocation should be ./optimise");

    /* Because series.gh is generated by another programme which is also used for other things
       check at run time that it has put enough digits in the prime series arrays for this programme.
    */
    for (unsigned base_index = 1; base_index  < ARRAYSIZE(db); base_index++)
    {
        if (db[base_index].requested_digits > db[base_index].max_digits)
            return status_message(7, "internal in one element of (odd prime base) db[] more digits requested than there are in the array for that base");
    }

    for (unsigned ii = 1; ii  < ARRAYSIZE(base2_db); ii++)
    {
        if (base2_db[ii].requested_digits > base2_db[ii].max_digits)
            return status_message(6, "internal in an element of base2_db[] more digits requested than there are in the array for the series");
    }

    even_prime_filter(filters, 0);
    return 0;
}