/* mbMod.cp */

#include "myTypedef.h"
#include "mbConv.h"
#include "mbMath.h"
#include "mbMod.h"
#include "mbRand.h"
#include <iostream>

using namespace std;

// returns theBit (0 or 1) at bitPosition in blockNum of x
static INT32 findBit(const mb& x, INT32 bitPosition, INT32 blockNum)
{
	UINT32				mask;
	
	// now construct mask for getting theBit, we want to AND with something like 00010000
	mask = 1;
	mask = (mask<<bitPosition);
	return (mask & x.b[blockNum]);
	
}/* findBit */


// Finds z = b^(m*2^k)(mod n) where k>=0
// Uses p(mod n)*q(mod n) = p*q(mod n)
void powerMod(mb& z, const mb& b, const mb& m, INT32 k, const mb& n)
{
	static mb		bmod, zt, zTemp, one, mt;
	static bool		arrayGood=false;
	INT32			i, blockNum, bitPosition;
	
	if(!arrayGood)
	{
		init(bmod);
		init(zt);
		init(zTemp);
		init(one, 1);
		init(mt);
		arrayGood = true;
	}
	
	if(m.n<0)
		return;
	
	if(n.n<=0)
	{
		z.n = 0;
		return;
	}
	
	if(!b.n	&& !m.n)
	{
		z.n = 0;
		return;
	}
	
	if(!compare(n, one))
	{
		z.n = 0;
		return;
	}
	
	if(!m.n || !compare(b, one))
	{
		equate(z, one);
		return;
	}
	   
	if(!b.n)
	{
		z.n = 0;
		return;
	}
	
	equate(mt, m);
	
	divRem(zTemp, bmod, b, n); // bmod = b(mod n)
	
	if(k<0)
	{
		mbShiftRight(mt, mt, -k); // x 2^(-|k|)
		if(!mt.n)
		{
			equate(z, one);
			return;
		}
	}
	else
	{
		// we build up to b^(2^k)(mod n)
		for(i=0;i<k;i++)
		{
			mul(bmod, bmod, bmod); // accumulates b^(2^i)(mod n)
			divRem(zTemp, bmod, bmod, n); // bmod = bmod(mod n)
		}
	}
	
	// now bmod = b^(2^k)(mod n), so need to build to b^(m*2^k)(mod n)
	// using mb form of m in binary
	equate(zt, 1.);  // value if m=0
	blockNum = 0;
	while(blockNum<m.n)
	{
		bitPosition = 0;
		while(bitPosition<32)
		{
			if(findBit(m, bitPosition, blockNum))
			{
				// bit = 1
				mul(zt, zt, bmod);
				divRem(zTemp, zt, zt, n); // reduce zt
			}
			mul(bmod, bmod, bmod);
			divRem(zTemp, bmod, bmod, n); // update and reduce bmod so k->k+1 etc.
			bitPosition++;
		}
		blockNum++;
	}
	
	equate(z, zt);

}/* powerMod */


// returns true if x is even, else false
static bool isItEven1(const mb& x)
{
	// test if low order bit is 1 or 0
	if(x.b[0] & 1)
		return false;
	return true;
	
}/* isItEven */


// if returns true then n is composite; if returns false it's inconclusive
bool Miller_Rabin(const mb& n)
{
	static mb		m, b, nn, zMod, one;
	static bool		arrayGood=false;
	INT32			i, q;
	
	if(!arrayGood)
	{
		init(m);
		init(b);
		init(nn);
		init(zMod);
		init(one, 1);
		arrayGood = true;
	}
	
	sub(nn, n, 1.); // nn = n-1 (mod value equivalent to -1)
	
	// get random 1 <= b <= n-1
	mbRandom(b, nn); // 0 <= b <= n-2
	add(b, b, 1.); // 1 <= b <= n-1
	
	// find q and m so n-1 = m*2^q, m an odd integer
	equate(m, nn);
	q = 0;
	while(isItEven1(m))
	{
		mbShiftRight(m, m, 1); // m=m/2
		q++;
	}
	
	// case a, if b^m = 1 (mod n) then return false
	powerMod(zMod, b, m, 0, n);
	if(!compare(zMod, one))
		return false;
	
	// case b, if there is an integer i in [0, q-1] such that b^(m*2^i) = nn (mod n) then return false
	for(i=0;i<q;i++)
	{
		powerMod(zMod, b, m, i, n);
		if(!compare(zMod, nn))
			return false;
	}
	
	return true;
	
}/* Miller_Rabin */


bool isItPrime(const mb& x, INT32 numTrials)
{
	INT32		i;
	
	for(i=0;i<numTrials;i++)  // For 30 trials that probability is less than about 8.7e-19 ((1/4)^30).
	{
		if(Miller_Rabin(x))
			return false;
	}
	return true;

}/* isItPrime */

/*
	z = x % n ; x can be negative, but let's always put out positive mods
	if z = x mod n then (x-z) is a multiple of n and so x mod n = z + n*k for all integers k
	If x>=0 just use divRem.
	if x<0 then
	z=-|x| is correct, find q so  -|x|+n*k is first >= 0. That's a suitable z. How to find k?
	Let q = |x|\n, then n*q is just <= |x|, so -|x|+n*q = z <= 0 and z = n*(q+1)-|x| >0 
*/
void myModulus(mb& z, const mb& x, const mb& n)
{
	static mb		temp, xt, q, nt;
	static bool		arrayGood=false;
		
	if(!arrayGood)
	{
		init(temp);
		init(xt);
		init(q);
		init(nt);
		arrayGood = true;
	}
	
	if(n.n<=0)
	{
		z.n = 0;
		return;
	}
	
	if(!x.n)
	{
		z.n = 0;
		return;
	}
	
	if(x.n>0)
	{
		divRem(temp, z, x, n);
		return;
	}
	
	// now x<0
	equate(xt, x);
	equate(nt, n);
	xt.n = -xt.n; // x>0
	
	div(q, xt, n);
	equate(z, q);
	add(z, z, 1.);
	mul(z, z, nt);
	sub(z, z, xt);
	if(!compare(z, nt))
		z.n = 0;
		
}/* myModulus */


// http://www.mathreference.com/num-mod,minv.html
// s = yt, t = xt
// returns false if x has no inverse mod, else z = x^-1 mod y
bool invMod(mb& z, const mb& x, const mb& y)
{
	static mb		s, t, q, a, b, m, r, temp, one; // xt = u  yt = v
	static bool		arrayGood=false;
	
	if(!arrayGood)
	{
		init(t);
		init(q);
		init(a);
		init(b);
		init(m);
		init(r);
		init(temp);
		init(one, 1);
		arrayGood = true;
	}
	
	equate(t, x);
	equate(s, y);
	equate(m, y);
	
	if(t.n<0)
		t.n = -t.n;
	
	if(s.n<0)
		s.n = -s.n;
	
	if(!s.n)
	{
		equate(z, t);
		return false;
	}
	
	a.n = 0;
	equate(b, 1.);
	
	while(t.n)
	{
		divRem(q, r, s, t); // r = s % t
		equate(s, t);
		equate(t, r);
		
		mul(temp, b, q);
		sub(temp, a, temp);
		myModulus(temp, temp, m); // temp = (a-b*q) % m;
		equate(a, b);
		equate(b, temp);
	}
	
	if(compare(s, one))
		return false;
	
	equate(z, a);
	return true;
	
}/* invMod */
