#!/usr/bin/python
import os
from operator import mod, sub
from uprms import radlst
from math import log as ilog
from random import random as rand, seed

def eea(a,b):
	v1 = [a,1,0]
	v2 = [b,0,1]
	while v2[0]<>0:
	   p = v1[0]//v2[0]
	   v2, v1 = map(sub,v1,[p*vi for vi in v2]), v2
	return v1

def minv(m,k):
	v = eea(m,k)
	if v[1]>=0:
		return v[1]
	else:
		return k+v[1]
		
def mdiv(dor, dend, rad):
    inv=minv(dend, rad)
    return mod(inv*dor, rad)

def euc(ar0, ar1):
    while ar0>0:
        rx=ar1%ar0
        ar1=ar0
        ar0=rx
    return ar1

class logs:
	def __init__(self, rad):
		self.rad=rad
		self.lgs=[0]*rad
		tot=1
		ind=0
		while ind<rad:
			self.lgs[tot]=ind
			tot<<=1
			if tot>rad:
				tot-=rad
			ind+=1
	def __getitem__(self, ind):
		return self.lgs[ind]
	def __str__(self):
		return '%d'%self.rad
		
class log:
	def __init__(self, dif, rut, exp):
		self.dif=dif
		self.rut=rut
		self.exp=exp
	def xp(self):
		return (cns.dexp(self.exp)*self.rut)%cns.mxm
	def __str__(self):
		return '%d,%d,%d'%(self.dif, self.rut, self.exp)

class ilog:
	def __init__(self, dif, irut, exp):
		self.dif=dif
		self.irut=irut
		self.exp=exp
	def xp(self):
		rut=rlst[self.irut]
		tlg=log(self.dif, rut, self.exp)
		return tlg.xp()
	def __str__(self):
		return '%d,%d,%d'%(self.dif, self.irut, self.exp)

class cons:
	def __init__(self, radx):
		self.radx=radx
		self.lradx=radx>>1
		self.rads=[]
		self.lrads=[]
		self.mxms=[]
		self.lcms=[]
		self.logss=[]
		self.mxm=1
		self.lcm=1
		self.rng=1
		self.logss=[]
		self.invs=[]
		self.linvs=[]
		for rad in radlst:
			lrad=rad>>1
			if euc(lrad, self.lcm)==1 and euc(lrad, self.lradx)==1:
				self.rads.append(rad)
				gcd=euc(rad-1, self.rng)
				self.rng=self.rng*(rad-1)/gcd
				self.mxms.append(self.mxm)
				self.lcms.append(self.lcm)
				self.logss.append(logs(rad))
				res=self.mxm%rad
				self.invs.append(minv(res, rad))
				self.mxm*=rad
				self.lrads.append(lrad)
				res=self.lcm%lrad
				self.linvs.append(minv(res, lrad))
				self.lcm*=lrad
				if self.lcm>self.lradx:
					break
		tot=2
		self.sqrs=[]
		bits=len(bin(self.mxm))
		ind=1
		while ind<bits:
			self.sqrs.append(tot)
			tot=(tot**2)%self.mxm
			ind+=1
		self.lglim=len(self.rads)**2
	def dexp(self, inum):
		str=bin(inum)[2:]
		dat=str[::-1]
		lim=len(dat)
		tot=1
		for ind in range(lim):
			if(dat[ind]=='1'):
				tot=(tot*self.sqrs[ind])%self.mxm
		return tot
	def dec2res(self, inum):
		lst=[]
		for rad in self.rads:
			lst.append(inum%rad)
		return lst
	def dec2lres(self, inum):
		lst=[]
		for lrad in self.lrads:
			lst.append(inum%lrad)
		return lst
	def lres2dec(self, ires):
		lim=len(ires)
		tot=ires[0]
		ind=1
		while ind<lim:
			lrad=self.lrads[ind]
			cry=mod(tot, lrad)
			dig=mod((ires[ind]-cry)*self.linvs[ind], lrad)
			tot+=dig*self.lcms[ind]
			ind+=1
		return tot
	def res2dec(self, ires):
		lim=len(ires)
		tot=ires[0]
		ind=1
		while ind<lim:
			rad=self.rads[ind]
			cry=mod(tot, rad)
			dig=mod((ires[ind]-cry)*self.invs[ind], rad)
			tot+=dig*self.mxms[ind]
			ind+=1
		return tot
	def vlg(self, inum):
		ind=0
		lim=len(self.rads)
		vlgs=[]
		while ind<lim:
			rad=self.rads[ind]
			lrad=self.rads[ind]>>1
			lgs=self.logss[ind]
			res=inum%rad
			lg=lgs[res]%lrad
			vlgs.append(lg%lrad)
			ind+=1
		return self.lres2dec(vlgs)
	def clg2dec(self, iclg):
		return (iclg.rut*self.dexp(iclg.exp))%self.mxm
	def clg(self, inum):
		vlg=self.vlg(inum)
		xnum=inum<<1
		if xnum>self.mxm:
			xnum-=self.mxm
		vlgx=self.vlg(xnum)
		dif=(vlgx-vlg)%self.lcm
		rng=self.rng
		gcd=euc(dif, rng)
		if gcd<1:
			dif/=gcd
			rng/=gcd
		inv=minv(dif, rng)
		rlg=(inv*vlg)%self.rng
		rut=mdiv(inum, self.dexp(rlg), self.mxm)
#		if rut>(self.mxm>>1):
#			rut=self.mxm-rut
#			rlg=(rlg+self.lcm)%self.rng
		return log(dif, rut, rlg)
	def rdx2dec(self, rut):
		ind=0
		lim=len(self.rads)
		res=[]
		msk=1
		while ind<lim:
			if msk&rut:
				res.append(self.rads[ind]-1)
			else:
				res.append(1)
			msk<<=1
			ind+=1
		return self.res2dec(res)
	def __str__(self):
		return '%d,%d,%d'%(self.mxm, self.rng, self.lcm)
		
def durdat():
	rdx=0
	dct={}
	lst=[]
	while rdx<cns.lglim:
		xp=cns.rdx2dec(rdx)
		dct[xp]=rdx
		lst.append(xp)
		rdx+=1
	return dct, lst
	
def clg2ilg(clg):
	iclg=ilog(clg.dif, 0, clg.exp)
	if clg.dif>1:
		return iclg
	if clg.rut in rdct:
		iclg.irut=rdct[clg.rut]
	return iclg
	
def ilg2clg(iclg):
	return log(iclg.dif, rlst[iclg.irut], iclg.exp)
	
def ilgmul(iclg0, iclg1):
	olg=ilog(1, 0, 0)
	if iclg0.dif>1 or iclg1.dif>1:
		return olg
	olg.irut=iclg0.irut^iclg1.irut
	olg.exp=iclg0.exp+iclg1.exp
	if olg.exp>cns.rng:
		olg.exp-=cns.rng
	return olg

def ilgdiv(dnd, dor):
	olg=ilog(1, 0, 0)
	if dnd.dif>1 or dor.dif>1:
		return olg
	olg.irut=dnd.irut^dor.irut
	olg.exp=dnd.exp-dor.exp
	if olg.exp<0:
		olg.exp+=cns.rng
	return olg
	
def duval():
	while 1:
		val=int(rand()*cns.mxm)
		clg=cns.clg(val)
		if clg.dif==1:
			return val, clg2ilg(clg)
			
def dumul():
	print 'multiply'
	seed(sed)
	cnt=0
	while cnt<10:
		val0=duval()
		val1=duval()
		olg=ilgmul(val0[1], val1[1])
		ans=(val0[0]*val1[0])%cns.mxm
		ansx=olg.xp()
		str='%d*%d=%d (%s)*(%s)=(%s)=%d'%(val0[0], val1[0], ans, val0[1].__str__(), val1[1].__str__(), olg.__str__(), ansx)
		if ansx==ans:
			print str
		else:
			print 'error', str
		cnt+=1
	
def dudiv():
	print 'divide'
	seed(sed)
	cnt=0
	while cnt<3:
		val0=duval()
		val1=duval()
		olg=ilgdiv(val0[1], val1[1])
		ans=mdiv(val0[0], val1[0], cns.mxm)
		ansx=olg.xp()
		str='%d/%d=%d (%s)/(%s)=(%s)=%d'%(val0[0], val1[0], ans, val0[1].__str__(), val1[1].__str__(), olg.__str__(), ansx)
		if ansx==ans:
			print str
		else:
			print 'error', str
		cnt+=1
		
def manx():
	print 'Residue range:', cns.mxm
	print 'Modules:', cns.rads
	print
	dumul()
	print
	dudiv()
	print

radx=101*199
sed=1000
cns=cons(radx)
rdat=durdat()
rdct=rdat[0]
rlst=rdat[1]
manx();