/*******************************************************************************
* A demo file to illustrate p-adic and 2-adic point counting on a surface of   *
* the shape                                                                    *
*                                                                              *
*     w^2 = x y z f_3(x,y,z) .                                                 *
*                                                                              *
* We use a Harvey style algorithm with a remainder tree and a moving simplex   *
* for the p-adic part.                                                         *
*                                                                              *
* The 2-adic part is initialised by naive point counting.                      *
*                                                                              *
* The code assumes that f_3 is a product of three integral linear forms.       *
*******************************************************************************/


/*******************************************************************************
*                                                                              *
* p-adic point counting. The Harvey style algorithm.                           *
*                                                                              *
*******************************************************************************/

// The example we work with.
r3<x,y,z> := PolynomialRing(Rationals(), 3);
lfl := [x, y, z, (x+y+z), (3*x+5*y+7*z), (-5*x + 11*y - 2*z)]; // Generic example
f6 := &*lfl;
f3 := f6 div (x*y*z);

/* The search bound */
lim := 10^5;

/* Stop magma before the machine runs out of memory. */
SetMemoryLimit(1*10^9);


// Lists to store the p-adic results
prime_list := [];
count_mod_p_list := [];

// Use (d/dx_i) (f^{n+1}) = (n+1) f^n (d/dx_i) f  to create relations between the coefficients.
/* Index ordering
   1st: offsets_2 of f^{n+1}
   2nd: offsets_3 \ offsets_2 of f^n
   3rd: offsets_2 of f^n
We have to treat n as a variable!
As everything is linear in n, we have a constant and a linear term. Each is represented by a matrix. */
function init_mov_simplex(f3)
 offsets_2 := [ [e1,e2,e3] : e1,e2,e3 in [-2..10] | e1+e2+e3 eq 0 and {e1,e2,e3} ne {0}  ];
 Append(~offsets_2,[0,0,0]);
 offsets_3 := [ [e1,e2,e3] : e1,e2,e3 in [-5..16] | e1+e2+e3 eq 0 and not [e1,e2,e3] in offsets_2 ];

 r_seq := [];
 qq := PolynomialRing(Rationals());
 n := qq.1;
 for i := 1 to 3 do
  col := Coefficients(Derivative(-f3, i));
  exl := [Exponents(a) : a in Monomials(Derivative(f3, i))];
  for j := 1 to #offsets_2 do
   act := [qq!0 : i in [1..2*#offsets_2 + #offsets_3]];
   ex := [(n+1) + offsets_2[j][1], (n+1) + offsets_2[j][2], (n+1) + offsets_2[j][3]];
   mul := ex[i];
   ex[i] := ex[i] - 1;
   act[j] := - mul;
   for k := 1 to #exl do
    exm := [ex[1] - exl[k][1] - n, ex[2] - exl[k][2] - n, ex[3] - exl[k][3] - n];
    ii := Index(offsets_2, exm);
    if ii eq 0 then
     ii := Index(offsets_3, exm);
     if ii eq 0 then continue j; end if; // Coefficient out of range. We can not use this relation.
     act[ii + #offsets_2] := (n+1) * col[k];
    else
     act[ii + #offsets_2 + #offsets_3] := (n+1) * col[k];
    end if;
   end for;
   Append(~r_seq, act);
  end for;
 end for;

 ff<u> := FunctionField(Rationals());
 mat0 := Matrix(ff, [[Coefficient(a,0) : a in b ] : b in r_seq]);
 mat1 := Matrix(ff, [[Coefficient(a,1) : a in b ] : b in r_seq]);

 time m2 := EchelonForm(mat0 + u * mat1);

 assert &and [m2[i,i] eq 1 : i in [1..#offsets_2]];
 assert &and [m2[i,j] eq 0 : i in [1..#offsets_2], j in [1..#offsets_2 + #offsets_3] | i ne j];

 r2_gen := Transpose(Submatrix(m2, 1, 1 + #offsets_2 + #offsets_3, #offsets_2, #offsets_2));
 r2_gen := RowSequence(r2_gen);

 ind := [i : i in [1..#offsets_2] | offsets_2[i] eq [0,0,0] or {r2_gen[i][j] : j in [1..#offsets_2] | j ne i} ne {0}];
 offsets_2 := offsets_2[ind];
 r2_gen := [r2_gen[i][ind] : i in ind];
 i_center := Index(offsets_2, [0,0,0]);

 return r2_gen, offsets_2, i_center;
end function;


// Initial coefficient vector containing the elements of the coefficient triangle (simplex) for f3^ex.
function init_iteration(f3, ex, r2_gen, offsets_2)
 f6 := f3^ex;
 co_act := Vector(Rationals(), [Min(a) lt -ex select 0 else MonomialCoefficient(f6, Monomial(Parent(f6), [ex+k : k in a]))
           : a in offsets_2]);
 init_bound := 1+2*Max(PrimeDivisors(LCM([2] cat [Denominator(a) : a in Coefficients(Numerator(b)), b in &cat r2_gen])));
 return co_act, init_bound;
end function;


/* Apply the iteration directly to compute successive coefficient triangles.
   Work with all exponents up to init_bound.
   Return the coefficient triangle and the corresponding exponent.
   As we work in rational arithmetic, there is no p-adic precision loss. */
function direct_iteration(co_act, ex, init_bound, r2_gen, i_center)
 prime_list := []; count_mod_p_list := [];

 while (ex le init_bound) do
  if IsPrime(2*ex+1) then
   Append(~prime_list, 2*ex+1);
   Append(~count_mod_p_list, Integers()!co_act[i_center] mod (2*ex+1));
  end if;
  r2 := Matrix([[Evaluate(f, ex) : f in r] : r in r2_gen]);
  co_act := co_act * r2;
  ex := ex + 1;
 end while;

 if IsPrime(2*ex+1) then
  Append(~prime_list, 2*ex+1);
  Append(~count_mod_p_list, Integers()!co_act[i_center] mod (2*ex+1));
 end if;

 return co_act, ex, prime_list, count_mod_p_list;
end function;


/* Build a binary tree of all the moduli for the next 2^level_count primes.
   The first one is the next prime after 2*ex+1. If 2*ex+1 is prime, it was already treated earlier. */
function modulus_tree_build(ex, level_count)
 cnt := 2^(level_count - 1);
 l1 := PrimesInInterval(2*ex+2, 2*ex+1+Round(cnt*Log(2*ex+2)));

 if #l1 lt cnt then
  ende := l1[#l1];
  l1 := l1 cat PrimesInInterval(ende+1, ende+Round(Log(ende)*(cnt-#l1)));
 end if;

 while #l1 lt cnt do
  Append(~l1, NextPrime(l1[#l1]));
 end while;

 tree := [l1[[1..cnt]]];
 while #tree[#tree] gt 2 do
  Append(~tree, [tree[#tree][i]*tree[#tree][i+1] : i in [2..#(tree[#tree])-1] | i mod 2 eq 0]);
 end while;
 d := #tree;
 Append(~tree, [tree[1][#tree[1]]]);

 return Reverse(tree);
end function;


/* Build the product tree of all matrices that match the modulus tree to continue the interation with starting exponent ex.

This function returns a tree of integer matrices and a second tree of denominators. The trees are stored as a sequence of sequences. The last sequence contains the leaves of the tree. These contain the matrices that describe the update data of the simplex to go from one prime to the next one. The sequence before contains products of consecutive pairs M_(2i-1) * M_(2i) in the level below. The top sequence contains only one matrix. It is the product of all the matices in the leaves. */
function matrix_tree_build(ex, mat_gen_list, gen_denominator_l, mod_tree)
 l1 := [];
 l2 := [];
 pzl := mod_tree[#mod_tree];

 for i := 1 to #pzl do
  ex_target := pzl[i] div 2;
  step := ex_target - ex;
  // i, pzl[i], step;
  assert step le #mat_gen_list; assert step gt 0;
  sp_z := Matrix([[Evaluate(a, ex) : a in b] : b in mat_gen_list[step]]);
  sp_n := Evaluate(gen_denominator_l[step], ex);
  dd := GCD([sp_n] cat Eltseq(sp_z));
  sp_z := sp_z div dd;
  sp_n := sp_n div dd;
  ex := ex_target;
  Append(~l1, sp_z);
  Append(~l2, sp_n);
 end for;

 mat_tree := [l1];
 denominator_tree := [l2];
 while #mat_tree[#mat_tree] gt 1 do
  d := #mat_tree;
  Append(~mat_tree, [mat_tree[d][i] * mat_tree[d][i+1] : i in [1..#mat_tree[d]] | i mod 2 eq 1]);
  Append(~denominator_tree, [denominator_tree[d][i] * denominator_tree[d][i+1] : i in [1..#denominator_tree[d]] | i mod 2 eq 1]);
  for i := 1 to #mat_tree[d+1] do
   dd := GCD([denominator_tree[d+1][i]] cat Eltseq(mat_tree[d+1][i]));
   mat_tree[d+1][i] := mat_tree[d+1][i] div dd;
   denominator_tree[d+1][i] := denominator_tree[d+1][i] div dd;
  end for;
 end while;

 return Reverse(mat_tree), Reverse(denominator_tree);
end function;


// Run through the tree recursively. Determine the point counts modulo all the primes p covered.
function tree_recursion(co_reds, mat_tree, denominator_tree, mod_tree, i_center);
 prime_list := []; count_mod_p_list := [];

 assert #mod_tree gt 0;

 d := #mod_tree;
 for z_act := d to 1 by -1 do
  modul := mod_tree[z_act][1];
  co_act := [Integers(modul)!a : a in co_reds];
  co_act := Eltseq(Vector(co_act) * ChangeRing(mat_tree[z_act][1], Integers(modul)) * (1/Integers(modul)!denominator_tree[z_act][1]));
  pz := mod_tree[d][2^(d-z_act)];

  Append(~prime_list, pz);
  Append(~count_mod_p_list, Integers()!co_act[i_center] mod pz);

  /* The recursive call */
  if z_act lt d then
   e1,e2 := tree_recursion(co_act, [mat_tree[i][1+2^(i-z_act)..2*2^(i-z_act)] : i in [z_act + 1..d]],
                                [denominator_tree[i][1+2^(i-z_act)..2*2^(i-z_act)] : i in [z_act + 1..d]],
                                [mod_tree[i][1+2^(i-z_act)..Min(2*2^(i - z_act), #mod_tree[i])] : i in [z_act+1..d]], i_center);
   prime_list := prime_list cat e1;
   count_mod_p_list := count_mod_p_list  cat e2;
  end if;
 end for;

 return prime_list, count_mod_p_list;
end function;


// Use the current matrix, denominator, and moduli trees to compute the points modulo all the primes covered.
function scan_tree(co_act, matrix_tree, denominator_tree, modul_tree, i_center);
 prime_list :=[]; count_mod_p_list := [];
 modul_list := [a[1] : a in modul_tree];
 modul := &*modul_list;
 co_null := [Integers(modul)!a : a in Eltseq(co_act)];

 for i := 1 to #modul_list do
  z_act := #matrix_tree + 1 - i;
  co_reds := [Integers(modul_list[z_act])!a : a in co_null];
  co_reds := Eltseq(Vector(co_reds) * ChangeRing(matrix_tree[z_act][1], Integers(modul_list[z_act]))*(1/Integers(modul)!denominator_tree[z_act][1]) );

  p_act := modul_tree[#modul_tree][2^(i-1)];
  assert modul_list[z_act] mod p_act eq 0;

  Append(~prime_list, p_act);
  Append(~count_mod_p_list, Integers()!co_reds[i_center] mod p_act);

  /* Scan the tree recursively. */
  if not i in [1,#modul_list] then
   e1,e2 := tree_recursion(co_reds, [matrix_tree[z_act + d][1+2^d..2*2^d] : d in [1..i-1] | d gt 0],
                                 [denominator_tree[z_act + d][1+2^d..2*2^d] : d in [1..i-1] | d gt 0],
                                 [modul_tree[z_act + d][1+2^d..Min(2*2^d, #modul_tree[z_act+d])] : d in [1..i-1] | d gt 0], i_center);
   prime_list := prime_list cat e1;
   count_mod_p_list := count_mod_p_list cat e2;
  end if;
 end for;

 return prime_list, count_mod_p_list;
end function;


// Initialise the computation.
// Start with p = 11 and the coefficients of f3^4.
r2_gen, offsets_2, i_center := init_mov_simplex(f3);
ex := 4;
co_act, init_bound := init_iteration(f3, ex, r2_gen, offsets_2);
init_bound := Max(init_bound, 10^3);

// dim := #Eltseq(co_act); // The number of coefficients in a simplex
time_0 := Cputime();

printf "Doing direct iteration for all primes up to ", 2*init_bound+1;
time co_act, ex, prime_list, count_mod_p_list := direct_iteration(co_act, ex, init_bound, r2_gen, i_center);
co_act := Vector([Integers()!a : a in Eltseq(co_act)]);
// As the code is designed for primes < 10^9, there is no gap larger than 300 between successive primes.

r2_gen_list := [Matrix(r2_gen)];

// step_limit of  60 is suitable for search bound 10^6.
// step_limit of  80 is suitable for search bound 10^7.
// step_limit of 115 is suitable for search bound 10^8.
step_limit := 60;

printf "Create list of matrices describing an increase of the exponent by 1 to ", step_limit;
time
for i := 1 to step_limit do
 printf "%o ", i;
 fac := [[Evaluate(a, Parent(a).1 + i)  : a in b] : b in r2_gen];
 Append(~r2_gen_list, Matrix(r2_gen_list[#r2_gen_list])*Matrix(fac));
end for;

// Store numerator and denominator separately.
rr_z := PolynomialRing(Integers());
qq_z := PolynomialRing(Rationals());

gen_denominator_l := [rr_z!Denominator(a) : a in r2_gen_list];
gen_mat_l := [[[qq_z!a : a in b] : b in RowSequence(r2_gen_list[i]*gen_denominator_l[i])] : i in [1..#r2_gen_list]];

for i := 1 to #gen_mat_l do
 mul := LCM([Denominator(a) : a in Coefficients(b), b in c, c in gen_mat_l[i]]);
 gen_denominator_l[i] := gen_denominator_l[i]*mul;
 gen_mat_l[i] := [[rr_z!(mul*a) : a in b] : b in gen_mat_l[i]];
end for;

gen_mat_l := [[[rr_z!a : a in b] : b in gen_mat_l[i]] : i in [1..#gen_mat_l]];


// SetProfile(true);

assert 2*step_limit ge Max({l[i+1] - l[i] : i in [1..#l-1]}) where l is PrimesInInterval(1, lim);
depth := 8;

printf "Starting remainder forest loop with prime limit ", lim, " and initial depth ", depth;
time
repeat
 t1 := Cputime();
 mod_tree := modulus_tree_build(ex, depth);
 matrix_tree, denominator_tree := matrix_tree_build(ex, gen_mat_l, gen_denominator_l, mod_tree);
 t2 := Cputime(t1);

 e1,e2 := scan_tree(co_act, matrix_tree, denominator_tree, mod_tree, i_center);
 /* Store the modulo p point counts found. */
 prime_list := prime_list cat e1;
 count_mod_p_list := count_mod_p_list cat e2;

 ex := a[#a] div 2 where a is mod_tree[#mod_tree];
 t1a := Cputime();
 co_act := co_act * matrix_tree[1][1] div denominator_tree[1][1];
 t3 := Cputime(t1a);

 t_abs := Round(Cputime());
 t_a1 := t_abs div 60;
 printf "%o %o  Total %o %o %o\n", 2*ex+1, Cputime(t1), t_a1 div 60, t_a1 mod 60, t_abs mod 60;

 if t3 gt t2 and depth le 16 then
  depth := depth + 1;
  printf "%o %o Increasing depth to %o\n", t2, t3, depth;
 end if;
until 2*ex+1 ge lim;


/************************************************************

Runtime and memory usage for the p-adic point counting

Example       Limit     Time (sec)   Memory usage (MB)
generic        10^4       14.96          150
               10^5       21.41          150
               10^6      222.46          450

************************************************************/


/*******************************************************************************
*                                                                              *
* 2-adic point counting. A modulo 16 point count.                              *
*                                                                              *
*******************************************************************************/

// A naive point counting method, used for initialisation.
function anz_deg_2(gl, p, e)
 k := GF(p, e);
 r_ku<u> := PolynomialRing(k);
 cnt := 0;

 // Points (1:a:u)
 for a in k do
  orbit := {a^(p^j) : j in [0..e]};
  if a eq Min(orbit) then
   mul := #orbit;
   spez := Evaluate(gl, [1,a,u]);
   for b in k do
    akt := Evaluate(spez, b);
    if akt eq 0 then cnt := cnt + 1*mul; else
     if IsSquare(akt) then cnt := cnt + 2*mul; end if;
    end if;
   end for;
  end if;
 end for;
 spez := Evaluate(gl, [0,1,u]);

 // Points (0:1;b)
 for b in k do
  akt := Evaluate(spez, b);
  if akt eq 0 then cnt := cnt + 1; else
   if IsSquare(akt) then cnt := cnt + 2; end if;
  end if;
 end for;

 // Point (0:0:1)
 akt := k!Evaluate(gl, [0,0,1]);
 if akt eq 0 then cnt := cnt + 1; else
  if IsSquare(akt) then cnt := cnt + 2; end if;
 end if;

 return cnt;
end function;


// Compute the bad primes of a surface of the shape w^2 = l[1]*...*l[6]. l is assumed to be a list of six linear forms.
function bad_primes(l)
 r := Parent(l[1]);
 mat := [[MonomialCoefficient(a, r.i) : i in [1..3]] : a in l];
 mino := Minors(Matrix(Integers(), mat), 3);
 bp := PrimeDivisors(2*LCM(mino));
 return [-1] cat bp;
end function;


// Run through the good primes until every vector of Jacobi symbols is realised.
function good_primes(bp)
 p := 2;
 res := [];
 list := [];
 repeat
  p := NextPrime(p);
  if not p in bp then
   pat := [JacobiSymbol(a, p) : a in bp];
   if not pat in list then
    Append(~res, p); Append(~list, pat);
   end if;
  end if;
 until #list eq 2^#bp;
 return res, list;
end function;


// Solving the linear system of equations encoding the traces.
// We work with the trace on T. I.e., all eigenvalues have absolute value 1.
function solve_trace_system(bp, tp, lfl)
 trl := [];

 /* Enlarge the system until the solution is unique in p-adic precision 1. */
 repeat
  p_akt := tp[#trl + 1];
  Append(~trl, anz_deg_2(&*lfl,p_akt,1) - (1+ p_akt + p_akt^2));
  tr_16 := [(trl[i] * Modinv(tp[i], 16) -6) mod 16 : i in [1..#trl]]; // TateTwist
  assert &and [a mod 2 eq 0 : a in tr_16];
  tr_rest := [(a div 2) mod 2 : a in tr_16];

  // Compute the first 2-adic digit of the trace on T \subset H^2_et(V, ZZ_2).
  mat := Matrix(GF(2), [[JacobiSymbol(g, p) eq -1 select 1 else 0 : g in bp] : p in tp[1..#tr_16]]);
  sol,ker := Solution(Transpose(mat), Vector(GF(2), tr_rest));
 until Dimension(ker) eq 0;
 d1 := [Integers()!sol[i] : i in [1..#bp]];

 /* Enlarge the system further until the solution is unique in p-adic precision 2. */
 repeat
  p_akt := tp[#trl + 1];
  Append(~trl, anz_deg_2(&*lfl,p_akt,1) - (1 + p_akt + p_akt^2));
  tr_16 := [(trl[i] * Modinv(tp[i], 16) -6) mod 16 : i in [1..#trl]];
  tr_rest := [tr_16[i] - 2*&+[(JacobiSymbol(bp[j], tp[i]) eq -1) select d1[j] else 0 : j in [1..#bp]] : i in [1..#tr_16]];
  assert {a mod 4 : a in tr_rest} eq {0};

  tr_rest := [(a div 4) mod 2 : a in tr_rest];
  // 2nd 2-adic digit of the traces and first 2-adic digit of the trace of products of pairs.
  paare := [[g1,g2] : g1,g2 in bp | g1 lt g2];
  mat := Matrix(GF(2), [[JacobiSymbol(g, p) eq -1 select 1 else 0 : g in bp] cat
                        [(JacobiSymbol(a[1], p) eq -1 and
                          JacobiSymbol(a[2], p) eq -1) select 1 else 0 : a in paare] : p in tp[1..#tr_16]]);
  sol,ker := Solution(Transpose(mat), Vector(GF(2), tr_rest));
 until Dimension(ker) eq 0;

 /* Enlarge even further, taking pairs and triples into consideration. */
 for i := 1 to #bp do
  if sol[i] eq 1 then d1[i] := d1[i] + 2; end if;
 end for;

 d2 := [0 : i in paare];
 for i := 1 to #paare do
  if sol[i+ #bp] eq 1 then d2[i] := 1; end if;
 end for;

 repeat
  p_akt := tp[#trl + 1];
  Append(~trl, anz_deg_2(&*lfl, p_akt, 1) - (1 + p_akt + p_akt^2));
  tr_16 := [(trl[i] * Modinv(tp[i], 16) - 6) mod 16 : i in [1..#trl]];

  tr_rest := [tr_16[i] - 2*&+[(JacobiSymbol(bp[j], tp[i]) eq -1) select d1[j] else 0 : j in [1..#bp]]
                       - 4*&+[(JacobiSymbol(paare[j][1], tp[i]) +
                              JacobiSymbol(paare[j][2], tp[i]) eq -2) select d2[j] else 0 : j in [1..#paare]] : i in [1..#tr_16]];

  assert {a mod 8 : a in tr_rest} eq {0};
  tr_rest := [(a div 8) mod 2 : a in tr_rest];
  tripel := [[g1,g2,g3] : g1,g2,g3 in bp | g1 lt g2 and g2 lt g3];

  mat := Matrix(GF(2), [[JacobiSymbol(g, p) eq -1 select 1 else 0 : g in bp] cat
                        [(JacobiSymbol(a[1], p) eq -1 and
                          JacobiSymbol(a[2], p) eq -1) select 1 else 0 : a in paare] cat
                        [(JacobiSymbol(a[1], p) eq -1 and
                          JacobiSymbol(a[2], p) eq -1 and
                          JacobiSymbol(a[3], p) eq -1) select 1 else 0 : a in tripel] : p in tp[1..#tr_16]]);

  sol,ker := Solution(Transpose(mat), Vector(GF(2), tr_rest));
 until Dimension(ker) eq 0;
 // Solution in principle found here.

 for i := 1 to #bp do
  if sol[i] eq 1 then d1[i] := d1[i] + 4; end if;
 end for;

 for i := 1 to #paare do
  if sol[i+ #bp] eq 1 then d2[i] := d2[i] + 2; end if;
 end for;

 d3 := [0 : i in tripel];
 for i := 1 to #tripel do
  if sol[i+ #bp + #paare] eq 1 then d3[i] := 1; end if;
 end for;

 res := [];
 for i := 1 to #tp do
  p_akt := tp[i];
  Append(~res, 6+&+([0] cat [2*d1[j] : j in [1..#d1] | JacobiSymbol(bp[j], p_akt) eq -1] cat
              [4*d2[j] : j in [1..#d2] | &and [JacobiSymbol(l, p_akt) eq -1 : l in  paare[j]]] cat
              [8*d3[j] : j in [1..#d3] | &and [JacobiSymbol(l, p_akt) eq -1 : l in tripel[j]]]));
 end for;

 res := [a mod 16 : a in res];
 res1 := [(res[i] * tp[i]) mod 16 : i in [1..#res]];

 printf "Initialisation finished. %o primes used.\n", #trl;
 return res;
end function;


// Compute the bad primes and a list of good primes that realises every vector of Jacobi symbols.
bp := bad_primes(lfl);
tp,pat := good_primes(bp);

// Compute the trace mod 16 for each prime in tp. This is the initialisation of the 2-adic method, in the more efficient form.
time trl_16 := solve_trace_system(bp, tp, lfl);

/* The actual determination of the point counts.
To find the trace modulo 16, only Jacobi symbols have to be calculated.
The exact count is determined from the modulo 16 count and the modulo p count by the Chinese remainder theorem. */
count_list := [];

for i := 1 to #prime_list do
 p := prime_list[i];

 /* Ignore about the bad primes. */
 if not p in bp then
  akt := [JacobiSymbol(a, p) : a in bp];
  ii := Index(pat, akt);
  /* Doing the CRT calculation. Use the p-adic information that was computed and stored above, as well as the 2-adic information just found. */
  tr := CRT([p*trl_16[ii], count_mod_p_list[i]], [16,p]); // Undo the Tate-Twist.
  count := (1 + p + p^2) + tr;

  /* Adaptation to Deligne's bound */
  if count gt (1 + 7*p + p^2) then count := count - 16 * p;  end if;
  Append(~count_list,<p,count>);
 end if;
end for;

// Output of the results.
count_list;

