# HG changeset patch
# User Chris Wuthrich <christian.wuthrich@gmail.com>
# Date 1300056078 0
# Node ID 9c43a1be2604bd263e2e839961dda17ef6306430
# Parent  361a4ad7d52c69b64ae2e658ffd0820af0d87e93
trac 10926: bug in is_gamma0_equiv

diff -r 361a4ad7d52c -r 9c43a1be2604 sage/modular/cusps.py
--- a/sage/modular/cusps.py	Fri Feb 25 18:56:01 2011 +0000
+++ b/sage/modular/cusps.py	Sun Mar 13 22:41:18 2011 +0000
@@ -1,3 +1,4 @@
+# -*- coding: utf-8 -*-
 r"""
 The set `\mathbb{P}^1(\QQ)` of cusps
 
@@ -606,16 +607,16 @@
            Gamma_0(N))
         
         -  ``transformation`` - bool (default: False), if True,
-           also return upper left entry of a matrix in Gamma_0(N) that sends
-           self to other.
+           also return a matrix in Gamma_0(N) that sends
+           self to other. The matrix is chosen such that the lower left entry is as small as possible in absolute value.
         
         
         OUTPUT:
         
         
-        -  ``bool`` - True if self and other are equivalent
+        -  a boolean - True if self and other are equivalent
         
-        -  ``integer`` - returned only if transformation is
+        -  a matrix - returned only if transformation is
            True
         
         
@@ -625,17 +626,23 @@
             sage: y = Cusp(4,5)
             sage: x.is_gamma0_equiv(y, 2)
             True
-            sage: x.is_gamma0_equiv(y, 2, True)
-            (True, 1)
+            sage: _, ga = x.is_gamma0_equiv(y, 2, True); ga
+            [-1  2]
+            [-2  3]
             sage: x.is_gamma0_equiv(y, 3)
             False
             sage: x.is_gamma0_equiv(y, 3, True)
             (False, None)
+            sage: Cusp(1/2).is_gamma0_equiv(1/3,11,True)
+            (True, [ -3   2]
+            [-11   7])
+
             sage: Cusp(1,0)
             Infinity
             sage: z = Cusp(1,0)
             sage: x.is_gamma0_equiv(z, 3, True)
-            (True, 2)
+            (True, [-1  1]
+            [-3  2])
         
         ALGORITHM: See Proposition 2.2.3 of Cremona's book "Algorithms for
         Modular Elliptic Curves", or Prop 2.27 of Stein's Ph.D. thesis.
@@ -674,6 +681,20 @@
 
         # Now we know the cusps are equivalent.  Use the proof of Prop 2.2.3
         # of Cremona to find a matrix in Gamma_0(N) relating them.
+        from sage.matrix.constructor import matrix
+        
+        if v1 == 0: # the first is oo
+            if v2 == 0: # both are oo
+                return (True, matrix(ZZ,[[1,0],[0,1]]))
+            else:
+                dum, s2, r2 = u2.xgcd(-v2)
+                assert dum.is_one()
+                return (True, matrix(ZZ, [[u2,r2],[v2,s2]]) ) 
+        elif v2 == 0: # the second is oo
+            dum, s1, r1 = u1.xgcd(-v1)
+            assert dum.is_one()
+            return (True, matrix(ZZ, [[s1,-r1],[-v1,u1]]) ) 
+        
         dum, s2, r2 = u2.xgcd(-v2)
         assert dum.is_one()
         dum, s1, r1 = u1.xgcd(-v1)
@@ -683,10 +704,34 @@
         # solve x*v1*v2 + a = 0 (mod N).
         d,x0,y0 = (v1*v2).xgcd(N)          # x0*v1*v2 + y0*N = d = g.
         # so x0*v1*v2 - g = 0 (mod N)
-        x = -x0 * (a/g)
+        x = -x0 * ZZ(a/g)
         # now  x*v1*v2 + a = 0 (mod N)
+        
+        # added in trac 10926
         s1p = s1+x*v1
-        return (True, (u2*s1p-r2*v1)%N)
+        
+        M = N//g
+        C = s1p*v2 - s2*v1
+        if C % (M*v1*v2) == 0 :
+            k = - C//(M*v1*v2)
+        else:    
+            k = - (C/(M*v1*v2)).round()
+            
+        s1pp = s1p + k *M* v1
+        # C += k*M*v1*v2  # is now the smallest in absolute value
+        C = s1pp*v2 - s2*v1
+        A = u2*s1pp - r2*v1
+        
+        r1pp = r1 + (x+k*M)*u1 
+        B = r2 * u1 - r1pp * u2
+        D = s2 * u1 - r1pp * v2
+        
+        ga = matrix(ZZ, [[A,B],[C,D]])
+        assert ga.det() == 1
+        assert C % N == 0
+        assert (A*u1 + B*v1)/(C*u1+D*v1) == u2/v2
+        
+        return (True, ga)
 
     def is_gamma1_equiv(self, other, N):
         """
diff -r 361a4ad7d52c -r 9c43a1be2604 sage/modular/modsym/boundary.py
--- a/sage/modular/modsym/boundary.py	Fri Feb 25 18:56:01 2011 +0000
+++ b/sage/modular/modsym/boundary.py	Sun Mar 13 22:41:18 2011 +0000
@@ -1262,9 +1262,13 @@
             sage: B._is_equiv(Cusp(2/3), Cusp(1/3))
             (True, 5)
             sage: B._is_equiv(Cusp(3/4), Cusp(1/4))
-            (True, 7)
+            (True, -5)
         """
-        return c1.is_gamma0_equiv(c2, self.level(), transformation=True)
+        res = c1.is_gamma0_equiv(c2, self.level(), transformation=True)
+        if res[1] is None:
+            return res
+        else:
+            return (res[0],res[1][0][0])
 
     def _cusp_index(self, cusp):
         """