def fpsp(n):
    #Si n es par acaba el programa
    if n%2==0:
        print str(n) + " es par."
        return (2,False)
    #Si n es impar
    else:
        #Calcula s y t tales que n-1=2^st con t impar
        s = mpot(2,n-1)
        t = (n-1)/pow(2,s)

        #Tomamos una base b aleatoria entre 1 y n-1
        base = random.randint(1,n-1)

        #Comprobamos si gcd(base,n) es 1
        #Si no lo es
        if not sp.gcd(base,n)==1:
            print str(sp.gcd(n,base)) + " es divisor de " + str(n) + "."
            return (base,False)
        #Si lo es
        else:
            #Comprueba ahora b^t=+-1mod n
            #Si lo es
            if pow(base,t,n) == 1 or pow(base,t,n) == (-1%n):
                return (base,True)
            #Si no lo es
            else:
                i=1
                while i<=s-1 :
                    if pow(base,t*pow(2,i),n) == (-1%n):
                        return (i,base,True)
                    i=i+1
            return (base,False)

def factordet(Mat):
    """
    Extract common factors of the determinant of the matrix.

    Examples
    ========

    >>> from sympy import *
    >>> x, y, z = symbols('x y z')
    >>> A = Matrix([[x,x*y],[y*z,z]])
    >>> factordet(A)
    x*z
    """
    fact = 1
    ncols = Mat.cols
    for i in range(ncols):
        col = Mat.col(i)
        common = gcd(list(col))
        if (common != 0)&(common != 1):
            fact *= common
            Mat[i] = Matrix(list(map(cancel, col/common)))
    for j in range(Mat.rows):
        row = Mat.row(j)
        common = gcd(list(row))
        if (common != 0)&(common != 1):
            fact *= common
            Mat[j*ncols] = Matrix([list(map(cancel, row/common))])
    return fact

 def sympy_gcd(self, p, q):
     sympy_p = self.sympy_from_upolynomial(p)
     sympy_q = self.sympy_from_upolynomial(q)
     if (p.ring().modulus() is None):
         return sympy.gcd(sympy_p, sympy_q)
     else:
         return sympy.gcd(sympy_p, sympy_q, modulus=p.ring().modulus())

	def __generate_polydesc(self):
		# Symbols
		a,x,y = sp.symbols('a x y')
		A,X,MAX = sp.symbols('A X MAX')
		
		# Attractor limits
		x_min = (4*a**2 - a**3) / 16
		x_max = a / 4
		
		# Functions
		#  f  : Attr  -> Attr
		#  g  : Attr  -> [0,1]
		#  gi : [0,1] -> Attr
		#  F  = g o f o gi : [0,1] -> [0,1]
		f  = lambda x: a*x*(1-x)
		g  = lambda x: (x - x_min) / (x_max - x_min)
		gi_expr = sp.solve(g(x) - y, x)[0]
		gi = lambda y_val: gi_expr.subs(y,y_val)
		F_expr = g(f(gi(x))).simplify()
		F = lambda x_val: F_expr.subs(x,x_val)
		
		# Parameters
		# From [Pisarchik]: 3.57 < a < 4
		a_min = sp.Rational(387,100) #357
		a_max = sp.Rational(400,100)
		
		# Map [0,1] -> [a_max,a_min]
		h  = lambda a: (a - a_min) / (a_max - a_min)
		hi_expr = sp.solve(h(a)-x,a)[0]
		hi = lambda x_val: hi_expr.subs(x,x_val)
		
		# Discretization of F
		FD = lambda X: (MAX * F(X/MAX)).subs(a, hi(A/MAX)).simplify()
		#FD(X).diff(C) = MAX*(1849*A**2 + 22102*A*MAX + 16049*MAX**2)/(1849*A**2 + 22102*A*MAX + 56049*MAX**2)

		# Taylor coefficients
		C  = [(FD(X).taylor_term(i,X) / X**i).simplify() for i in range(3)]
		
		# Polynomial descriptor
		num = [c.as_numer_denom()[0] for c in C]
		den = [c.as_numer_denom()[1] for c in C]
		ngcd = sp.gcd(sp.gcd(num[0],num[1]),num[2])
		dgcd = sp.gcd(sp.gcd(den[0],den[1]),den[2])

		# Descriptor
		discrete.__polydesc = {
			'num' : tuple([n/ngcd for n in num]),
			'den' : tuple([d/dgcd for d in den]),
			'N'   : ngcd,
			'D'   : dgcd
		}
		discrete.__polydesc = {
			'num' : tuple([n/ngcd for n in num]),
			'den' : tuple([d/dgcd for d in den]),
			'N'   : ngcd,
			'D'   : dgcd
		}

def psp(n):
    #Tomamos una base aleatoria entre 1 y n-1
    base=random.randint(1,n-1)
    #Comprobamos si gcd(b,n)=1
    #Si es false
    if not sp.gcd(base,n)==1:
        print str(sp.gcd(n,base)) + " es divisor de " + str(n) + "."
        return (base,False)
    #Si es true
    else:
        #comprueba si b(n−1)≡1 mod n
        if not pow(base,n-1,n)==1:
            return (base,False)
        else:
            return (base,True)

def smith_column_step(col, t, var):

    nr = len(col)
    L0 = sp.eye(nr)
    col = col.expand()
    at = col[t]

    for k, ak in enumerate(col):
        if k == t or ak == 0:
            continue
        GCD = sp.gcd(at, ak)
        alpha_t = sp.simplify(at/GCD)
        gamma_k = sp.simplify(ak/GCD)

        sigma, tau = solve_bezout_eq(alpha_t, gamma_k, var)

        L0[t, t] = sigma
        L0[t, k] = tau
        L0[k, t] = -gamma_k
        L0[k, k] = alpha_t

        new_col = sp.expand(L0*col)
        # Linksmultiplikation der Spalte mit L0 liefert eine neue Spalte
        # mit Einträgen beta bei t und 0 bei k
        break

    return new_col, L0

def decodersa(n, e, c):
    
    # using sympy factor int function obtain p,q.
    # takes significant time for long numbers
    primes = factorint(n)
    
    p, _ = dict.popitem(primes)
    q, _ = dict.popitem(primes)

    p1, q1 = p-1, q-1

    # check p-1 and q-1
    if not(gcd(p1 * q1, n) == 1):
        raise Exception('Incorrect p-1 and q-1, their GCD is not 1')
  
    p1q1 = p1 * q1

    # now we solve e*d = 1 (mod (p-1)(q-1))
    # e^-1 (mod (p-1)(q-1))
    e1 = modinv(e, p1q1)

    # check that e1 is e's modular inverse
    if not(pow(e * e1, 1, p1q1) == 1):
        raise Exception('Incorrect e^-1 and e, e*e^-1 not 1')

    # solve for d
    d = e1 % p1q1

    # solve c^d (mod n)
    decoded = pow(long(c), d, n)

    return num2alph(str(decoded))

    def order(self, strategy="relator_based"):
        """
        Returns the order of the finitely presented group ``self``. It uses
        the coset enumeration with identity group as subgroup, i.e ``H=[]``.

        Examples
        ========

        >>> from sympy.combinatorics.free_groups import free_group
        >>> from sympy.combinatorics.fp_groups import FpGroup
        >>> F, x, y = free_group("x, y")
        >>> f = FpGroup(F, [x, y**2])
        >>> f.order(strategy="coset_table_based")
        2

        """
        from sympy import S, gcd
        if self._order != None:
            return self._order
        if self._coset_table != None:
            self._order = len(self._coset_table.table)
        elif len(self.generators) == 1:
            self._order = gcd([r.array_form[0][1] for r in self.relators])
        elif self._is_infinite():
            self._order = S.Infinity
        else:
            gens, C = self._finite_index_subgroup()
            if C:
                ind = len(C.table)
                self._order = ind*self.subgroup(gens, C=C).order()
            else:
                self._order = self.index([])
        return self._order

def regular(S,X,Y,nterms,degree_bound):
    """
    INPUT:

    -- ``S``: a finite set of pairs `\{(\pi_k,F_k)\}_{1 \leq k \leq B}` 
         where `\pi_k` is a finite `\mathbb{K}`-expansion, `F_k
         \in \bar{\mathbb{K}}[X,Y]` with `F_k(0,0) = 0`, `\partial_Y
         F_k(0,0) \neq 0`, and `F_k(X,0) \neq 0`.

    -- ``H``: a positive integer

    OUTPUT:
    
    -- ``(list)``: a set `\{ \pi_k' \}_{1 \leq k \leq B}` of finite
         `\mathbb{K}` expansions such that `\pi_k'` begins with
         `\pi_k` and contains at least `H` `\mathbb{K}`-terms.

    """
    R = []
    for (pi,F) in S:        
        # grow each expansion to the number of desired terms
        Npi = len(pi)

        # if a degree bound is specified, get the degree of the
        # singular part of the Puiseux series. We also want to compute
        # enough terms to distinguish the puiseux series.
        q,mu,m,beta,eta = pi[-1]
        e = reduce( lambda q1,q2: q1*q1, (tau[0] for tau in pi) )
        ydeg = sum( tau[0]*tau[2] for tau in pi )
        need_more_terms = True

        while ((Npi < nterms) and (ydeg < e*degree_bound)) or need_more_terms:
            # if the set of all (0,j), j!=0 is empty, then we've 
            # encountered a finite puiseux expansion
            a = dict(F.terms())
            ms = [j for (j,i) in a.keys() if i==0 and j!=0]
            if ms == []: break
            else:        m = min(ms)

            # if a degree bound is specified, break pi-series
            # construction once we break the bound.
            ydeg += m
            Npi += 1
                
            beta = sympy.together(-a[(m,0)]/a[(0,1)])
            tau = (1,1,m,beta,1)
            pi.append(tau)
            F = _new_polynomial(F,X,Y,tau,m)

            # if the degree in the y-series isn't divisible by the
            # ramification index then we have enough terms to
            # distinguish the puiseux series.
            if sympy.gcd(ydeg,e) == 1:
                need_more_terms = False

        R.append(tuple(pi))

    return tuple(R)

def epsp(n):
    #Primer paso, si n es par devuelve (2,false)
    if(n%2==0):
        print str(n)+" es par"
        return (2,False)
    #Elegimos una base b al azar
    base=random.randint(2,n-1)
    #Comprobamos si gcd(b,n)=1
    #Si no es cierto
    if not sp.gcd(base,n)==1:
        print str(sp.gcd(n,base)) + " es divisor de " + str(n) + "."
        return (base,False)
    #Si es cierto
    else:
        if pow(base,(n-1)/2,n) == sp.jacobi_symbol(base,n)%n :
            return (base,True)
        else:
            return (base,False)

    def apply(self, args, evaluation):
        'CoprimeQ[args__]'

        py_args = [arg.to_python() for arg in args.get_sequence()]
        if not all(isinstance(i, int) or isinstance(i, complex) for i in py_args):
            return Symbol('False')

        if all(sympy.gcd(n,m) == 1 for (n,m) in combinations(py_args, 2)):
            return Symbol('True')
        else:
            return Symbol('False')

    def apply(self, ns, evaluation):
        'GCD[ns___Integer]'

        ns = ns.get_sequence()
        result = 0
        for n in ns:
            value = n.get_int_value()
            if value is None:
                return
            result = sympy.gcd(result, value)
        return Integer(result)

def smith_step(A, t, var):
    # erste Spalte (Index: j), die nicht komplett 0 ist
    # j soll größer als Schrittzähler sein
    nr, nc = A.shape

    row_op_list = []

    cols = st.col_split(A)

    # erste nicht verschwindende Spalte finden
    for j, c in enumerate(cols):
        if j < t:
            continue
        if not c == c*0:
            break
    # Eintrag mit Index t soll ungleich 0 sein, ggf. Zeilen tauschen
    if c[t] == 0:
        i, elt = first_nonzero_element(c)
        ro = row_swap(nr, t, i)
        c = ro*c

        row_op_list.append(ro)

    col = c.expand()
    while True:
        new_col, L0 = smith_column_step(col, t, var)
        if L0 == sp.eye(nr):
            # nothing has changed
            break
        row_op_list.append(L0)
        col = new_col

    # jetzt teilt col[t] alle Einträge in col
    # Probe in der nächsten Schleife

    col.simplify()
    col = col.expand()
    for i,a in enumerate(col):
        if i == t:
            continue
        if not sp.simplify(sp.gcd(a, col[t]) - col[t]) == 0:
            IPS()
            raise ValueError, "col[t] should divide all entries in col"
        quotient = sp.simplify(a/col[t])
        if a == 0:
            continue

        # eliminiere a
        ro = row_op(nr, i, t, -1, quotient)
        row_op_list.append(ro)


    return row_op_list

def keycreation(p, q, e):
    n = p * q

    # check p, q are primes
    if not(isprime(p) and isprime(q)):
        raise Exception('p, q are not primes')

    # check gcd(e, (p-1)(q-1))=1
    if (gcd(e, (p-1) * (q-1)) == 1): 
        return (n, e)
    else:
        raise Exception('Improper e')

def test_sincos_rewrite_sqrt():
    # equivalent to testing rewrite(pow)
    for p in [1, 3, 5, 17, 3*5*17]:
        for t in [1, 8]:
            n = t*p
            for i in xrange(1, (n + 1)//2 + 1):
                if 1 == gcd(i, n):
                    x = i*pi/n
                    s1 = sin(x).rewrite(sqrt)
                    c1 = cos(x).rewrite(sqrt)
                    assert not s1.has(cos, sin), "fails for %d*pi/%d" % (i, n)
                    assert not c1.has(cos, sin), "fails for %d*pi/%d" % (i, n)
                    assert 1e-10 > abs( sin(float(x)) - float(s1) )
                    assert 1e-10 > abs( cos(float(x)) - float(c1) )

def test_sincos_rewrite_sqrt():
    # equivalent to testing rewrite(pow)
    for p in [1, 3, 5, 17]:
        for t in [1, 8]:
            n = t*p
            for i in xrange(1, (n + 1)//2 + 1):
                if 1 == gcd(i, n):
                    x = i*pi/n
                    s1 = sin(x).rewrite(sqrt)
                    c1 = cos(x).rewrite(sqrt)
                    assert not s1.has(cos, sin), "fails for %d*pi/%d" % (i, n)
                    assert not c1.has(cos, sin), "fails for %d*pi/%d" % (i, n)
                    assert 1e-3 > abs(sin(x.evalf(5)) - s1.evalf(2)), "fails for %d*pi/%d" % (i, n)
                    assert 1e-3 > abs(cos(x.evalf(5)) - c1.evalf(2)), "fails for %d*pi/%d" % (i, n)
    assert cos(pi/14).rewrite(sqrt) == sqrt(cos(pi/7)/2 + S.Half)

def test_tancot_rewrite_sqrt():
    # equivalent to testing rewrite(pow)
    for p in [1, 3, 5, 17, 3*5*17]:
        for t in [1, 8]:
            n = t*p
            for i in xrange(1, (n + 1)//2 + 1):
                if 1 == gcd(i, n):
                    x = i*pi/n
                    if  2*i != n and 3*i != 2*n:
                        t1 = tan(x).rewrite(sqrt)
                        assert not t1.has(cot, tan), "fails for %d*pi/%d" % (i, n)
                        assert 1e-10 > abs( tan(float(x)) - float(t1) )
                    if  i != 0 and i != n:
                        c1 = cot(x).rewrite(sqrt)
                        assert not c1.has(cot, tan), "fails for %d*pi/%d" % (i, n)
                        assert 1e-10 > abs( cot(float(x)) - float(c1) )

def desingularize(f,x,y):
    """
    If f is singular, it is desginularized. Outputs new f and Puiseux 
    series expansion data.
    """

    coeffs = _coefficient(f)
    c = coeffs.pop((0,0))
    
    # for each monomial c x**j y**i find the dominant term: that is
    # the one that cancels out the constant term c. This is done
    # by substituting x = T**q, y = eta T**m giving the monomial 
    # c eta**iT**(qj+mi). To balance the equation (kill the constant term) 
    # we need qj+mi=0. q = i, m = -j satisfies this equation.
    #
    # Finally, we need to check that this choice of q,m doesn't introduce 
    # terms with negative exponent in the curve.
    q,m = (1,1)
    for (i,j),aij in coeffs.iteritems():
        # compute q,m
        g = sympy.gcd(i,j)
        q = sympy.Rational(i,g)
        m = -sympy.Rational(j,g)

        # check if the other terms remain positive.
        if all(q*jj+m*ii>=0 for ii,jj in coeffs.keys()):
            break

    if (q,m) == (1,1):
        raise ValueError("Unable to compute singular term.")            
        
    # now compute the values of eta that cancel the constant term c
    Phi = [aij*_Z**sympy.Rational(i,q) for (i,j) in coeffs.keys() 
           if q*j+m*i == 0]
    Phi = sympy.Poly(sum(Phi) + c, _Z)

    return [(q,m,0,Phi)]

def test_H5():
    assert gcd(p1, p2, x) == 1

import socket
import fractions
import sympy

x = sympy.symbols('x')
HOST = "54.64.40.172"
PORT = 5454

s = socket.socket(socket.AF_INET, socket.SOCK_STREAM);
s.connect((HOST, PORT))
while 1:
    n = int(s.recv(20000))
    all = s.recv(20000)
    all = all.split("\n")
    c1 = int(all[1])
    c2 = int(all[2])
    f = x**3-int(c1)
    g = (x+1)**3-int(c2)
    q, r = sympy.div(f, g, x)
    print q
    print r
    print sympy.gcd(g, r)
    print f
    print g
    print x
    #gcd = x - M (find m from expression above) send m
    s.send(str(int(1)) + "\n")
    print s.recv(10000)