-- BIGNUM by rnd v05122018b
-- functions: 
--	new, tostring, rnd, import/exportdec,import/exporthex _add, _sub, mul, div2, div, is_larger, is_equal, add, sub,
--	binary2base, base2binary
--	bignum.barrett, bignum.mod, bignum.modpow

bignum = {};
bignum.new = function(base,sgn, digits)
	local ret = {};
	ret.base = base -- base of digit system
	ret.digits = {};
	ret.sgn = sgn -- sign of number,+1 or -1
	local data = ret.digits;
	local m = #digits;
	ret.digits = digits; -- THIS SEEMS TO MAKE A NEW COPY! if you work on this original wont change
	--for i=1,m do data[i] = digits[m-i+1] end -- copy
	return ret
end

bignum.rnd = function(base,sgn, length) -- random number
	local ret = {};
	for i =1,length do ret[#ret+1] = math.random(base)-1 end
	return bignum.new(base,sgn,ret)
end

bignum.tostring = function(n)
	local ret =  {};
	for i = #n.digits,1,-1 do ret[#ret+1] = n.digits[i] end
	return (n.sgn>0 and "" or "-") .. table.concat(ret,"'") .. "_" ..n.base
end

--n1 = bignum.new(10,-1,{5,7,3,1})
--say(bignum.tostring(n1))

bignum.importdec = function(ndec)
	local ret = {};
	local sgn = ndec>0 and 1 or -1;
	local base = 10;
	local n = ndec*sgn;
	local data = {};
	while n>0 do
		local r = n%base
		data[#data+1] = r;
		n=(n-r)/base
	end
	ret.base = base; ret.sgn = sgn; ret.digits = data;
	return ret
end

local importdec_test = function()
	local ndec = math.random(10^9);
	local n = bignum.importdec(ndec)
	say("importdec_test : " .. ndec .. " -> " .. bignum.tostring(n))
end
--importdec_test()

bignum.exportdec = function(n) -- warning: can cause overflow if number larger than 2^52 ~ 4.5*10^15
	local ndec = 0;
	for i = #n.digits,1,-1 do ndec = 10*ndec + n.digits[i] end
	return ndec*n.sgn
end

bignum.importhex = function(hex) -- nhex is string with characters 0-9(48-57) and a-f(97-102)
	local ret = {sgn=1,base = 16, digits = {}};
	local data = ret.digits;
	local length = string.len(hex);
	for i = length,1,-1 do
		local c = string.byte(hex,i)
		if c>=48 and c<=57 then
			data[length-i+1] = c-48
		elseif c>=97 and c<=102 then
			data[length-i+1]=c-97+10
		end
	end
	return ret
end

local importhex_test = function()
	local hex = "deadbeef";
	local n = bignum.importhex(hex)
	say(bignum.tostring(n))
end
--importhex_test()

bignum.exporthex = function(nhex) -- returns string with hex
	if nhex.base~=16 then return end
	local data = nhex.digits;
	local ret = {};
	for i = #data,1,-1 do
		local c = data[i];
		if c<10 then ret[#ret+1] =  string.char(48+c) else ret[#ret+1] =  string.char(97+c-10) end
	end
	return table.concat(ret,"")
end

local exporthex_test = function()
	local n = {sgn=1,base=16,digits = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15}}
	say(bignum.exporthex(n))
end
--exporthex_test()

-----------------------------------------------
--	ADDITION
-----------------------------------------------


bignum._add = function(n1,n2,res) -- assume both >0, same base: n1+n2 -> res
	local b = n1.base;
	local m1 = #n1.digits;
	local m2 = #n2.digits
	local m = m1; if m2<m then m = m2 end
	local M = m1;if m2>M then M = m2 end
	
	local data1 = n1.digits; local data2 = n2.digits;
	res.digits = {} -- expensive?
	local data = res.digits; local carry = 0;
	for i = 1,M do
		local j = (data1[i] or 0) +(data2[i] or 0) + carry;
		if j >=b then carry = 1; j = j-b else carry = 0 end
		data[i] = j
	end
	if carry== 1 then data[M+1] = 1 end
	res.base = n1.base
end

local _add_test = function()
	local n1 = bignum.rnd(10,1,5)
	local n2 = bignum.rnd(10,1,5)
	local res = bignum.new(10,1,{})
	bignum._add(n1,n2,res)
	say("_add_test: " .. bignum.tostring(n1) .. " + " .. bignum.tostring(n2) .. " = " .. bignum.tostring(res))
end
--_add_test()

-----------------------------------------------
--	SUBTRACTION
-----------------------------------------------

bignum._sub = function(n1,n2,res) -- assume n1>n2>0, same base: n1-n2 -> res
	local b = n1.base;
	local m1 = #n1.digits;
	local m2 = #n2.digits
	local m = m1; if m2<m then m = m2 end
	local M = m1;if m2>M then M = m2 end
	
	local data1 = n1.digits; local data2 = n2.digits;
	res.digits = {};
	local data = res.digits; local carry = 0;
	local maxi = 0;
	for i = 1,M do
		local j = (data1[i] or 0) - (data2[i] or 0) + carry;
		if j < 0 then carry = -1; j = j+b else carry = 0 end
		if j~=0 then maxi = i end -- max nonzero digit
		data[i] = j
	end
	
	for i = maxi+1,M do	data[i] = nil end -- remove trailing zero digits if any
	res.base = n1.base
end

local _sub_test = function()
	local n1 = bignum.rnd(10,1,5)
	local n2 = bignum.rnd(10,1,5)
	local res = bignum.new(10,1,{})
	bignum._sub(n1,n2,res)
	say("_sub_test: " .. bignum.tostring(n1) .. " - " .. bignum.tostring(n2) .. " = " .. bignum.tostring(res))
end
--_sub_test()

bignum.is_equal = function(n1,n2) -- assume both >0, same base. return true if n1==n2
	local b = n1.base;
	local data1 = n1.digits; local data2 = n2.digits;
	if #data1~=#data2 then return false end
	for i =#data1,1,-1 do -- from high bits
		local d1 = data1[i];
		local d2 = data2[i];
		if d1~=d2 then return false end
	end
	return true -- all digits were ==
end

bignum.is_larger = function(n1,n2) -- assume both >0, same base. return true if n1>=n2
	local b = n1.base;
	local data1 = n1.digits; local data2 = n2.digits;
	if #data1>#data2 then return true elseif #data1<#data2 then return false end
	--remains when both same lentgth
	for i =#data1,1,-1 do -- from high bits
		local d1 = data1[i];
		local d2 = data2[i];
		if d1>d2 then return true elseif d1<d2 then return false end
	end
	return true -- all digits were >=, still larger
end

local is_larger_test = function()
	local n1 = bignum.rnd(10,1,5)
	local n2 = bignum.rnd(10,1,5)
	local res = bignum.is_larger(n1,n2);
	if res then res = "larger" else res = "smaller" end
	say("is_larger_test : " .. bignum.tostring(n1) .. " is ".. res .. " than "  .. bignum.tostring(n2))
end
--is_larger_test()

bignum.add = function(n1,n2,res) -- handle all cases, >0 or <0
	local sgn1 = n1.sgn;
	local sgn2 = n2.sgn;
	if sgn1*sgn2>0 then bignum._add(n1,n2,res); res.sgn = sgn1; return end -- simple case
	
	local is_larger = bignum.is_larger(n1,n2) -- is abs(n1)>abs(n2) ?
	local sgn = 1;
	if is_larger then sgn = sgn1 else sgn = sgn2 end
	
	if is_larger then 
		bignum._sub(n1,n2,res);
	else
		bignum._sub(n2,n1,res);
	end
	res.sgn = sgn
end

local add_test = function()
	local ndec1 = math.random(10^5) * (2*math.random(2)-3);
	local ndec2 = math.random(10^5) * (2*math.random(2)-3);
	
	local n1 = bignum.importdec(ndec1)
	local n2 = bignum.importdec(ndec2)
	local res = bignum.new(10,1,{})
	bignum.add(n1,n2,res)
	
	local resdec = bignum.exportdec(res);
	say("add_test: " .. bignum.tostring(n1) .. " + " .. bignum.tostring(n2) .. " = " .. bignum.tostring(res) .. " CHECK : " .. resdec-(ndec1+ndec2))
end
--add_test()

bignum.sub = function(n1,n2,res) -- handle all cases, >0 or <0
	--just add(n1,-n2)
	local sgn1 = n1.sgn;
	local sgn2 = -n2.sgn;
	if sgn1*sgn2>0 then bignum._add(n1,n2,res); res.sgn = sgn1; return end -- simple case
	
	local is_larger = bignum.is_larger(n1,n2) -- is abs(n1)>abs(n2) ?
	local sgn = 1;
	if is_larger then sgn = sgn1 else sgn = sgn2 end
	
	if is_larger then 
		bignum._sub(n1,n2,res);
	else
		bignum._sub(n2,n1,res);
	end
	res.sgn = sgn
end

local sub_test = function()
	local ndec1 = math.random(10^5) * (2*math.random(2)-3);
	local ndec2 = math.random(10^5) * (2*math.random(2)-3);
	
	local n1 = bignum.importdec(ndec1)
	local n2 = bignum.importdec(ndec2)
	local res = bignum.new(10,1,{})
	bignum.sub(n1,n2,res)
	local resdec = bignum.exportdec(res);
	say("sub_test: " .. bignum.tostring(n1) .. " - " .. bignum.tostring(n2) .. " = " .. bignum.tostring(res) .. " CHECK : " .. resdec-(ndec1-ndec2))
end
--sub_test()

-----------------------------------------------
--	MULTIPLY 
-----------------------------------------------

bignum.mul = function(n1,n2,res)
	
	local base = n1.base
	local sgn = n1.sgn*n2.sgn;
	
	local data1 = n1.digits; local m1 = #data1;
	local data2 = n2.digits; local m2 = #data2;
	
	res.digits = {}; res.base = base
	local data = res.digits; local m = m1+m2;
	
	local carry = 0
	for i = 1, m1 do
		-- multiply i-th digit of data1 and add to res
		local d1 = data1[i];
		carry  = 0
		for j = 1,m2 do
			local d2 = data2[j];
			local d = carry + d1*d2;
			local r =  (data[i+j-1] or 0) + d
			if r>=base then 
				data[i+j-1] = r % base; carry = (r - (r%base))/base 
			else 
				data[i+j-1] = r; carry = 0
			end
		end
		if carry>0 then data[i+m2] = carry % base end
	end
end

local mul_test = function()
	local ndec1 = math.random(10^8) 
	local ndec2 = math.random(10^8)
	
	local n1 = bignum.importdec(ndec1)
	local n2 = bignum.importdec(ndec2)
	local res = bignum.new(10,1,{})
	bignum.mul(n1,n2,res)
	local resdec = bignum.exportdec(res);
	say("mul_test: " .. bignum.tostring(n1) .. "*" .. bignum.tostring(n2) .. " = " .. bignum.tostring(res) .. " CHECK : " .. resdec-(ndec1*ndec2))
end
--mul_test()

-- m = 300, base 2^26, 100 repeats: amd ryzen 1200: 0.1s, amd-e350 apu 1.6ghz (2010) : 5.15s
mul_bench = function()
	local m = 300;
	local base = 2^26
	local r = 100
	
	local n1 = bignum.rnd(base, 1, m)
	local n2 = bignum.rnd(base, 1, m)
	local res = {digits = {}};
	local t = os.clock()
	for i = 1, r do	bignum.mul(n1,n2,res) end
	local elapsed = os.clock() - t;
	--say("n1 = " .. bignum.tostring(n1) .. ", n2 = " .. bignum.tostring(n2))
	say("mul benchmark. ".. m .. " digits, base " .. base .. ", repeats " .. r ..  " -> time " .. elapsed)
end
--mul_bench()

local exp_test = function()
	local n1 = bignum.importdec(2);
	local res1 = bignum.importdec(2);
	local res2 = bignum.importdec(1);
	
	local m=128
	for i = 1, m do
		bignum.mul(n1,res1, res2) -- n1*res1 = res2
		bignum.mul(n1,res2, res1) -- n1*res1 = res2
	end
	
	say("2^" .. (2*m) .. " = " .. bignum.tostring(res2) .. " CHECK " ..2^(2*m))

end
--exp_test()

-----------------------------------------------
--	DIVIDE
-----------------------------------------------

bignum.div2 = function(n,res) -- res = n/2, return n % 2. note: its safe to do: bignum.div2(res,res);
	
	local base = n.base;
	local data = n.digits; local m = #data;
	
	res.digits = {};
	local rdata = res.digits;
	local carry = 0
	
	local q = data[m]/2; 
	local fq = math.floor(q);
	if q~=fq then carry = base end
	if fq>0 then rdata[m] = fq else rdata[m]=nil end -- maybe digits shrink by 1?
	
	for i = m-1,1,-1 do
		local q = (data[i]+carry)/2;
		local fq = math.floor(q)
		if q~= fq then carry = base else carry = 0 end
		rdata[i] = fq;
	end
	if carry ~= 0 then return 1 else return 0 end
end
	
local div2_test = function()
	local ndec1 = math.random(10^8) 
	local n1 = bignum.importdec(ndec1)
	local res = bignum.new(10,1,n1.digits)
	bignum.div2(res,res) -- res = res/2
	
	say("div2_test: n1/2 = " .. bignum.tostring(n1) .. "/2 = " .. bignum.tostring(res) .. " = res")
	local rescheck = bignum.new(10,1,{})
	bignum._add(res,res,res);bignum._sub(n1,res, res);
	
	say("CHECK: n1-2*res = " .. bignum.tostring(res))
	
end

--div2_test()

--[[
	very simple division that works reasonably well (we only need 1 division for barrett reduction anyway, could use precomputed too)
	
	strategy: bisection for f(x) = x*D + comparison with N, takes around Log_2(initial range) steps (sums+mults), 
		so ~O(D^2*log base^(n2-n1))
		low, mid, high. pick reasonably good initial range guess, like near order of magnitude close.
		mid =  (low+high)/2
		compute: compare N and mid*D, if N bigger then low = mid else high = mid..
	BENCHMARKS: (amd ryzen 1200)
		HUGE N=10k bit number/ D=5k bit number : 1.5 s
		N = 8k bit number / D = 4k bit number : 0.7s
		N = 1040 bit, D = 520 bit: 0.0042 s ( typical application srp,diffie-hellman in Z_~2^512 group)
		if D is 3900 bits it takes around 3900 steps of iteration
			amd-e350 apu 1.6ghz (sep 2010)
				N = 8k bit, D = 4k bit, divide takes 44s ( 60x slower than ryzen)
			AMD Dual Core E2-1800 (july 2012)
				N = 8k bit, D = 4k bit, 5.25s
	TODO: 
		possible speed improve ?: after there are some digits correct
		reduce N by N = N-q0*D which will effectively decrease N and multiplies of mid ( smaller numbers ). then keep adding 
		obtained q's together to get final quotient.
		
	--]]

bignum.div = function(N,D, res) -- res = [N/D]
	
	local base = N.base;
	res.base = base
	res.digits = {}; 
	local data = res.digits; 
	
	local n1 = #N.digits;local n2 = #D.digits;
	-- trivial cases, prevent wasting time here
	if n1<n2 then res.digits = {0}; return end -- clearly N<D
	if n2 == 1 and D.digits[1] == 1 then res.digits = N.digits return end -- division by 1!
	
	local low = bignum.new(base,1,{})
	local high = bignum.new(base,1,{})
	-- better initial range for less needed iterations
	local ldigits = low.digits;local hdigits = high.digits;
	for i = 1,n1-n2 do ldigits[i]=0;hdigits[i]=0 end
	ldigits[n1-n2]=N.digits[n1];hdigits[n1-n2+1] = ldigits[n1-n2];
	--say("low " .. bignum.tostring(low) .. " high " .. bignum.tostring(high))
	
	local mid = bignum.new(base,1,{});
	local temp = bignum.new(base,1,{});
	local step = 0;
	
	while step < 100000 do -- in practice this uses around log_2 (base^(n2-n1)) iterations, for example dividing 8192 bit number by 4096 takes ~4000 iterations..
		step = step + 1
		bignum._add(low,high,mid); bignum.div2(mid,mid); -- mid = (low+high)/2
		
		if bignum.is_equal(low,mid) then 
			if DEBUG then say("DONE. step  " .. step) end-- .. " low = " .. bignum.tostring(low) .. " high = " .. bignum.tostring(high) .. " mid = " .. bignum.tostring(mid))
			res.digits = mid.digits
			return
		end

		bignum.mul(D, mid, temp) -- temp = D*mid
		if bignum.is_larger(N,temp) then low.digits = mid.digits else high.digits = mid.digits end
	end
end


local div_test = function()
	local ndec1 = math.random(10^8) 
	local n1 = bignum.importdec(ndec1)
	local ndec2 = math.random(10^6)
	local n2 = bignum.importdec(ndec2)
	
	local res = bignum.new(10,1,{})
	bignum.div(n1,n2,res)
	
	local temp = bignum.new(10,1,{})
	bignum.mul(n2,res,temp);bignum._sub(n1,temp,temp) -- temp = n1 - n2*res
	
	say(ndec1/ndec2)
	say("n1/n2 =  " .. bignum.tostring(n1) .. " / " .. bignum.tostring(n2) .. " = res = " .. bignum.tostring(res) .. ", residue n1-n2*res = " .. bignum.tostring(temp) .. (bignum.is_larger(n2,temp) and " (IS SMALLER THAN n2) " or " FAIL."))
end
--div_test()

local divbignum_test = function()
	local m = 300;
	local base = 2^26
	local n1 = bignum.rnd(base, 1, m)
	local n2 = bignum.rnd(base, 1, m/2)
	local res = {sgn=1, digits = {}};
	DEBUG = true -- to display how many steps were needed
	local t = os.clock();bignum.div(n1,n2,res); local elapsed = os.clock() - t;
	DEBUG = false
	local temp = {sgn=1, digits = {}};
	bignum.mul(n2,res,temp);bignum._sub(n1,temp, res); -- res = n1-n2*res
	if bignum.is_larger(n2, res) then 
		say("divbignum_test : residue n1 - n2*res is smaller than n2. OK.") 
	else
		say("divbignum_test : residue n1 - n2*res is NOT smaller than n2. FAIL.") 
	end
end

--divbignum_test()

div_bench = function()
	local m = 300;
	local base = 2^26
	local r = 1
	
	local n1 = bignum.rnd(base, 1, m)
	local n2 = bignum.rnd(base, 1, m/2)
	local res = {sgn=1, digits = {}};
	local t = os.clock()
	for i = 1, r do	bignum.div(n1,n2,res) end
	local elapsed = os.clock() - t;
	say("n1 = " .. bignum.tostring(n1) .. "\nn2 = " .. bignum.tostring(n2) .. "\nn1/n2 = "  .. bignum.tostring(res))
	say("div benchmark. n1 (".. m .. " digits ( " .. 26*m .." bits)), n2 (" .. m/2 .. " digits), base " .. base .. ", repeats " .. r ..  " -> time " .. elapsed)
end

--div_bench()

bignum.base2binary = function(_n)
	local base = _n.base;
	local n = {sgn = 1, base = base, digits = _n.digits }
	local data = n.digits;
	local i = 0;
	local out =  {};
	while (#data > 1 or (#data==1 and data[1] > 0)) do -- n>0
		i=i+1
		out[i]  = bignum.div2(n,n); data = n.digits;
	end
	return {sgn=1,base = 2, digits = out}
end

local base2binary_test = function()
	local base = 10;
	local m = 2;
	local n = bignum.rnd(base, 1, m)
	local nb = bignum.base2binary(n);
	say(bignum.tostring(n) .. " -> " .. bignum.tostring(nb))
end
--base2binary_test()


bignum.binary2base = function(n, newbase) -- newbase must be even
	local base = n.base;
	local ret = {sgn=1,base=newbase, digits = {0}}
	local out = ret.digits
	local data = n.digits
	for i = #data,1,-1 do
		bignum._add(ret,ret,ret) -- ret = 2*ret
		out = ret.digits
		out[1]=out[1]+ data[i]; -- if newbase is even no carry, else more complication here
	end
	return ret
end

local binary2base_test = function()
	local base = 10;
	local m = 4;
	local nb = bignum.rnd(2, 1, m)
	local n = bignum.binary2base(nb,base);
	say(bignum.tostring(nb) .. " -> " .. bignum.tostring(n))
end
--binary2base_test()


-----------------------------------------------
--	MODULAR MULTIPLY
-----------------------------------------------

-- a,b in Z_n -> a*b mod n = ?
-- how to compute a % n efficiently? We can use barrett reduction trick.
-- normally: a%n = a - [a/n]*n. Instead of division we compute [a/n] with multiply and shift ( base = b)
-- [a/n] = [a*(B^k/n)/B^k] = [a*m/B^k]. Here integer m is [B^k/n] for some k, where B^k>=n. since
-- a*(m/B^k-1/n) < 1 we get a*(m-B^k/n) < B^k or m-B^k/n < B^k/a. since left side is always <1 this will be true if 
-- 1 < B^k/a or a < B^k. note since a*m/B^k - a/n < 1 after applying [ ] we can still get difference = 1 (but not more),
-- so need to check if  a - [a*m/B^k]*n is smaller than n. If not additional -n is needed.
-- so REQUIREMENTS: n<=B^k, a< B^k.
-- if we need a<n^2 (like in modulo multiply in Z_n) then this means: (n-1)^2 < B^k. So if n<B^N then k should be 2N.

-- barret = {n = bignum,  m =  from barrett.., k= .., }
-- returns barrett data. used to cheaply compute a/n without division as bitshiftright_k(a*m) + possible 1, where a <= (n-1)^2
bignum.get_barrett = function(n) 
	local base = n.base;
	local k = 2*#n.digits+2; -- n<B^(n1+1) -> k = 2*(n1+1)
	local Bk = bignum.new(base,1,{})
	local res = bignum.new(base,1,{})
	local data = Bk.digits;
	for i =1,k do data[i]= 0 end; data[k+1]=1; -- this is B^k
	bignum.div(Bk, n,res);
	return {n=n, m=res, k=k};
end

local get_barrett_test = function()
	local d=4
	local ndec2 = math.random(10^d) 
	local n2 = bignum.importdec(ndec2)
	local barrett = bignum.get_barrett(n2)
	local barrettm = math.floor(10^(2*d+2)/ndec2)
	say(bignum.tostring(barrett.m) .. "(CHECK: " .. barrettm .. ")")
end
--get_barrett_test()

-- mod using barrett. possible improvement: montgomery.
bignum.mod = function(a,barrett,res) -- a should be less or equal (n-1)^2, stores a%n into res
	
	local k = barrett.k;
	local n = barrett.n;
	local m = barrett.m;
	local base = a.base;
	
	bignum.mul(a,m,res); -- large multiply 1: res = a*m

	local data = res.digits;local n1 = #data;  --res = res / B^k
	for i = 1, n1-k do data[i]=data[i+k] end; for i = n1-k+1,n1 do data[i] = nil end -- bitshift
	
	local temp = bignum.new(base,1,{});
	bignum.mul(res,n, temp); -- multiply 2: res*n
	bignum._sub(a,temp,res); -- subtract: res = a - res*n
	if bignum.is_larger(res,n) then bignum._sub(res,n,res) end
end

local mod_test = function()
	local m = 3;
	local base = 10;
	local n1 = bignum.rnd(base, 1, 2*m)
	local n2 = bignum.rnd(base, 1, m)
	local barrett = bignum.get_barrett(n2)
	local res = bignum.new(base,1,{});
	bignum.mod(n1, barrett, res);
	local is_larger = bignum.is_larger(n2,res);
	
	say("barrett mod_test: n1 " .. bignum.tostring(n1) .. " n2 " .. bignum.tostring(n2) .. " res = n1 % n2 = " .. bignum.tostring(res) .. " CHECK: res<n2 " .. (is_larger and "OK" or "FAIL") )
	
end
--mod_test()

--  computing a^b mod n using modular exponentiation : 
--  1. suppose we have b = sum_i b_i*2^i ( binary expansion). Then a^b = prod_i a^(2^i) a^b_i.
-- 		we start from low bits and multiply with new term each step. we compute a^(2^i) mod n iteratively, with a = a^2 mod n.
--		In each step if b_i == 1 we multiply with a. Work: 2-3 multiplies each step
--  2. we could first compute all b_i and then start from high end ( like efficient polynomial evaluate):
--  	res = res^2*a^(b_i) mod n . work (1-2 multiplies per step + all b_i).
--	3. sliding window. works like 2 but it uses binary representation of b and does larger parts (window size k) at once. Window it
--		uses is of the form 1xxxxxxxx1, where total length less than prescribed k. Potentially k* faster but with precompute of 2^k-1 values

-- method 1 
bignum.modpow = function(a_,b_,barrett) -- efficiently calculate a^b mod n, need log_2 b steps
	local base = a_.base
	local a = bignum.new(base,1,a_.digits); -- base
	local b = bignum.new(base,1,b_.digits); -- exponent
	local bdata = b.digits;
	
	local ret = bignum.new(base,1,{1});
	local temp = bignum.new(base,1,{});
	
	while (#bdata > 1 or (#bdata==1 and bdata[1] > 0)) do -- b>0
		if bdata[1] % 2 == 1 then
			bignum.mul(ret, a, temp) 
			bignum.mod(temp,barrett, ret) -- ret = a*ret % n
		end
		bignum.div2(b,b); bdata = b.digits -- b = b/2
		
		bignum.mul(a,a,temp);bignum.mod(temp,barrett, a) -- a=a^2 % n
	end
	return ret
end

-- times: base 2^26. 
-- 	ryzen 3 1200: m=150 (4000 bits): 7.2s, m=80(2080bits): 1.1s,  m=40: 0.14, m = 30: 0.07, m = 20: 0.028
-- 	amd-e350 apu 1.6ghz (2010), m = 20: 0.59s, m = 40: 6.28s
DH_bench = function() -- test for modpow with diffie-hellman key exchange
	local m = 40;
	local bits = 26
	local base = 2^bits;
	local a = bignum.rnd(base, 1, m)
	local b = bignum.rnd(base, 1, m)
	local c = bignum.rnd(base, 1, m)
	local n = bignum.rnd(base, 1, m)
	
	local t = os.clock()
	local barrett = bignum.get_barrett(n) -- precompute this from modulo n
	local t1 = os.clock(); say("barret pre computation : " .. t1-t);t=t1
	
	local resb = bignum.modpow(a,b, barrett); -- efficient modular a^b using barrett reductions
	t1 = os.clock();say("resb = modpow(a,b) : " .. t1-t);t=t1
	
	local resc = bignum.modpow(a,c, barrett); -- a^c
	t1 = os.clock();say("resc = modpow(a,c) : " .. t1-t);t=t1
	
	local resbc = bignum.modpow(resb,c, barrett); -- (a^b)^c
	t1 = os.clock();say("resbc = modpow(resb,c) : " .. t1-t);t=t1
	
	local rescb = bignum.modpow(resc,b, barrett); -- (a^c)^b
	t1 = os.clock();say("rescb = modpow(resc,b) : " .. t1-t);t=t1
	
	say("DH benchmark: (a^b)^c = " .. bignum.tostring(resbc) .. "\nCHECK equality (a^b)^c == (a^c)^b : " .. (bignum.is_equal(resbc,rescb) and "OK" or "FAIL"))
	say("settings: base " .. base .. ", digits " .. m .. " (" ..  m*bits .." bits )")
end
DH_bench()

-- note: DH security ( log problem ) can be solved for 512 bit by expert team, 1024 prime modulo for state actors ( logjam attack, "Imperfect Forward Secrecy:
--	How Diffie-Hellman Fails in Practice", Oct 2015

-- rnd key exchange v2 (improved key hide when doing verification):
-- auth + key exchange in one idea: A alice (client), B bob (server)
-- key generation: A generates its shared secret as v = g^x for secret random x and exchanges it with B. 
-- this initial exchange can use DH with very large group to prevent any snooping.

-- verification: 
-- 		1. B picks random r and tells A y=g^r+v. 
--			Note here that listeners only see g^r+v so they have no clue what v is or what g^r is. Even from multiple sessions they
--			learn nothing if g is generator of whole group and r truly random, since then g^r can be anything with same probability.
-- 		2. A confirms identity by sending back hash(g^rx) = hash(v^r) = hash((y-v)^x) to prove you know a. 
-- 		3. shared secret is then another hash of K=g^rx = (g^r)^x =  (g^x)^r.

--
-- basic idea is this: i tell you y=g^r for random r, you need to tell me y^x, so i believe you know x. But hide y by adding v and with hash so
-- listeners dont get any clues.

--observations: 
-- 0. note that K lies in potentially smaller group generated by g^x, srp 6a is better here since for srp K = g^(b(a + ux)).
-- 		where a,b are random
-- 		The order of Z_p* is 2q for prime q, where p-1 = 2q ( p safe prime ). Order of g^x is either 2  or q or 2q (if (g^x)^2~= 1 then its not 2)
-- 		So order of g^x is still (p-1)/2 in this case.
-- thm: if G is finite group then order of any element divides |G|. Furthermore, if p is prime dividing |G| there exists element
-- in G of order p (cauchy thm)

-- 1. Suppose there is observer C. first time of interaction he can see v=g^x sent from  A to B, but getting x is not possible ( log problem )
-- 2. later C can see B send v+g^r, and if he doesnt know v already this does not reveal it to him.. If we encrypt this this might introduce
-- 		weakness if not properly done, cause only good password would encrypt to number, bad password would be random mess.

-- WEAKNESS: if C knows v ...?
-- idea: first time exchange uses larger DH group for more safety to exchange: v = g^x to prevent C from learning x.

--srp v6a : http://srp.stanford.edu/design.html
--NOTE: srp v1 http://srp.stanford.edu/design1.html used weak approach to compute secret: S = (Wp*Xp)^Ys, where Xp = could be leaked 
-- pass. verifier, Wp = controlled by client, Ys = random server secret
-- if rogue client know Xp he could select Wp so that Wp*Xp becomes value of his choice, hence controlling what the final shared secret
-- will be.
-- rnd key exchange v2 does not have that weakness. even if rogue learns v = g^x by other means he still needs to solve log problem for x.

-- prime generation: openssl prime -generate -bits 2048 -hex

self.remove() -- prevent loop in robot csm