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

   This file is a programme to demonstrate Fermat factorisation.

   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 file must be linked with GMP, the GNU Multiple Precision Arithmetic Library.
   If both the header file from the library, gmp.h and the static version of the library, libgmp.a
   are in the current directory then an executable fermat1 can be built as follows:

   g++ fermat1.cc -l:libgmp.a -o fermat1
*/

#include <iostream>
#include <fstream>
#include <sstream>
#include <csignal>
#include <ctime>

#include "gmp.h"

typedef mpz_t Bignum;

static int status_message(int status, const char* message)
{
    std::cerr << "fermat1: error " << status << ' ';
    std::cerr << message << '\n';
    return status;
}

static void printbig(const Bignum& number, const char* message, int base=10)
{
    std::cout << message << ' ';
    mpz_out_str(0, base, number);
    std::cout.put('\n');
}

// the following flag gives factorise() ability to be safely called multiple times, without danger of multiple mpz_init()s

static bool static_variables_are_initialised_f = false;
static Bignum delta;
static Bignum candidate;

static time_t initial_seconds;
static Bignum initial_delta;
static Bignum initial_candidate;
static Bignum modulus;

void static initialise_static_variables()
{
    mpz_init(candidate);
    mpz_init(delta);
    mpz_init(initial_delta);
    mpz_init(initial_candidate);
    mpz_init(modulus);

    static_variables_are_initialised_f = true;
}

void report(int i)
{
    unsigned long seconds_elapsed;
    mpz_t temp;
    mpz_init(temp);
    mpz_t iterations;
    mpz_init(iterations);
    mpz_t delta_increase;
    mpz_init(delta_increase);
    Bignum check;
    mpz_init(check);


    printf("--------- Report on ");

    if (i == -1 )
        printf("successful completion\n");
    else
    {
        printf("signal %d\n", i);
        printf("Note there is no protection against data races on this asynchronous signal; variables could be in an inconsistent state\n");
    }

    printbig(candidate, "current small square candidate", 10);

    printbig(modulus, "modulus");
    mpz_add(check, candidate, modulus);
    printbig(check, "big square back-calculated from small square + modulus");
    if (!mpz_perfect_square_p(check))
        status_message(133, "internal problem - back-calculated big square is not square");

    printbig(initial_delta, "initial delta", 10);
    printbig(delta,         "current delta", 10);

    mpz_sub(delta_increase, delta, initial_delta);
    printbig(delta_increase, "delta increase", 10);

    mpz_tdiv_q_2exp(iterations, delta_increase, 3); // since delta += 8 per iteration
    printbig(iterations, "iterations", 10);

    mpz_tdiv_r_2exp(temp, delta_increase, 3);
    if (mpz_sgn(temp))
       printf("Note: inconsistent value for delta, increase is not a multiple of 8\n");

    time_t now = time(0);
    seconds_elapsed = now - initial_seconds;
    printf("seconds elapsed %lu\n", seconds_elapsed);

    if (seconds_elapsed > 0)
    {
        mpz_tdiv_q_ui(temp, iterations, seconds_elapsed);
        printbig(temp, "iterations per second", 10);
    }

    fflush(stdout);
}

int factorise(unsigned base, const char* modulus_string)
{
    unsigned long long search_count = 0;

    /* following two variables have been made global for benefit of signal handler which wants to read them
       so let's initialise at the earliest possible point
    */
    if (!static_variables_are_initialised_f)
        initialise_static_variables();

    if (mpz_set_str(modulus, modulus_string, base))
        return 29;

    if (mpz_perfect_square_p(modulus))
    {
        return status_message(85, "fermat: modulus is a perfect square.");
    }

    int probab_prime = mpz_probab_prime_p(modulus, 50);
    if (probab_prime == 2)
    {
        return status_message(77, "fermat: modulus is certainly prime says gmp");
    }
    std::cout << "fermat: mpz_probab_prime: " << (probab_prime ? "probably": "certainly not") << '\n';
    printbig(modulus, "modulus");

    Bignum sqrt_N;
    mpz_init(sqrt_N);

    mpz_sqrt(sqrt_N, modulus);
    printbig(sqrt_N, "sqrt_N debug"); // debug


    // signal(SIGINT, report);
    signal(SIGQUIT, report);

    mpz_t initial_big;
    mpz_init_set(initial_big, sqrt_N);

    /*
       It is only necessary to consider either odd or even big presquares,
       whether odd or even depending on the modulus.
       The following line adjusts the initial guess onto an odd or even value depending on the modulus.
       See the explanation in https://www.cambridgeclarion.org/41.html
    */

    if (mpz_tstbit(modulus, 1) == mpz_tstbit(initial_big, 0))
        mpz_add_ui(initial_big, initial_big, 1); // change from odd to even or vice versa

    printbig(initial_big, "first big presquare"); // debug

    mpz_mul(candidate, initial_big, initial_big);
    printbig(candidate, "first big square"); // debug
    mpz_sub(candidate, candidate, modulus);

    mpz_set(initial_candidate, candidate);
    printbig(initial_candidate, "first small square candidate"); // debug

    mpz_mul_2exp(delta, initial_big, 2);
    mpz_add_ui(delta, delta, 4); // so delta = 4b + 4, where b is the big presquare

    mpz_set(initial_delta, delta);
    printbig(delta, "initial delta"); // debug

    initial_seconds = time(0);
    printf("initial seconds %lu\n", initial_seconds);

    mpz_t debug_big_sq;
    mpz_init(debug_big_sq);

    while (!mpz_perfect_square_p(candidate))
    {
        mpz_add(candidate, candidate, delta);
        mpz_add_ui(delta, delta, 8);
    }
    report(-1);

    Bignum check;
    mpz_init(check);

    printbig(modulus, "modulus, ie number factorised");
    mpz_add(check, candidate, modulus);
    printbig(check, "big square back-calculated from small square + modulus");
    if (!mpz_perfect_square_p(check))
        status_message(175, "big internal problem - back-calculated big square is not square");

    Bignum big;
    mpz_init(big);

    mpz_sqrt(big, check);
    printbig(big, "big back-calculated from big square");

    /* now check based on delta = 4b + 4 */
    mpz_sub_ui(check, delta, 4);
    mpz_tdiv_q_2exp(check, check, 2);
    printbig(check, "big back-calculated from delta");

    if (mpz_cmp(check, big))
    {
        status_message(177, "big internal problem - big back-calculated from delta is not consistent with small square + modulus");
    }

    mpz_sub(check, big, initial_big);
    printbig(check, "Absolute increase in big presquare from initial value");

    Bignum pp;
    mpz_init(pp);
    Bignum qq;
    mpz_init(qq);

    Bignum small;
    mpz_init(small);

    mpz_sqrt(small, candidate);
    printbig(small, "small");

    mpz_add(pp, big, small);
    mpz_sub(qq, big, small);

    printbig(pp, "prime factor p");
    printbig(qq, "prime factor q");

    mpz_mul(check, pp, qq);
    if (mpz_cmp(check, modulus))
    {
        return status_message(75, "pq is not equal to modulus");
    }
    printf("good for you pq = modulus\n");
    return 0;
}

static const char* examples[] =
{
#include "modulus1.h"
};

int main(int argc, char* argv[])
{
    if (argc != 1)
        return status_message(45, "this simple demonstration programme does not take any command line arguments");
    printf("Fermat factorisation demonstration\n");
    printf("Copyright 2022 Stephen Hewitt\n");
    printf("This is free software with ABSOLUTELY NO WARRANTY.\n");
    printf("Licence: GNU GPL version 3\n");
    return factorise(16, examples[0]);
}