Changeset 8354:99efc33c8b7f

Ignore:
Timestamp:
11/03/07 10:56:51 (6 years ago)
Branch:
default
Message:

Totally real field enumeration, significant changes and upgrades.

Location:
sage
Files:
6 edited

Unmodified
Removed
• sage/combinat/combinat.py

 r6995 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

 r6974 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

 r8353 from totallyreal import enumerate_totallyreal_fields from totallyreal_data import hermite_constant from totallyreal_dsage import totallyreal_dsage
• sage/rings/number_field/totallyreal.py

 r8353 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 \$\$ 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: # 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: 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: 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 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 """ 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

 r8353 import math, numpy # Other global variables ZZx = PolynomialRing(IntegerRing(), 'x') #*********************************************************************************************** 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) #*********************************************************************************************** # (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: 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 # 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. 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 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: # 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

 r8352 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 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 []
Note: See TracChangeset for help on using the changeset viewer.