#include <stdio.h>
#include <stdlib.h>
#include "../ecc/ecc.h"
#include "../inv32.h"
#include "../crypto.h"
#include "sieve.h"

	/*  test for primality (probably); 1 if x is prime. test done n times	*
	 *  ; chance of wrong identification << (1/4)^n. Note however that		*
	 *  this is an extreme upper bound. For example for a 100 digit "prime"	*
	 *  the chances of false witness are actually < (.00000003)^n			*
	 *  See Kim & Pomerance "The probability that a random probable prime	*
	 *  is Composite", Math. Comp. October 1989 pp.721-741					*
	 *  The value of n is now adjusted internally to account for this.		*/

unsigned long is_prime (unsigned long *x, unsigned long small_prime)
{
	unsigned long	j, k, ndash = inv32 (0-x[1]);
	unsigned long	w1[ECC_WORDS + 2], w2[ECC_WORDS + 2], w8[ECC_WORDS + 2], w9[ECC_WORDS + 2];
	
	if (x[0] < 1) return 0;
	
/* Miller-Rabin */
	
    big_psub (x, ECC_1, w1); /* calculate k and mr_w8 ...    */
    k = 0;
    while (big_sdiv (w1, 2, w1) == 0)
    {
        k++;
        big_copy (w1, w8);
    }	// so that x = w8 * 2^k + 1
	
	j = 0;
	nres_powltr (small_prime, w8, w9, x, ndash);
	big_redc (w9, w9, x, ndash);
	big_psub (x, ECC_1, w2);
	while ((j > 0 || w9[0] != 1) && big_compare (w9, w2) != 0)
	{
		j++;
		if ((j > 1 && w9[0] == 1) || j == k) return 0;	// definitely not prime
		big_mad (w9, w9, w9, x, x, w9);
	}
    return 1;	// probably prime
}

//  find next highest prime from w using probabilistic primality test

void next_prime (unsigned long *w, unsigned long *x)
{
	unsigned long			w1[ECC_WORDS + 2];
	
	if (x[0] < 2)
	{
		big_convert (2, x);
		return;
	}
	big_copy (w, x);
	if (big_sdiv (x, 2, w1) == 0)	big_padd (x, ECC_1, x);
	else							big_padd (x, ECC_2, x);
	while (!is_prime (x, 2)) big_padd (x, ECC_2, x);
}

  /* If type=0 finds next highest "safe" prime p >= w *
   * A safe prime is one for which q=(p-1)/2 is also  *
   * prime, and will be congruent to 3 mod 4          *
   * However if type=1 finds a prime p for which      *
   * q=(p+1)/2 is prime, congruent to 1 mod 4         *
   * Set subset=1 for q=1 mod 4, subset=3 for         *
   * q=3 mod 4, subset=0 if you don't care which      */

static __inline void big_vadd (unsigned long *x, unsigned long n, unsigned long *z)
{
	unsigned long	y[3] = { 1, n, 0 };
	
	big_padd (x, y, z);
}

unsigned long next_safe_prime (int type, int subset, unsigned long *w, unsigned long *p)
{
	unsigned long		rem, increment, test0, test1, test2, test3;
	unsigned long		w10[ECC_WORDS + 2];
	
	big_copy (w, p);
	
    rem = p[1] & 7;
    if (subset == 0)
    {
        rem &= 3;
        if (type == 0)		big_vadd (p, 3 - rem, p);
        else
        {
            if (rem > 1)	big_vadd (p, 5 - rem, p);
            else			big_vadd (p, 1 - rem, p);
        }
        increment = 4;
    }
    else
    {
        if (subset == 1) 
        {
            if (type == 0)
            {
                if (rem > 3) big_vadd (p, 11 - rem, p);
                else         big_vadd (p, 3 - rem, p);
            }
            else
            {
                if (rem > 1) big_vadd(p, 9 - rem, p);
                else         big_vadd(p, 1 - rem, p);
            }
        }
        else 
        { 
            if (type == 0) big_vadd (p, 7 - rem, p);
            else
            {
                if (rem > 5) big_vadd (p, 13 - rem, p);
                else         big_vadd (p, 5 - rem, p);
            }
        }
        increment = 8;
    }
    if (type == 0) big_psub (p, ECC_1, w10);
    else           big_padd (p, ECC_1, w10);
    big_sdiv (w10, 2, w10);
	
    for (test1 = 0, test2 = 1, test3 = 1; (test1 < 64*64);)
    {
		test0 = 0;
        do
		{
            big_vadd (p, increment, p);
            big_vadd (w10, increment >> 1, w10);
			test0++;
			test1++;
			if (test0 > 64) goto nf;	// gave up on it
        } while (sieve_prime (w10) || sieve_prime (p));
		test2++;
        if (!is_prime (p, 2)) continue;
		test3++;
		// just return; allow for interruption and output;
		// if (is_prime (w10, 2))
		{
		//	*tries1 = test1;
		//	*tries2 = test2 - 1;
		//	*tries3 = test3 - 1;
			return 1;
		}
    }
nf:	// *tries1 = test1;
	// *tries2 = test2 - 1;
	// *tries3 = test3 - 1;
	return 0;
}