/*
   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.

    It calculates and prints the first terms in series based on a prime number and
    related to reduction ratios in filters used to speed up the search in Fermat factorisation.

    The data generated by this programme is used in another programme "optimise", which generates
    the set of possible optimum filter values that are in the Appendix of
    article www.cambridgeclarion.org/61.html

    With further processing, the information can also be used for other purposes, including
    to generate the tables of reduction ratios published in article /42.html

    It prints series for the first few primes.
    For base 2 it prints 3 series.
    For an odd prime base it prints a single series.

    Each series represents the number of LSD patterns passing a filter based on that
    prime as increasing numbers of LSDs are considered in the filter.

    The filter is for candidates for b where we are trying to factorise N by finding
    b^2 = N + s^2

    Each series is printed in the form of a C/C++ declaration and initialisation of an
    array but this format can easily be changed.

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

*/

#include <memory>
#include <cstdio>

#define ARRAYSIZE(a) (sizeof a/sizeof a[0])
#include "stop2.h"

static const unsigned primes[] = { 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,};

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

/*
    prints a C array declaration for the requested RR series initialised with correct values.
    The array is named according to passed @series_designator

    Eg  unsigned long tr_base3_series_A[] = {2, 2, 4, 8, 22, };

    prints the number of pass patterns (patterns that pass the filter, in distinction to stop patterns)
    but NB it is the number of pass patterns out of the pattern space at that number of digits
    so the implied denominator for the RR (which is not printed) is different at each digit
    Eg for base 5 it prints 3, 7, 31 ...  meaning:
    for k = 1 (1 LSD) RR = 3/5
    for k = 2 (2 LSDs) RR = 7/25
    for k = 3 (3 LSDs) RR = 31/125

    One motivation for this is that the pattern space is then implicit so it avoids having to also specify it
*/

void print_reduction_series(const char* series_designator, unsigned print_digits_n, unsigned long *base_powers, unsigned long *culm_stop_at_digits, unsigned pattern_digits_n)
{
    printf("unsigned long tr_base%lu%s[] = {", base_powers[1], series_designator);
    for (unsigned i = 0; i < print_digits_n; i++)
    {
        /* The following division makes the number of stops be relative to the LSD pattern space at the associated number of digits
           which is (i + 1) digits.
           Before the division it is relative to the total LSD pattern space that was used in the counting
           For example at the second digit of base 7 the LSD pattern space is 49
           but the total patterns used in the counting might have been (say) 823543 (7^7)
           so it would divide the count down by 7^5.  The division result is always an integer.
        */
        printf("%lu, ", (base_powers[pattern_digits_n] - culm_stop_at_digits[i])/base_powers[pattern_digits_n - i - 1]);
    }
    printf("};\n");
}

void print_first_digit_tr(const char* name_stem, const char* suffix, unsigned long patterns_n, unsigned long *culm_stop_at_digits)
{
    printf("unsigned long %s_%s = %lu;\n", name_stem, suffix, (patterns_n - culm_stop_at_digits[0]));
}

/* For a given base and a given number of digits and a given stop table passed to it:
   generates the reduction ratio series for increasing number of digits from 1 up and including digits_n

   Note this will be for a specific LSB pattern (aka modular residue of N) but this pattern is implicit in the stop table
   - it does not need to be passed to this function.

   output is put into passed array culm_stopped.
   It is total stop patterns for each number of digits, first in array is 1 digit
   Note this is the total (cumulative) number of stop patterns for k digits, not the increase from k-1 digits
*/

int generate_reduction_series(char* stop_table, unsigned digits_n, unsigned long patterns_n, unsigned long *culm_stopped)
{
    // for convenience, the number of digits is the index into this array, but note stop_at_digits[0] is used and must be 0
    unsigned long stop_at_digits[digits_n + 1];

    for (unsigned i = 0; i < ARRAYSIZE(stop_at_digits); i++)
        stop_at_digits[i] = 0;

    for (unsigned long i = 0; i < patterns_n; i++)
    {
        if (stop_table[i] < ARRAYSIZE(stop_at_digits))
            stop_at_digits[stop_table[i]]++;
    }

    unsigned long culm_stop_at_digits = 0;
    for (unsigned i = 1; i < ARRAYSIZE(stop_at_digits); i++)
    {
        culm_stop_at_digits += stop_at_digits[i];
        culm_stopped[i - 1] = culm_stop_at_digits;
    }
    return 0;
}

/* prime_index is index into array of primes where 2 is at index 0, 3 at index 1, 5 at index 2, etc
   Returns 0 if all OK, error code otherwise
*/

int tabulate(unsigned prime_index, unsigned requested_digits_n)
{
    unsigned long base_powers[requested_digits_n + 1];

    if (prime_index >= ARRAYSIZE(primes))
        exit(status_message(36, "requested prime index for base is too big for prime list"));
    unsigned base = primes[prime_index];

    /* check requested_digits_n is not too many and fill in base_powers[] array at the same time.  */
    base_powers[0] = 1;
    for (unsigned ii = 1; ii <= requested_digits_n; ii++)
    {
        unsigned long new_powers = base_powers[ii - 1] * base;

        if (new_powers < base_powers[ii - 1])
                exit(status_message(34, "integer overflow too many digits (for given base) requested"));
        base_powers[ii] = new_powers;
    }

    unsigned pattern_digits_n = requested_digits_n;
    unsigned long patterns_n = base_powers[pattern_digits_n];

    auto uniq_ptr_valid_squares = std::make_unique<char[]>(patterns_n);
    char* valid_squares = uniq_ptr_valid_squares.get();

    for (unsigned i = 0; i < patterns_n; i++)
        valid_squares[i] = 0;

    fill_valid_squares(valid_squares, patterns_n, 100); // 100 is an arbitrary large number of digits, more than will ever be used

    /* the following will contain the number of digits before presquare becomes invalid
       so value of zero means that it is never valid (ie always stopped)
    */
    auto uniq_ptr_stop_table = std::make_unique<char[]>(patterns_n);
    char* stop_table = uniq_ptr_stop_table.get();

    if (2 == base)
    {
        unsigned long culm_stops[pattern_digits_n];

        /* For base 2  there are three cases defined as modulus ending in binary 11, 001 or 101 */

        fill_stop_table(stop_table, 1, base_powers, pattern_digits_n, valid_squares);
        generate_reduction_series(stop_table, pattern_digits_n, patterns_n, culm_stops);
        print_reduction_series("_series_001", pattern_digits_n, base_powers, culm_stops, pattern_digits_n);

        fill_stop_table(stop_table, 3, base_powers, pattern_digits_n, valid_squares);
        generate_reduction_series(stop_table, pattern_digits_n, patterns_n, culm_stops);

        unsigned digits_n = pattern_digits_n < 3? pattern_digits_n: 3;
        print_reduction_series("_series_11", digits_n, base_powers, culm_stops, pattern_digits_n);

        fill_stop_table(stop_table, 5, base_powers, pattern_digits_n, valid_squares);
        generate_reduction_series(stop_table, pattern_digits_n, patterns_n, culm_stops);

        digits_n = pattern_digits_n < 5? pattern_digits_n: 5;
        print_reduction_series("_series_101", digits_n, base_powers, culm_stops, pattern_digits_n);
    }
    else
    {
        char tag[2];
        unsigned long series_culm_stops[pattern_digits_n];

        tag[1] = 0; // string terminator
        tag[0] = 'A';

        /* lsd of 1 always generates series A */
        fill_stop_table(stop_table, 1, base_powers, pattern_digits_n, valid_squares);
        generate_reduction_series(stop_table, pattern_digits_n, patterns_n, series_culm_stops);
        print_reduction_series("_series_A", pattern_digits_n, base_powers, series_culm_stops, pattern_digits_n);
    }
    return 0;
}

int main(int argc, char* argv[])
{
    int ret = 0;
    unsigned j;

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

    tabulate(0, 20);

    for (j = 1; j < 5; j++)
        tabulate(j, 7);

    for (; j < 9; j++)
        tabulate(j, 4);

    for (; j < 19; j++)
        tabulate(j, 3);

    return 0;
}