/*
   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 information on series based on a prime number and
    related to reduction ratios in a filter used to speed up the search in Fermat factorisation.

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

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

    In the terms used in the articles, this programme generates certain information for each (odd prime) base.
    At the same time it empirically checks several conjectures, made in Part 4 and Part 6 of the series of articles.

    For the filter, the information this programme generates shows
    for each possible value of the single LSD of N in this base whether it:
    * causes the lower reduction ratio (series A) OR
    * causes the higher level reduction ratio (series B)

    It also checks the following conjectures and halts with an error message if one fails:

    1. For each prime base (except 2), there are exactly two different reduction ratios for a filter using one LSD.
    2. There are equal numbers of LSD values of N causing each series.
    (Value 0 is not considered a possible lsd, since that would mean N was a multiple of the small prime base)
    3. The number of stop patterns for a filter with 1 LSD is
       stops = (base - 1)/2 for series A
       stops = (base + 1)/2 for series B
       So for example for base 7, there are:
           3 stop patterns (ie big presquare LSD values) in series A
           4 stop patterns (ie big presquare LSD values) in series B
       - note this means they always sum to base.
    4. The pattern of series A/B is identical to the Legendre symbol (lsd/base)
        with the following mapping:

    Legendre(lsd/base) = 1 ie a square, corresponds to series A
    Legendre(lsd/base) = -1 ie not a square, corresponds to series B

    Files needed are this one (lsd2.cc), stop2.cc and stop2.h

    On GNU/Linux it has been built with

    g++  lsd2.cc stop2.cc -o lsd2

    and run with no arguments.
*/

#include <cstdio>

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

static const char 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, "lsd: error %d %s\n", status, message);
    return status;
}

/* base must be prime
   returns non zero on error
*/

int tabulate(unsigned base)
{
    char stop_table[base];
    char valid_squares[base];
    unsigned long base_powers[2]; // needed by fill_stop_table()

    /*
        Every value of lsd is mapped to one of two series.
        (Here "series" refers to the series of reduction ratios that will be obtained as the number of ls digits is increased above 1)
        The function of the following array is to record this mapping, as it is discovered by the code in this function.
        The two series always have a different initial reduction ratio, (initial meaning when considering only the 1st ls digit.)
        This array maps value of lsd (1 to base - 1) (as array position) to a series designation (1 or 2) (as element value).
        Note position 0 is unused, since 0 is not a possible lsd value of N (which is the product of a pair of large primes).
        We will designate as series 1 whichever series the first lsd tested belongs to.
        The other series will be series 2.
        After we have found both series we will see which one has the higher reduction ratio on single lsd.
    */
    unsigned lsd_mapping[base];

    base_powers[0] = 1;
    base_powers[1] = base;

    if (2 == base)
        return status_message(75, "internal, base 2 requested, for which this module is not suitable");

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

    fill_valid_squares(valid_squares, base, 1);

    /*
        this code will go through all digit values of N to demonstrate the conjecture
        that both series always have equal numbers of digit values and that all digit values belong to one of the two series.
    */

    for (unsigned i = 0; i < ARRAYSIZE(lsd_mapping); i++)
        lsd_mapping[i] = 0; // 0 means unknown, 'A' means series A. 'B' means series B

    unsigned long seriesA_stopped_n = (base - 1)/2; // conjecture in Part 4 article
    unsigned long seriesB_stopped_n = base - seriesA_stopped_n; // ditto

    unsigned seriesA_lsd_count = 0; // count of distinct values actually seen going through all lsd values
    unsigned seriesB_lsd_count = 0;

    // conjecture that lsd=1 is always series A:
    fill_stop_table(stop_table, 1, base_powers, 1, valid_squares);
    if (seriesA_stopped_n != number_stopped_at_first_digit(stop_table, base))
        return status_message(210, "conjecture failure (See Part 6, cambridgeclarion.org/62.html) stop count for LSD=1 is not (base + 1)/2");
    seriesA_lsd_count = 1;
    lsd_mapping[1] = 'A';

    for (unsigned long lsd = 2; lsd < base; lsd++)
    {
        fill_stop_table(stop_table, lsd, base_powers, 1, valid_squares);
        unsigned long stopped = number_stopped_at_first_digit(stop_table, base);

        if (stopped == seriesA_stopped_n)
        {
            lsd_mapping[lsd] = 'A';
            seriesA_lsd_count++;
            if (!valid_squares[lsd])
                return status_message(203, "conjecture failure: (See Part 6, cambridgeclarion.org/62.html) Legendre symbol incorrect: series A lsd value is not a square");
        }
        else if (stopped == seriesB_stopped_n)
        {
            lsd_mapping[lsd] = 'B';
            seriesB_lsd_count++;
            if (valid_squares[lsd])
                return status_message(204, "conjecture failure: (See Part 6, cambridgeclarion.org/62.html) Legendre symbol incorrect: series B lsd value is a square");
        }
        else
        {
            return status_message(200, "conjecture failure: (See Part 6, cambridgeclarion.org/62.html) stopped count value not equal to series A or series B");
        }
    }

    if (seriesA_lsd_count != seriesB_lsd_count)
        return status_message(213, "conjecture failure: (See Part 4, cambridgeclarion.org/60.html) unequal numbers of LSD values in each series");

    putchar('"');
    putchar('x');
    for (unsigned i = 1; i < base; i++)
        putchar(lsd_mapping[i]);
    putchar('"');
    printf(", // %u\n", base);

    return 0;
}

int main(int argc, char* argv[])
{
    if (argc != 1)
        return status_message(9, "bad invocation should be ./lsd");

    printf("\"xx\", // 2\n"); // dummy line for base 2
    for (unsigned i = 1; i < ARRAYSIZE(primes); i++)
    {
        int ret;

        ret = tabulate(primes[i]);
        if (ret)
            return ret;
    }
    return 0;
}