Ticket #3967: trac-3967-rebase.patch

File trac-3967-rebase.patch, 5.6 kB (added by craigcitro, 4 months ago)

Rebased John's patch for 3.1.2, depends on trac #4155 (both patches)

  • a/sage/rings/number_field/totallyreal.pyx

    old new  
    102102        dB = 33.9508 
    103103    return dB 
    104104 
    105 def enumerate_totallyreal_fields_prim(n, B, a = [], verbose=0, return_seqs=False, \ 
    106                                       phc=False, keep_fields=False, t_2=False): 
     105def enumerate_totallyreal_fields_prim(n, B, a = [], verbose=0, return_seqs=False, 
     106                                      phc=False, keep_fields=False, t_2=False, 
     107                                      just_print=False): 
    107108    r""" 
    108109    This function enumerates primitive totally real fields of degree 
    109110    $n>1$ with discriminant $d \leq B$; optionally one can specify the 
    …… 
    124125     
    125126    If t_2 = T, then keep only polynomials with t_2 norm >= T. 
    126127 
     128    If just_print is not False, instead of creating a sorted list of 
     129    totally real number fileds, we simply write each totally real 
     130    field we find to the file whose filename is given by 
     131    just_print. In this case, we don't return anything. 
     132 
    127133    NOTE: 
    128     This is guaranteed to give all primitive such fields, and 
    129     seems in practice to give many imprimitive ones. 
     134        This is guaranteed to give all primitive such fields, and 
     135        seems in practice to give many imprimitive ones. 
    130136 
    131137    INPUT: 
    132138        n -- integer, the degree 
    …… 
    222228    cdef pari_gen B_pari, d, d_poly, keepB, nf, t2val, ngt2, ng 
    223229    cdef int *f_out 
    224230    cdef int counts[4] 
    225     cdef int i, n_int, j 
     231    cdef int i, n_int, j, skip_jp 
    226232    cdef bint found, use_t2, phc_flag, verb_int, temp_bint 
    227233    cdef Py_ssize_t k0, ind, lenS 
    228234    cdef tr_data T 
    …… 
    295301    else: 
    296302        phc_flag = 0 
    297303 
     304    if just_print: 
     305        skip_jp = 0 
     306        jp_file = open(just_print, "w") 
     307    else: 
     308        skip_jp = 1 
     309 
    298310    # Trivial case 
    299311    if n == 1: 
     312        sage_free(f_out) 
    300313        if return_seqs: 
    301314            return [[0,0,0,0],[[1,[-1,1]]]] 
    302315        else: 
    …… 
    345358 
    346359                        dng = [d, ng] 
    347360 
    348                         # Check if K is contained in the list. 
    349                         found = 0 
    350                         ind = bisect.bisect_left(S, dng) 
    351                         while ind < lenS: 
    352                             if S[ind][0] != d: 
    353                                 break 
    354                             if S[ind][1] == ng: 
    355                                 if verb_int: 
    356                                     print "but is not new" 
    357                                 found = 1 
    358                                 break 
    359                             ind += 1 
    360                              
    361                         ngt2 = <pari_gen>(ng[n_int-1]**2-2*ng[n_int-2]) 
    362                         if not found: 
    363                             temp_bint = ngt2 >= t2val 
    364                             if ((not use_t2) or temp_bint): 
    365                                 if verb_int: 
    366                                     print "and is new!" 
    367                                 S.insert(ind, dng) 
    368                                 lenS += 1 
     361                        if skip_jp: 
     362                            # Check if K is contained in the list. 
     363                            found = 0 
     364                            ind = bisect.bisect_left(S, dng) 
     365                            while ind < lenS: 
     366                                if S[ind][0] != d: 
     367                                    break 
     368                                if S[ind][1] == ng: 
     369                                    if verb_int: 
     370                                        print "but is not new" 
     371                                    found = 1 
     372                                    break 
     373                                ind += 1 
     374 
     375                            ngt2 = <pari_gen>(ng[n_int-1]**2-2*ng[n_int-2]) 
     376                            if not found: 
     377                                temp_bint = ngt2 >= t2val 
     378                                if ((not use_t2) or temp_bint): 
     379                                    if verb_int: 
     380                                        print "and is new!" 
     381                                    S.insert(ind, dng) 
     382                                    lenS += 1 
     383                        else: 
     384                            if ((not use_t2) or ngt2 >= t2val): 
     385                                jp_file.write(str([d, ng.Vecrev()]) + "\n") 
    369386 
    370387                    else: 
    371388                        if verb_int: 
    …… 
    387404            T.incr(f_out,verb_int,0,phc_flag) 
    388405        else: 
    389406            T.incr(f_out,0,0,phc_flag) 
     407 
     408    if not skip_jp: 
     409        if n_int == 2 and B >= 5 and ((not use_t2) or t2val <= 5): 
     410            jp_file.write(str([2,[-1,-1,1]]) + "\n") 
     411            if B >= 8 and B < 32: 
     412                jp_file.write(str([2,[-2,0,1]]) + "\n") 
     413        elif n_int == 3 and B >= 49 and ((not use_t2) or 5 >= t2val): 
     414            jp_file.write(str([3,[1,-2,-1,1]]) + "\n") 
     415        jp_file.close() 
     416        sage_free(f_out) 
     417        return 
     418             
    390419 
    391420    # In the application of Smyth's theorem above (and easy 
    392421    # irreducibility test), we exclude finitely many possibilities 
    …… 
    418447            fsock.close() 
    419448        sys.stdout = saveout 
    420449 
     450    sage_free(f_out) 
    421451    if return_seqs: 
    422452        return [[ counts[i] for i in range(4) ], 
    423453                [[s[0],s[1].reverse().Vec()] for s in S]]