Changeset 8376:fe09ee4c519e

Timestamp:

11/18/07 13:26:19 (6 years ago)

Author:

jvoight@…

Branch:

default

Message:

Many bug fixes to totallyreal.

Location:

sage/rings/number_field

Files:

• all.py (modified) (1 diff)
• totallyreal.py (modified) (1 diff)
• totallyreal_data.pyx (modified) (1 diff)
• totallyreal_dsage.py (modified) (1 diff)
• totallyreal_rel.py (modified) (19 diffs)

sage/rings/number_field/all.py

@@ -14 +14 @@
 from totallyreal_data import hermite_constant
 from totallyreal_dsage import totallyreal_dsage
+from totallyreal_rel import enumerate_totallyreal_fields_imprim, enumerate_totallyreal_fields_rel
+from hot import my_matrix

sage/rings/number_field/totallyreal.py

@@ -225 +225 @@
         d = nf.poldisc()
         counts[0] += 1
-        if d > 0 and nf.polsturm_full() == n:
+        if d > 0:
             da = int_has_small_square_divisor(Integer(d))
             if d > dB or d <= B*da:

sage/rings/number_field/totallyreal_data.pyx

@@ -630 +630 @@
                         print ""
 
-                    # Check for double roots
                     for i from 0 <= i < n-k-1:
                         if abs(self.beta[k*np1+i] 
-                                 - self.beta[k*np1+(i+1)]) < 2*eps_global:
+                                 - self.beta[k*np1+(i+1)]) < 10*eps_global:
                             # This happens reasonably infrequently, so calling
                             # the Python routine should be sufficiently fast...
                             f = ZZx([self.gnk[(k+1)*n+i] for i in range(n-(k+1)+1)])
-                            df = ZZx([self.gnk[(k+2)*n+i] for i in range(n-(k+2)+1)])
+                            # Could just take self.gnk(k+2)*n+i, but this may not be initialized...
+                            df = ZZx([i*self.gnk[(k+1)*n+i] for i in range(1,n-(k+1)+1)])
                             if gcd(f,df) <> 1:
                                 if verbose:

sage/rings/number_field/totallyreal_dsage.py

@@ -316 +316 @@
         while True:
             time.sleep(self.timeout)
-            self.split_timed_out_jobs()
             self.compile_fields()
+            self.split_all_jobs()
 
     def num_jobs(self):

sage/rings/number_field/totallyreal_rel.py

@@ -22 +22 @@
 #***********************************************************************************************
 
-from sage.rings.arith import binomial, gcd
-from sage.rings.number_field.totallyreal_data import ZZx, lagrange_degree_3, int_has_small_square_divisor
+from sage.rings.arith import binomial, gcd, divisors
+from sage.rings.integer import Integer
+from sage.rings.integer_ring import IntegerRing
+from sage.rings.number_field.totallyreal_data import ZZx, lagrange_degree_3, int_has_small_square_divisor, hermite_constant
 from sage.rings.number_field.number_field import NumberField
-
-from sage.rings.number_field.totallyreal import weed_fields, odlyzko_bound_totallyreal
-
-import math, numpy, bisect
+from sage.rings.polynomial.polynomial_ring import PolynomialRing
+from sage.rings.number_field.totallyreal import weed_fields, odlyzko_bound_totallyreal, enumerate_totallyreal_fields
+from sage.libs.pari.gen import pari
+
+import math, numpy, bisect, sys
+from numpy.linalg import inv
 
 #***********************************************************************************************
@@ -42 +46 @@
     B = Z_F.reduced_basis
 
-    L = matrix([ [v(b) for b in B] for v in Foo])
-    Linv = L.inverse()
+    L = numpy.array([ [v(b) for b in B] for v in Foo])
+    Linv = numpy.linalg.inv(L)
     Vi = [[C[0][0]],[C[0][1]]]
     for i in range(1,d):
         Vi = sum([ [v + [C[i][0]], v + [C[i][1]]] for v in Vi], [])
-    V = Linv*(matrix(Vi).transpose())
+    V = numpy.matrix(Linv)*(numpy.matrix(Vi).transpose())
     j = 0
     while j < 2**d:
         for i in range(d):
             if V[i,j] < V[i,j+1]:
-                V[i,j] = floor(V[i,j])
-                V[i,j+1] = ceil(V[i,j+1])
+                V[i,j] = math.floor(V[i,j])
+                V[i,j+1] = math.ceil(V[i,j+1])
             else:
-                V[i,j] = ceil(V[i,j])
-                V[i,j+1] = floor(V[i,j+1])
+                V[i,j] = math.ceil(V[i,j])
+                V[i,j+1] = math.floor(V[i,j+1])
         j += 2
-    W0 = Linv*matrix([Vi[0]]*d).transpose()
-    W = Linv*matrix([Vi[2**i] for i in range(d)]).transpose()
+    W0 = (Linv*numpy.array([Vi[0]]*d)).transpose()
+    W = (Linv*numpy.array([Vi[2**i] for i in range(d)])).transpose()
     for j in range(d):
         for i in range(d):
             if W[i,j] < W0[i,j]:
-                W[i,j] = floor(W[i,j])
-                W0[i,j] = ceil(W0[i,j])
+                W[i,j] = math.floor(W[i,j])
+                W0[i,j] = math.ceil(W0[i,j])
             else:
-                W[i,j] = ceil(W[i,j])
-                W0[i,j] = floor(W0[i,j])
-    M = matrix(ZZ,d,2**d,[floor(v) for v in V.list()]).columns()
-    M += matrix(ZZ,d,d,[floor(w) for w in W0.list()]).columns()
-    M += matrix(ZZ,d,d,[floor(w) for w in W.list()]).columns()
+                W[i,j] = math.ceil(W[i,j])
+                W0[i,j] = math.floor(W0[i,j])
+    M = [[int(V[i,j]) for i in range(V.shape[0])] for j in range(V.shape[1])]
+    M += [[int(W0[i,j]) for j in range(W0.shape[0])] for i in range(W0.shape[0])]
+    M += [[int(W[i,j]) for j in range(W.shape[1])] for i in range(W.shape[0])]
+
+    from sage.matrix.constructor import matrix
+    M = (matrix(IntegerRing(),len(M),len(M[0]), M).transpose()).columns()
+
     i = 0
     while i < len(M):
@@ -81 +89 @@
         i += 1
 
+    from sage.geometry.lattice_polytope import LatticePolytope
     P = LatticePolytope(matrix(M).transpose())
     S = []
@@ -91 +100 @@
     for p in pts:
         theta = sum([ p.list()[i]*B[i] for i in range(d)])
-        if prod([Foo[i](theta) >= C[i][0] and Foo[i](theta) <= C[i][1] for i in range(d)]):
+        inbounds = True
+        for i in range(d):
+            inbounds = inbounds and Foo[i](theta) >= C[i][0] and Foo[i](theta) <= C[i][1]
+        if inbounds:
             S.append(theta)
     return S
@@ -98 +110 @@
     r"""
     """
+    Z_F = F.maximal_order()
+    B = Z_F.reduced_basis
+    T = Z_F.T
+    P = T.qfminim((C[1]**2)*(1./2), 10**6)[2]
+
+    S = []
+    for p in P:
+        theta = sum([ p.list()[i]*B[i] for i in range(F.degree())])
+        if theta.trace() >= C[0] and theta.trace() <= C[1]:
+            inbounds = True
+            for v in F.real_embeddings():
+                inbounds = inbounds and v(theta) > 0
+            if inbounds:
+                S.append(F(theta))
+    return S
+
+r"""
+def integral_elements_with_trace(F, C):
     if C[0] == C[1]:
         C = [C[0]-0.5, C[1]+0.5]
+
+    print F, C
 
     d = F.degree()
@@ -106 +138 @@
     B = Z_F.reduced_basis
 
-    L = matrix([ [v(b) for b in B] for v in Foo])
-    Linv = L.inverse()
+    L = numpy.array([ [v(b) for b in B] for v in Foo])
+    Linv = numpy.linalg.inv(L)
     V = [ [[0]*i + [C[0]] + [0]*(d-1-i), [0]*i + [C[1]] + [0]*(d-1-i) ] for i in range(d) ]
     V = sum(V, [])
-    V = Linv*(matrix(V).transpose())
+    V = numpy.matrix(Linv)*(numpy.matrix(V).transpose())
     j = 0
     while j < 2*d:
         for i in range(d):
             if V[i,j] < V[i,j+1]:
-                V[i,j] = floor(V[i,j])
-                V[i,j+1] = ceil(V[i,j+1])
+                V[i,j] = math.floor(V[i,j])
+                V[i,j+1] = math.ceil(V[i,j+1])
             else:
-                V[i,j] = ceil(V[i,j])
-                V[i,j+1] = floor(V[i,j+1])
+                V[i,j] = math.ceil(V[i,j])
+                V[i,j+1] = math.floor(V[i,j+1])
         j += 2
-    M = matrix(ZZ,d,2*d,[floor(v) for v in V.list()]).columns()
+    M = [[int(V[i,j]) for j in range(V.shape[1])] for i in range(V.shape[0])]
+
+    M = (matrix(IntegerRing(),d,2*d, M).transpose()).columns()
     i = 0
     while i < len(M):
@@ -139 +173 @@
             S.append(theta)
     return S
+"""
 
 #***********************************************************************************************
@@ -181 +216 @@
         # Initialize constants.
         self.m = m
-        self.d = F.degree()
-        self.n = m*self.d
+        d = F.degree()
+        self.d = d
+        self.n = m*d
         self.B = B
         self.gamma = hermite_constant(self.n-self.d)
@@ -198 +234 @@
 
         Z_Fbasis = self.Z_F.basis()
-        T = pari(matrix([[(x*y).trace() for x in Z_Fbasis] for y in Z_Fbasis])).qflllgram()
+        from sage.matrix.constructor import matrix
+        M = matrix(IntegerRing(),self.d,self.d, [[(x*y).trace() for x in Z_Fbasis] for y in Z_Fbasis])
+        T = pari(M).qflllgram()
         self.Z_F.reduced_basis = [sum([T[i][j].__int__()*Z_Fbasis[j] for j in range(self.d)]) for i in range(self.d)]
+        self.Z_F.T = pari(matrix(IntegerRing(),self.d,self.d, [[(x*y).trace() for x in self.Z_F.reduced_basis] for y in self.Z_F.reduced_basis]))
 
         # Initialize variables.
@@ -215 +254 @@
             for i in range(len(anm1s)):
                 Q = [ [ v(m*x) for v in self.Foo] + [0] for x in self.Z_F.basis()] + [[v(anm1s[i]) for v in self.Foo] + [10**6]]
-                adj = pari(matrix(Q).transpose()).qflll()[self.d]
+                Q = str(numpy.array(Q).transpose())
+                Q = Q.replace(']\n [',';').replace('\n ', '').replace(' ', ',')
+                Q = Q[1:len(Q)-1]
+                adj = pari(Q).qflll()[self.d]
                 anm1s[i] += sum([m*self.Z_F.basis()[i]*adj[i].__int__()//adj[self.d].__int__() for i in range(self.d)])
                 
@@ -221 +263 @@
             self.a[m-1] = self.amaxvals[m-1].pop()
             self.k = m-2
+
+            bl = math.ceil(1.7719*self.n)
+            br = max([1./m*(am1**2).trace() + \
+                            self.gamma*(1./(m**d)*self.B/self.dF)**(1./(self.n-d)) for am1 in anm1s])
+            br = math.floor(br)
+            T2s = integral_elements_with_trace(self.F, [bl,br])
+            self.trace_elts.append([bl,br,T2s])
 
         elif len(a) <= m+1:
@@ -241 +290 @@
             # Bounds come from an application of Lagrange multipliers in degrees 2,3.
             self.b_lower = [-1./m*(v(self.a[m-1]) +  
-                              (m-1.)*sqrt(v(self.a[m-1])**2 - 2.*(1+1./(m-1))*v(self.a[m-2]))) for v in self.Foo]
+                              (m-1.)*math.sqrt(v(self.a[m-1])**2 - 2.*(1+1./(m-1))*v(self.a[m-2]))) for v in self.Foo]
             self.b_upper = [-1./m*(v(self.a[m-1]) -  
-                              (m-1.)*sqrt(v(self.a[m-1])**2 - 2.*(1+1./(m-1))*v(self.a[m-2]))) for v in self.Foo]
+                              (m-1.)*math.sqrt(v(self.a[m-1])**2 - 2.*(1+1./(m-1))*v(self.a[m-2]))) for v in self.Foo]
             if k < m-2:
                 bminmax = [lagrange_degree_3(n,v(self.a[m-1]),v(self.a[m-2]),v(self.a[m-3])) for v in self.Foo]
@@ -332 +381 @@
 
                 if k == m-2:
-                    # We only know the value of a[n-1], the trace.  Need to apply
-                    # basic bounds from Hunter's theorem and Siegel's theorem, with
-                    # improvements due to Smyth, to get bounds on T_2.
+                    # We only know the value of a[n-1], the trace.  
                     bl = max(math.ceil(1.7719*self.n), ((self.a[m-1]**2).trace()*1./m))
-                    br = 1./m*(self.a[m-1].trace())**2 + \
+                    br = 1./m*(self.a[m-1]**2).trace() + \
                            self.gamma*(1./(m**d)*self.B/self.dF)**(1./(self.n-d))
                     br = math.floor(br)
@@ -371 +418 @@
                     for t2 in T2s:
                         am2 = (self.a[m-1]**2-t2)/2
-                        if am2.is_integral() and prod([v((m-1)*self.a[m-1]**2-2*m*am2) > 0 for v in self.Foo]):
-                            am2s.append(am2)
+                        if am2.is_integral():
+                            ispositive = True
+                            for v in self.Foo:
+                                ispositive = ispositive and v((m-1)*self.a[m-1]**2-2*m*am2) > 0
+                            if ispositive:
+                                am2s.append(am2)
 
                     if verbose >= 2:
@@ -389 +440 @@
                     # Initialize the second derivative.
                     self.b_lower = [-1./m*(v(self.a[m-1]) +  
-                                      (m-1.)*sqrt(v(self.a[m-1])**2 - 2.*(1+1./(m-1))*v(self.a[m-2]))) for v in self.Foo]
+                                      (m-1.)*math.sqrt(v(self.a[m-1])**2 - 2.*(1+1./(m-1))*v(self.a[m-2]))) for v in self.Foo]
                     self.b_upper = [-1./m*(v(self.a[m-1]) -  
-                                      (m-1.)*sqrt(v(self.a[m-1])**2 - 2.*(1+1./(m-1))*v(self.a[m-2]))) for v in self.Foo]
+                                      (m-1.)*math.sqrt(v(self.a[m-1])**2 - 2.*(1+1./(m-1))*v(self.a[m-2]))) for v in self.Foo]
                     self.beta[k] = [[self.b_lower[i], -self.Foo[i](self.a[m-1])/m, self.b_upper[i]] for i in range(d)]
                     self.gnk[k] = [0, (m-1)*self.a[m-1], m*(m-1)/2]
@@ -429 +480 @@
                     if k == m-3:
                         self.b_lower = [-1./m*(v(self.a[m-1]) +  
-                                          (m-1.)*sqrt(v(self.a[m-1])**2 - 2.*(1+1./(m-1))*v(self.a[m-2]))) for v in self.Foo]
+                                          (m-1.)*math.sqrt(v(self.a[m-1])**2 - 2.*(1+1./(m-1))*v(self.a[m-2]))) for v in self.Foo]
                         self.b_upper = [-1./m*(v(self.a[m-1]) -  
-                                          (m-1.)*sqrt(v(self.a[m-1])**2 - 2.*(1+1./(m-1))*v(self.a[m-2]))) for v in self.Foo]
+                                          (m-1.)*math.sqrt(v(self.a[m-1])**2 - 2.*(1+1./(m-1))*v(self.a[m-2]))) for v in self.Foo]
                     elif k == m-4:
                         # New bounds from Lagrange multiplier in degree 3.
@@ -768 +819 @@
                 S.insert(ind, [d,ng])
                 Srel.insert(ind, Fx([-1,1,1]))
-        else:
+        elif d == 2:
             for ff in [[1,-7,13,-7,1],[1,-8,14,-7,1]]:
                 f = Fx(ff).factor()[0][0]
@@ -823 +874 @@
     for d in divisors(n):
         if d > 1 and d < n:
-            Sds = enumerate_totallyreal_fields(d, floor((1.*B)**(1.*d/n)), verbose=verbose)
+            Sds = enumerate_totallyreal_fields(d, int(math.floor((1.*B)**(1.*d/n))), verbose=verbose)
             for i in range(len(Sds)):
                 if verbose:
@@ -837 +888 @@
                     S += [[t[0],t[1]] for t in T]
                 j = i+1
-                for E in enumerate_totallyreal_fields(n/d, floor((1.*B)**(1./d)/(1.*Sds[i][0])**(n*1./d**2))):
+                for E in enumerate_totallyreal_fields(n/d, int(math.floor((1.*B)**(1./d)/(1.*Sds[i][0])**(n*1./d**2)))):
                     for EF in F.composite_fields(NumberField(ZZx(E[1]), 'u')):
-                        if EF.degree() == n:
+                        if EF.degree() == n and EF.disc() <= B:
                             S.append([EF.disc(), pari(EF.absolute_polynomial())])
     S.sort()