Ticket #6037: patch-2__QF_local_densities__3.4.1.patch

File patch-2__QF_local_densities__3.4.1.patch, 105.2 KB (added by jonhanke, 12 years ago)
  • sage/quadratic_forms/count_local_2.pyx

    # HG changeset patch
    # User Jonathan Hanke <jonhanke@gmail.com>
    # Date 1242337470 25200
    # Node ID b3d0bbb703abcd636299538541c6706d6fd19b94
    # Parent  af0f19cc490e35ac3b250e6365ad6d504f7d71e2
    Major rewrite/reorganization/documentation of all QuadraticForm local density routines.
    
    diff -r af0f19cc490e -r b3d0bbb703ab sage/quadratic_forms/count_local_2.pyx
    a b  
    22include "../ext/cdefs.pxi"
    33include "../ext/gmp.pxi"
    44
    5 from sage.rings.arith import valuation
     5from sage.rings.arith import valuation, kronecker_symbol, is_prime
    66from sage.rings.integer_mod import IntegerMod, Mod
    77from sage.rings.integer_mod_ring import IntegerModRing
    88
    99from sage.rings.integer_ring import ZZ
    1010
    1111from sage.rings.integer_mod cimport IntegerMod_gmp
     12from sage.sets.set import Set
     13
     14
     15
     16
     17def extract_sublist_indices(Biglist, Smalllist):
     18    """
     19    Returns the indices of Biglist which index the entries of
     20    Smalllist appearing in Biglist.  (Note that Smalllist may not be a
     21    sublist of Biglist.)
     22
     23    NOTE 1: This is an internal routine which deals with reindexing
     24    lists, and is not exported to the QuadraticForm namespace!
     25
     26    NOTE 2: This should really by applied only when BigList has no
     27    repeated entries.
     28
     29    TO DO: *** Please revisit this routine, and eliminate it! ***
     30
     31
     32
     33    INPUT:
     34        Biglist, Smalllist -- two lists of a common type, where
     35                              Biglist has no repeated entries.
     36
     37    OUTPUT:
     38        a list of integers >= 0
     39
     40    EXAMPLES:
     41        sage: from sage.quadratic_forms.quadratic_form__local_density_congruence import extract_sublist_indices
     42
     43        sage: biglist = [1,3,5,7,8,2,4]
     44        sage: sublist = [5,3,2]
     45        sage: sublist == [biglist[i]  for i in extract_sublist_indices(biglist, sublist)]  ## Ok whenever Smalllist is a sublist of Biglist
     46        True
     47
     48        sage: extract_sublist_indices([1,2,3,6,9,11], [1,3,2,9])
     49        [0, 2, 1, 4]
     50
     51        sage: extract_sublist_indices([1,2,3,6,9,11], [1,3,10,2,9,0])
     52        [0, 2, 1, 4]
     53
     54        sage: extract_sublist_indices([1,3,5,3,8], [1,5])
     55        Traceback (most recent call last):
     56        ...
     57        TypeError: Biglist must not have repeated entries!
     58    """ 
     59    ## Check that Biglist has no repeated entries
     60    Big_set = Set(Biglist)
     61    if len(Set(Biglist)) != len(Biglist):
     62        raise TypeError, "Biglist must not have repeated entries!"
     63
     64    ## Extract the indices of Biglist needed to make Sublist
     65    index_list = []
     66    for x in Smalllist:
     67        try:
     68            index_list.append(Biglist.index(x))
     69        except ValueError:                                   ## This happens when an entry of Smalllist is not contained in Biglist
     70            None
     71
     72    ## Return the list if indices
     73    return index_list
     74
     75
     76
     77
     78def count_modp__by_gauss_sum(n, p, m, Qdet):
     79    """
     80    Returns the number of solutions of Q(x) = m over the finite field
     81    Z/pZ, where p is a prime number > 2 and Q is a non-degenerate
     82    quadratic form of dimension n >= 1 and has Gram determinant Qdet.
     83
     84    REFERENCE:
     85               These are defined in Table 1 on p363 of Hanke's "Local
     86        Densities..." paper.
     87
     88    INPUT:
     89        n -- an integer >= 1
     90        p -- a prime number > 2
     91        m -- an integer
     92        Qdet -- a integer which is non-zero mod p
     93
     94    OUTPUT:
     95        an integer >= 0
     96
     97    EXAMPLES:
     98        sage: from sage.quadratic_forms.count_local_2 import count_modp__by_gauss_sum
     99
     100        sage: count_modp__by_gauss_sum(3, 3, 0, 1)    ## for Q = x^2 + y^2 + z^2  => Gram Det = 1 (mod 3)
     101        9
     102        sage: count_modp__by_gauss_sum(3, 3, 1, 1)    ## for Q = x^2 + y^2 + z^2  => Gram Det = 1 (mod 3)
     103        6
     104        sage: count_modp__by_gauss_sum(3, 3, 2, 1)    ## for Q = x^2 + y^2 + z^2  => Gram Det = 1 (mod 3)
     105        12
     106
     107        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1])
     108        sage: [Q.count_congruence_solutions(3, 1, m, None, None) == count_modp__by_gauss_sum(3, 3, m, 1)  for m in range(3)]
     109        [True, True, True]
     110
     111
     112        sage: count_modp__by_gauss_sum(3, 3, 0, 2)    ## for Q = x^2 + y^2 + 2*z^2  => Gram Det = 2 (mod 3)
     113        9
     114        sage: count_modp__by_gauss_sum(3, 3, 1, 2)    ## for Q = x^2 + y^2 + 2*z^2  => Gram Det = 2 (mod 3)
     115        12
     116        sage: count_modp__by_gauss_sum(3, 3, 2, 2)    ## for Q = x^2 + y^2 + 2*z^2  => Gram Det = 2 (mod 3)
     117        6
     118
     119        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,2])
     120        sage: [Q.count_congruence_solutions(3, 1, m, None, None) == count_modp__by_gauss_sum(3, 3, m, 2)  for m in range(3)]
     121        [True, True, True]
     122
     123       
     124    """
     125    ## Check that Qdet is non-degenerate
     126    if Qdet % p == 0:
     127        raise RuntimeError, "Qdet must be non-zero."
     128
     129    ## Check that p is prime > 0
     130    if not is_prime(p) or p == 2:
     131        raise RuntimeError, "p must be a prime number > 2."
     132
     133    ## Check that n >= 1
     134    if n < 1:
     135        raise RuntimeError, "the dimension n must be >= 1."
     136
     137    ## Compute the Gauss sum
     138    neg1 = -1   
     139    if (m % p == 0):
     140        if (n % 2 != 0):
     141            count = (p**(n-1))
     142        else:   
     143            count = (p**(n-1)) + (p-1) * (p**((n-2)/2)) * kronecker_symbol(((neg1**(n/2)) * Qdet) % p, p)
     144    else:
     145        if (n % 2 != 0):
     146            count = (p**(n-1)) + (p**((n-1)/2)) * kronecker_symbol(((neg1**((n-1)/2)) * Qdet * m) % p, p)
     147        else:
     148            count = (p**(n-1)) - (p**((n-2)/2)) * kronecker_symbol(((neg1**(n/2)) * Qdet) % p, p)
     149
     150    ## Return the result
     151    return count
     152
    12153
    13154
    14155
     
    16157
    17158cdef CountAllLocalTypesNaive_cdef(Q, p, k, m, zvec, nzvec):
    18159    """
    19     ///////////////////////////////////////////////////////////////////
    20     /// Naively counts the number of solutions of Q(x) = m mod p^k   //
    21     /// of type solntype, satisfying the mod p congruence conditions //
    22     /// at the indices of the vectors "zero" and "nonzero"           //
    23     ///////////////////////////////////////////////////////////////////
    24 
    25     valarray <mpz_class> Matrix_mpz::CountAllLocalTypesNaive(const mpz_class & p, unsigned long k, const mpz_class & m,
    26                                              const valarray<size_t> & zero, const valarray<size_t> & nonzero) const
     160    This cython routine is documented in its Python wrapper method
     161    QuadraticForm.count_congruence_solutions_by_type().
    27162    """
    28 
    29     ## DIAGNOSTIC
    30     #print "   --> CountAllLocalTypesNaive is using the form Q \n" + str(Q)
    31     #print "       p = " + str(p) + "  and   m = " + str(m)
    32 
    33     #cdef mpz_t* v
    34163    cdef long n, i
    35164    cdef long a, b    ## Used to quickly evaluate Q(v)
    36165    cdef long ptr     ## Used to increment the vector
    37166    cdef long solntype    ## Used to store the kind of solution we find
    38167
    39168
    40 
    41 
     169    ## Some shortcuts and definitions
    42170    n = Q.dim()
    43171    R = p ** k
     172    Q1 = Q.base_change_to(IntegerModRing(R))
     173
    44174
    45175    ## Cython Variables
    46176    cdef IntegerMod_gmp zero, one
     
    48178    one = IntegerMod_gmp(IntegerModRing(R), 1)
    49179
    50180
    51     Q1 = Q.base_change_to(IntegerModRing(R))
    52 
    53 
    54 
    55     ###########################################
    56     #v = <mpz_t*> sage_malloc(sizeof(mpz_t)*n)
    57     #for i from 0 <= i < n:
    58     #    mpz_init(v[i])       
    59     ###########################################
    60181
    61182    ## Initialize the counting vector
    62183    count_vector = [0  for i in range(6)]
    63184
    64185    ## Initialize v = (0, ... , 0)
    65186    v = [Mod(0, R)  for i in range(n)]
    66     #v = []
    67     #for i in 1 <= i < n:
    68     #    v.append(zero)
    69187
    70188
    71189    ## Some declarations to speed up the loop
     
    100218                count_vector[solntype] += 1
    101219
    102220
    103 
    104221    ## Generate the Bad-type and Total counts
    105222    count_vector[3] = count_vector[4] + count_vector[5]
    106223    count_vector[0] = count_vector[1] + count_vector[2] + count_vector[3]
    107224
     225    ## Return the solution counts
     226    return count_vector   
    108227
    109     ## DIAGNOSTIC
    110     #cout << "R = " << R << "\n";
    111     #cout << "n = " << n << "\n";
    112     #cout << "R_n = " << R_n << "\n";
    113     #
    114     #for(i=1; i<=25; i++) {
    115     #    cout << "v = " << v << "\n";
    116     #    Increment(v,R);
    117     #}
    118     #
    119     #cout << "Q1 = " << Q1 << "\n";
    120     #cout << "v = " << v << "\n";
    121     #cout << "Q1 * v = " << Q1 * v<< "\n";
    122 
    123     ###########################
    124     #for i from 0 <= i < n:
    125     #    mpz_clear(v[i])
    126     #sage_free(v)
    127     ###########################
    128 
    129     return count_vector   
    130228
    131229
    132230
    133231def CountAllLocalTypesNaive(Q, p, k, m, zvec, nzvec):
    134232    """
    135     Python wrapper function for this (now) cython function.
     233    This is an internal routine, which is called by
     234    QuadraticForm.count_congruence_solutions_by_type().  See that
     235    documentation for more details.
     236
     237    INPUT:
     238        Q -- quadratic form over ZZ
     239        p -- prime number > 0
     240        k -- an integer > 0
     241        m -- an integer (depending only on mod p^k)
     242        zvec, nzvec -- a list of integers in range(Q.dim()), or None
     243
     244    OUTPUT:
     245        a list of six integers >= 0 representing the solution types:
     246            [All, Good, Zero, Bad, BadI, BadII]
     247
     248
     249    EXAMPLES:
     250        sage: from sage.quadratic_forms.count_local_2 import CountAllLocalTypesNaive
     251        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
     252        sage: CountAllLocalTypesNaive(Q, 3, 1, 1, None, None)
     253        [6, 6, 0, 0, 0, 0]
     254        sage: CountAllLocalTypesNaive(Q, 3, 1, 2, None, None)
     255        [6, 6, 0, 0, 0, 0]
     256        sage: CountAllLocalTypesNaive(Q, 3, 1, 0, None, None)
     257        [15, 12, 1, 2, 0, 2]
     258
    136259    """
    137260    return CountAllLocalTypesNaive_cdef(Q, p, k, m, zvec, nzvec)
    138261
     
    144267
    145268cdef local_solution_type_cdef(Q, p, w, zvec, nzvec):
    146269    """
    147     ////////////////////////////////////////////////////////////////////////////////////
    148     /// Private routine to check if a given solution vector w (of Q(w) = m mod p^k)   //
    149     /// is of a certain local type and satisfies certain congruence conditions mod p. //
    150     ///   (Personal Note: For p=2, we should still use p=2 and not p=8.)              //
    151     ////////////////////////////////////////////////////////////////////////////////////
     270    Internal routine to check if a given solution vector w (of Q(w) =
     271    m mod p^k) is of a certain local type and satisfies certain
     272    congruence conditions mod p.
    152273
    153     size_t Matrix_mpz::local_solution_type(const mpz_class & p, const valarray<mpz_class> & w,
    154                const valarray<size_t> & zero, const valarray<size_t> & nonzero) const
     274    NOTE: No internal checking is done to test if p is a prime >=2, or
     275    that Q has the same size as w.
    155276
    156277    """
    157 
    158     ## Note: Here p is assumed to be a prime >= 2, though the routine still works if not...
    159 
    160     ## ToDo?: Add a check that Q is square and has the same size as w.
    161 
    162278    cdef long i
    163279    cdef long n 
    164280
    165281    n = Q.dim()
    166282
    167     zero_flag = False        ## Tests the zero mod p congruence conditions
    168     nonzero_flag = False     ## Tests the nonzero congruence conditions
    169 
    170283
    171284    ## Check if the solution satisfies the zvec "zero" congruence conditions
    172285    ## (either zvec is empty or its components index the zero vector mod p)
    173     if (len(zvec) == 0):
     286    if (zvec == None) or (len(zvec) == 0):
    174287        zero_flag = True
    175288    else:
     289        zero_flag = False
    176290        i = 0
    177         while ( (i < len(zvec)) and ((w[zvec[i]] % p) == 0) ):
     291        while ( (i < len(zvec)) and ((w[zvec[i]] % p) == 0) ):  ## Increment so long as our entry is zero (mod p)
    178292            i += 1
    179         if (i == len(zvec)):
     293        if (i == len(zvec)):      ## If we make it through all entries then the solution is zero (mod p)
    180294            zero_flag = True
    181295
    182296
     
    184298    #print "IsLocalSolutionType: Finished the Zero congruence condition test \n"
    185299
    186300    if (zero_flag == False):
    187         return 0
     301        return <long> 0
    188302
    189303    ## DIAGNOSTIC
    190304    #print "IsLocalSolutionType: Passed the Zero congruence condition test \n"
    191305
    192306
    193307    ## Check if the solution satisfies the nzvec "nonzero" congruence conditions
    194     ## (either nzvec is empty or its components index a non-zero vector mod p)
    195     if (len(nzvec) == 0):
     308    ## (nzvec is non-empty and its components index a non-zero vector mod p)
     309    if (nzvec == None):
    196310        nonzero_flag = True
     311    elif (len(nzvec) == 0):
     312        nonzero_flag = False           ## Trivially no solutions in this case!
    197313    else:
     314        nonzero_flag = False 
    198315        i = 0
    199316        while ((nonzero_flag == False) and (i < len(nzvec))):
    200317            if ((w[nzvec[i]] % p) != 0):
    201                 nonzero_flag = True
     318                nonzero_flag = True           ## The non-zero condition is satisfied when we find one non-zero entry
    202319            i += 1
    203320 
    204321    if (nonzero_flag == False):
    205322        return <long> 0
    206323 
    207324 
    208     ## Check if the solution has the appropriate (local) type
    209  
     325    ## Check if the solution has the appropriate (local) type:
     326    ## -------------------------------------------------------
    210327 
    211328    ## 1: Check Good-type
    212329    for i from 0 <= i < n:
     
    253370        return <long> 4
    254371
    255372
    256     ##    cout << " Bad II Soln :  " << w << "  wS1_nonzero_flag = " << wS1_nonzero_flag << endl;
    257 
    258373    ## 5: Check Bad-type II
    259374    if (wS1_nonzero_flag == False):
    260375        #print " Bad II Soln :  " + str(w)
  • sage/quadratic_forms/quadratic_form.py

    diff -r af0f19cc490e -r b3d0bbb703ab sage/quadratic_forms/quadratic_form.py
    a b  
    9090    ## Import specialized methods:
    9191    ## ---------------------------
    9292
    93     from sage.quadratic_forms.quadratic_form__local_normal_form \
    94     import find_entry_with_minimal_scale_at_prime, local_normal_form, jordan_blocks_by_scale_and_unimodular, \
    95         jordan_blocks_in_unimodular_list_by_scale_power              ## Routines to compute the p-adic local normal form.
     93    ## Routines to compute the p-adic local normal form
     94    from sage.quadratic_forms.quadratic_form__local_normal_form import \
     95            find_entry_with_minimal_scale_at_prime, \
     96            local_normal_form, \
     97            jordan_blocks_by_scale_and_unimodular, \
     98            jordan_blocks_in_unimodular_list_by_scale_power
    9699
    97     from sage.quadratic_forms.quadratic_form__variable_substitutions \
    98     import swap_variables, multiply_variable, divide_variable, scale_by_factor, \
    99         extract_variables, elementary_substitution, add_symmetric    ## Routines to perform elementary variable substitutions.
     100    ## Routines to perform elementary variable substitutions
     101    from sage.quadratic_forms.quadratic_form__variable_substitutions import \
     102            swap_variables, \
     103            multiply_variable, \
     104            divide_variable, \
     105            scale_by_factor, \
     106            extract_variables, \
     107            elementary_substitution, \
     108            add_symmetric   
    100109
    101     from sage.quadratic_forms.quadratic_form__local_field_invariants \
    102     import rational_diagonal_form, signature_vector, signature, local_diagonal, hasse_invariant, hasse_invariant__OMeara, \
    103         is_hyperbolic, is_anisotropic, is_isotropic, anisotropic_primes, compute_definiteness, \
    104         is_positive_definite, is_negative_definite, is_indefinite, is_definite    ## Routines to compute p-adic field invariants
     110    ## Routines to compute p-adic field invariants
     111    from sage.quadratic_forms.quadratic_form__local_field_invariants import \
     112            rational_diagonal_form, \
     113            signature_vector, \
     114            signature, \
     115            local_diagonal, \
     116            hasse_invariant, \
     117            hasse_invariant__OMeara, \
     118            is_hyperbolic, \
     119            is_anisotropic, \
     120            is_isotropic, \
     121            anisotropic_primes, \
     122            compute_definiteness, \
     123            is_positive_definite, \
     124            is_negative_definite, \
     125            is_indefinite, \
     126            is_definite
    105127
    106128    ## Routines to compute local densities by the reduction procedure
    107129    from sage.quadratic_forms.quadratic_form__local_density_congruence import \
    108             reindex_vector_from_extraction, \
    109             count_modp__by_gauss_sum, \
     130            count_modp_solutions__by_Gauss_sum, \
    110131            local_good_density_congruence_odd, \
    111132            local_good_density_congruence_even, \
    112             local_good_density_congruence_even__experimental, \
    113133            local_good_density_congruence, \
    114134            local_zero_density_congruence, \
    115135            local_badI_density_congruence, \
     
    120140
    121141    ## Routines to compute local densities by counting solutions of various types
    122142    from sage.quadratic_forms.quadratic_form__count_local_2 import \
    123             VecIncrement__deprecated, \
    124             local_solution_type__deprecated, \
    125             CountAllLocalTypesNaive__deprecated, \
    126             count_local_type, \
    127             count_local_good_type, \
    128             count_local_zero_type, \
    129             count_local_bad_type, \
    130             count_local_bad_typeI, \
    131             count_local_bad_typeII
     143            count_congruence_solutions_as_vector, \
     144            count_congruence_solutions, \
     145            count_congruence_solutions__good_type, \
     146            count_congruence_solutions__zero_type, \
     147            count_congruence_solutions__bad_type, \
     148            count_congruence_solutions__bad_type_I, \
     149            count_congruence_solutions__bad_type_II
    132150       
    133151    ## Routines to be called by the user to compute local densities
    134152    from sage.quadratic_forms.quadratic_form__local_density_interfaces import \
    135             local_good_density, \
    136             local_zero_density, \
    137             local_bad_density, \
    138             local_badI_density, \
    139             local_badII_density, \
    140153            local_density, \
    141154            local_primitive_density
    142155
     
    166179            is_zero_nonsingular, \
    167180            is_zero_singular
    168181
    169     from sage.quadratic_forms.quadratic_form__theta \
    170     import theta_series, theta_by_pari, theta_by_cholesky         ## Routines to compute the theta function
     182    ## Routines to compute the theta function
     183    from sage.quadratic_forms.quadratic_form__theta import \
     184            theta_series, \
     185            theta_by_pari, \
     186            theta_by_cholesky
    171187
    172188    ## Routines to compute the product of all local densities
    173189    from sage.quadratic_forms.quadratic_form__siegel_product import \
    174190            siegel_product
    175191
     192    ## Routines to compute p-neighbors
    176193    from sage.quadratic_forms.quadratic_form__neighbors import \
    177     find_primitive_p_divisible_vector__random, \
    178     find_primitive_p_divisible_vector__all, \
    179     find_primitive_p_divisible_vector__next, \
    180     find_p_neighbor_from_vec                       ## Routines to compute p-neighbors
     194            find_primitive_p_divisible_vector__random, \
     195            find_primitive_p_divisible_vector__all, \
     196            find_primitive_p_divisible_vector__next, \
     197            find_p_neighbor_from_vec
    181198
    182199    ## Routines to reduce a given quadratic form
    183200    from sage.quadratic_forms.quadratic_form__reduction_theory import \
     
    187204            minkowski_reduction, \
    188205            minkowski_reduction_for_4vars__SP
    189206 
    190     from sage.quadratic_forms.quadratic_form__genus \
    191     import global_genus_symbol, local_genus_symbol, CS_genus_symbol_list   ## Wrappers for Conway-Sloane genus routines (in ./genera/)
     207    ## Wrappers for Conway-Sloane genus routines (in ./genera/)
     208    from sage.quadratic_forms.quadratic_form__genus import \
     209            global_genus_symbol, \
     210            local_genus_symbol, \
     211            CS_genus_symbol_list
     212
    192213
    193214    ## Routines to compute local masses for ZZ.
    194215    from sage.quadratic_forms.quadratic_form__mass import \
     
    227248            is_locally_represented_number_at_place, \
    228249            is_locally_represented_number
    229250
     251    ## Routines to make a split local covering of the given quadratic form.
    230252    from sage.quadratic_forms.quadratic_form__split_local_covering import \
    231     cholesky_decomposition, vectors_by_length, \
    232     complementary_subform_to_vector, split_local_cover     ## Routines to make a split local covering of the given quadratic form.
     253            cholesky_decomposition, \
     254            vectors_by_length, \
     255            complementary_subform_to_vector, \
     256            split_local_cover
    233257
    234258    ## Routines to make automorphisms of the given quadratic form.
    235259    from sage.quadratic_forms.quadratic_form__automorphisms import \
  • sage/quadratic_forms/quadratic_form__count_local_2.py

    diff -r af0f19cc490e -r b3d0bbb703ab sage/quadratic_forms/quadratic_form__count_local_2.py
    a b  
     1##################################################################
     2## Methods for counting/computing the number of representations ##
     3## of a number by a quadratic form in Z/(p^k)Z of various types ##
     4##################################################################
    15
    26
    37from sage.quadratic_forms.count_local_2 import CountAllLocalTypesNaive
    48
    5 from sage.rings.arith import valuation
    69
    710
     11def count_congruence_solutions_as_vector(self, p, k, m, zvec, nzvec):
     12    """
     13    Gives the number of integer solution vectors x satisfying the
     14    congruence Q(x) = m (mod p^k) of each solution type (i.e. All,
     15    Good, Zero, Bad, BadI, BadII) which satisfy the additional
     16    congruence conditions of having certain coefficients = 0 (mod p)
     17    and certain collections of coefficients not congruent to the zero
     18    vector (mod p).
    819
     20    A solution vector x satisfies the additional congruence conditions
     21    specified by zvec and nzvec (and therefore is counted) iff both of
     22    the following conditions hold:
    923
     24        1) x[i] == 0 (mod p) for all i in zvec
    1025
    11 def VecIncrement__deprecated(self, v, R):
    12     """
    13     Performs an in-place imcrement operation on the vector v, whose
    14     entries are satisfy 0 <= v[i] <= R-1.  No values are returned.
     26        2) x[i] != 0 (mod p) for all i in nzvec
    1527
    16 void Matrix_mpz::Increment(valarray<mpz_class> & v, const mpz_class & R) const
    17 00008 {
    18 00009     size_t i;
    19 00010     i = v.size();
    20 00011
    21 00012     // Do the carry operations
    22 00013     while ((i > 0) && (v[i-1] == R-1))   // Assumes that all components satisfy 0 <= v[i] <= R-1
    23 00014       {
    24 00015         v[i-1] = 0;
    25 00016         i--;
    26 00017       }
    27 00018
    28 00019     // Only increment if we're not already at the zero vector =)
    29 00020     if (i > 0)
    30 00021       v[i-1]++;
    31 00022 }
    32     """
    33     i = len(v)
    3428
    35     ## Do the carry operations
    36     while ((i > 0) and (v[i-1] == R-1)):   ## Assumes that all components satisfy 0 <= v[i] <= R-1
    37         v[i-1] = 0
    38         i += -1
     29    REFERENCES: See Hanke's (????) paper "Local Densities and explicit
     30        bounds...", p??? for the definitions of the solution types and
     31        congruence conditions.
    3932
    40     ## Only increment if we're not already at the zero vector =)
    41     if (i > 0):
    42         v[i-1] += 1
     33    INPUT:
     34        p -- prime number > 0
     35        k -- an integer > 0
     36        m -- an integer (depending only on mod p^k)
     37        zvec, nzvec -- a list of integers in range(self.dim()), or None
    4338
     39    OUTPUT:
     40        a list of six integers >= 0 representing the solution types:
     41            [All, Good, Zero, Bad, BadI, BadII]
    4442
    4543
    46 
    47 def local_solution_type__deprecated(self, p, w, zvec, nzvec):
    48     """
    49     ////////////////////////////////////////////////////////////////////////////////////
    50     /// Private routine to check if a given solution vector w (of Q(w) = m mod p^k)   //
    51     /// is of a certain local type and satisfies certain congruence conditions mod p. //
    52     ///   (Personal Note: For p=2, we should still use p=2 and not p=8.)              //
    53     ////////////////////////////////////////////////////////////////////////////////////
    54 
    55     size_t Matrix_mpz::local_solution_type(const mpz_class & p, const valarray<mpz_class> & w,
    56                const valarray<size_t> & zero, const valarray<size_t> & nonzero) const
     44    EXAMPLES:
     45        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
     46        sage: Q.count_congruence_solutions_as_vector(3, 1, 1, [], [])
     47        [0, 0, 0, 0, 0, 0]
     48        sage: Q.count_congruence_solutions_as_vector(3, 1, 1, None, [])
     49        [0, 0, 0, 0, 0, 0]
     50        sage: Q.count_congruence_solutions_as_vector(3, 1, 1, [], None)
     51        [6, 6, 0, 0, 0, 0]
     52        sage: Q.count_congruence_solutions_as_vector(3, 1, 1, None, None)
     53        [6, 6, 0, 0, 0, 0]
     54        sage: Q.count_congruence_solutions_as_vector(3, 1, 2, None, None)
     55        [6, 6, 0, 0, 0, 0]
     56        sage: Q.count_congruence_solutions_as_vector(3, 1, 0, None, None)
     57        [15, 12, 1, 2, 0, 2]
    5758
    5859    """
     60    return CountAllLocalTypesNaive(self, p, k, m, zvec, nzvec)
    5961
    60     ## Note: Here p is assumed to be a prime >= 2, though the routine still works if not...
    61 
    62     ## ToDo?: Add a check that Q is square and has the same size as w.
    63  
    64 
    65     zero_flag = False        ## Tests the zero mod p congruence conditions
    66     nonzero_flag = False     ## Tests the nonzero congruence conditions
    67 
    68 
    69     ## Check if the solution satisfies the zvec "zero" congruence conditions
    70     ## (either zvec is empty or its components index the zero vector mod p)
    71     if (len(zvec) == 0):
    72         zero_flag = True
    73     else:
    74         i = 0
    75         while ( (i < len(zvec)) and ((w[zvec[i]] % p) == 0) ):
    76             i += 1
    77         if (i == len(zvec)):
    78             zero_flag = True
    79 
    80 
    81     ## DIAGNOSTIC
    82     #print "IsLocalSolutionType: Finished the Zero congruence condition test \n"
    83 
    84     if (zero_flag == False):
    85         return 0
    86 
    87     ## DIAGNOSTIC
    88     #print "IsLocalSolutionType: Passed the Zero congruence condition test \n"
    89 
    90 
    91     ## Check if the solution satisfies the nzvec "nonzero" congruence conditions
    92     ## (either nzvec is empty or its components index a non-zero vector mod p)
    93     if (len(nzvec) == 0):
    94         nonzero_flag = True
    95     else:
    96         i = 0
    97         while ((nonzero_flag == False) and (i < len(nzvec))):
    98             if ((w[nzvec[i]] % p) != 0):
    99                 nonzero_flag = True
    100             i += 1
    101  
    102     if (nonzero_flag == False):
    103         return 0
    104  
    105  
    106     ## Check if the solution has the appropriate (local) type
    107  
    108  
    109     ## 1: Check Good-type
    110     for i in range(len(w)):
    111         if (((w[i] % p) != 0)  and ((self[i,i] % p) != 0)):
    112             return 1
    113     if (p == 2):
    114         for i in range(len(w) - 1):
    115             if (((self[i,i+1] % p) != 0) and (((w[i] % p) != 0) or ((w[i+1] % p) != 0))):
    116                 return 1 
    117 
    118 
    119     ## 2: Check Zero-type
    120     Zero_flag = True
    121     for i in range(len(w)):
    122         if ((w[i] % p) != 0):
    123             Zero_flag = False
    124     if (Zero_flag == True):
    125         return 2
    126 
    127 
    128     ## Check if wS1 is zero or not
    129     wS1_nonzero_flag = False
    130     for i in range(self.dim()):     
    131  
    132         ## Compute the valuation of each index, allowing for off-diagonal terms
    133         if (self[i,i] == 0):
    134             if (i == 0):
    135                 val = valuation(self[i,i+1], p)    ## Look at the term to the right
    136             elif (i == self.dim() - 1):
    137                 val = valuation(self[i-1,i], p)    ## Look at the term above
    138             else:
    139                 val = valuation(self[i,i+1] + self[i-1,i], p)    ## Finds the valuation of the off-diagonal term since only one isn't zero
    140         else:
    141             val = valuation(self[i,i], p) 
    142 
    143         ## Test each index
    144         if ((val == 1) and ((w[i] % p) != 0)):
    145             wS1_nonzero_flag = True
    146 
    147      
    148     ## 4: Check Bad-type I
    149     if (wS1_nonzero_flag == True):
    150         #print " Bad I Soln :  " + str(w)
    151         return 4
    152 
    153 
    154     ##    cout << " Bad II Soln :  " << w << "  wS1_nonzero_flag = " << wS1_nonzero_flag << endl;
    155 
    156     ## 5: Check Bad-type II
    157     if (wS1_nonzero_flag == False):
    158         #print " Bad II Soln :  " + str(w)
    159         return 5
    160 
    161 
    162     ## Error if we get here! =o
    163     print "   Solution vector is " + str(w)
    164     print "   and Q is \n" + str(self) + "\n"
    165     raise RuntimeError, "Error in IsLocalSolutionType: Should not execute this line... =( \n"
    166 
    167 
    168 
    169 
    170 
    171 
    172 def CountAllLocalTypesNaive__deprecated(self, p, k, m, zvec, nzvec):
    173     """
    174     ///////////////////////////////////////////////////////////////////
    175     /// Naively counts the number of solutions of Q(x) = m mod p^k   //
    176     /// of type solntype, satisfying the mod p congruence conditions //
    177     /// at the indices of the vectors "zero" and "nonzero"           //
    178     ///////////////////////////////////////////////////////////////////
    179 
    180     valarray <mpz_class> Matrix_mpz::CountAllLocalTypesNaive(const mpz_class & p, unsigned long k, const mpz_class & m,
    181                                              const valarray<size_t> & zero, const valarray<size_t> & nonzero) const
    182     """
    183 
    184     ## DIAGNOSTIC
    185     #print "   --> CountAllLocalTypesNaive is using the form Q \n" + str(Q)
    186     #print "       p = " + str(p) + "  and   m = " + str(m)
    187 
    188 
    189     if (True or (p == 2)):          ##  <----------------- THIS IS WIERD...
    190 
    191         n = self.dim()
    192         R = p ** k
    193 
    194         ## Initialize the counting vector
    195         count_vector = [0  for i in range(6)]
    196 
    197         ## Initialize v = (0, ... , 0)
    198         v = [0  for i in range(n)]
    199 
    200 
    201         ## Some declarations to speed up the loop
    202         R_n = R ** n
    203         m1 = m % R
    204 
    205         ## Count the local solutions
    206         #for(size_t i=1; i<=(R_n).get_ui(); i++):
    207         for i in range(R_n):
    208 #            self.VecIncrement(v, R)                          ## Increments the vector v in-place.
    209             VecIncrement_cdef(v, R)                          ## Increments the vector v in-place.
    210             if (self(v) % R  == m1):                        ## Evaluates the quadratic form (mod R) at the vector v -- can be sped up!
    211                 solntype = self.local_solution_type(p, v, zvec, nzvec)
    212                 if (solntype != 0):
    213                     count_vector[solntype] += 1
    214 
    215 
    216 
    217         ## Generate the Bad-type and Total counts
    218         count_vector[3] = count_vector[4] + count_vector[5]
    219         count_vector[0] = count_vector[1] + count_vector[2] + count_vector[3]
    220 
    221 
    222         ## DIAGNOSTIC
    223         #cout << "R = " << R << "\n";
    224         #cout << "n = " << n << "\n";
    225         #cout << "R_n = " << R_n << "\n";
    226         #
    227         #for(i=1; i<=25; i++) {
    228         #    cout << "v = " << v << "\n";
    229         #    Increment(v,R);
    230         #}
    231         #
    232         #cout << "Q1 = " << Q1 << "\n";
    233         #cout << "v = " << v << "\n";
    234         #cout << "Q1 * v = " << Q1 * v<< "\n";
    235        
    236 
    237         return count_vector   
    238 
    239     else:
    240         raise RuntimeError, "Error in count_local_typeNaive: Matrix \n" + str(self) + "\n is not symmetric!"
    24162
    24263
    24364
     
    24768##/// Front-ends for our counting routines //
    24869##///////////////////////////////////////////
    24970
    250 def count_local_type(self, p, k, m, solntype, zvec, nzvec):
     71def count_congruence_solutions(self, p, k, m, zvec, nzvec):
    25172    """
    252     mpz_class Matrix_mpz::count_local_type(const mpz_class & p, long k, const mpz_class & m, size_t solntype,
    253                          const valarray<size_t> & zero, const valarray<size_t> & nonzero) const
     73    Counts all solutions of Q(x) = m (mod p^k) satisfying the
     74    additional congruence conditions described in
     75    QuadraticForm.count_congruence_solutions_as_vector().
    25476
    255     // Ideally this would use the count_local_typeWithSymmetry routine, but this is fine for now. =)
     77    INPUT:
     78        p -- prime number > 0
     79        k -- an integer > 0
     80        m -- an integer (depending only on mod p^k)
     81        zvec, nzvec -- a list of integers in range(self.dim()), or None
     82
     83    EXAMPLES:
     84        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
     85        sage: Q.count_congruence_solutions(3, 1, 0, None, None)
     86        15
     87
    25688    """
    257     #return self.CountAllLocalTypesNaive(p, k, m, zvec, nzvec)[0]
    25889    return CountAllLocalTypesNaive(self, p, k, m, zvec, nzvec)[0]
    25990
    26091
    261 def count_local_good_type(self, p, k, m, zvec, nzvec):
     92
     93def count_congruence_solutions__good_type(self, p, k, m, zvec, nzvec):
    26294    """
    263     mpz_class Matrix_mpz::count_local_good_type(const mpz_class & p, long k, const mpz_class & m,
    264                              const valarray<size_t> & zero, const valarray<size_t> & nonzero) const
     95    Counts the good-type solutions of Q(x) = m (mod p^k) satisfying the
     96    additional congruence conditions described in
     97    QuadraticForm.count_congruence_solutions_as_vector().
     98
     99    INPUT:
     100        p -- prime number > 0
     101        k -- an integer > 0
     102        m -- an integer (depending only on mod p^k)
     103        zvec, nzvec -- a list of integers up to dim(Q)
     104
     105    EXAMPLES:
     106        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
     107        sage: Q.count_congruence_solutions__good_type(3, 1, 0, None, None)
     108        12
    265109
    266110    """
    267     #return self.CountAllLocalTypesNaive(p, k, m, zvec, nzvec)[1]
    268111    return CountAllLocalTypesNaive(self, p, k, m, zvec, nzvec)[1]
    269112
    270113
    271 def count_local_zero_type(self, p, k, m, zvec, nzvec):
     114
     115def count_congruence_solutions__zero_type(self, p, k, m, zvec, nzvec):
    272116    """
    273     mpz_class Matrix_mpz::count_local_zero_type(const mpz_class & p, long k, const mpz_class & m,
    274                              const valarray<size_t> & zero, const valarray<size_t> & nonzero) const
     117    Counts the zero-type solutions of Q(x) = m (mod p^k) satisfying the
     118    additional congruence conditions described in
     119    QuadraticForm.count_congruence_solutions_as_vector().
     120
     121    INPUT:
     122        p -- prime number > 0
     123        k -- an integer > 0
     124        m -- an integer (depending only on mod p^k)
     125        zvec, nzvec -- a list of integers up to dim(Q)
     126
     127    EXAMPLES:
     128        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
     129        sage: Q.count_congruence_solutions__zero_type(3, 1, 0, None, None)
     130        1
    275131
    276132    """
    277     #return self.CountAllLocalTypesNaive(p, k, m, zvec, nzvec)[2]
    278133    return CountAllLocalTypesNaive(self, p, k, m, zvec, nzvec)[2]
    279134
    280135
    281 def count_local_bad_type(self, p, k, m, zvec, nzvec):
     136def count_congruence_solutions__bad_type(self, p, k, m, zvec, nzvec):
    282137    """
    283     mpz_class Matrix_mpz::count_local_bad_type(const mpz_class & p, long k, const mpz_class & m,
    284                              const valarray<size_t> & zero, const valarray<size_t> & nonzero) const
     138    Counts the bad-type solutions of Q(x) = m (mod p^k) satisfying the
     139    additional congruence conditions described in
     140    QuadraticForm.count_congruence_solutions_as_vector().
     141
     142    INPUT:
     143        p -- prime number > 0
     144        k -- an integer > 0
     145        m -- an integer (depending only on mod p^k)
     146        zvec, nzvec -- a list of integers up to dim(Q)
     147
     148    EXAMPLES:
     149        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
     150        sage: Q.count_congruence_solutions__bad_type(3, 1, 0, None, None)
     151        2
    285152
    286153    """
    287     #return self.CountAllLocalTypesNaive(p, k, m, zvec, nzvec)[3]
    288154    return CountAllLocalTypesNaive(self, p, k, m, zvec, nzvec)[3]
    289155
    290156
    291 def count_local_bad_typeI(self, p, k, m, zvec, nzvec):
     157def count_congruence_solutions__bad_type_I(self, p, k, m, zvec, nzvec):
    292158    """
    293     mpz_class Matrix_mpz::count_local_bad_typeI(const mpz_class & p, long k, const mpz_class & m,
    294                              const valarray<size_t> & zero, const valarray<size_t> & nonzero) const
     159    Counts the bad-typeI solutions of Q(x) = m (mod p^k) satisfying
     160    the additional congruence conditions described in
     161    QuadraticForm.count_congruence_solutions_as_vector().
     162
     163    INPUT:
     164        p -- prime number > 0
     165        k -- an integer > 0
     166        m -- an integer (depending only on mod p^k)
     167        zvec, nzvec -- a list of integers up to dim(Q)
     168
     169    EXAMPLES:
     170        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
     171        sage: Q.count_congruence_solutions__bad_type_I(3, 1, 0, None, None)
     172        0
    295173
    296174    """
    297     #return self.CountAllLocalTypesNaive(p, k, m, zvec, nzvec)[4]
    298175    return CountAllLocalTypesNaive(self, p, k, m, zvec, nzvec)[4]
    299176
    300177
    301 def count_local_bad_typeII(self, p, k, m, zvec, nzvec):
     178def count_congruence_solutions__bad_type_II(self, p, k, m, zvec, nzvec):
    302179    """
    303     mpz_class Matrix_mpz::count_local_bad_typeII(const mpz_class & p, long k, const mpz_class & m,
    304                              const valarray<size_t> & zero, const valarray<size_t> & nonzero) const
     180    Counts the bad-typeII solutions of Q(x) = m (mod p^k) satisfying
     181    the additional congruence conditions described in
     182    QuadraticForm.count_congruence_solutions_as_vector().
     183
     184    INPUT:
     185        p -- prime number > 0
     186        k -- an integer > 0
     187        m -- an integer (depending only on mod p^k)
     188        zvec, nzvec -- a list of integers up to dim(Q)
     189
     190    EXAMPLES:
     191        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
     192        sage: Q.count_congruence_solutions__bad_type_II(3, 1, 0, None, None)
     193        2
    305194
    306195    """
    307     #return self.CountAllLocalTypesNaive(p, k, m, zvec, nzvec)[5]
    308196    return CountAllLocalTypesNaive(self, p, k, m, zvec, nzvec)[5]
  • sage/quadratic_forms/quadratic_form__local_density_congruence.py

    diff -r af0f19cc490e -r b3d0bbb703ab sage/quadratic_forms/quadratic_form__local_density_congruence.py
    a b  
    1 
    2 ## Include headers from the front-end (repeated execution) routines.
    3 #include <assert.h>
    4 
    5 
    6 
    7 #####################################################
    8 ## Historical Note:
    9 ## ----------------
    10 ## (Much) Older versions of the NTL code can be found at
    11 ## Old-Laptop-Odyssues_Extended-Home-directory_8-28-2003/home/jonhanke/TEMP/C++/Modular_Project/
    12 ##
    13 ##  Updated for GMP libraries on 1/21/04
    14 #####################################################
     1##########################################################################
     2## Methods which compute the local densities for representing a number
     3## by a quadratic form at a prime (possibly subject to additional
     4## congruence conditions).
     5##########################################################################
    156
    167from copy import deepcopy
    178
     
    2112from sage.rings.integer_ring import ZZ
    2213from sage.misc.misc import prod, verbose
    2314
    24 ##############################################################################
    25 ### Creates a new index vector which points to the extracted rows/columns
    26 ### in the extracted matrix where the ...?
    27 ###
    28 ### This is used internally in our reduction procedure which allows us
    29 ### to keep track of a list of (row/column) indices through the
    30 ### process of extracting a submatrix corresponding to some other set
    31 ### of (row/column) indices.
    32 ###
    33 ### (Note: This is probably not very efficient,
    34 ### and could be improved by using vectors,
    35 ### but we're using valarrays for now...) 
    36 #############################################################################
     15from sage.quadratic_forms.count_local_2 import count_modp__by_gauss_sum, extract_sublist_indices
    3716
    3817
    39 def reindex_vector_from_extraction(self, Original, Extracted):
     18
     19
     20def count_modp_solutions__by_Gauss_sum(self, p, m):
    4021    """
    41     This takes a vector of indices, and applies a  truncated permutation to them,
    42     defined by Extracted.
     22    Returns the number of solutions of Q(x) = m (mod p) of a
     23    non-degenerate quadratic form over the finite field Z/pZ,
     24    where p is a prime number > 2. 
    4325
    44     TO DO: Please revisit this routine, and eliminate it!
     26    Note: We adopt the useful convention that a zero-dimensional
     27    quadratic form has exactly one solution always (i.e. the empty
     28    vector).
    4529
    46     valarray<size_t> Matrix_mpz::reindex_vector_from_extraction(const valarray<size_t> & Original, const valarray<size_t> & Extracted) const
     30    These are defined in Table 1 on p363 of Hanke's "Local
     31    Densities..." paper.
     32
     33    INPUT:
     34        p -- a prime number > 2
     35        m -- an integer
     36
     37    OUTPUT:
     38        an integer >= 0
     39
     40    EXAMPLES:
     41        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1])
     42        sage: [Q.count_modp_solutions__by_Gauss_sum(3, m)  for m in range(3)]
     43        [9, 6, 12]
     44
     45        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,2])
     46        sage: [Q.count_modp_solutions__by_Gauss_sum(3, m)  for m in range(3)]
     47        [9, 12, 6]
     48
    4749    """
    48     ## Set some default values (all to be overwritten)
    49     Reindex = [-1  for i in range(len(Original))]
    50  
    51     ## Replace all Original indices entries with the position of the matching Extraced index
    52     ind = 0
    53     for i in range(len(Original)):
    54         for i in range(len(Extracted)):
    55             if (Original[i] == Extracted[j]):
    56                 Reindex[ind] = j
    57                 ind += 1
    58  
    59     ## Copy these to a new vecarray of the appropriate length -- Since Reindex may be too big
    60     Final = [Reindex[i]  for i in range(ind)]
    61  
    62     return Final
     50    if self.dim() == 0:
     51        return 1
     52    else:
     53        return count_modp__by_gauss_sum(self.dim(), p, m, self.Gram_det())
    6354
    6455
    6556
    6657
    67 def count_modp__by_gauss_sum(self, n, p, m, Qdet):
    68     """
    69     Returns the number of solutions of Q(x) = m over the finite field F_p,
    70     where p is a prime number > 2 and Q has determinant Qdet.
    71 
    72     These are defined in Table 1 on p363 of Hanke's "Local Densities..." paper.
    73     """
    74     neg1 = -1
    75 
    76     ## DIAGNOSTIC
    77     verbose(" n = " + str(n))
    78     verbose(" neg1 = " + str(neg1))
    79     verbose(" p**(n-1) = " + str(p**(n-1)))
    80    
    81     if ((p != 2) and (n >= 1)):   ## TO DO: Check that p is an odd prime and n >= 1  --  To Do: Need to check if p is a prime...
    82 
    83         if (m % p == 0):
    84             if (n % 2 != 0):
    85                 count = (p**(n-1))
    86             else:
    87                 count = (p**(n-1)) + (p-1) * (p**((n-2)/2)) * kronecker_symbol(((neg1**(n/2)) * Qdet) % p, p)
    88         else:
    89             if (n % 2 != 0):
    90                 count = (p**(n-1)) + (p**((n-1)/2)) * kronecker_symbol(((neg1**((n-1)/2)) * Qdet * m) % p, p)
    91             else:
    92                 count = (p**(n-1)) - (p**((n-2)/2)) * kronecker_symbol(((neg1**(n/2)) * Qdet) % p, p)
    93 
    94     else:
    95         raise RuntimeError, "\n Error in count_modp__by_gauss_sum: Either p is not prime or n<1. \n"
    96  
    97     return count
    98 
    99 
    100 
    101 
    102 
     58 
    10359def local_good_density_congruence_odd(self, p, m, Zvec, NZvec):
    10460    """
    10561    Finds the Good-type local density of Q representing m at p. 
    10662    (Assuming that p > 2 and Q is given in local diagonal form.)
    10763
    108     mpq_class Matrix_mpz::local_good_density_congruence_odd(const mpz_class & p, const mpz_class & m,
    109                                             const valarray<size_t> & Zvec, const valarray<size_t> & NZvec) const
     64    The additional congruence condition arguments Zvec and NZvec can
     65    be either a list of indices or None.  Zvec = [] is equivalent to
     66    Zvec = None which both impose no additional conditions, but NZvec
     67    = [] returns no solutions always while NZvec = None imposes no
     68    additional condition.
     69
     70    TO DO: Add type checking for Zvec, NZvec, and that Q is in local
     71    normal form.
     72
     73    INPUT:
     74        Q -- quadratic form assumed to be diagonal and p-integral
     75        p -- a prime number
     76        m -- an integer
     77        Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
     78
     79    OUTPUT:
     80        a rational number
     81
     82    EXAMPLES:
     83        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
     84        sage: Q.local_good_density_congruence_odd(3, 1, None, None)
     85        2/3
     86
     87        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
     88        sage: Q.local_good_density_congruence_odd(3, 1, None, None)
     89        8/9
     90
    11091    """
     92    n = self.dim()
    11193
    112     ## To Do: Check to see if Q_int is diagonal (only ok for p>2)
    113  
     94    ## Put the Zvec congruence condition in a standard form
     95    if Zvec == None:
     96        Zvec = []
    11497
    115     n = self.dim()
    116  
    117  
     98
     99    ## Sanity Check on Zvec and NZvec:
     100    ## -------------------------------
     101    Sn = Set(range(n))
     102    if (Zvec != None) and (len(Set(Zvec) + Sn) > n):
     103        raise RuntimeError, "Zvec must be a subset of {0, ..., n-1}."
     104    if (NZvec != None) and (len(Set(NZvec) + Sn) > n):
     105        raise RuntimeError, "NZvec must be a subset of {0, ..., n-1}."
     106
     107
     108
    118109    ## Assuming Q is diagonal, find the indices of the p-unit (diagonal) entries
    119     Qtrimvec = [i  for i in range(n)  if (self.__getitem__((i,i)) % p) != 0]
    120    
     110    UnitVec = [i  for i in range(n)  if (self[i,i] % p) != 0]
     111    NonUnitVec = list(Set(range(n)) - Set(UnitVec))
    121112
    122113
    123     ## DEBUGGING:  Check the matrix is primitive --
    124     ##   (WARNING: We may be passed imprimitive form from the BI-reduction...)  <==  This should be fixed now. =)
    125     if (Qtrimvec == []):
    126         raise RuntimeError, "Uh oh.  There should be some p-unit diagonal elements at this stage...  but Qtrimvec is empty. =("
     114    ## Take cases on the existence of additional non-zero congruence conditions (mod p)
     115    UnitVec_minus_Zvec = list(Set(UnitVec) - Set(Zvec))
     116    NonUnitVec_minus_Zvec = list(Set(NonUnitVec) - Set(Zvec))
     117    Q_Unit_minus_Zvec = self.extract_variables(UnitVec_minus_Zvec)
    127118
    128 
    129 
    130     ## DIAGNOSTIC
    131     #print " Stage 1"
    132     #print " Q = " + str(Q)
    133     #print " n = " + str(n) + "  and  Qtrimvec = " + str(Qtrimvec)
    134 
    135 
    136     Qtrim = self.extract_variables(Qtrimvec)
    137 
    138 
    139     ## print "We're here now..." << endl;
    140 
    141  
    142     ## Construct the big and small matrices:
    143     ## -------------------------------------
    144 
    145     ## Construct new congruence condition indices for the trimmed matrix
    146     trimZvec = self.reindex_vector_from_extraction(Zvec, Qtrimvec)
    147 
    148 
    149     ## DIAGNOSTIC
    150     #print "\n Zvec = " + str(Zvec)
    151     #print " NZvec = " + str(NZvec)
    152     #print " Qtrim = " + str(Qtrim)
    153     #print " Qtrimvec = " + str(Qtrimvec)
    154     #print " trimZvec = " + str(trimZvec)
    155     #####  print " Vector indexing Qtrim = " << range(len(Qtrimvec))
    156     #####  print " and its complement by trimZvec  = " << VectorComplement(range(len(Qtrimvec)), trimZvec)
    157 
    158 
    159 
    160     ## print "Getting closer..."
    161 
    162 
    163     ## Make the big trim vector (all trim indices not somewhere in Zvec)
    164     ##   and the big free vector (all non-trim indices not somewhere in Zvec)
    165     ##   [To Do: This could probably be faster if we assumed Zvec was ordered.]
    166     bigtrimvec = [ i  for i in range(len(Qtrimvec))  if i not in trimZvec]
    167     bigfreelength = (n - len(Qtrimvec)) - (len(Zvec) - len(trimZvec))
    168 
    169 
    170 
    171 
    172     ## Make the small vector from the big one (all indices not somewhere in Zvec or NZvec)
    173     new_vec = list(Set(Zvec + NZvec))
    174     trim_new_vec = self.reindex_vector_from_extraction(new_vec, Qtrimvec)
    175 
    176 
    177     smalltrimvec = [ i  for i in range(len(Qtrimvec))  if i not in new_vec]
    178     smallfreelength = (n - len(Qtrimvec)) - (len(new_vec) - len(trim_new_vec))
    179 
    180 
    181 
    182     ## Make the big and small matrices
    183     bigmatrix = Qtrim.extract_variables(bigtrimvec)
    184     smallmatrix = Qtrim.extract_variables(smalltrimvec)
    185 
    186 
    187     ## DIAGNOSTIC
    188     #print "\n Q is : " + str(self)
    189     #print " m is : " + str(m)
    190     #print " p is : " + str(p)
    191     #print " Qtrim is : " + str(Qtrim)
    192     #print " bigtrimvec is : " + str(bigtrimvec)
    193     #print " bigfreelength is : " + str(bigfreelength)
    194     #print " smalltrimvec is : " + str(smalltrimvec)
    195     #print " smallfreelength is : " + str(smallfreelength)
    196     #print " bigmatrix is : \n" + str(bigmatrix)
    197     #print " smallmatrix is : \n" + str(smallmatrix)
    198 
    199 
    200     ## print "And closer still..."
    201 
    202  
    203     ## Number of representations
    204  
    205 
    206     ## Adjust the local solutions to count only the good-type ones
    207     if (bigtrimvec == []):    ## Check if we have the empty matrix...
    208         big_factor = 0
     119    if (NZvec == None):
     120        if m % p != 0:
     121            total = Q_Unit_minus_Zvec.count_modp_solutions__by_Gauss_sum(p, m) * p**len(NonUnitVec_minus_Zvec)          ## m != 0 (mod p)
     122        else:
     123            total = (Q_Unit_minus_Zvec.count_modp_solutions__by_Gauss_sum(p, m) - 1) * p**len(NonUnitVec_minus_Zvec)    ## m == 0 (mod p)
    209124
    210125    else:
    211         big_dim = bigmatrix.dim()
    212         big_det = prod([bigmatrix[i,i]  for i in range(bigmatrix.dim())])   ## This gives the product of the diagonal entries, not the Gram matrix det
     126        UnitVec_minus_ZNZvec = list(Set(UnitVec) - (Set(Zvec) + Set(NZvec)))
     127        NonUnitVec_minus_ZNZvec = list(Set(NonUnitVec) - (Set(Zvec) + Set(NZvec)))
     128        Q_Unit_minus_ZNZvec = self.extract_variables(UnitVec_minus_ZNZvec)
    213129
    214         if (m % p != 0):
    215             big_factor = (p**bigfreelength) \
    216                 * self.count_modp__by_gauss_sum(big_dim, p, m, big_det)
    217         else:
    218             big_factor = (p**bigfreelength) \
    219                 * (self.count_modp__by_gauss_sum(big_dim, p, m, big_det) - 1)
     130        if m % p != 0:         ## m != 0 (mod p)
     131            total = Q_Unit_minus_Zvec.count_modp_solutions__by_Gauss_sum(p, m) * p**len(NonUnitVec_minus_Zvec) \
     132                    - Q_Unit_minus_ZNZvec.count_modp_solutions__by_Gauss_sum(p, m) * p**len(NonUnitVec_minus_ZNZvec)
     133        else:                  ## m == 0 (mod p)
     134            total = (Q_Unit_minus_Zvec.count_modp_solutions__by_Gauss_sum(p, m) - 1) * p**len(NonUnitVec_minus_Zvec) \
     135                    - (Q_Unit_minus_ZNZvec.count_modp_solutions__by_Gauss_sum(p, m) - 1) * p**len(NonUnitVec_minus_ZNZvec)
    220136
    221 
    222     ## Similarly for the smallmatrix if it exists
    223     if (smalltrimvec == []):    ## Check if we have the empty matrix...
    224         small_factor = 0
    225 
    226     else:
    227         if (len(NZvec) > 0):
    228             small_dim = smallmatrix.dim()
    229             small_det = prod([smallmatrix[i,i]  for i in range(smallmatrix.dim())])   ## This gives the product of the diagonal entries, not the Gram matrix det
    230 
    231             if (m % p != 0):
    232                 small_factor = (p**smallfreelength) \
    233                     * self.count_modp__by_gauss_sum(small_dim, p, m, small_det)
    234             else:
    235                 small_factor = (p**smallfreelength) \
    236                     * (self.count_modp__by_gauss_sum(small_dim, p, m, small_det) - 1)
    237  
    238         else:
    239             small_factor = 0
    240  
    241     total = big_factor - small_factor
    242     good_density = QQ(total / p**(n-1))
    243  
    244  
    245     ## DIAGNOSTIC
    246     #print "\n big factor = " + str(big_factor)
    247     #print " small factor = " + str(small_factor)
    248     #print " big_det = " + str(prod([bigmatrix[i,i]  for i in range(bigmatrix.dim())]))
    249     #print " smal_det = " + str(prod([smallmatrix[i,i]  for i in range(smallmatrix.dim())]))
    250     #print " total = " + str(total)
    251     #print " denominator = " + str(p**(n-1))
    252     #print " Good Density = " + str(good_density)
    253     #print " count_modp__by_gauss_sum output = " + str(self.count_modp__by_gauss_sum(big_dim, p, m, prod([bigmatrix[i,i]  for i in range(bigmatrix.dim())])))
    254 
     137    ## Return the Good-type representation density
     138    good_density = QQ(total) / p**(n-1)
    255139    return good_density
    256140
    257141
    258142
    259143
     144def local_good_density_congruence_even(self, m, Zvec, NZvec):
     145    """
     146    Finds the Good-type local density of Q representing m at p=2. 
     147    (Assuming Q is given in local diagonal form.)
    260148
    261 def local_good_density_congruence_even(self, p, m, Zvec, NZvec):
     149
     150    The additional congruence condition arguments Zvec and NZvec can
     151    be either a list of indices or None.  Zvec = [] is equivalent to
     152    Zvec = None which both impose no additional conditions, but NZvec
     153    = [] returns no solutions always while NZvec = None imposes no
     154    additional condition.
     155
     156    WARNING: Here the indices passed in Zvec and NZvec represent
     157    indices of the solution vector x of Q(x) = m (mod p^k), and *not*
     158    the Jordan components of Q.  They therefore are required (and
     159    assumed) to include either all or none of the indices of a given
     160    Jordan component of Q.  This is only important when p=2 since
     161    otherwise all Jordan blocks are 1x1, and so there the indices and
     162    Jordan blocks coincide.
     163
     164    TO DO: Add type checking for Zvec, NZvec, and that Q is in local
     165    normal form.
     166
     167
     168    INPUT:
     169        Q -- quadratic form assumed to be block diagonal and 2-integral
     170        p -- a prime number
     171        m -- an integer
     172        Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
     173
     174    OUTPUT:
     175        a rational number
     176
     177    EXAMPLES:
     178        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
     179        sage: Q.local_good_density_congruence_even(1, None, None)
     180        1
     181
     182        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
     183        sage: Q.local_good_density_congruence_even(1, None, None)
     184        1
     185        sage: Q.local_good_density_congruence_even(2, None, None)
     186        3/2
     187        sage: Q.local_good_density_congruence_even(3, None, None)
     188        1
     189        sage: Q.local_good_density_congruence_even(4, None, None)
     190        1/2
     191
     192        sage: Q = QuadraticForm(ZZ, 4, range(10))
     193        sage: Q[0,0] = 5
     194        sage: Q[1,1] = 10
     195        sage: Q[2,2] = 15
     196        sage: Q[3,3] = 20
     197        sage: Q
     198        Quadratic form in 4 variables over Integer Ring with coefficients:
     199        [ 5 1 2 3 ]
     200        [ * 10 5 6 ]
     201        [ * * 15 8 ]
     202        [ * * * 20 ]
     203        sage: Q.theta_series(20)
     204        1 + 2*q^5 + 2*q^10 + 2*q^14 + 2*q^15 + 2*q^16 + 2*q^18 + O(q^20)
     205        sage: Q.local_normal_form(2)
     206        Quadratic form in 4 variables over Integer Ring with coefficients:
     207        [ 0 1 0 0 ]
     208        [ * 0 0 0 ]
     209        [ * * 0 1 ]
     210        [ * * * 0 ]
     211        sage: Q.local_good_density_congruence_even(1, None, None)
     212        3/4
     213        sage: Q.local_good_density_congruence_even(2, None, None)
     214        1
     215        sage: Q.local_good_density_congruence_even(5, None, None)
     216        3/4
     217
    262218    """
    263     Finds the Good-type local density of Q representing m at p. 
    264     (Assuming that p = 2 and Q is given in local normal form.)
     219    n = self.dim()
     220   
     221    ## Put the Zvec congruence condition in a standard form
     222    if Zvec == None:
     223        Zvec = []
    265224
    266     mpq_class Matrix_mpz::local_good_density_congruence_even(const mpz_class & p, const mpz_class & m,
    267                                              const valarray<size_t> & Zvec, const valarray<size_t> & NZvec) const
    268     """
    269     #print " Break 0"
    270     #print "\n Q is : " + str(Q)
    271225
    272     n = self.dim()
    273  
    274  
    275     ## Assuming Q is diagonal, trim it to ensure it's non-degenerate mod 8      <--- ??? Do we really assume that Q is diagonal?!?!?  =(
    276     Qtrimvec = []
     226    ## Sanity Check on Zvec and NZvec:
     227    ## -------------------------------
     228    Sn = Set(range(n))
     229    if (Zvec != None) and (len(Set(Zvec) + Sn) > n):
     230        raise RuntimeError, "Zvec must be a subset of {0, ..., n-1}."
     231    if (NZvec != None) and (len(Set(NZvec) + Sn) > n):
     232        raise RuntimeError, "NZvec must be a subset of {0, ..., n-1}."
    277233
    278     #print " Break 0.1"
    279234
    280     ## Find the indices of the non-zero blocks mod 8
     235
     236    ## Find the indices of x for which the associated Jordan blocks are non-zero mod 8    TO DO: Move this to special Jordan block code separately!
     237    ## -------------------------------------------------------------------------------
     238    Not8vec = []
    281239    for i in range(n):
    282240 
    283241        ## DIAGNOSTIC
    284242        verbose(" i = " + str(i))
    285243        verbose(" n = " + str(n))
    286         verbose(" Qtrimvec = " + str(Qtrimvec))
     244        verbose(" Not8vec = " + str(Not8vec))
    287245
    288246        nz_flag = False
    289247
     248        ## Check if the diagonal entry isn't divisible 8
    290249        if  ((self[i,i] % 8) != 0): 
    291250            nz_flag = True
     251
     252        ## Check appropriate off-diagonal entries aren't divisible by 8
    292253        else:
    293             #print " here 1"
     254
     255            ## Special check for first off-diagonal entry
    294256            if ((i == 0) and ((self[i,i+1] % 8) != 0)):
    295257                nz_flag = True
     258
     259            ## Special check for last off-diagonal entry
     260            elif ((i == n-1) and ((self[i-1,i] % 8) != 0)):
     261                nz_flag = True
     262
     263            ## Check for the middle off-diagonal entries
    296264            else:
    297                 #print " here 2" << endl;
    298                 if ((i == n-1) and ((self[i-1,i] % 8) != 0)):
     265                if ( (i > 0)  and  (i < n-1)  and  (((self[i,i+1] % 8) != 0) or ((self[i-1,i] % 8) != 0)) ):
    299266                    nz_flag = True
    300                 else:
    301                     #print " here 3" << endl;
    302                     if ( (i > 0)  and  (i < n-1)  and  (((self[i,i+1] % 8) != 0) or ((self[i-1,i] % 8) != 0)) ):
    303                         nz_flag = True
     267                   
     268        ## Remember the (vector) index if it's not part of a Jordan block of norm divisible by 8
     269        if (nz_flag == True):
     270            Not8vec += [i]
    304271
    305      
    306         if (nz_flag == True):
    307             Qtrimvec += [i]
    308272
    309273
    310     #print " Break 1"
     274    ## Compute the number of Good-type solutions mod 8:
     275    ## ------------------------------------------------
    311276
     277    ## Setup the indexing sets for additional zero congruence solutions
     278    Q_Not8 = self.extract_variables(Not8vec)
     279    Not8 = Set(Not8vec)
     280    Is8 = Set(range(n)) - Not8
    312281
    313     ## DEBUGGING: Quick tests to make sure the form isn't zero mod 8
    314     if (Qtrimvec == []):
    315         raise RuntimeError, "Oh no!  The quadratic form should always be non-degenerate mod 8... =("
    316  
     282    Z = Set(Zvec)
     283    Z_Not8 = Not8.intersection(Z)
     284    Z_Is8 = Is8.intersection(Z)
     285    Is8_minus_Z = Is8 - Z_Is8
    317286
    318     Qtrim = self.extract_variables(Qtrimvec)
    319 
    320 
    321     ## DEBUGGING: Quick tests to make sure the extracted form isn't zero mod 8
    322     assert(Qtrim.dim() > 0);
    323 
    324     #print " Break 2"
    325 
    326 
    327     ## Construct the big and small matrices:
    328     ## -------------------------------------
    329 
    330     ## Construct new congruence condition indices for the trimmed matrix
    331     trimZvec = self.reindex_vector_from_extraction(Zvec, Qtrimvec)
    332     trimNZvec = self.reindex_vector_from_extraction(NZvec, Qtrimvec)
    333 
    334     ## Make the trimmed congruence vector
    335     new_vec = list(Set(Zvec + NZvec))
    336     trim_vec = self.reindex_vector_from_extraction(new_vec, Qtrimvec)
    337287
    338288    ## DIAGNOSTIC
    339     verbose("")
    340     verbose("Zvec = " + str(Zvec))
    341     verbose("NZvec = " + str(NZvec))
    342     verbose("new_vec = " + str(new_vec))
    343     verbose("")
     289    verbose("Z = " + str(Z))
     290    verbose("Z_Not8 = " + str(Z_Not8))
     291    verbose("Z_Is8 = " + str(Z_Is8))
     292    verbose("Is8_minus_Z = " + str(Is8_minus_Z))
     293
     294
     295    ## Take cases on the existence of additional non-zero congruence conditions (mod 2)   
     296    if NZvec == None:
     297        total = (4 ** len(Z_Is8)) * (8 ** len(Is8_minus_Z)) \
     298            * Q_Not8.count_congruence_solutions__good_type(2, 3, m, list(Z_Not8), None)
     299    else:
     300        ZNZ = Z + Set(NZvec)
     301        ZNZ_Not8 = Not8.intersection(ZNZ)
     302        ZNZ_Is8 = Is8.intersection(ZNZ)
     303        Is8_minus_ZNZ = Is8 - ZNZ_Is8
     304
     305        ## DIAGNOSTIC
     306        verbose("ZNZ = " + str(ZNZ))
     307        verbose("ZNZ_Not8 = " + str(ZNZ_Not8))
     308        verbose("ZNZ_Is8 = " + str(ZNZ_Is8))
     309        verbose("Is8_minus_ZNZ = " + str(Is8_minus_ZNZ))
     310
     311        total = (4 ** len(Z_Is8)) * (8 ** len(Is8_minus_Z)) \
     312            * Q_Not8.count_congruence_solutions__good_type(2, 3, m, list(Z_Not8), None) \
     313            - (4 ** len(ZNZ_Is8)) * (8 ** len(Is8_minus_ZNZ)) \
     314            * Q_Not8.count_congruence_solutions__good_type(2, 3, m, list(ZNZ_Not8), None)
    344315
    345316
    346317    ## DIAGNOSTIC
    347     verbose("\n Q is : " + str(self))
    348     verbose(" m is : " + str(m))
    349     verbose(" Qtrim is : " + str(Qtrim))
    350     verbose(" Qtrimvec is : " + str(Qtrimvec))
    351     verbose(" trimZvec is : " + str(trimZvec))
    352     verbose(" trimNZvec is : " + str(trimNZvec))
     318    verbose("total = " + str(total))
    353319
    354320
    355     ## DEBUGGING: Check that partlyfreenum is in range...
    356     if not (len(new_vec) >= len(trim_vec)):
    357         raise RuntimeError, "Oh no!  Some variable called 'partlyfreenum' is out of range. =("
    358 
    359  
    360     ## Compute the number of different free components
    361     partlyfreenum = len(new_vec) - len(trim_vec)
    362     veryfreenum = (n - len(Qtrimvec)) - partlyfreenum
    363 
    364 
    365     ## In the free part, each component with a congruence condition contrubite a factor of 4,
    366     ## while components with no congruence conditions contribute a factor of 8.
    367 
    368     total = (4 ** partlyfreenum) * (8 ** veryfreenum) \
    369         * Qtrim.count_local_good_type(2, 3, m, trimZvec, trimNZvec)
    370     good_density = ZZ(total) / ZZ(8**(n-1))
    371  
    372  
    373     ## DIAGNOSTIC
    374     verbose(" partlyfreenum = " + str(partlyfreenum))
    375     verbose(" veryfreenum = " + str(veryfreenum))
    376     verbose(" Qtrim.count_local_good_type(2, 3, m, trimZvec, trimNZvec) = " + str(Qtrim.count_local_good_type(2, 3, m, trimZvec, trimNZvec)))
    377     verbose(" total = " + str(total))
    378     verbose("    total has type ", type(total))
    379     verbose(" denominator = " + str(8**(n-1)))
    380     verbose(" Good Density = " + str(good_density))
    381 
     321    ## Return the associated Good-type representation density
     322    good_density = QQ(total) / 8**(n-1)
    382323    return good_density
    383324
    384325
    385326
    386327
    387 def local_good_density_congruence_even__experimental(self, p, m, Zvec, NZvec):
    388     """
    389     Finds the Good-type local density of Q representing m at p. 
    390     (Assuming that p = 2 and Q is given in local normal form.)
    391328
    392     mpq_class Matrix_mpz::local_good_density_congruence_even(const mpz_class & p, const mpz_class & m,
    393                                              const valarray<size_t> & Zvec, const valarray<size_t> & NZvec) const
    394     """
    395     #print " Break 0"
    396     #print "\n Q is : " + str(Q)
    397329
    398     n = self.dim()
    399  
    400  
    401     ## Assuming Q is diagonal, trim it to ensure it's non-degenerate mod 8      <--- ??? Do we really assume that Q is diagonal?!?!?  =(
    402     Qtrimvec = []
    403330
    404     #print " Break 0.1"
    405331
    406     ## Find the indices of the non-zero blocks mod 8
    407     for i in range(n):
    408  
    409         ## DIAGNOSTIC
    410         verbose(" i = " + str(i))
    411         verbose(" n = " + str(n))
    412         verbose(" Qtrimvec = " + str(Qtrimvec))
    413 
    414         nz_flag = False
    415 
    416         if  ((self[i,i] % 8) != 0): 
    417             nz_flag = True
    418         else:
    419             #print " here 1"
    420             if ((i == 0) and ((self[i,i+1] % 8) != 0)):
    421                 nz_flag = True
    422             else:
    423                 #print " here 2" << endl;
    424                 if ((i == n-1) and ((self[i-1,i] % 8) != 0)):
    425                     nz_flag = True
    426                 else:
    427                     #print " here 3" << endl;
    428                     if ( (i > 0)  and  (i < n-1)  and  (((self[i,i+1] % 8) != 0) or ((self[i-1,i] % 8) != 0)) ):
    429                         nz_flag = True
    430 
    431      
    432         if (nz_flag == True):
    433             Qtrimvec += [i]
    434 
    435 
    436     #print " Break 1"
    437 
    438 
    439     ## DEBUGGING: Quick tests to make sure the form isn't zero mod 8
    440     if (Qtrimvec == []):
    441         raise RuntimeError, "Oh no!  The quadratic form should always be non-degenerate mod 8... =("
    442  
    443 
    444     Qtrim = self.extract_variables(Qtrimvec)
    445 
    446 
    447     ## DEBUGGING: Quick tests to make sure the extracted form isn't zero mod 8
    448     assert(Qtrim.dim() > 0);
    449 
    450     #print " Break 2"
    451 
    452 
    453     ## Construct the big and small matrices:
    454     ## -------------------------------------
    455 
    456     ## Construct new congruence condition indices for the trimmed matrix
    457     trimZvec = self.reindex_vector_from_extraction(Zvec, Qtrimvec)
    458     trimNZvec = self.reindex_vector_from_extraction(NZvec, Qtrimvec)
    459 
    460     ## Make the trimmed congruence vector
    461     new_vec = list(Set(Zvec + NZvec))
    462     trim_vec = self.reindex_vector_from_extraction(new_vec, Qtrimvec)
    463 
    464     ## DIAGNOSTIC
    465     verbose("")
    466     verbose("Zvec = " + str(Zvec))
    467     verbose("NZvec = " + str(NZvec))
    468     verbose("new_vec = " + str(new_vec))
    469     verbose("")
    470 
    471 
    472     ## DIAGNOSTIC
    473     verbose("\n Q is : " + str(self))
    474     verbose(" m is : " + str(m))
    475     verbose(" Qtrim is : " + str(Qtrim))
    476     verbose(" Qtrimvec is : " + str(Qtrimvec))
    477     verbose(" trimZvec is : " + str(trimZvec))
    478     verbose(" trimNZvec is : " + str(trimNZvec))
    479 
    480 
    481     ## DEBUGGING: Check that partlyfreenum is in range...
    482     if not (len(new_vec) >= len(trim_vec)):
    483         raise RuntimeError, "Oh no!  Some variable called 'partlyfreenum' is out of range. =("
    484 
    485  
    486     ## Compute the number of different free components
    487     partlyfreenum = len(new_vec) - len(trim_vec)
    488     veryfreenum = (n - len(Qtrimvec)) - partlyfreenum
    489 
    490 
    491     ## In the free part, each component with a congruence condition contrubite a factor of 4,
    492     ## while components with no congruence conditions contribute a factor of 8.
    493 
    494     total = (4 ** partlyfreenum) * (8 ** veryfreenum) \
    495         * Qtrim.count_local_good_type(2, 3, m, trimZvec, trimNZvec)
    496     good_density = ZZ(total) / ZZ(8**(n-1))
    497  
    498  
    499     ## DIAGNOSTIC
    500     verbose(" partlyfreenum = " + str(partlyfreenum))
    501     verbose(" veryfreenum = " + str(veryfreenum))
    502     verbose(" Qtrim.count_local_good_type(2, 3, m, trimZvec, trimNZvec) = " + str(Qtrim.count_local_good_type(2, 3, m, trimZvec, trimNZvec)))
    503     verbose(" total = " + str(total))
    504     verbose("    total has type ", type(total))
    505     verbose(" denominator = " + str(8**(n-1)))
    506     verbose(" Good Density = " + str(good_density))
    507 
    508     return good_density
    509 
    510 
    511 
    512     ############################################
    513      ## Note: Assumes all forms are primitive       <== What does this mean?!?
    514     ############################################
    515 
    516 
    517 
    518 def local_good_density_congruence(self, p, m, Zvec, NZvec):
     332def local_good_density_congruence(self, p, m, Zvec=None, NZvec=None):
    519333    """
    520334    Finds the Good-type local density of Q representing m at p. 
    521335    (Front end routine for parity specific routines for p.)     
    522336
    523     mpq_class Matrix_mpz::local_good_density_congruence(const mpz_class & p, const mpz_class & m,
    524                                         const valarray<size_t> & Zvec, const valarray<size_t> & NZvec) const
     337    TO DO: Add Documentation about the additional congruence
     338    conditions Zvec and NZvec.
     339
     340
     341
     342    INPUT:
     343        Q -- quadratic form assumed to be block diagonal and p-integral
     344        p -- a prime number
     345        m -- an integer
     346        Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
     347
     348    OUTPUT:
     349        a rational number
     350
     351    EXAMPLES:
     352        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
     353        sage: Q.local_good_density_congruence(2, 1, None, None)
     354        1
     355        sage: Q.local_good_density_congruence(3, 1, None, None)
     356        2/3
     357
     358        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
     359        sage: Q.local_good_density_congruence(2, 1, None, None)
     360        1
     361        sage: Q.local_good_density_congruence(3, 1, None, None)
     362        8/9
     363
    525364    """
    526365    ## DIAGNOSTIC
    527366    verbose(" In local_good_density_congruence with ")
     
    531370    verbose(" Zvec = " + str(Zvec))
    532371    verbose(" NZvec = " + str(NZvec))
    533372
     373    ## Put the Zvec congruence condition in a standard form
     374    if Zvec == None:
     375        Zvec = []
     376
     377
     378    n = self.dim()
     379
     380    ## Sanity Check on Zvec and NZvec:
     381    ## -------------------------------
     382    Sn = Set(range(n))
     383    if (Zvec != None) and (len(Set(Zvec) + Sn) > n):
     384        raise RuntimeError, "Zvec must be a subset of {0, ..., n-1}."
     385    if (NZvec != None) and (len(Set(NZvec) + Sn) > n):
     386        raise RuntimeError, "NZvec must be a subset of {0, ..., n-1}."
     387
     388
    534389
    535390    ## Check that Q is in local normal form -- should replace this with a diagonalization check?
    536391    ##   (it often may not be since the reduction procedure
     
    540395    #    print "Warning in local_good_density_congruence: Q is not in local normal form! \n";
    541396
    542397
    543 
    544     ## Check that the congruence conditions don't overlap
    545     if (len(Set(Zvec).intersection(Set(NZvec))) != 0):
    546         return QQ(0)
    547 
    548     #print "We're here!"
    549398 
    550399
    551400    ## Decide which routine to use to compute the Good-type density
     
    554403
    555404    if (p == 2):
    556405        #print "\n Using the (p=2) Local_Good_Density_Even routine! \n"
    557         return self.local_good_density_congruence_even(p, m, Zvec, NZvec)
     406        return self.local_good_density_congruence_even(m, Zvec, NZvec)
    558407
    559408    raise RuntimeError, "\n Error in Local_Good_Density: The 'prime' p = " + str(p) + " is < 2. \n"
    560409
     
    564413
    565414
    566415
    567 def local_zero_density_congruence(self, p, m, Zvec, NZvec):
     416def local_zero_density_congruence(self, p, m, Zvec=None, NZvec=None):
    568417    """
    569418    Finds the Zero-type local density of Q representing m at p,
    570419    allowing certain congruence conditions mod p.
    571420
    572     mpq_class Matrix_mpz::local_zero_density_congruence(const mpz_class & p, const mpz_class & m,
    573                                         const valarray<size_t> & Zvec, const valarray<size_t> & NZvec) const
     421
     422    INPUT:
     423        Q -- quadratic form assumed to be block diagonal and p-integral
     424        p -- a prime number
     425        m -- an integer
     426        Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
     427
     428    OUTPUT:
     429        a rational number
     430
     431    EXAMPLES:
     432        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
     433        sage: Q.local_zero_density_congruence(2, 2, None, None)
     434        0
     435        sage: Q.local_zero_density_congruence(2, 4, None, None)
     436        1/2
     437        sage: Q.local_zero_density_congruence(3, 6, None, None)
     438        0
     439        sage: Q.local_zero_density_congruence(3, 9, None, None)
     440        2/9
     441
     442        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
     443        sage: Q.local_zero_density_congruence(2, 2, None, None)
     444        0
     445        sage: Q.local_zero_density_congruence(2, 4, None, None)
     446        1/4
     447        sage: Q.local_zero_density_congruence(3, 6, None, None)
     448        0
     449        sage: Q.local_zero_density_congruence(3, 9, None, None)
     450        8/81
     451
    574452    """
    575453    ## DIAGNOSTIC
    576     #print " In local_zero_density_congruence with "
    577     #print " Q is: \n" + str(self)
    578     #print " p = " + str(p)
    579     #print " m = " + str(m)
    580     #print " Zvec = " + str(Zvec)
    581     #print " NZvec = " + str(NZvec)
     454    verbose(" In local_zero_density_congruence with ")
     455    verbose(" Q is: \n" + str(self))
     456    verbose(" p = " + str(p))
     457    verbose(" m = " + str(m))
     458    verbose(" Zvec = " + str(Zvec))
     459    verbose(" NZvec = " + str(NZvec))
     460
     461    ## Put the Zvec congruence condition in a standard form
     462    if Zvec == None:
     463        Zvec = []
     464
     465
     466    n = self.dim()
     467
     468    ## Sanity Check on Zvec and NZvec:
     469    ## -------------------------------
     470    Sn = Set(range(n))
     471    if (Zvec != None) and (len(Set(Zvec) + Sn) > n):
     472        raise RuntimeError, "Zvec must be a subset of {0, ..., n-1}."
     473    if (NZvec != None) and (len(Set(NZvec) + Sn) > n):
     474        raise RuntimeError, "NZvec must be a subset of {0, ..., n-1}."
     475
    582476
    583477    p2 = p * p
    584478
    585     if ((m % (p2) != 0) or (len(NZvec) > 0)):
     479    ## Check some conditions for no zero-type solutions to exist
     480    if ((m % (p2) != 0) or (NZvec != None)):
    586481        return 0
    587     else:
    588         ## Need 2 steps since we can't take negative powers... =|
    589         EmptyVec = []
    590         if (self.dim() > 2):
    591             return QQ(1 / p**(self.dim() - 2)) \
    592                 * self.local_density_congruence(p, m / (p2), EmptyVec, EmptyVec);
    593         else:
    594             return QQ(p**(2 - self.dim())) \
    595                 * self.local_density_congruence(p, m / (p2), EmptyVec, EmptyVec);
    596482
     483    ## Use the reduction procedure to return the result
     484    return self.local_density_congruence(p, m / p2, None, None) / p**(self.dim() - 2)
    597485
    598486
    599487
    600488
    601489
    602 def local_badI_density_congruence(self, p, m, Zvec, NZvec):
     490
     491
     492def local_badI_density_congruence(self, p, m, Zvec=None, NZvec=None):
    603493    """
    604494    Finds the Bad-type I local density of Q representing m at p.
    605495    (Assuming that p > 2 and Q is given in local diagonal form.)
    606496
    607     mpq_class Matrix_mpz::local_badI_density_congruence(const mpz_class & p, const mpz_class & m,
    608                                         const valarray<size_t> & Zvec, const valarray<size_t> & NZvec) const
     497
     498    INPUT:
     499        Q -- quadratic form assumed to be block diagonal and p-integral
     500        p -- a prime number
     501        m -- an integer
     502        Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
     503
     504    OUTPUT:
     505        a rational number
     506
     507    EXAMPLES:
     508        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
     509        sage: Q.local_badI_density_congruence(2, 1, None, None)
     510        0
     511        sage: Q.local_badI_density_congruence(2, 2, None, None)
     512        1
     513        sage: Q.local_badI_density_congruence(2, 4, None, None)
     514        0
     515        sage: Q.local_badI_density_congruence(3, 1, None, None)
     516        0
     517        sage: Q.local_badI_density_congruence(3, 6, None, None)
     518        0
     519        sage: Q.local_badI_density_congruence(3, 9, None, None)
     520        0
     521
     522        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
     523        sage: Q.local_badI_density_congruence(2, 1, None, None)
     524        0
     525        sage: Q.local_badI_density_congruence(2, 2, None, None)
     526        0
     527        sage: Q.local_badI_density_congruence(2, 4, None, None)
     528        0
     529        sage: Q.local_badI_density_congruence(3, 2, None, None)
     530        0
     531        sage: Q.local_badI_density_congruence(3, 6, None, None)
     532        0
     533        sage: Q.local_badI_density_congruence(3, 9, None, None)
     534        0
     535
     536        sage: Q = DiagonalQuadraticForm(ZZ, [1,3,3,9])
     537        sage: Q.local_badI_density_congruence(3, 1, None, None)
     538        0
     539        sage: Q.local_badI_density_congruence(3, 3, None, None)
     540        4/3
     541        sage: Q.local_badI_density_congruence(3, 6, None, None)
     542        4/3
     543        sage: Q.local_badI_density_congruence(3, 9, None, None)
     544        0
     545        sage: Q.local_badI_density_congruence(3, 18, None, None)
     546        0
     547
     548
    609549    """
    610550    ## DIAGNOSTIC
    611     #print " In local_badI_density_congruence with "
    612     #print " Q is: \n" + str(self)
    613     #print " p = " + str(p)
    614     #print " m = " + str(m)
    615     #print " Zvec = " + str(Zvec)
    616     #print " NZvec = " + str(NZvec)
     551    verbose(" In local_badI_density_congruence with ")
     552    verbose(" Q is: \n" + str(self))
     553    verbose(" p = " + str(p))
     554    verbose(" m = " + str(m))
     555    verbose(" Zvec = " + str(Zvec))
     556    verbose(" NZvec = " + str(NZvec))
     557
     558    ## Put the Zvec congruence condition in a standard form
     559    if Zvec == None:
     560        Zvec = []
    617561
    618562
    619563    n = self.dim()
    620     if n == 0:
    621         raise RuntimeException, "Oops!  We don't currently allow 0-dim'l forms... =( "
    622564
    623565
    624     ## Define the indexing sets S_i
     566
     567    ## Sanity Check on Zvec and NZvec:
     568    ## -------------------------------
     569    Sn = Set(range(n))
     570    if (Zvec != None) and (len(Set(Zvec) + Sn) > n):
     571        raise RuntimeError, "Zvec must be a subset of {0, ..., n-1}."
     572    if (NZvec != None) and (len(Set(NZvec) + Sn) > n):
     573        raise RuntimeError, "NZvec must be a subset of {0, ..., n-1}."
     574
     575
     576
     577    ## Define the indexing set S_0, and determine if S_1 is empty:
     578    ## -----------------------------------------------------------
    625579    S0 = []
    626580    S1_empty_flag = True    ## This is used to check if we should be computing BI solutions at all!
    627581                            ## (We should really to this earlier, but S1 must be non-zero to proceed.)
    628582
    629 
    630583    ## Find the valuation of each variable (which will be the same over 2x2 blocks),
    631584    ## remembering those of valuation 0 and if an entry of valuation 1 exists.
    632585    for i in range(n):
     
    647600            S0 += [i]       
    648601        elif (val == 1):
    649602            S1_empty_flag = False    ## Need to have a non-empty S1 set to proceed with Bad-type I reduction...
     603
     604
     605
     606
    650607     
    651     ## Check that S1 is non-empty to proceed, otherwise return no solutions.
    652     if (S1_empty_flag == True):
     608    ## Check that S1 is non-empty and p|m to proceed, otherwise return no solutions.
     609    if (S1_empty_flag == True) or (m % p != 0):
    653610        return 0
    654611
     612    ## Check some conditions for no bad-type I solutions to exist
     613    if (NZvec != None) and (len(Set(S0).intersection(Set(NZvec))) != 0):
     614        return 0
    655615
    656     ## Check that the form is primitive...
     616
     617
     618    ## Check that the form is primitive...                     WHY DO WE NEED TO DO THIS?!?
    657619    if (S0 == []):
    658620        print " Using Q = " + str(self)
    659621        print " and p = " + str(p)
    660622        raise RuntimeError, "Oops! The form is not primitive!"
    661623
    662624
    663     S0New = S0    ## REDUNDANT ALERT!!
    664 
    665 
    666 
    667     ## Note: The following lines assume that NZvec is an ordered vector.  Should check this here...
    668625
    669626    ## DIAGNOSTIC
    670     #print " m = " + str(m) + "   p = " + str(p)
    671     #print " S0 = " + str(S0) + "   NZvec = " str(NZvec) "   IsDisjoint = " + str(IsDisjointOrdered(S0, NZvec))
    672     #print " len(S0) = " + str(len(S0))
     627    verbose(" m = " + str(m) + "   p = " + str(p))
     628    verbose(" S0 = " + str(S0))
     629    verbose(" len(S0) = " + str(len(S0)))
    673630
    674631
    675     ## Make the form Qnew for the reduction procedure
    676     if ((m % p == 0) and (len(Set(S0New).intersection(Set(NZvec))) == 0) and (len(S0New) != n)):
    677         #print "Test 1.1 "
    678         Qnew = deepcopy(self)
    679         #Qnew = QuadraticForm(self.base_ring(), self.dim(), self.entries)
    680632
     633    ## Make the form Qnew for the reduction procedure:
     634    ## -----------------------------------------------
     635    Qnew = deepcopy(self)                                          ## TO DO:  DO THIS WITHOUT A copy(). =)
     636    for i in range(n):
     637        if i in S0:
     638            Qnew[i,i] = p * Qnew[i,i]
     639            if ((p == 2) and (i < n-1)):
     640                Qnew[i,i+1] = p * Qnew[i,i+1]
     641        else:
     642            Qnew[i,i] = Qnew[i,i] / p
     643            if ((p == 2) and (i < n-1)):
     644                Qnew[i,i+1] = Qnew[i,i+1] / p
    681645
    682         ## Run over all indices with i, while keeping track of the next index in S0 as S0[j].
    683         j = 0       ## The index for following the increasing set (of indices) S0
    684         for i in range(n):
    685             ## Note: Short circuit && is necessary here:
    686             if ((j < len(S0New)) and (i == S0New[j])):                 ## i is in S0
    687                 j += 1
    688                 Qnew[i,i] = p * Qnew[i,i]
    689                 if ((p == 2) and (i < n)):
    690                     Qnew[i+1,i] = p * Qnew[i+1,i]
    691                     Qnew[i,i+1] = p * Qnew[i,i+1]
    692             else:                             ## i not in S0;  Since S0 is increasing, we know i < S0(j)
    693                 Qnew[i,i] = Qnew[i,i] / p
    694646
    695                 #print "  dividing in row " + str(i)
    696647
    697                 if ((p == 2) and (i < n-1)):
    698                     Qnew[i+1,i] = Qnew[i+1,i] / p
    699                     Qnew[i,i+1] = Qnew[i,i+1] / p
     648    ## DIAGNOSTIC   
     649    verbose("\n\n Check of Bad-type I reduction: \n")
     650    verbose(" Q is " + str(self))
     651    verbose(" Qnew is " + str(Qnew))
     652    verbose(" p = " + str(p))
     653    verbose(" m / p = " + str(m/p))
     654    verbose(" NZvec " + str(NZvec))
    700655
    701         ## DIAGNOSTIC   
    702         #print "\n\n Check of Bad-type I reduction: \n";
    703         #print " Q is " + str(self)
    704         #print " Qnew is " + str(Qnew)
    705         #print " p = " + str(p)
    706         #print " m / p = " + str(m/p)
    707         #print " VectorComplement(Zvec, S0) is " + str(VectorComplement(Zvec, S0))       <------  FIX THIS!
    708         #print " NZvec " << + str(NZvec)
    709656
    710657   
    711         ## Do the reduction
    712         ## (Need 2 steps since we can't take negative powers... =| )
    713         VC = list(Set([i  for i in Zvec  if i not in S0New]))
    714         if (len(S0New) > 1):
    715             return QQ(1 / p**(len(S0New) - 1)) * Qnew.local_good_density_congruence(p, m / p, VC, NZvec)
    716         else:
    717             return QQ(p**(1 - len(S0New))) * Qnew.local_good_density_congruence(p, m / p, VC, NZvec)
     658    ## Do the reduction
     659    Zvec_geq_1 = list(Set([i  for i in Zvec  if i not in S0]))
     660    if NZvec == None:
     661        NZvec_geq_1 = NZvec
     662    else:
     663        NZvec_geq_1 = list(Set([i  for i in NZvec  if i not in S0]))
    718664
    719     else:
    720         return 0
     665    return QQ(p**(1 - len(S0))) * Qnew.local_good_density_congruence(p, m / p, Zvec_geq_1, NZvec_geq_1)
    721666
    722667
    723668
    724 def local_badII_density_congruence(self, p, m, Zvec, NZvec):
     669
     670
     671
     672
     673def local_badII_density_congruence(self, p, m, Zvec=None, NZvec=None):
    725674    """
    726675    Finds the Bad-type II local density of Q representing m at p.
    727676    (Assuming that p > 2 and Q is given in local diagonal form.)
    728677
    729     mpq_class Matrix_mpz::local_badII_density_congruence(const mpz_class & p, const mpz_class & m,
    730                                         const valarray<size_t> & Zvec, const valarray<size_t> & NZvec) const
     678
     679     INPUT:
     680        Q -- quadratic form assumed to be block diagonal and p-integral
     681        p -- a prime number
     682        m -- an integer
     683        Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
     684
     685    OUTPUT:
     686        a rational number
     687
     688    EXAMPLES:
     689        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
     690        sage: Q.local_badII_density_congruence(2, 1, None, None)
     691        0
     692        sage: Q.local_badII_density_congruence(2, 2, None, None)
     693        0
     694        sage: Q.local_badII_density_congruence(2, 4, None, None)
     695        0
     696        sage: Q.local_badII_density_congruence(3, 1, None, None)
     697        0
     698        sage: Q.local_badII_density_congruence(3, 6, None, None)
     699        0
     700        sage: Q.local_badII_density_congruence(3, 9, None, None)
     701        0
     702        sage: Q.local_badII_density_congruence(3, 27, None, None)
     703        0
     704
     705        sage: Q = DiagonalQuadraticForm(ZZ, [1,3,3,9,9])
     706        sage: Q.local_badII_density_congruence(3, 1, None, None)
     707        0
     708        sage: Q.local_badII_density_congruence(3, 3, None, None)
     709        0
     710        sage: Q.local_badII_density_congruence(3, 6, None, None)
     711        0
     712        sage: Q.local_badII_density_congruence(3, 9, None, None)
     713        4/27
     714        sage: Q.local_badII_density_congruence(3, 18, None, None)
     715        4/9
     716
     717
    731718    """
    732719    ## DIAGNOSTIC
    733     #print " In local_badII_density_congruence with "
    734     #print " Q is: \n" + str(self)
    735     #print " p = " + str(p)
    736     #print " m = " + str(m)
    737     #print " Zvec = " + str(Zvec)
    738     #print " NZvec = " + str(NZvec)
     720    verbose(" In local_badII_density_congruence with ")
     721    verbose(" Q is: \n" + str(self))
     722    verbose(" p = " + str(p))
     723    verbose(" m = " + str(m))
     724    verbose(" Zvec = " + str(Zvec))
     725    verbose(" NZvec = " + str(NZvec))
     726
     727    ## Put the Zvec congruence condition in a standard form
     728    if Zvec == None:
     729        Zvec = []
    739730
    740731
    741732    n = self.dim()
    742733
    743     ## Define the indexing sets S_i
     734
     735    ## Sanity Check on Zvec and NZvec:
     736    ## -------------------------------
     737    Sn = Set(range(n))
     738    if (Zvec != None) and (len(Set(Zvec) + Sn) > n):
     739        raise RuntimeError, "Zvec must be a subset of {0, ..., n-1}."
     740    if (NZvec != None) and (len(Set(NZvec) + Sn) > n):
     741        raise RuntimeError, "NZvec must be a subset of {0, ..., n-1}."
     742
     743
     744    ## Define the indexing sets S_i:
     745    ## -----------------------------
    744746    S0 = []
    745747    S1 = []
    746     S2 = []
     748    S2plus = []
    747749
    748750    for i in range(n):
    749751
     
    764766        elif (val == 1):
    765767            S1 += [i]
    766768        elif (val >= 2):
    767             S2 += [i]
     769            S2plus += [i]
    768770
    769771
     772     
     773    ## Check that S2 is non-empty and p^2 divides m to proceed, otherwise return no solutions.
     774    p2 = p * p
     775    if (S2plus == []) or (m % p2 != 0):
     776        return 0
    770777
    771     ## Check that the form is primitive...
     778
     779
     780
     781    ## Check some conditions for no bad-type II solutions to exist
     782    if (NZvec != None) and (len(Set(S2plus).intersection(Set(NZvec))) == 0):
     783        return 0
     784
     785
     786
     787    ## Check that the form is primitive...                     WHY IS THIS NECESSARY?
    772788    if (S0 == []):
    773789        print " Using Q = " + str(self)
    774790        print " and p = " + str(p)
    775791        raise RuntimeError, "Oops! The form is not primitive!"
    776792
    777793
    778     ## DIAGNOSTIC
    779     #print "\n Entering BII routine "
    780     #print " S0 is " + str(S0)
    781     #print " S1 is " + str(S1)
    782     #print " S2 is " + str(S2)
    783794
    784795
    785 
    786     ## Note: The following lines assume that NZvec is an ordered vector.  Should check this here...
    787 
    788796    ## DIAGNOSTIC
    789     #print " m = " + str(m) + "   p = " + str(p)
    790     #print " S0 = " + str(S0) + "   NZvec = " str(NZvec) "   IsDisjoint = " + str(len(Set(S0 + S1).intersection(Set(NZvec))) == 0)
    791     #print " len(S0) = " + str(len(S0))
     797    verbose("\n Entering BII routine ")
     798    verbose(" S0 is " + str(S0))
     799    verbose(" S1 is " + str(S1))
     800    verbose(" S2plus is " + str(S2plus))
     801    verbose(" m = " + str(m) + "   p = " + str(p))
    792802
    793803
    794     p2 = p * p
    795804
    796     ## Make the form Qnew for the reduction procedure
    797     if ((m % (p2) == 0) and (S2 != []) and (len(Set(S0 + S1).intersection(Set(NZvec))) == 0)):   ##  <=====  CHECK THIS!!! ****  Should this be (S0 U S1) instead of S0 ???
    798         #print "Test 1.1 "
    799805
    800         Qnew = deepcopy(self)
    801         #Qnew = QuadraticForm(self.base_ring(), self.dim(), self.entries)
    802806
    803         ## Run over all indices with i, while keeping track of the next index in S2 as S2[j].
    804         j = 0       ## The index for following the increasing set (of indices) S2
    805         for i in range(n):
    806             ## Note: Short circuit && is necessary here:
    807             if ((j < len(S2)) and (i == S2[j])):                  ## i is in S2
    808                 j += 1
    809                 Qnew[i,i] = Qnew[i,i] / p2
    810                 if ((p == 2) and (i < n-1)):
    811                     Qnew[i+1,i] = Qnew[i+1,i] / p2
    812                     Qnew[i,i+1] = Qnew[i,i+1] / p2
     807    ## Make the form Qnew for the reduction procedure:
     808    ## -----------------------------------------------
     809    Qnew = deepcopy(self)                                          ## TO DO:  DO THIS WITHOUT A copy(). =)
     810    for i in range(n):
     811        if i in S2plus:
     812            Qnew[i,i] = Qnew[i,i] / p2
     813            if (p == 2) and (i < n-1):
     814                Qnew[i,i+1] = Qnew[i,i+1] / p2
    813815
    814         ## DIAGNOSTIC
    815         #print "\n\n Check of Bad-type II reduction: \n";
    816         #print " Q is " + str(self)
    817         #print " Qnew is " + str(Qnew)
     816    ## DIAGNOSTIC
     817    verbose("\n\n Check of Bad-type II reduction: \n")
     818    verbose(" Q is " + str(self))
     819    verbose(" Qnew is " + str(Qnew))
    818820
    819821
    820         ## Do the reduction
    821         ## (Need 2 steps for each case since we can't take negative powers... =| )
    822         ## ------------------------------------------------------------------------
    823         ##new_Zvec = VectorComplement(Zvec, list(Set(S0 + S1))) 
    824         new_Zvec = list(Set([i  for i in Zvec  if i not in S0 + S1]))
    825822
     823    ## Perform the reduction formula
     824    Zvec_geq_2 = list(Set([i  for i in Zvec  if i in S2plus]))
     825    if NZvec == None:
     826        NZvec_geq_2 = NZvec
     827    else:
     828        NZvec_geq_2 = list(Set([i  for i in NZvec  if i in S2plus]))
    826829
    827         ## DIAGNOSTIC
    828         #print "  m = " << m << ",  m/p2 = " + str(m / p2)
    829         #print "  new_Zvec = " + str(new_Zvec)
    830         #print "  NZvec = " + str(NZvec)
    831         #print "  local_density_congruence(Qnew, p, m / p2, new_Zvec, NZvec) = " \
    832         #   + str(local_density_congruence(Qnew, p, m / p2, new_Zvec, NZvec))
    833         #print "  local_density_congruence(Qnew, p, m / p2, List(Set(S2 + new_Zvec)), NZvec)) = " \
    834         #   + str(local_density_congruence(Qnew, p, m / p2, List(Set(S2 + new_Zvec)), NZvec))
    835         #print "  count_local_type(Qnew, 2, 2, m=1, 0, trimZvec, trimNZvec) = " \
    836         #   + str(count_local_type(Qnew, 2, 2, 1, 0, Zvec, NZvec))
    837         #print "  count_local_type(Qnew, 2, 3, m=1, 0, trimZvec, trimNZvec) = " \
    838         #   + str(count_local_type(Qnew, 2, 3, 1, 0, Zvec, NZvec))
    839         #print "  count_local_type(Qnew, 2, 4, m=1, 0, trimZvec, trimNZvec) = " \
    840         #   + str(count_local_type(Qnew, 2, 4, 1, 0, Zvec, NZvec))
    841         #print "  count_local_type(Qnew, 2, 5, m=1, 0, trimZvec, trimNZvec) = " \
    842         #   + str(count_local_type(Qnew, 2, 5, 1, 0, Zvec, NZvec))
    843         #print "  count_local_type(Q, 2, 2, m=2, 0, trimZvec, trimNZvec) = " \
    844         #   + str(count_local_type(Q, 2, 2, 2, 0, Zvec, NZvec))
    845         #print "  count_local_type(Q, 2, 3, m=2, 0, trimZvec, trimNZvec) = " \
    846         #   + str(count_local_type(Q, 2, 3, 2, 0, Zvec, NZvec))
    847         #print "  count_local_type(Q, 2, 4, m=2, 0, trimZvec, trimNZvec) = " \
    848         #   + str(count_local_type(Q, 2, 4, 2, 0, Zvec, NZvec))
    849         #print "  count_local_type(Q, 2, 5, m=2, 0, trimZvec, trimNZvec) = " \
    850         #   + str(count_local_type(Q, 2, 5, 2, 0, Zvec, NZvec))
     830    return QQ(p**(len(S2plus) + 2 - n)) \
     831        * (Qnew.local_density_congruence(p, m / p2, Zvec_geq_2, NZvec_geq_2) \
     832        - Qnew.local_density_congruence(p, m / p2, S2plus , NZvec_geq_2))
    851833
    852834
    853         if (n > len(S2) + 2):
    854             return QQ(1 / p**(n - len(S2) - 2)) \
    855                 * (Qnew.local_density_congruence(p, m / p2, new_Zvec, NZvec) \
    856                 - Qnew.local_density_congruence(p, m / p2, List(Set(S2 + new_Zvec)), NZvec))
    857         else:
    858             return QQ(p**(len(S2) + 2 - n)) \
    859                 * (Qnew.local_density_congruence(p, m / p2, new_Zvec, NZvec) \
    860                 - Qnew.local_density_congruence(p, m / p2, List(Set(S2 + new_Zvec)), NZvec))
    861835
    862     else:
    863         return 0
    864836
    865837
    866838
    867 
    868 
    869 def local_bad_density_congruence(self, p, m, Zvec, NZvec):
     839def local_bad_density_congruence(self, p, m, Zvec=None, NZvec=None):
    870840    """
    871841    Finds the Bad-type local density of Q representing
    872842    m at p, allowing certain congruence conditions mod p.
     843
     844    INPUT:
     845        Q -- quadratic form assumed to be block diagonal and p-integral
     846        p -- a prime number
     847        m -- an integer
     848        Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
     849
     850    OUTPUT:
     851        a rational number
     852
     853    EXAMPLES:
     854        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
     855        sage: Q.local_bad_density_congruence(2, 1, None, None)
     856        0
     857        sage: Q.local_bad_density_congruence(2, 2, None, None)
     858        1
     859        sage: Q.local_bad_density_congruence(2, 4, None, None)
     860        0
     861        sage: Q.local_bad_density_congruence(3, 1, None, None)
     862        0
     863        sage: Q.local_bad_density_congruence(3, 6, None, None)
     864        0
     865        sage: Q.local_bad_density_congruence(3, 9, None, None)
     866        0
     867        sage: Q.local_bad_density_congruence(3, 27, None, None)
     868        0
     869
     870        sage: Q = DiagonalQuadraticForm(ZZ, [1,3,3,9,9])
     871        sage: Q.local_bad_density_congruence(3, 1, None, None)
     872        0
     873        sage: Q.local_bad_density_congruence(3, 3, None, None)
     874        4/3
     875        sage: Q.local_bad_density_congruence(3, 6, None, None)
     876        4/3
     877        sage: Q.local_bad_density_congruence(3, 9, None, None)
     878        4/27
     879        sage: Q.local_bad_density_congruence(3, 18, None, None)
     880        4/9
     881        sage: Q.local_bad_density_congruence(3, 27, None, None)
     882        8/27
     883
    873884    """
    874885    return self.local_badI_density_congruence(p, m, Zvec, NZvec) + self.local_badII_density_congruence(p, m, Zvec, NZvec)
    875886 
    876887
    877888
    878 
    879889#########################################################
    880890## local_density and local_density_congruence routines ##
    881891#########################################################
    882892
    883 def local_density_congruence(self, p, m, Zvec, NZvec):
     893def local_density_congruence(self, p, m, Zvec=None, NZvec=None):
    884894    """
    885895    Finds the local density of Q representing m at p,
    886896    allowing certain congruence conditions mod p.
    887897
     898    INPUT:
     899        Q -- quadratic form assumed to be block diagonal and p-integral
     900        p -- a prime number
     901        m -- an integer
     902        Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
     903
     904    OUTPUT:
     905        a rational number
     906
    888907    EXAMPLES:
    889908        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
    890         sage: Q.local_density_congruence(p=2, m=1, Zvec=[], NZvec=[])
     909        sage: Q.local_density_congruence(p=2, m=1, Zvec=None, NZvec=None)
    891910        1
    892         sage: Q.local_density_congruence(p=3, m=1, Zvec=[], NZvec=[])
     911        sage: Q.local_density_congruence(p=3, m=1, Zvec=None, NZvec=None)
    893912        8/9
    894         sage: Q.local_density_congruence(p=5, m=1, Zvec=[], NZvec=[])
     913        sage: Q.local_density_congruence(p=5, m=1, Zvec=None, NZvec=None)
    895914        24/25
    896         sage: Q.local_density_congruence(p=7, m=1, Zvec=[], NZvec=[])
     915        sage: Q.local_density_congruence(p=7, m=1, Zvec=None, NZvec=None)
    897916        48/49
    898         sage: Q.local_density_congruence(p=11, m=1, Zvec=[], NZvec=[])
     917        sage: Q.local_density_congruence(p=11, m=1, Zvec=None, NZvec=None)
    899918        120/121
    900919
     920        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
     921        sage: Q.local_density_congruence(2, 1, None, None)
     922        1
     923        sage: Q.local_density_congruence(2, 2, None, None)
     924        1
     925        sage: Q.local_density_congruence(2, 4, None, None)
     926        3/2
     927        sage: Q.local_density_congruence(3, 1, None, None)
     928        2/3
     929        sage: Q.local_density_congruence(3, 6, None, None)
     930        4/3
     931        sage: Q.local_density_congruence(3, 9, None, None)
     932        14/9
     933        sage: Q.local_density_congruence(3, 27, None, None)
     934        2
     935
     936        sage: Q = DiagonalQuadraticForm(ZZ, [1,3,3,9,9])
     937        sage: Q.local_density_congruence(3, 1, None, None)
     938        2
     939        sage: Q.local_density_congruence(3, 3, None, None)
     940        4/3
     941        sage: Q.local_density_congruence(3, 6, None, None)
     942        4/3
     943        sage: Q.local_density_congruence(3, 9, None, None)
     944        2/9
     945        sage: Q.local_density_congruence(3, 18, None, None)
     946        4/9
     947
    901948    """
    902949    return self.local_good_density_congruence(p, m, Zvec, NZvec) \
    903950                + self.local_zero_density_congruence(p, m, Zvec, NZvec) \
     
    905952
    906953
    907954
    908 def local_primitive_density_congruence(self, p, m, Zvec, NZvec):
     955def local_primitive_density_congruence(self, p, m, Zvec=None, NZvec=None):
    909956    """
    910957    Finds the primitive local density of Q representing
    911958    m at p, allowing certain congruence conditions mod p.
    912959
    913960    Note: The following routine is not used internally, but is included for consistency.
     961
     962    INPUT:
     963        Q -- quadratic form assumed to be block diagonal and p-integral
     964        p -- a prime number
     965        m -- an integer
     966        Zvec, NZvec -- non-repeating lists of integers in range(self.dim()) or None
     967
     968    OUTPUT:
     969        a rational number
     970
     971    EXAMPLES:
     972        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])
     973        sage: Q.local_primitive_density_congruence(p=2, m=1, Zvec=None, NZvec=None)
     974        1
     975        sage: Q.local_primitive_density_congruence(p=3, m=1, Zvec=None, NZvec=None)
     976        8/9
     977        sage: Q.local_primitive_density_congruence(p=5, m=1, Zvec=None, NZvec=None)
     978        24/25
     979        sage: Q.local_primitive_density_congruence(p=7, m=1, Zvec=None, NZvec=None)
     980        48/49
     981        sage: Q.local_primitive_density_congruence(p=11, m=1, Zvec=None, NZvec=None)
     982        120/121
     983
     984        sage: Q = DiagonalQuadraticForm(ZZ, [1,2,3])
     985        sage: Q.local_primitive_density_congruence(2, 1, None, None)
     986        1
     987        sage: Q.local_primitive_density_congruence(2, 2, None, None)
     988        1
     989        sage: Q.local_primitive_density_congruence(2, 4, None, None)
     990        1
     991        sage: Q.local_primitive_density_congruence(3, 1, None, None)
     992        2/3
     993        sage: Q.local_primitive_density_congruence(3, 6, None, None)
     994        4/3
     995        sage: Q.local_primitive_density_congruence(3, 9, None, None)
     996        4/3
     997        sage: Q.local_primitive_density_congruence(3, 27, None, None)
     998        4/3
     999
     1000        sage: Q = DiagonalQuadraticForm(ZZ, [1,3,3,9,9])
     1001        sage: Q.local_primitive_density_congruence(3, 1, None, None)
     1002        2
     1003        sage: Q.local_primitive_density_congruence(3, 3, None, None)
     1004        4/3
     1005        sage: Q.local_primitive_density_congruence(3, 6, None, None)
     1006        4/3
     1007        sage: Q.local_primitive_density_congruence(3, 9, None, None)
     1008        4/27
     1009        sage: Q.local_primitive_density_congruence(3, 18, None, None)
     1010        4/9
     1011        sage: Q.local_primitive_density_congruence(3, 27, None, None)
     1012        8/27
     1013        sage: Q.local_primitive_density_congruence(3, 81, None, None)
     1014        8/27
     1015        sage: Q.local_primitive_density_congruence(3, 243, None, None)
     1016        8/27
     1017
    9141018    """
    9151019    return self.local_good_density_congruence(p, m, Zvec, NZvec) \
    9161020                + self.local_bad_density_congruence(p, m, Zvec, NZvec)
  • sage/quadratic_forms/quadratic_form__local_density_interfaces.py

    diff -r af0f19cc490e -r b3d0bbb703ab sage/quadratic_forms/quadratic_form__local_density_interfaces.py
    a b  
    88from sage.rings.rational_field import QQ, RationalField
    99
    1010
    11 #  ////////////////////////////////
    12 #  // Private Front-end Routines //
    13 #/////////////////////////////////////////////////////////////////////////////////////////////
    1411
    15 def local_good_density(self, p, m):
    16     """
    17     Finds the Good-type local density of Q representing m at p. 
    18     (Front end routine for its congruence counterpart.)         
    19     """
    20     #print "Doing Good Density with p = " + str(p) +  " and m = " + str(m)
    21     return self.local_good_density_congruence(p, m, [], [])
    22 
    23 
    24 def local_zero_density(self, p, m):
    25     """
    26     Finds the Zero-type local density of Q representing m at p. 
    27     (Front end routine for its congruence counterpart.)         
    28     """
    29     #print "Doing Good Density with p = " + str(p) +  " and m = " + str(m)
    30     return self.local_zero_density_congruence(p, m, [], [])
    31 
    32 
    33 def local_bad_density(self, p, m):
    34     """
    35     Finds the Bad-type local density of Q representing m at p. 
    36     (Front end routine for its congruence counterpart.)         
    37     """
    38     #print "Doing Bad Density with p = " + str(p) +  " and m = " + str(m)
    39     return self.local_bad_density_congruence(p, m, [], [])
    40 
    41 
    42 def local_badI_density(self, p, m):
    43     """
    44     Finds the Bad-type I local density of Q representing m at p. 
    45     (Front end routine for its congruence counterpart.)         
    46     """
    47     #print "Doing Bad I Density with p = " + str(p) +  " and m = " + str(m)
    48     return self.local_badI_density_congruence(p, m, [], [])
    49 
    50 
    51 def local_badII_density(self, p, m):
    52     """
    53     Finds the Bad-type II local density of Q representing m at p. 
    54     (Front end routine for its congruence counterpart.)         
    55     """
    56     #print "Doing Bad II Density with p = " + str(p) +  " and m = " + str(m)
    57     return self.local_badII_density_congruence(p, m, [], [])
    58 
    59 
    60 ## ---------------  These are the important ones, which we'll filter for primitive forms!!!  ------------------
    61 
    62 
    63 ### TODO:  THESE TWO ROUTINES HAVE **A LOT** OF CODE IN COMMON, BUT NEITHER PUTS THE FORM IN LOCAL NORMAL FORM FIRST!
    6412
    6513
    6614def local_density(self, p, m):
    6715    """
    6816    Gives the local density -- should be called by the user. =)
    6917
    70     NOTE: This screens for imprimitive forms, but *doesn't* put the
    71     quadratic form in local normal form, which is a *requirement* of
    72     the routines performing the computations!
     18    NOTE: This screens for imprimitive forms, and puts the quadratic
     19    form in local normal form, which is a *requirement* of the
     20    routines performing the computations!
    7321
    74     mpq_class Matrix_mpz::local_density(const mpz_class & p, const mpz_class & m) const {
     22    INPUT:
     23        p -- a prime number > 0
     24        m -- an integer
     25
     26    OUTPUT:
     27        a rational number
    7528
    7629    EXAMPLES:
    7730        sage: Q = DiagonalQuadraticForm(ZZ, [1,1,1,1])   ## NOTE: This is already in local normal form for *all* primes p!
     
    9144    if (n == 0):
    9245        raise TypeError, "Oops!  We currently don't handle 0-dim'l forms. =("
    9346
    94 
    95     ## Set the modulus to check for imprimitive forms at p
    96     no_val_flag = True
    97 
    98     ## Check for imprimitive forms at p -- ASSUMES THE FORM IS UPPER TRIANGULAR!
    99     ## (NOTE: We could do better if we know it's normalized!)
    100     for i in range(n):
    101         for j in range(i, n):
    102             if (self[i,j] != 0):
    103                 if (no_val_flag == True):
    104                     no_val_flag = False
    105                     p_valuation = valuation(self[i,j], p)
    106                 else:
    107                     p_valuation = min(p_valuation, valuation(self[i,j], p))
    108 
    109 
    110     ## DIAGNOSTIC
    111     #cout << " Using the matrix: \n " << (*this) << endl;
    112     #cout << "Valuation(m,p) = " << Valuation(m,p) << endl;
    113     #cout << "p_valuation = " << p_valuation << endl;
    114 
     47    ## Find the local normal form and p-scale of Q     --  Note: This uses the valuation ordering of local_normal_form.
     48    ##                                                     TO DO:  Write a separate p-scale and p-norm routines!
     49    Q_local = self.local_normal_form(p)
     50    if n == 1:                 
     51        p_valuation = valuation(Q_local[0,0], p)
     52    else:
     53        p_valuation = min(valuation(Q_local[0,0], p), valuation(Q_local[0,1], p))
    11554
    11655    ## If m is less p-divisible than the matrix, return zero 
    11756    if ((m != 0) and (valuation(m,p) < p_valuation)):   ## Note: The (m != 0) condition protects taking the valuation of zero.
    11857        return QQ(0)
    119  
    120     ## If the form is imprimitive, divide it (and m) by p-powers to get a primitive form
    121     else:
    122         if (p_valuation > 0):
    123      
    124             ## Make a new (primitive) matrix
    125             Q1 = deepcopy(self)
    126             #Q1 = QuadraticForm(self.base_ring(), self.dim())
    12758
    128             ## DIAGNOSTIC
    129             #print " p = " << p
    130             #print " p_valuation = " << p_valuation
    131             #print " p_mod = " << (p ** p_valuation)
    132      
    133             p_mod = p ** p_valuation      ## This should give a power...
    134             for i in range(n):
    135                 for j in range(i, n):
    136                     Q1[i,j] = self[i,j] / p_mod
    137      
    138             ## Make a new number mm
    139             mm = m / p_mod
    140      
    141             ## Then return the densities for the reduced problem
    142             return Q1.local_good_density(p, mm) + Q1.local_zero_density(p, mm) + Q1.local_bad_density(p, mm)
    14359
    144    
    145         ## Otherwise, proceed as usual... =)
    146         else:
    147             return self.local_good_density(p, m) + self.local_zero_density(p, m) + self.local_bad_density(p, m)
     60    ## If the form is imprimitive, rescale it and call the local density routine
     61    p_adjustment = QQ(1) / p**p_valuation
     62    m_prim = QQ(m) / p**p_valuation
     63    Q_prim = Q_local.scale_by_factor(p_adjustment)
     64
     65    ## Return the densities for the reduced problem
     66    return Q_prim.local_density_congruence(p, m_prim)
    14867
    14968
    15069
     
    15372    """
    15473    Gives the local primitive density -- should be called by the user. =)
    15574
    156     NOTE: This screens for imprimitive forms, but *doesn't* put the
     75    NOTE: This screens for imprimitive forms, and puts the
    15776    quadratic form in local normal form, which is a *requirement* of
    15877    the routines performing the computations!
    15978
    160     mpq_class Matrix_mpz::local_density(const mpz_class & p, const mpz_class & m) const {
     79    INPUT:
     80        p -- a prime number > 0
     81        m -- an integer
     82
     83    OUTPUT:
     84        a rational number
     85
     86    EXAMPLES:
     87        sage: Q = QuadraticForm(ZZ, 4, range(10))
     88        sage: Q[0,0] = 5
     89        sage: Q[1,1] = 10
     90        sage: Q[2,2] = 15
     91        sage: Q[3,3] = 20
     92        sage: Q
     93        Quadratic form in 4 variables over Integer Ring with coefficients:
     94        [ 5 1 2 3 ]
     95        [ * 10 5 6 ]
     96        [ * * 15 8 ]
     97        [ * * * 20 ]
     98        sage: Q.theta_series(20)
     99        1 + 2*q^5 + 2*q^10 + 2*q^14 + 2*q^15 + 2*q^16 + 2*q^18 + O(q^20)
     100        sage: Q.local_normal_form(2)
     101        Quadratic form in 4 variables over Integer Ring with coefficients:
     102        [ 0 1 0 0 ]
     103        [ * 0 0 0 ]
     104        [ * * 0 1 ]
     105        [ * * * 0 ]
     106       
     107        sage: Q.local_primitive_density(2, 1)
     108        3/4
     109        sage: Q.local_primitive_density(5, 1)
     110        24/25
     111
     112        sage: Q.local_primitive_density(2, 5)
     113        3/4
     114        sage: Q.local_density(2, 5)
     115        3/4
     116
    161117    """
    162 
    163118    n = self.dim()
    164119    if (n == 0):
    165120        raise TypeError, "Oops!  We currently don't handle 0-dim'l forms. =("
    166121
    167 
    168     ## Set the modulus to check for imprimitive forms at p
    169     no_val_flag = True
    170 
    171     ## Check for imprimitive forms at p -- ASSUMES THE FORM IS UPPER TRIANGULAR!
    172     ## (NOTE: We could do better if we know it's normalized!)
    173     for i in range(n):
    174         for j in range(i, n):
    175             if (self[i,j] != 0):
    176                 if (no_val_flag == True):
    177                     no_val_flag = False
    178                     p_valuation = valuation(self[i,j], p)
    179                 else:
    180                     p_valuation = min(p_valuation, valuation(self[i,j], p))
    181 
    182 
    183     ## DIAGNOSTIC
    184     #cout << " Using the matrix: \n " << (*this) << endl;
    185     #cout << "Valuation(m,p) = " << Valuation(m,p) << endl;
    186     #cout << "p_valuation = " << p_valuation << endl;
     122    ## Find the local normal form and p-scale of Q     --  Note: This uses the valuation ordering of local_normal_form.
     123    ##                                                     TO DO:  Write a separate p-scale and p-norm routines!
     124    Q_local = self.local_normal_form(p)
     125    if n == 1:                 
     126        p_valuation = valuation(Q_local[0,0], p)
     127    else:
     128        p_valuation = min(valuation(Q_local[0,0], p), valuation(Q_local[0,1], p))
    187129
    188130
    189131    ## If m is less p-divisible than the matrix, return zero 
    190132    if ((m != 0) and (valuation(m,p) < p_valuation)):   ## Note: The (m != 0) condition protects taking the valuation of zero.
    191133        return QQ(0)
    192  
    193     ## If the form is imprimitive, divide it (and m) by p-powers to get a primitive form
    194     else:
    195         if (p_valuation > 0):
    196      
    197             ## Make a new (primitive) matrix
    198             Q1 = deepcopy(self) 
    199             #Q1 = QuadraticForm(self.base_ring(), self.dim())
    200134
    201             ## DIAGNOSTIC
    202             #print " p = " << p
    203             #print " p_valuation = " << p_valuation
    204             #print " p_mod = " << (p ** p_valuation)
    205      
    206             p_mod = p ** p_valuation      ## This should give a power...
    207             for i in range(n):
    208                 for j in range(i, n):
    209                     Q1[i,j] = self[i,j] / p_mod
    210      
    211             ## Make a new number mm
    212             mm = m / p_mod
    213      
    214             ## Then return the densities for the reduced problem
    215             return Q1.local_good_density(p, mm) + Q1.local_bad_density(p, mm)
    216135
    217    
    218         ## Otherwise, proceed as usual... =)
    219         else:
    220             return self.local_good_density(p, m) + self.local_bad_density(p, m)
     136    ## If the form is imprimitive, rescale it and call the local density routine
     137    p_adjustment = QQ(1) / p**p_valuation
     138    m_prim = QQ(m) / p**p_valuation
     139    Q_prim = Q_local.scale_by_factor(p_adjustment)
    221140
     141    ## Return the densities for the reduced problem
     142    return Q_prim.local_primitive_density_congruence(p, m_prim)
     143