# -*- coding: cp1252 -*-
# -*- coding: ISO-8859-1 -*-
# Python Math Lab library
# Version 0.9.4
# This software is under GPL
# (c) 2007 Stefan Mueller-Stach 
# and Michael Mardaus and Tobias Nagel
# uses Python for Nokia Series 60

from __future__ import generators
import math, operator, random, sys, os

def bsgs(base,arg,mod):
    if not isprime(mod):
        return 'mod is not prime'
    #if not is_primitive_root(mod,base):
    #    return 'base is not primitive root'
    Q = int(math.floor(math.sqrt(mod)))
    while Q*Q+Q+1 < mod:
        Q = Q+1
    base_inv = inverse(base,mod)
    babysteps, giantsteps = [],[]
    for i in range(Q+1):
        babysteps.append(powermod(base,Q*i,mod))
    for l in range(Q):
        giantsteps.append((arg*powermod(base_inv,l,mod))%mod)
        if giantsteps[l] in babysteps:
            k = babysteps.index(giantsteps[l])
            return k*Q+l
        
def prho(base,arg,mod):
    b, y, z = [], [], []
    b.append(arg)
    y.append(0)
    z.append(1)
    found = k = 0
    while not found:
        if b[k]%3 == 0:
            b.append( (b[k]*b[k])%mod )
            y.append( (2*y[k])%(mod-1) )
            z.append( (2*z[k])%(mod-1) )
        elif b[k]%3 == 1:
            b.append( (b[k]*arg)%mod )
            y.append( (y[k])%(mod-1) )
            z.append( (z[k]+1)%(mod-1) )
        else:  # entspricht b[k]= 2 mod 3
            b.append( (b[k]*base)%mod )
            y.append( (y[k]+1)%(mod-1) )
            z.append( (z[k])%(mod-1) )
        if b[k+1] in b[:k+1]:
            j = k+1
            i = b.index(b[j])
            found = 1
        k = k + 1
    
    g, c, _ = gcd(z[j]-z[i],mod-1) 
    mod2 = (mod-1)/g                
    x = (((y[i]-y[j])//g)*c) % mod2
    for l in range(g):
        if powermod(base,(x+mod2*l),mod)==arg:
            return x+mod2*l
        
def fastpow(base,expo):
    binary = []
    while expo > 0:
        binary.append(expo%2)
        expo = expo / 2
    erg = 1
    for i in range(len(binary)):
        erg = erg*erg
        if binary[-(i+1)] == 1:
            erg = erg * base
    return erg

def pohell(base,arg,mod):
    factors = factorlist(mod-1)
    #fset = set(factors)
    factorbase, factorexp = [], []
    #for i in range(len(fset)):
    #    factorbase.append(fset.pop())                
    for i in factors:
        if i not in factorbase:
            factorbase.append(i)
    for i in range(len(factorbase)): 
        factorexp.append(factors.count(factorbase[i]))
    n , gamma, alpha, x, faktoren = [], [], [], [], []
    for i in range(len(factorbase)):
        faktoren.append(fastpow(factorbase[i],factorexp[i]))
        n.append((mod-1)/faktoren[i])
        gamma.append(powermod(base,n[i],mod))
        alpha.append(powermod(arg,n[i],mod))
        x.append(bsgs(gamma[i],alpha[i],mod))
        if i == 0:
            chinesenliste = [x[i],faktoren[i]]
        elif i > 0:
            chinesenliste = chinese(x[i],faktoren[i],chinesenliste[0],chinesenliste[1])
    return chinesenliste[0]
    
            
def dl(base,n,modulus,limit):
    if n >= modulus:
        n = n % modulus
    if n <= 0:
        return 'impossible'
    for i in range(limit):
        if powermod(base, i, modulus)==n:
            return i
    return 'nothing found'

def Z_sqrt_d_gcd(a1,a2,b1,b2,d):
    a, b = Zsqrtd(a1,a2,d), Zsqrtd(b1,b2,d)
    x1,x2,y1,y2=1,0,0,1                 
    while b!=0:                            
        q=a//b                          
        x1,x2,y1,y2=x2,x1-q*x2,y2,y1-q*y2  #x1->x2,x2->x3,y1->y2,y2->y3
        a,b=b,a%b
    return a,x1,y1 # ggt>0 and  ax1+bx2=ggt

def rsaKeymaker(bits=768):
    bits = bits//2 - 1
    p = nextprime(bigrand(int(bits*math.log10(2))))
    q = nextprime(bigrand(int(bits*math.log10(2))))    
    n = p * q
    phi = (p-1)*(q-1)
    e = nextprime(bigrand(int(bits*2*math.log10(2))))
    d = inverse(e,phi)
    publicKey = (n,e)
    privateKey = (n,d)
    return (publicKey, privateKey)

def rsaKey_to_file(name,length):
    dir = 'c:'
    if not os.path.isdir(dir+'\\RSA'):
        os.mkdir(dir+'\\RSA')
    if os.path.isfile(dir+'\\RSA\\keys\\prv\\own.prv'):
        return 2
    if not os.path.isdir(dir+'\\RSA\\keys'):
        pubKey, prvKey = rsaKeymaker(int(length))
        os.mkdir(dir+'\\RSA\\keys')
        pub = open(dir+'\\RSA\\keys\\'+name+'.pub','w')
        pub.write(str(pubKey[0])+'\n'+str(pubKey[1])+'\n')
        pub.close()
        os.mkdir(dir+'\\RSA\\keys\\prv')
        prv = open(dir+'\\RSA\\keys\\prv\\own.prv','w')
        prv.write(str(prvKey[0])+'\n'+str(prvKey[1])+'\n')
        prv.close()
        return 1
    else:
        return 0

def rsaEncrypt(msg, key):
    a=[]
    if len(msg) % 16 != 0:
            msg += " "*(16-len(msg)%16)
    while len(msg) != 0:
        tmp= ""
        for i in range(16):
            char = ord(msg[i])
            if char>=100:
                tmp += str(char)
            else:
                tmp += "0"+str(char)
        base=long(tmp)    
        a.append(to_base62(pow(base,key[1],key[0])))
        msg = msg[16:]
    return a

def rsaEn_to_file(msg, recp):
    dir = 'c:\\RSA\\keys'
    if os.path.isfile(dir+'\\'+str(recp)+'.pub'):
        file = open(dir+'\\'+str(recp)+'.pub','r')
        lines = file.readlines()
        file.close()
        N = long(lines[0].replace('\n',''))
        e = long(lines[1].replace('\n',''))
        cipher = rsaEncrypt(msg,(N,e))
        if not os.path.isdir('c:\\RSA\\msg'):
            os.mkdir('c:\\RSA\\msg')
        if os.path.isfile('c:\\RSA\\msg\\msg.py'):
            os.remove('c:\\RSA\\msg\\msg.py')
        file = open('c:\\RSA\\msg\\msg.py','w')
        for i in range(len(cipher)):
            file.write(cipher[i]+'\n')
        file.close()
        return 1
    else:
        return 0
            

def rsaDecrypt(msg, key):
    a = ""
    for i in range(1,len(msg)+1):
        tmp = pow(from_base62(msg[-i]),key[1],key[0])
        while tmp > 0:
            a = chr(tmp%1000) + a
            tmp //= 1000
    return a

def rsaDe_from_file():
    if os.path.isfile('c:\\RSA\\keys\\prv\\own.prv') and os.path.isfile(sys.path[0]+'\\my\\msg.py'):
        homedir = sys.path[0]+'\\my'
        key = open('c:\\RSA\\keys\\prv\\own.prv','r')
        lines = key.readlines()
        key.close()
        N = long(lines[0].replace('\n',''))
        d = long(lines[1].replace('\n',''))
        msg = open(homedir+'\\msg.py','r')
        lines2 = msg.readlines()
        msg.close()
        for i in range(len(lines2)):
            lines2[i]=lines2[i].replace('\n','')
        return rsaDecrypt(lines2,(N,d))
    else:
        return None
    

def to_base62(n):
    alphabet=['0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F','G','H',
              'I','J','K','L','M','N','O','P','Q','R','S','T','U','V','W','X','Y','Z',
              'a','b','c','d','e','f','g','h','i','j','k','l','m','n','o','p','q','r',
              's','t','u','v','w','x','y','z']
    v=""
    while n > 0:
        v = (alphabet[n%62])+v
        n //= 62
    return v

def from_base62(v):
    alphabet=['0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F','G','H',
              'I','J','K','L','M','N','O','P','Q','R','S','T','U','V','W','X','Y','Z',
              'a','b','c','d','e','f','g','h','i','j','k','l','m','n','o','p','q','r',
              's','t','u','v','w','x','y','z']
    n = 0
    for i in range(len(v)):
        n = 62 * n + alphabet.index(v[i])
    return n
    
def bigrand(n):
    erg = str(random.randint(1,9))
    for i in range(n-1):
        erg += str(random.randint(0,9))
    return long(erg)    
    
def factorlist(m):
    list_of_factors = []
    sq = math.sqrt(m)
    while m%2 == 0:
        list_of_factors.append(2)
        m = m//2
        sq = math.sqrt(m)
    d = 3
    while  d <= sq:
        if m % d == 0:
            list_of_factors.append(d)
            m = m//d
            sq = math.sqrt(m)
        else:
            d = d+2
    if m > 1:
            list_of_factors.append(m)
    return list(list_of_factors)

def primes(n):
    if n <= 1: return []
    X = [i for i in range(3,n+1) if i%2 != 0]     
    P = [2]                                       
    sqrt_n = math.sqrt(n)                         
    while len(X) > 0 and X[0] <= sqrt_n:          
        p = X[0]                                  
        P.append(p)                               
        X = [a for a in X if a%p != 0]            
    return P + X                                  

def isprime(n):
         if n % 2 == 0 and not n == 2: return (0)
         else:
                 limit = long(math.sqrt(n)) + 1
                 for l in range (3, limit, 2):
                         if n % l == 0: return (0)
                 return (1)
 
def nextprime(n):
         np = n + 1
         while np:
                 #if np%2!=0 or np%3!=0 or np%5!=0: 
                 #if miller_rabin(np): return (np)#if isprime(np): return (np)
                 if miller_rabin(np): return (np)
                 np = np + 1
def sgn(n):
     if n > 0:
         return 1
     elif n < 0:
         return -1
     else:
         return 0
    
def gcd(a,b):
    x1,x2,y1,y2=1,0,0,1                 
    sgna=sgn(a)  #unnoetig? 
    sgnb=sgn(b)
    a=abs(a)    #unnoetig? 
    b=abs(b)
    while b!=0:                            
        q=a//b                          
        x1,x2,y1,y2=x2,x1-q*x2,y2,y1-q*y2  #x1->x2,x2->x3,y1->y2,y2->y3
        a,b=b,a%b
    return a,x1*sgna,y1*sgnb # ggt>0 and  ax1+bx2=ggt

def lcm(a,b): 
    return a//gcd(a,b)[0]*b

def binomial(n,k):
    out = 1
    for i in range(1,k+1):
        out = (out * (n-i+1))//i
    return out

def factorial(n):
    f = 1
    for i in range(1, n+1):
        f = f*i
    return f

def fibonacci(n):
    a,b = 0,1
    list=[1,1]
    for x in range(1,n-1):
        a,b=b,a+b
        list.append(a+b)
    return list
    
def eulerphi(n):
    negative_list=[]
    for p in factorlist(n):
        if not p in negative_list:
            n = n*(p-1)//p
            negative_list.append(p) 
    return n

def contfrac(number, steps=8):
    l = []
    for i in range(steps):
        i = long(math.floor(number))
        l.append(i)
        if abs(number-i) < 0.00000001:
             break
        number = 1/(number-i)
    return l

def legendre(a,m): 
    a=a%m
    t=1
    while a!=0:
        while a%2==0:
             a=a//2
             if m%8==3 or m%8==5:
                 t=-t
        a,m=m,a
        if a%4==3 and m%4==3:
             t=-t
        a=a%m
    if m==1:
        return t
    return 0

def is_primitive_root( modulus, b = 2L ):
    less_one = modulus - 1
    list_of_factors = factorlist( modulus - 1 )
    for factor in list_of_factors:
        if pow(b, less_one//factor, modulus) == 1:
            return 0
    return 1
 
def primitive_root(modulus):
    b = 2
    while ( is_primitive_root( modulus, b) == 0):
        b = b+1
    return b

def inverse(a,b):
    if a<0:
        while a<0:
             a += b
    y = gcd(a,b)
    if y[0] != 1:
        return None
    return y[1]%b

def chinese(a,m,b,n):
    modulus= lcm(m,n)
    d, u, v = gcd(m,n) #um+vn=d > 0 
    if (b-a)%d == 0:  # d must divide b-a, otherwise there is no solution
        c= (b-a)//d    # c(um+vn)=cd=b-a, hence x=b-vcn=a+ucm is the solution 
        return (a+u*c*m)%modulus, modulus
    else:
        pass

def solve_linear(a,b,n):
    g, c, _ = gcd(a,n)                 
    if b%g != 0: return None
    return ((b//g)*c) % n            

def powermod(a, m, n):
    ans = 1
    apow = a
    if (m < 0 or n < 1): return None
    while m != 0:
        if m%2 != 0:
            ans = (ans * apow) % n            
        apow = (apow * apow) % n              
        m //= 2
    return ans % n

def sqrtmod(a, p): #p Primzahl
    a %= p
    if p == 2: return a 
    if legendre(a, p) != 1: return None
    if p%4 == 3: return powermod(a, (p+1)//4, p)

    def mul(x, y):   
        return ((x[0]*y[0] + a*y[1]*x[1]) % p, (x[0]*y[1] + x[1]*y[0]) % p)
    
    def pow(x, n):  
        ans = (1,0)
        xpow = x
        while n != 0:
           if n%2 != 0: ans = mul(ans, xpow)
           xpow = mul(xpow, xpow)
           n //= 2
        return ans

    while True:
        z = random.randrange(2,p)
        u, v = pow((1,z), (p-1)//2)
        if v != 0:
            vinv = inverse(v, p)
            for x in [-u*vinv, (1-u)*vinv, (-1-u)*vinv]:
                if (x*x)%p == a: return x%p
#           assert False, "Bug in sqrtmod."
            return None
            

def real_class_number(d):
    #class number of real quadratic field with discriminant d =0,1 mod 4
    counter=0
    arrayB=[] # sammelt b's
    array2C=[] # sammelt 2c's
    class_number=0
    f=long(math.floor(math.sqrt(d)))
    t=f*f
    if t==d:
        #print "d ist Quadrat"
        return
    u=d%4
    if u<>0 and u<>1:
	#print "d ist nicht 0,1 mod 4"
        return
    if u==1:
        e=1
    else:
        e=2
    g = f//2;
    for a in range(1,f+1):      
        for b in range(e,f+1,2):
            h = b*b - d
            i = 2*a
            j = 4*a
            if (h%j == 0) and (a<= g and f-i < b): #test reduced ! 
                c = -h//j
                done=0
                for i in range(0,counter):
                    if (b==arrayB[i]) and (2*c==array2C[i]):
                        done=1
                if (gcd(gcd(a, b)[0],c)[0] == 1) and done==0:              
                    u,v=b,2*c
                    s,r=v,u
                    k=counter 
                    while (counter==k) or (s<>v or u<>r): 
                        a=(f+u)//v
                        u=a*v-u
                        v=(d-u*u)//v 		
                        arrayB.append(u)
                        array2C.append(v)
                        counter=counter+1	    	
                    class_number = class_number + 1
                    #print "Reduzierte Form [", class_number, "]: (", a, ",", b, ",", -c, ")"
                    #print "Zykellaenge=", counter-k
            else:
                pass
    return class_number
    

def imaginary_class_number(d):
# class number of imaginary quadratic field with discriminant -d =0,1 mod 4
#   d = abs(d)
    f=long(math.floor(math.sqrt(d)))
    t=f*f
    if t==d:
        #print "d ist Quadrat"
        return
    u=-d%4
    if u<>0 and u<>1:
	#print "d ist nicht 0,1 mod 4"
        return
    h = 0
    g = 1
    if d%4 == 0:
        b = 0
    else:
        b = 1
    bound = math.sqrt(d//3)
    while b <= bound:
        q =(b**2+d)//4
        a = b
        if a <= 1:
            a = 1
        while a**2 <= q:
            if q%a == 0:
                t=q//a
                ggt=gcd(a,b)[0]
                ggt=gcd(ggt,t)[0]
                if ggt > 1:
                    g = 0
                if g == 1:
                    if a == b or a**2 == q or b == 0:
                        h = h+1
                    else:
                        h = h+2
                else:
                    g = 1
            a = a+1
        b = b+2
    return h

def pollard(N, m):
    for a in [2, 3]:
        x = powermod(a, m, N) - 1
        g = gcd(x, N)[0]
        if g != 1 and g != N:
            return g
    return N

def randcurve(m):
    if m<2: return None
#   assert m > 2, "m must be > 2."
    a = random.randrange(m)
    while gcd(4*a**3 + 27, m)[0] != 1:
        a = random.randrange(m)
    return (a, 1, m), (0,1)

def elliptic_curve_method(N, m, tries=5):
    for _ in range(tries):                     
        E, P = randcurve(N)                    
        try:                                   
            Q = ellcurve_mul(E, m, P)          
        except ZeroDivisionError, x:           
            g = gcd(x[0],N)[0]                    
            if g != 1 or g != N: return g
    return N            

def ellcurve_add(E, P1, P2):
    a, b, p = E
    if p <=2: return None
#   assert p > 2, "p must be odd."
    if P1 == "Identity": return P2
    if P2 == "Identity": return P1
    x1, y1 = P1; x2, y2 = P2
    x1 %= p; y1 %= p; x2 %= p; y2 %= p
    if x1 == x2 and y1 == p-y2: return "Identity"
    if P1 == P2:
        if y1 == 0: return "Identity"
        lam = (3*x1**2+a) * inverse(2*y1,p)
    else:
        lam = (y1 - y2) * inverse(x1 - x2, p)
    x3 = lam**2 - x1 - x2
    y3 = -lam*x3 - y1 + lam*x1
    return (x3%p, y3%p)

def ellcurve_mul(E, m, P):
    if m<0: return None
#   assert m >= 0, "m must be nonnegative."
    power = P
    mP = "Identity"
    while m != 0:
        if m%2 != 0: mP = ellcurve_add(E, mP, power)
        power = ellcurve_add(E, power, power)
        m /= 2
    return mP

def lcm_to(B):
    ans = 1
    logB = math.log(B)
    for p in primes(B):
        ans *= p**long(logB/math.log(p))
    return ans

def miller_rabin(n, num_trials=4):
    if n < 0: n = -n
    if n in [2,3]: return True
    if n <= 4: return False
    m = n - 1
    k = 0
    while m%2 == 0:
        k += 1; m /= 2
    # Now n - 1 = (2**k) * m with m odd
    for i in range(num_trials):
        #ab hier
        if n > 1000000000:
            laenge = len(str(n-1))
            endlaenge = random.randint(2,laenge-1)
            a = bigrand(endlaenge)
        else:
        #fuer diese zeile
            a = random.randrange(2,n-1)                  
        apow = powermod(a, m, n)
        if not (apow in [1, n-1]):            
            some_minus_one = False
            for r in range(k-1):              
                apow = (apow**2)%n
                if apow == n-1:
                    some_minus_one = True
                    break                     
        if (apow in [1, n-1]) or some_minus_one:
            prob_prime = True
        else:
            return False
    return True    

def gen(n,c=1):
    x = 1
    while True:
        x = (x**2 + c) % n
        yield x

def rho(n, max=500, maxc=10):  # Pollard's Rho Method
    seqslow = gen(n) 
    seqfast = gen(n)
    trials = 0
    c = 1   
    while c < maxc:
        while trials < max:        
            xb = seqslow.next() 
            seqfast.next()
            xk = seqfast.next() 
            trials += 1
            diff = abs(xk-xb)
            if not diff:  continue
            d = gcd(diff,n)[0]
            if n>d>1:                
                return rho(d, max, maxc) + rho(n//d, max, maxc)
        c += 1
        seqslow = gen(n,c) 
        seqfast = gen(n,c)
        trials = 0
    return [n]    # failure to factor


def sum(v):
    if len(v) == 0:
        return 0
    return reduce(operator.add, v)

def ggt(a,b):
    while b !=0:
        a,b=b,a%b
    return a

def isInteger(n):
    return type(n) in (int, long)

class Zsqrtd:
    def __init__(self,a,b=0,d=1):
        self.a,self.b,self.d = int(a),int(b),int(d)
        if self.d >= 0 and math.ceil(math.sqrt(self.d))==math.floor(math.sqrt(self.d)):
            self.z = (a+b*math.sqrt(self.d),0)
            self.d = 1
        else:
            self.z,self.d = (a,b),d
            
    def __ne__(self,v):
        if not isinstance(v, Zsqrtd):
            self,v = self._unify(v)
        if self.z[0]!=v.z[0] or self.z[1]!=v.z[1]:
            return True
        else:
            return False
        
    def _unify(self, other):
        return (self,Zsqrtd(other,d=self.d))    
        
    def __repr__(self):
        if self.z[1] < 0:
            return '%d - %d{%d}' %(self.z[0],int(math.fabs(self.z[1])),self.d)
        elif self.z[1] > 0:
            return '%d + %d{%d}' %(self.z[0],self.z[1],self.d)
        else:
            return '%d'%self.z[0]
    
    def __neg__(self):
        return Zsqrtd(-self.z[0],-self.z[1],self.d)
    
    def conj(self):
        return Zsqrtd(self.z[0],-self.z[1],self.d)
        
    def __add__(self,v):
        if not isinstance(v, Zsqrtd):
            self,v = self._unify(v)
        if self.d == v.d or v.d==1 or self.d==1:
            return Zsqrtd(self.z[0] + v.z[0],self.z[1] + v.z[1],self.d)
        else:
            return 'Fehler'
        
    def __radd__(self,v):
        return self + v
    
    def __sub__(self,v):
        return self + (-v)
        
    def __rsub__(self,v):
        return v + (-self)
    
    def __mul__(self,v):
        if not isinstance(v, Zsqrtd):
            self,v = self._unify(v)
        if self.d == v.d or v.d==1 or self.d==1:
            return Zsqrtd(self.z[0]*v.z[0]+self.z[1]*v.z[1]*self.d,self.z[0]*v.z[1]+self.z[1]*v.z[0] ,self.d)
        else:
            return 'Fehler'
        
    def __rmul__(self,v):
        return self * v
    
    def __div__(self,v):
        if not isinstance(v, Zsqrtd):
            self,v = self._unify(v)
        if self.d == v.d or v.d==1 or self.d==1:
            numerator = self * v.conj()
            denominator = v.z[0]**2 - v.z[1]**2*v.d
            erg = numerator.z[0]//denominator
            if math.fabs(numerator.z[0]%denominator) > math.fabs(denominator//2):
                erg += 1
            erg2 = numerator.z[1]//denominator
            if math.fabs(numerator.z[1]%denominator) > math.fabs(denominator//2):
                erg2 += 1
            return Zsqrtd(erg, erg2 ,self.d)
        else:
            return 'Fehler'
        
    def __rdiv__(self,v):
        if not isinstance(v, Zsqrtd):
            self,v = self._unify(v)
        if v.d == self.d or v.d==1 or self.d==1:
            numerator = v * self.conj()
            denominator = self.z[0]**2 - self.z[1]**2*self.d
            erg = numerator.z[0]//denominator
            if math.fabs(numerator.z[0]%denominator) > math.fabs(denominator//2):
                erg += 1
            erg2 = numerator.z[1]//denominator
            if math.fabs(numerator.z[1]%denominator) > math.fabs(denominator//2):
                erg2 += 1
            return Zsqrtd(erg, erg2 ,v.d)
        else:
            return 'Fehler'
        
    def __mod__(self,v):
        if not isinstance(v, Zsqrtd):
            self,v = self._unify(v)
        if self.d == v.d or v.d==1 or self.d==1:
            return self - (self/v)*v
        else:
            return 'Fehler'
        
    def __rmod__(self,v):
        if not isinstance(v, Zsqrtd):
            self,v = self._unify(v)
        if self.d == v.d or v.d==1 or self.d==1:
            return v - (v/self)*self
        else:
            return 'Fehler'
    
    def __truediv__(self, b):   return self.__div__(b)
    def __rtruediv__(self, b):  return self.__rdiv__(b)
    def __floordiv__(self, b):  return self.__div__(b)
    def __rfloordiv__(self, b): return self.__rdiv__(b)
        
class Rational (object): 
    def operators(self):
        pass

    def __init__(self, p, q=1):
        if type(p) == str:   #  z.B. '3/17'
            l = p.split('/')
            p = long(l[0])
            if len(l) > 1:
                q = long(l[1])
        p=long(p)        
        if q < 0:
            p,q = -p, -q  
        d = ggt(abs(p), q)     
        if d > 1:         
            p //= d
            q //= d
        self.p = p
        self.q = q

    def __abs__(self):
        return Rational(abs(self.p), self.q)

 #   def abs(self):
 #       return A - B

    def __float__(self):
        return float(self.p) / float(self.q)

    def __int__(self):
#       assert self.q > 0
        return (self.p + self.q/2) // self.q

    def __long__(self):
#       assert self.q > 0
        return (self.p + self.q/2) // self.q

    def _unify(self, other):
        return (self, Rational(other))

    def __repr__(self):
        if self.q == 1:
            return str(self.p)
        else:
            return '%d/%d' % (self.p, self.q)

    def __cmp__(self, b):
        if type(b) == float:
            return cmp(float(self), b)
        if not isinstance(b, Rational):
            self,b = self._unify(b)
        return cmp(self.p*b.q, b.p*self.q)

    def __neg__(self):
        return Rational(-self.p, self.q)

    def neg(self):
        return -self

    def __add__(self, b):
        if type(b) == float:
            return float(self) + b
        if not isinstance(b, Rational):
            self,b = self._unify(b)
        return Rational(self.p*b.q + b.p*self.q, self.q*b.q)

    def add(self,b):
        return self + b

    def __radd__(self, b):
        return self + b

    def __sub__(self, b):
        if type(b) == float:
            return float(self) - b
        if not isinstance(b, Rational):
            self,b = self._unify(b)
        return Rational(self.p*b.q - b.p*self.q, self.q*b.q)

    def sub(self,b):
        return self - b

    def __rsub__(self, a):
        if type(a) == float:
            return a - float(self)
        if not isinstance(a, Rational):
            self,a = self._unify(a)
        return a - self

    def __mul__(self, b):
        if type(b) == float:
            return float(self) * b
        if not isinstance(b, Rational):
            self,b = self._unify(b)
        return Rational(self.p*b.p, self.q*b.q)

    def mul(self,b):
        return self * b

    def __rmul__(self, a):
        return self * a
    
    def __div__(self, b):
        if type(b) == float:
            return float(self) / b
        if not isinstance(b, Rational):
            self,b = self._unify(b)
        return Rational(self.p*b.q, self.q*b.p)

    def div(self,b):
        return self / b

    def __rdiv__(self, a):
        if type(a) == float:
            return a / float(self)
        if not isinstance(a, Rational):
            self,a = self._unify(a)
        return Rational(a.p*self.q, a.q*self.p)

    def __truediv__(self, b):   return self.__div__(b)
    def __rtruediv__(self, b):  return self.__rdiv__(b)
    def __floordiv__(self, b):  return self.__div__(b)
    def __rfloordiv__(self, b): return self.__rdiv__(b)

def toNumber(self):
    if '.' in self or 'e' in self or 'E' in self:
        return float(self)
    else:
        return Rational(self)


class Vector (list):
    
    def operators(self):
        pass
    
    def __init__(self, n):
        if isinstance(n, list):
            list.__init__(self, n)
        elif isinstance(n, tuple):
            list.__init__(self, list(n))
            
    def __repr__(self):
        """ Darstellung des Objekts als String
        """
        return 'Vector(%s)' % list.__repr__(self)

    def __invert__(self):
        """ <b>v.__invert__()</b><br/>
            <i>Transposition</i><br/>
            berechnet die transponierte Matrix (Spaltenvektor)<br/>
            Auch Operatorschreibweise: <b>~v</b><br/>
        """
        return Matrix([[a] for a in self])

    def transp(self):
        """ <b>v.transp()</b><br/>
            <i>berechnet die transponierte Matrix (Spaltenvektor)</i><br/>
            kuerzer: ~v
        """
        return Matrix([[a] for a in self])

    def concat(self, y):
        """ <b>.concat(v)</b><br/>
            <i>Verkettung mit dem Vektor v</i>
        """
        return Vector(list(self)+list(y))
    def __neg__(self):
        return Vector([-x for x in self])

    def __add__(self, y):
        n = len(self)
#       assert n == len(y)
        return Vector([self[i]+y[i] for i in range(n)])

    def __radd__(self, b):
        return self + b

    def __sub__(self, y):
        return self + (-y)

    def __rsub__(self, x):
        return x - self

    def __mul__(self, y):
        if isinstance(y, Vector): 
            n = len(self)
#           assert n == len(y)
            return sum([self[i]*y[i] for i in range(n)])
        else:
            return Vector([y*self[i] for i in range(len(self))])

    def __rmul__(self, c):
        return self * c

    def __imul__(self, c):
        for i in range(len(self)):
            self[i] *= c
        return self



class Matrix (Vector): 
    def operators(self):
        pass
    
    def __init__(self, v):
        Vector.__init__(self, [Vector(v[i]) for i in range(len(v))])
        self._makeRational()

    def null(self, n=0):
        if n == 0:
            n = self
        return Matrix([ [0 for k in range(n)]
                        for i in range(self)])

#    @staticmethod
    def fromFunction(self, n, fn, offset=0):
        return Matrix([ [fn(i,k) for k in range(n)]
                        for i in range(self)])

#    @staticmethod
    def fromString(self):
        return Matrix([ [toNumber(a) for a in row.split(',')]
                            for row in self.split(';')])

    def _makeRational(self):
        for i in self.rowRange():
            for k in self.colRange():
                if isInteger(self[i,k]):
                    self[i,k] = Rational(self[i,k])

    def copy(self):
        return Matrix([ [self[i,k] for k in self.colRange()]
                        for i in self.rowRange()])

    def isSquare(self):
        return self.height() == self.width()

    def str(self, i,k):
        if type(self[i,k]) == float:
            return '%g' % self[i,k]
        else:
            return str(self[i,k])

    def _getColWidth(self, k):
        return max([len(self.str(i,k)) for i in self.rowRange()])

    def _pp(self, name=''):
        w = [self._getColWidth(k) for k in self.colRange()]
        strW = sum(w) + 3 * (self.width()-1)
        if self.height() == 1:
            left, right = '[', ']'
        else:
            mid = (2*self.height()-3) * '|'
            left, right = '/'+mid+'\\',  '\\'+mid+'/'
        s = ''
        fill = len(name) * ' '

        n = 2*self.height()-1
        for j in range(n):
            if j == self.height()-1:
                s += name
            else:
                s += fill
            i = j // 2
            if j % 2 == 0:
                s += left[j] + ' '
                for k in self.colRange():
                    s += self.str(i,k).center(w[k])
                    if k < self.width()-1:
                        s += '   '
                s += ' ' + right[j] +  '\n'
            else:
                s += '|' + (strW+2) * ' ' +  '|\n'
        return s

    def __repr__(self):
        return self._pp()

    def __getitem__(self, i):
        if isinstance(i, tuple):
            return self[i[0]][i[1]] 
        else:
            return Vector.__getitem__(self, i)

    def __setitem__(self, i, x):
        if isinstance(i, tuple):
            self[i[0]][i[1]] = x
        else:
            Vector.__setitem__(self, i, x)

    def height(self): #Zeilenzahl
        return len(self)

    def width(self): #Spaltenzahl
        return len(self[0])

    def rowRange(self):
        return range(self.height())

    def colRange(self):
        return range(self.width())

    def __neg__(self):
        return Matrix([ [-self[i,k] for k in self.colRange()]
                        for i in self.rowRange()])

    def __cmp__(self, B):
        if self.heigth() != B.heigth():
            return self.heigth() - B.heigth()
        for i in self.rowRange():
            if self[i] != self[i]:
                return cmp(self[i], B[i])
        return 0

    def concat(self, B):
#       assert self.height() == B.height()
        return Matrix([ list(self[i])+list(B[i]) for i in self.rowRange()])

    def __add__(self, B):
#       assert self.height() == B.height() and self.width() == B.width()
        return Matrix([ [self[i,k]+B[i,k] for k in self.colRange()]
                        for i in self.rowRange()
                      ])

    def __sub__(self, B):
        return self + (-B)

    def __mul__(self, B):
        # Wenn beide Operanden Matrizen sind wird Matrizenmultiplikation durchgefuehrt,
        # sonst Skalarmultiplikation
        if isinstance(B, Matrix): # Matrizenmultiplikation
#           assert self.width() == B.height()
            return Matrix(
                [ [ sum([self[i,j] * B[j,k] for j in self.colRange()])
                    for k in B.colRange()]
                  for i in self.rowRange()])
        else: # Skalarmultiplikation
            return Matrix([ [B*self[i,k] for k in self.colRange()]
                            for i in self.rowRange()])

    def __rmul__(self, B):
#       assert not isinstance(B, Matrix)
        return self * B

    def transp(self):
        return Matrix([ [self[i,k] for i in self.rowRange()]
                        for k in self.colRange()])

    def __invert__(self):
        return self.transp()

    def submatrix(self, i, k, m, n):
#       assert i+m-1 < self.height() and k+n-1 < self.width()
        return Matrix([ [self[i1,k1] for k1 in range(k,k+n)]
                        for i1 in range(i,i+m)])

    def complement(self, i, k):
        return Matrix([ [self[i1,k1] for k1 in self.colRange() if k1 != k]
                        for i1 in self.rowRange() if i1 != i])

    def minor(self, i, k):
        return self.complement(i,k).det()

    def cofactor(self, i, k):
        return (-1)**(i+k) * self.minor(i,k)

    def adjoint(self):
#       assert A.isSquare()
        return Matrix([ [self.cofactor(k,i) for k in self.rowRange()]
                        for i in self.colRange()])
                        
    def det(self):
        if not self.isSquare():
            raise 'Matrix ist nicht quadratisch'
        n = self.height()
        if n == 1:
            return self[0,0]
        else:
            return sum([self[0,k] * self.cofactor(0,k)
                        for k in self.colRange()])

    def _gaussElim0(self, jordan=True):
        m, n = self.height(), self.width()
        k = 0
        for i0 in self.rowRange():
            k = i0
            # oberste Zeile fuer fuehrende Eins passend durchmultiplizieren
            self[i0] *= 1/self[i0,k]

            # Passende Vielfache zu anderen Zeilen addieren, so dass
            # unterhalb der fuehrenden Eins Nullen entstehen
            for i in range(i0+1, m):
                self[i] -= self[i0] * self[i,k]
            if jordan:
                # Gauss-Jordan: dto. auch oberhalb der fuehrenden Eins
                for i in range(i0):
                    self[i] -= self[i0] * self[i,k]

    def gaussElim(self, jordan=True):
        m, n = self.height(), self.width()
        k = -1
        for i0 in self.rowRange():
            # Bestimme die am weitesten links stehende Spalte k, die
            # (ab Zeile i0) von Null verschiedene Elemente enthaelt:
            k += 1
            aMax, iMax = 0, i0
            while aMax == 0 and k < n:
                for i in range(i0,m):
                    if abs(self[i,k]) >aMax: 
                        aMax, iMax =  abs(self[i,k]), i
                if aMax == 0:
                    k += 1
            if k < n:
                # Oberste Zeile mit der vertauschen, die das (dem Betrag nach)
                # groesste Element in Spalte i0 enthaelt
                self[iMax], self[i0] = self[i0], self[iMax]

                # oberste Zeile fuer fuehrende Eins passend durchmultiplizieren
                self[i0] *= 1/self[i0,k]

                # Passende Vielfache zu anderen Zeilen addieren, so dass
                # unterhalb der fuehrenden Eins Nullen entstehen
                for i in range(i0+1, m):
                    self[i] -= self[i0] * self[i,k]
                if jordan:
                    # Gauss-Jordan: dto. auch oberhalb der fuehrenden Eins
                    for i in range(i0):
                        self[i] -= self[i0] * self[i,k]

    def inverse(self):
        if not self.isSquare():
            raise 'not a square matrix'
        n = self.height()
        I=Matrix([ [long(i==k) for k in range(n)] for i in range(n)])
        B = self.concat(I)
        B.gaussElim()
        if B.submatrix(0,0,n,n) != I:
            raise 'matrix not invertible'
        return B.submatrix(0,n,n,n)

    def solve(self, b):
        if not isinstance(b, Matrix):
#           assert isinstance(b, Vector)
            b = b.transp()
        v=self.inverse() * b
        return v.transp()

def test(n):
    a, bs, rho,poh = [],[],[],[]
    import timeit
    for i in range(n):
        a.append(random.randint(1,15485862))
        teststr = 'print bsgs(6,' + str(a[i]) + ',15485863)'
        teststr1= 'print prho(6,' + str(a[i]) + ',15485863)'
        teststr2= 'print pohell(6,' + str(a[i]) + ',15485863)'
        t = timeit.Timer(teststr,"from mathe import bsgs")
        s = timeit.Timer(teststr1,"from mathe import prho")
        r = timeit.Timer(teststr2,"from mathe import pohell")
        bs.append(t.timeit(1))
        rho.append(s.timeit(1))
        poh.append(r.timeit(1))
    print 'as:',a
    print 'bs:',bs
    print 'pr:',rho 
    print 'ph:',poh
    
#rsaKey_to_file('123')
#rsaEn_to_file('das ist easdasdasdasdasdasdasdasdasdasdasdasdasdasdasdasdasdin test tst',123)
#print rsaDe_from_file()