Changeset 8354:99efc33c8b7f

Timestamp:

11/03/07 10:56:51 (6 years ago)

Author:

jvoight@…

Branch:

default

Message:

Totally real field enumeration, significant changes and upgrades.

Location:

sage

Files:

• combinat/combinat.py (modified) (1 diff)
• libs/pari/gen.pyx (modified) (1 diff)
• rings/number_field/all.py (modified) (1 diff)
• rings/number_field/totallyreal.py (modified) (7 diffs)
• rings/number_field/totallyreal_data.pyx (modified) (8 diffs)
• rings/number_field/totallyreal_phc.py (modified) (2 diffs)

sage/combinat/combinat.py

@@ -193 +193 @@
 
 from sage.interfaces.all import gap, maxima
-from sage.rings.all import QQ, RR, ZZ
+from sage.rings.all import QQ, ZZ
 from sage.rings.arith import binomial
 from sage.misc.sage_eval import sage_eval

sage/libs/pari/gen.pyx

@@ -4705 +4705 @@
         return n
 
+    def polsturm_full(self):
+        _sig_on
+        n = sturmpart(self.g, NULL, NULL)
+        _sig_off
+        return n
+
     def polsylvestermatrix(self, g):
         t0GEN(g)

sage/rings/number_field/all.py

@@ -12 +12 @@
 
 from totallyreal import enumerate_totallyreal_fields
+from totallyreal_data import hermite_constant
+from totallyreal_dsage import totallyreal_dsage

sage/rings/number_field/totallyreal.py

@@ -40 +40 @@
 from sage.rings.number_field.totallyreal_data import tr_data, int_has_small_square_divisor
 
-def enumerate_totallyreal_fields(n, B, a = [], verbose = False, return_seqs=False):
+def enumerate_totallyreal_fields(n, B, a = [], verbose=0, return_seqs=False, phc=False):
     r"""
     This function enumerates (primitve) totally real fields of 
@@ -48 +48 @@
         $$ a[d]*x^n + ... + a[0]*x^(n-d) $$
     if length(a) = d+1, so in particular always a[d] = 1.  
-    If verbose == True, then print to the screen verbosely; if
-    verbose is a string, then print to the file specified by verbose
-    verbosely.
+    If verbose == 1 (or 2), then print to the screen (really) verbosely; if
+    verbose is a string, then print verbosely to the file specified by verbose.
+    If return_seqs, then return the polynomials as sequences (for easier
+    exporting to a file).
 
     NOTE: 
@@ -158 +159 @@
         # Else, print to screen
     f_out = [0]*n + [1]
-    T.incr(f_out)
+    if verbose == 2:
+        T.incr(f_out,verbose,phc=phc)
+    else:
+        T.incr(f_out,phc=phc)
  
     while f_out[n] <> 0:
@@ -167 +171 @@
         d = nf.poldisc()
         counts[0] += 1
-        if d > 0:
+        if d > 0 and nf.polsturm_full() == n:
             da = int_has_small_square_divisor(Integer(d))
             if d > dB or d <= B*da:
@@ -214 +218 @@
                     print "is not totally real"
 
-        T.incr(f_out)
+        if verbose == 2:
+            T.incr(f_out,verbose=verbose,phc=phc)
+        else:
+            T.incr(f_out,phc=phc)
 
     # In the application of Smyth's theorem above, we exclude finitely
@@ -241 +248 @@
 
     if return_seqs:
-        return [counts,[[s[0],s[1].Vec()] for s in S]]
+        return [counts,[[s[0],s[1].reverse().Vec()] for s in S]]
     else:
         return S
@@ -324 +331 @@
     """
 
-    ppi = float(pi)
-    return ((16./3*ppi)*(g+1))**(2./(3*n))*(2*ppi)**(4./3)
+    ppi = 3.1415926535897931
+    return ((16./3)*(g+1))**(2./(3*n))*(2*ppi)**(4./3)

sage/rings/number_field/totallyreal_data.pyx

@@ -34 +34 @@
 
 import math, numpy
+
+# Other global variables
+ZZx = PolynomialRing(IntegerRing(), 'x')
 
 #***********************************************************************************************
@@ -242 +245 @@
     return asq
 
+cdef eval_seq_as_poly_int(int *f, int n, int x):
+    r"""
+    Evaluates the sequence a, thought of as a polynomial with
+        $$ f[n]*x^n + f[n-1]*x^(n-1) + ... + f[0]. $$
+    """
+    cdef int s, xp
+
+    s = f[n]
+    for i from n > i >= 0:
+        s = s*x+f[i]
+    return s
+
+cdef double eps_abs, phi, sqrt2
+eps_abs = 10.**(-12)
+phi = 0.618033988749895
+sqrt2 = 1.41421356237310
+
 cdef easy_is_irreducible(int *a, int n):
     r"""
-    Very often, polynomials have roots in {-2,-1,1,2}, so we rule
+    Very often, polynomials have roots in {+/-1, +/-2, +/-phi, sqrt2}, so we rule
     these out quickly.  Returns 0 if reducible, 1 if inconclusive.
     """
-    cdef s, sgn
-
-    # Check if x-1 is a factor.
-    s = 0
-    for i from 0 <= i <= n:
-        s += a[i]
-    if s == 0:
+    cdef int s, t, st, sgn
+
+    # Check if a has a root in {1,-1,2,-2}.
+    if eval_seq_as_poly_int(a,n,1) == 0 or eval_seq_as_poly_int(a,n,-1) == 0 or eval_seq_as_poly_int(a,n,2) == 0 or eval_seq_as_poly_int(a,n,-2) == 0:
         return 0
 
-    # Check if x+1 is a factor.
-    s = 0
-    sgn = 1
-    for i from 0 <= i <= n:
-        s += sgn*a[i]
-        sgn *= -1
-    if s == 0:
-        return 0
-
-    # Check if x-2 is a factor.
-    s = 0
-    for i from 0 <= i <= n:
-        s += (2**i)*a[i]
-    if s == 0:
-        return 0
-
-    # Check if x+2 is a factor.
-    s = 0
-    sgn = 1
-    for i from 0 <= i <= n:
-        s += sgn*(2**i)*a[i]
-        sgn *= -1
-    if s == 0:
-        return 0
+    # Check if f has factors x^2-x-1, x^2+x-1, x^2-2, respectively.
+    # Note we only call the ZZx constructor if we're almost certain to reject.
+    if abs(eval_seq_as_poly(a,n,-phi)) < eps_abs:
+        s = 2*a[n]
+        t = 0
+        for i from n > i >= 0:
+            st = (s+t)//2
+            s = 2*t+st+2*a[i]
+            t = st
+        if s == 0 and t == 0:
+            return 0
+    if abs(eval_seq_as_poly(a,n,phi)) < eps_abs:
+        s = 2*a[n]
+        t = 0
+        for i from n > i >= 0:
+            st = (s-t)//2
+            s = 2*t-st+2*a[i]
+            t = st
+        if s == 0 and t == 0:
+            return 0
+    if abs(eval_seq_as_poly(a,n,sqrt2)) < eps_abs:
+        s = a[n]
+        t = 0
+        for i from n > i >= 0:
+            st = s
+            s = 2*t+a[i]
+            t = st
+        if s == 0 and t == 0:
+            return 0
 
     return 1
 
+def easy_is_irreducible_py(f):
+    cdef int a[10]
+
+    for i from 0 <= i < len(f):
+        a[i] = f[i]
+    return easy_is_irreducible(a, len(f)-1)
 
 #***********************************************************************************************
@@ -293 +321 @@
 # (which gives trivial bounds on coefficients) nor too small (which
 # spends needless time computing higher precision on the roots).
+cdef double eps_global
 eps_global = 10.**(-4)
 
-# Other global variables
-ZZx = PolynomialRing(IntegerRing(), 'x')
+from totallyreal_phc import lagrange_bounds_phc
 
 cdef class tr_data:
@@ -363 +391 @@
                 self.a[i] = a[i]
                 self.amax[i] = a[i]
-            self.a[n-1] = -(n/2)
+            self.a[n-1] = -(n//2)
             self.amax[n-1] = 0
             self.k = n-2
@@ -382 +410 @@
 
             # Bounds come from an application of Lagrange multipliers in degrees 2,3.
-            if k == n-2:
-                self.b_lower = 1./n*(-self.a[n-1] - 
-                                 (n-1)*sqrt((1.*self.a[n-1])**2 - 2.*(1-1./n)*self.a[n-2]))
-                self.b_upper = -(2.*self.a[n-1]/n + self.b_lower)
-            else:
+            self.b_lower = 1./n*(-self.a[n-1] - 
+                             (n-1)*sqrt((1.*self.a[n-1])**2 - 2.*(1-1./n)*self.a[n-2]))
+            self.b_upper = -(2.*self.a[n-1]/n + self.b_lower)
+            if k < n-2:
                 bminmax = __lagrange_degree_3(n,a[n-1],a[n-2],a[n-3])
-                self.b_lower = bminmax[0]
-                self.b_upper = bminmax[1]
+                self.b_lower = max(bminmax[0],self.b_lower)
+                self.b_upper = min(bminmax[1],self.b_upper)
 
             # Annoying, but must reverse coefficients for numpy.
@@ -419 +446 @@
         sage_free(self.gnk)
 
-    def incr(self, f_out, verbose=False, haltk=0):
+    def incr(self, f_out, verbose=False, haltk=0, phc=False):
         r"""
         This function 'increments' the totally real data to the next value
@@ -437 +464 @@
         haltk -- integer, the level at which to halt the inductive 
             coefficient bounds
+        phc -- boolean, if PHCPACK is available, use it when k == n-5 to
+            compute an improved Lagrange multiplier bound
 
         OUTPUT:
@@ -578 +607 @@
                         # New bounds from Lagrange multiplier in degree 3.
                         bminmax = __lagrange_degree_3(n,self.a[n-1],self.a[n-2],self.a[n-3])
-                        self.b_lower = bminmax[0]
-                        self.b_upper = bminmax[1]
+                        self.b_lower = max(self.b_lower,bminmax[0])
+                        self.b_upper = min(self.b_upper,bminmax[1])
+                    elif k == n-5 and phc:
+                        # New bounds using phc/Lagrange multiplier in degree 4.
+                        bminmax = lagrange_bounds_phc(n, 4, [self.a[i] for i from 0 <= i <= n])
+                        if len(bminmax) > 0:
+                            self.b_lower = max(self.b_lower,bminmax[0])
+                            self.b_upper = min(self.b_upper,bminmax[1])
+                        else:
+                            maxoutflag = 1
+                            break
 
                     if verbose:

sage/rings/number_field/totallyreal_phc.py

@@ -48 +48 @@
 
 import os, math
+from sage.combinat.combinat import partitions_list
 
-def lagrange_bounds(n, m, a):
+def lagrange_bounds_phc(n, m, a):
     r"""
     This function determines the bounds on the roots in
@@ -118 +119 @@
             output_data += [[P, min(crits), max(crits)]]
 
-    return output_data
+    if len(output_data) > 0:
+        return [min([v[1] for v in output_data]), max([v[2] for v in output_data])]
+    else:
+        return []