// SPDX-License-Identifier: GPL-3.0-only
// Copyright (c) 2022, Sylvain Huet, Ambermind
// Minimacy (r) System

// Primal generation : bigPrimal(nbits, nloops)
// using Miller–Rabin primality test

const _BIG1= \big 1;;

const PRIMAL_PRODUCTS_23_743= \big 0x1bf05ca9ec4207d5e82daf4b48410614b1b2e9b0313303bf2daace79df7a0c71938fb835261e2c88fa1a16c8919ab02f9f8e97396caf94f8d05d6f6614e4cc94f2b81668783646bf1e315a76b81921fc13519e80f2d34d95e81eba938f0009eddd2bb579d2ebcb9aa8e529279ad85ace2408e46657b3d7a9574e20ab8dd3;;

// Miller-Rabin test
// https://defeo.lu/in420/DM3%20-%20Test%20de%20Miller-Rabin
// n prime => for any random number p, p^(n-1)=1[n] (Fermat theorem)
// n prime => the only square roots are 1 and -1
fun _randomMillerRabin(n, Rinv, one, nMinusOne, s, d)=	// nMinusOne= d.2^s, with s>0
	// select a random number y in [2, n-2]
	let bigMod(bigRand(bigNbits(n), false), n) -> y in
	if bigIsNull(y) || y==one || y==nMinusOne then _randomMillerRabin(n, Rinv, one, nMinusOne, s, d)
	else // compute y^d [n]
	let bigExpModMontgomeryFast(y, d, n, Rinv) -> x in
	if x==one || x==nMinusOne then true	// x=1 or x=-1[n] => n might be prime
	else (
		for r = 1; r < s do (
			set x=bigSquareModMontgomery(x, n, Rinv);	// compute x^2[n]
			if x==one then return false	// x=1 => n can't be prime
			else if x==nMinusOne then return true	// x=-1 => n might be prime
		);
		false	// => n can't be prime
	);;

fun _millerRabinLoop(n, nloops)=
//	bigDump ("n",n);
	if bigGcd(PRIMAL_PRODUCTS_23_743,n)==_BIG1 then
	let bigSub(n, _BIG1) -> nMinusOne in
	let bigLowestBit(nMinusOne) -> s in
	let bigASR(nMinusOne, s) -> d in	// nMinusOne= d.2^s, with s>0
	let bigMontgomery(n) -> Rinv in
	let bigMontgomeryForm(\big 1, n) -> one in
	let bigMontgomeryForm(nMinusOne, n) -> nMinusOne in
	(
		for i=0;i<nloops do if !_randomMillerRabin(n, Rinv, one, nMinusOne, s, d) then return false;
		true
	);;

fun _primalCheckLoop(n, nloops)=
	threadResign();	// on slow platforms
	// we check n with Miller-Rabin. If it fails we continue with n+2
	let bigCheckFirstFactors(n) -> delta in
	let if delta==0 then n else bigAdd(n,bigFromInt(delta)) -> n in
	if !_millerRabinLoop(n, nloops) then _primalCheckLoop(\big n+2, nloops)
	else n;;

fun _primalFirstCandidate(nbits, squareMin)=
	let bigRand(max(2, nbits), true) -> b in
	// we check that b^2 has 2*nbits (useful for RSA p and q: p and q have nbits, and p*q has 2*nbits, not 2*nbits-1)
	if \big b*b<squareMin then _primalFirstCandidate(nbits, squareMin)
	else
	if bigIsEven(b) then bigAdd(b, _BIG1)
	else b;;

export fun bigPrimal(nbits, nloops)=
	let bigPower2(nbits*2-1) -> squareMin in
	let _primalFirstCandidate(nbits, squareMin) -> n in
	_primalCheckLoop(n, nloops);;