From 9bfc685d2bc5fb92857edb00c56dddd8d761afc4 Mon Sep 17 00:00:00 2001
From: Aurore Guillevic <aurore.guillevic@inria.fr>
Date: Thu, 5 Oct 2023 18:09:08 +0200
Subject: [PATCH] include better [Barbulescu-Duquesne JoC 2019] model

---
 sage/test_curves.sage              |  28 ++-
 sage/tnfs/simul/simulation_tnfs.py | 345 ++++++++++++++++++-----------
 2 files changed, 234 insertions(+), 139 deletions(-)

diff --git a/sage/test_curves.sage b/sage/test_curves.sage
index a74441f..d91545e 100644
--- a/sage/test_curves.sage
+++ b/sage/test_curves.sage
@@ -45,7 +45,7 @@ from tnfs.alpha.alpha3d import alpha3d
 # k=16 D=1 deg(px) = 16
 
 args=sys.argv
-print("usage: sage {} [--kss | --bls | --bn | --freeman | --cyclo | --auri | --fm | --fst62 | --fst63 | --fst64 | --fst66] [-choice_kss <choice>] [-code <code FM>] -k <k> -D <D> [-e0 <e0>] [-a <auri_a>] [--special | --conj | --jouxlercier | --sarkarsingh | --basem] [--lower_degree_px] [-deg_px_dec <deg_px_dec>] [-deg_h <deg_h>] [-cost <cost>] [-sievingdim <sieving_dim>] [--no-alpha] -i <idx> [-samples <samples>]".format(args[0]))
+print("usage: sage {} [--kss | --bls | --bn | --freeman | --cyclo | --auri | --fm | --fst62 | --fst63 | --fst64 | --fst66] [-choice_kss <choice>] [-code <code FM>] -k <k> -D <D> [-e0 <e0>] [-a <auri_a>] [--special | --conj | --jouxlercier | --sarkarsingh | --basem] [--lower_degree_px] [-deg_px_dec <deg_px_dec>] [-deg_h <deg_h>] [-cost <cost>] [-sievingdim <sieving_dim>] [--no-alpha] [--barbulescu_duquesne] -i <idx> [-samples <samples>]".format(args[0]))
 
 curve = None # KSS or BLS or BN or Aurifeuillean or FotiadisMartindale
 k = None
@@ -59,6 +59,7 @@ idx = None
 samples = None
 cost = None
 sieving_dim = None
+barbulescu_duquesne=False
 
 Special = True; Conj = False; JL = False; JLSV1 = False; SSingh = False; Base_m = False
 deg_h = None
@@ -67,7 +68,7 @@ deg_px_dec = 2 # decrease the degree of px by a factor deg_px_dec -> P(z^deg_px_
 compute_alpha = True
 
 i=1
-cmd_options = ["--kss", "--bls", "--bn", "--freeman", "--cyclo", "--auri", "--fm", "--fst62", "--fst63", "--fst64", "--fst66", "--fst67", "-choice_kss", "-code", "-k", "-D", "-e0", "-a", "--special", "--conj", "--jouxlercier", "--jlsv1", "--sarkarsingh", "--basem", "-deg_h", "-cost", "--lower_degree_px", "-deg_px_dec", "--no-alpha", "-sievingdim"]
+cmd_options = ["--kss", "--bls", "--bn", "--freeman", "--cyclo", "--auri", "--fm", "--fst62", "--fst63", "--fst64", "--fst66", "--fst67", "-choice_kss", "-code", "-k", "-D", "-e0", "-a", "--special", "--conj", "--jouxlercier", "--jlsv1", "--sarkarsingh", "--basem", "-deg_h", "-cost", "--lower_degree_px", "-deg_px_dec", "--no-alpha", "-sievingdim", "--barbulescu_duquesne"]
 
 while i < len(args):
     if args[i] == "--bls":
@@ -127,6 +128,8 @@ while i < len(args):
         lower_degree_px = True
     elif args[i] == "--no-alpha":
         compute_alpha = False
+    elif args[i] == "--barbulescu_duquesne":
+        barbulescu_duquesne = True
     elif args[i] == "-sievingdim" and (i+1) < len(args) and args[i+1] not in cmd_options:
         sieving_dim = Integer(args[i+1])
         i += 1
@@ -619,8 +622,12 @@ if deg_h > 1:
         print("inv_zeta_Kh, w = {:.6f},{:2d}".format(inv_zeta_Kh, w))
         print("h = {} # {} class number Kh = {}".format(h, h.list(), clKh))
         if not compute_alpha:
-            alpha_f = float(2.0)
-            alpha_g = float(2.0)
+            if barbulescu_duquesne:
+                alpha_f = float(0.0)
+                alpha_g = float(0.0)
+            else:
+                alpha_f = float(2.0)
+                alpha_g = float(2.0)
         elif clKh == 1:
             alpha_f = float(alpha_TNFS_2d(f,h,1000,test_principal=True))
             alpha_g = float(alpha_TNFS_2d(g,h,1000,test_principal=True))
@@ -635,7 +642,18 @@ if deg_h > 1:
         print("aut = {}".format(aut))
         print("alpha_f = {:.4f} alpha_g = {:.4f}".format(alpha_f,alpha_g))
         print("    ({:.8f},{:2d}, {:40s}, {:.4f}, {:.4f}, {:.4f}),".format(float(inv_zeta_Kh),int(w),str(h),float(alpha_f),float(alpha_g),float(sum_alpha)))
-        simul = Simulation_TNFS(p,ell,Fp,Fpz,h,f,g,Rxy,cost,aut,inv_zeta_Kh,count_sieving=True,alpha_f=alpha_f, alpha_g=alpha_g)
+        if barbulescu_duquesne:
+            print("barbulescu_duquesne")
+            count_sieving = False
+            weight_per_row = 128
+            cst_filtering=1.0
+            label="model of cost BD19"
+        else:
+            count_sieving = True
+            weight_per_row = 200
+            cst_filtering=20
+            label="model of cost GS21"
+        simul = Simulation_TNFS(p,ell,Fp,Fpz,h,f,g,Rxy,cost,aut,inv_zeta_Kh,None,None,alpha_f,alpha_g,weight_per_row,cst_filtering,count_sieving,label,barbulescu_duquesne)
         simul.print_params()
         simul.simulation(samples=samples)
         simul.print_results()
diff --git a/sage/tnfs/simul/simulation_tnfs.py b/sage/tnfs/simul/simulation_tnfs.py
index 0601185..f54903a 100644
--- a/sage/tnfs/simul/simulation_tnfs.py
+++ b/sage/tnfs/simul/simulation_tnfs.py
@@ -33,6 +33,7 @@ from sage.rings.integer import Integer
 from sage.rings.integer_ring import Z, ZZ
 from sage.rings.real_mpfr import RR
 from sage.misc.prandom import randint
+from sage.arith.misc import gcd
 
 from sage.functions.transcendental import dickman_rho
 from sage.symbolic.constants import euler_gamma
@@ -45,13 +46,17 @@ from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
 from tnfs.simul.polyselect import Polyselect
 
 default_weight_per_row=200
+barbulescu_duquesne_weight_per_row=128 # Barbulescu Duquesne assumes that the weight per row is 2^7 = 128
+barbulescu_duquesne_cst_linalg = float(0.25) # Barbulescu-Duquesne paper assumes that there is a constant factor c_linalg = 1/4
+machine_wordsize = 64 # assume the linear algebra step is run on a 64-bit architecture platform
 # min and max according to appendix A in https://eprint.iacr.org/2019/431
 cst_filtering_min = 9
 cst_filtering_max = 382
 default_cst_filtering = 20
+# in Barbulescu-Duquesne, the cost of sieving is modeled as corr_sieve = 1
+log2 = float(log(2.0))
 
 def log2_L_N_alpha_c(log2N, alpha, c):
-    log2 = float(log(2.0))
     logN = float(log2N)*log2
     return float(float(c)*(log2N)**(float(alpha))*(log(logN,2))**(float(1.0-alpha)))
 
@@ -87,7 +92,7 @@ def core_volume_sieving_space(A, h, inv_zeta_Kh=1.0, compute_zeta_with_pari=Fals
     :returns: volume_sieving_space(A,deg_h)*inv_zeta_Kh/w, w, inv_zeta_Kh \
     where w is the number of roots of unity up to multiplication by -1.
 
-    Thanks to Paul Zimmerman (Inria Nancy) on how to compute zeta_Kh with PARI
+    Thanks to Paul Zimmermann (Inria Nancy) on how to compute zeta_Kh with PARI
     See also http://doc.sagemath.org/html/en/reference/number_fields/sage/rings/number_field/unit_group.html
     """
     deg_h = h.degree()
@@ -112,57 +117,66 @@ def lg2_core_volume_sieving_space(A, h, inv_zeta_Kh=1.0, compute_zeta_with_pari=
     w = get_w(Kh)
     return float(RR(lg2_volume_sieving_space(A,deg_h) + log(inv_zeta_Kh,2.0) - log(w,2.0))),w,inv_zeta_Kh
 
-def find_A(deg_h, sieving_cost, start_A,log2_B=0):
+def find_A(deg_h, target_cost, start_A, count_sieving_cost=False, log2_B=None):
     """ 
     Binary search on the value of A such that given h and an expected sieving 
     cost, the sieving volume will be minimal.
-    :param deg_h: degree of polynomial h. Elements are considered modulo h. 
+    :param deg_h: int, degree of polynomial h
+    :param target_cost: target cost, formula: volume_sieving_space(A, deg_h)*(cost of sieving ~= log log B)
+    :param start_A: int, a starting value for A
+    :param log2_B: float
+
+    Elements of Kh are considered as polynomials modulo h.
+    The cost of sieving is estimated to be log(log(B)) in basis e (not in basis 2)
+    where log log B = log (log_2 B * log(2))
     """
-    print("find_A, log2_B={}".format(log2_B))
+    print("find_A, target cost {}, count_sieving_cost is {}, log2_B={}".format(target_cost, count_sieving_cost, log2_B))
     A = start_A
-    v = volume_sieving_space(A, deg_h)
-    if log2_B > 1:
-        corr_log2_B = log(log(log2_B*log(2.0)),2.0)
+    lg2_v = lg2_volume_sieving_space(A, deg_h)
+    # if count_sieving_cost then assume that the cost of sieving is log log B [Guillevic Singh 21, Section 6.1]
+    if count_sieving_cost and log2_B is not None and log2_B > 1.0:
+        # log log B = log (log_2 B * log(2))
+        corr_log2_B = float(log(log(log2_B*log2),2.0))
     else:
-        corr_log2_B = 0.0
-    s = float(log(v,2.0) + corr_log2_B)
+        corr_log2_B = float(0.0)
+    s = lg2_v + corr_log2_B
     A_min = A
     A_max = A
-    while s < sieving_cost: # increase A
+    while s < target_cost: # increase A
         A = max(A+1, floor(1.1*A_min)) # if A < 9 it will loop if there is no max
-        v = volume_sieving_space(A, deg_h)
-        s = float(log(v,2.0) + corr_log2_B)
-        if s < sieving_cost:
+        lg2_v = lg2_volume_sieving_space(A, deg_h)
+        s = lg2_v + corr_log2_B
+        if s < target_cost:
             A_min = A
         else:
             A_max = A
-    while s > sieving_cost: # decrease A
+    while s > target_cost: # decrease A
         A = min(A-1, floor(0.9*A_max))
-        v = volume_sieving_space(A, deg_h)
-        s = float(log(v,2.0) + corr_log2_B)
-        if s > sieving_cost:
+        lg2_v = lg2_volume_sieving_space(A, deg_h)
+        s = lg2_v + corr_log2_B
+        if s > target_cost:
             A_max = A
         else:
             A_min = A
     while ((A_max - A_min) > 1):
         A_mid = (A_min + A_max)//2
-        v = volume_sieving_space(A_mid, deg_h)
-        s = float(log(v,2.0) + corr_log2_B)
-        if s < sieving_cost:
+        lg2_v = lg2_volume_sieving_space(A_mid, deg_h)
+        s = lg2_v + corr_log2_B
+        if s < target_cost:
             A_min = A_mid
         else:
             A_max = A_mid
-    v1 = volume_sieving_space(A_min, deg_h)
-    s1 = float(log(v1,2.0) + corr_log2_B)
-    v2 = volume_sieving_space(A_max, deg_h)
-    s2 = float(log(v2,2.0) + corr_log2_B)
+    lg2_v1 = lg2_volume_sieving_space(A_min, deg_h)
+    s1 = lg2_v1 + corr_log2_B
+    lg2_v2 = lg2_volume_sieving_space(A_max, deg_h)
+    s2 = lg2_v2 + corr_log2_B
     return A_min, s1, A_max, s2
 
-def print_A(s_min, s_max, deg_h, start_A):
+def print_A(s_min, s_max, deg_h, start_A, count_sieving_cost=False, log2_B=None):
     print("params_A = {")
     A1 = start_A
     for vol in range(s_max-1, s_min-2, -1):
-        A0, s0, A1, s1 = find_A(deg_h, vol, A1)
+        A0, s0, A1, s1 = find_A(deg_h, vol, A1, count_sieving_cost, log2_B)
         print("    ({},{}): {},".format(deg_h,(vol+1), A1))
     print("}")
 
@@ -181,32 +195,46 @@ def get_w(Kh):
     # count those which have a small size of coefficients, for example when the constant coefficient of h(y) is +/-1, the element 'y' has a unit, it is an algebraic integer and a unit.
     return w
 
-def time_linalg(log2_B, ell_word_size, cst_filtering=cst_filtering_min, weight_per_row=default_weight_per_row):
+def time_linalg(log2_B, ell_word_size, weight_per_row=default_weight_per_row, cst_filtering=default_cst_filtering, barbulescu_duquesne=False, aut=None):
     """ Computes an estimate of the cost of the linear algebra step.
 
     :param log2_B: the log in basis 2 of the factor bound (relations involve primes <= B). Assume the same bound on both sides.
     :param ell_word_size: the number of machine-words (64-bit words) to write ell. The linear algebra operations are performed modulo ell.
-    :param cst_filtering: a constant s.t. after the filtering step, the number of rows of the matrix is divided by this constant. Experimentally, 7.43 <= cst_filtering <= 386 in recent computations (< 10 years)
     :param weight_per_row: the number of non-zero entries per row, arbitrarily 200. It is from 125 to 200 in recent computations.
+    :param cst_filtering: a constant s.t. after the filtering step, the number of rows of the matrix is divided by this constant. Experimentally, 7.43 <= cst_filtering <= 386 in recent computations (< 10 years)
+    :param barbulescu_duquesne: if True, apply Conjecture 1 in [BD19] claiming a variable cst_filtering = (log B)^(1+o(1))
     :returns log2_time_linalg, log2_size_FB, log2_size_matrix, logB (all float)
 
     Thanks to Emmanuel Thome (Inria Nancy).
     The size of the factor basis is estimated to be 2*log_integral(B)
     The number of rows of the matrix is the above size divided by a constant factor 'cst_filtering'
+    or a variable factor log(B) if barbulescu_duquesne is set to True
+    [BD19] Barbulescu, Duquesne, Updating key size estimations for pairings, Journal of Cryptology
+    https://doi.org/10.1007/s00145-018-9280-5
+    Conjecture 1:
+        In the filtering step of NFS, one reduces the matrix by a
+        factor (log B)^(1+o(1)), where B is the smoothness bound.
+
     The cost is estimated to be
     word_size(ell) * weight_per_row * (matrix size (=number of rows) after filtering step)^2
     """
-    logB = RR(log2_B*log(2.0))
-    #log2_size_FB = float(1+log2_B - log(logB,2))
-    log2_size_FB = RR(1.0+log(log_integral(2.0**log2_B),2.0))
-    # assume that the number of rows in the matrix is divided by a factor 'cst_filtering'
-    log2_size_matrix = RR(log2_size_FB - log(cst_filtering,2.0))
-    # assume that the time is (word_size(ell) * weight_per_row * (number_of_rows_matrix)^2)
-    log2_time_linalg = RR(log(ell_word_size,2.0) + log(weight_per_row,2.0) + 2.0*log2_size_matrix)
-    return float(log2_time_linalg), float(log2_size_FB), float(log2_size_matrix), float(logB)
-
-
-def find_log2_B(log2_target_cost, start_log2B, prec, fct_time_linalg, ell_word_size, cst_filtering=cst_filtering_min, weight_per_row=default_weight_per_row):
+    logB = float(log2_B*log2)
+    if not barbulescu_duquesne:
+        #log2_size_FB = log2(2*log_integral(B)) = 1 + log2(log_integral(B))
+        log2_size_FB = float(1.0+log(log_integral(2.0**log2_B),2.0))
+        # assume that the number of rows in the matrix is divided by a constant factor 'cst_filtering'
+        log2_size_matrix = float(log2_size_FB - log(cst_filtering, 2.0))
+        # assume that the time is (word_size(ell) * weight_per_row * (number_of_rows_matrix)^2)
+        log2_time_linalg = float(log(ell_word_size,2.0) + log(weight_per_row,2.0) + 2.0*log2_size_matrix)
+    else:
+        log2_size_FB = float(1.0 + log2_B - log(logB,2.0))
+        # assume that the number of rows in the matrix is divided by a variable factor log(B) according to Conjecture 1
+        log2_size_matrix = float(log2_size_FB - log(logB, 2.0) - log(aut, 2.0))
+        # assume that the time is (weight_per_row * (number_of_rows_matrix)^2)
+        log2_time_linalg = float(log(barbulescu_duquesne_weight_per_row,2.0) + log(barbulescu_duquesne_cst_linalg,2.0) + 2.0*log2_size_matrix)
+    return log2_time_linalg, log2_size_FB, log2_size_matrix, logB
+
+def find_log2_B(log2_target_cost, start_log2B, prec, fct_time_linalg, ell_word_size, weight_per_row=default_weight_per_row, cst_filtering=default_cst_filtering, barbulescu_duquesne=False, aut=None):
     """ 
     Binary search on the value of log2(B) such that given a function on input B and word_size(ell) that computes the cost (as 2^c) of the linear algebra step, 
     find B such that the estimated cost is 2^log2_target_cost
@@ -214,45 +242,47 @@ def find_log2_B(log2_target_cost, start_log2B, prec, fct_time_linalg, ell_word_s
     :param start_log2B: a hint on the starting value of log2(B)
     :param prec: the precision, stops when (B_max-B_min) <= prec
     :param fct_time_linalg: a function whose inputs are log2_B and ell_word_size and 1st output is the estimated cost of linear algebra (in bits)
-    :param ell_word_size: the number
     :param ell_word_size: the number of machine-words (64-bit words) to write ell. The linear algebra operations are performed modulo ell.
+    :param weight_per_row: the number of non-zero entries per row, arbitrarily 200. It is from 125 to 200 in recent computations.
+    :param cst_filtering: a constant s.t. after the filtering step, the number of rows of the matrix is divided by this constant. Experimentally, 7.43 <= cst_filtering <= 386 in recent computations (< 10 years)
+    :param barbulescu_duquesne: if True, apply Conjecture 1 in [BD19] claiming a variable cst_filtering = (log B)^(1+o(1)), and a linear algebra cost of barbulescu_duquesne_weight_per_row*barbulescu_duquesne_cst_linalg = 128/4 = 2^5
     """
     B = start_log2B
-    tt, FB, M, logB = fct_time_linalg(B,ell_word_size, cst_filtering=cst_filtering, weight_per_row=weight_per_row)
-    tt = float(tt)
+    tt, FB, M, logB = fct_time_linalg(B, ell_word_size, weight_per_row, cst_filtering, barbulescu_duquesne, aut)
+    #tt = float(tt)
     B_min = B
     B_max = B
     while tt < log2_target_cost:
         B = floor(1.1*B_min)
-        tt, FB, M, logB = fct_time_linalg(B,ell_word_size)
-        tt = float(tt)
+        tt, FB, M, logB = fct_time_linalg(B, ell_word_size, weight_per_row, cst_filtering, barbulescu_duquesne, aut)
+        #tt = float(tt)
         if tt < log2_target_cost:
             B_min = B
         else:
             B_max = B
     while tt > log2_target_cost:
         B = floor(0.9*B_max)
-        tt, FB, M, logB = fct_time_linalg(B,ell_word_size)
-        tt = float(tt)
+        tt, FB, M, logB = fct_time_linalg(B, ell_word_size, weight_per_row, cst_filtering, barbulescu_duquesne, aut)
+        #tt = float(tt)
         if tt > log2_target_cost:
             B_max = B
         else:
             B_min = B
     while ((B_max - B_min) > prec):
         B_mid = float((B_min + B_max)/2)
-        tt, FB, M, logB = fct_time_linalg(B_mid,ell_word_size)
-        tt = float(tt)
+        tt, FB, M, logB = fct_time_linalg(B_mid, ell_word_size, weight_per_row, cst_filtering, barbulescu_duquesne, aut)
+        #tt = float(tt)
         if tt < log2_target_cost:
             B_min = B_mid
         else:
             B_max = B_mid
-    t1, FB, M, logB = fct_time_linalg(B_min,ell_word_size)
-    t2, FB, M, logB = fct_time_linalg(B_max,ell_word_size)
-    t1 = float(t1)
-    t2 = float(t2)
+    t1, FB, M, logB = fct_time_linalg(B_min, ell_word_size, weight_per_row, cst_filtering, barbulescu_duquesne, aut)
+    t2, FB, M, logB = fct_time_linalg(B_max, ell_word_size, weight_per_row, cst_filtering, barbulescu_duquesne, aut)
+    #t1 = float(t1)
+    #t2 = float(t2)
     return B_min, t1, B_max, t2
     
-def print_B(s_min, s_max, start_log2B, fct_time_linalg, ell_word_size):
+def print_B(s_min, s_max, start_log2B, fct_time_linalg, ell_word_size, weight_per_row=default_weight_per_row, cst_filtering=default_cst_filtering, barbulescu_duquesne=False, aut=None):
     """
     how to call this function:
     print_B(100, 150, 83.5, time_linalg, 4)
@@ -261,11 +291,11 @@ def print_B(s_min, s_max, start_log2B, fct_time_linalg, ell_word_size):
     print("params_log2B_ell_w = {")
     B1 = start_log2B
     for t in range(s_max-1, s_min-2, -1):
-        B0, t0, B1, t1 = find_log2_B(t, B1, 0.001, fct_time_linalg, ell_word_size)
+        B0, t0, B1, t1 = find_log2_B(t, B1, 0.001, fct_time_linalg, ell_word_size, weight_per_row, cst_filtering, barbulescu_duquesne, aut)
         print("    ({0},{1}): {2:.3f},".format(ell_word_size, (t+1), float(B1)))
     print("}")
 
-def get_parameters_A_log2B(cost, ell_wordsize, aut, deg_h, cst_filtering=cst_filtering_min, weight_per_row=default_weight_per_row,count_sieving_cost=False):
+def get_parameters_A_log2B(cost, ell_wordsize, aut, deg_h, weight_per_row=default_weight_per_row, cst_filtering=default_cst_filtering, barbulescu_duquesne=False, count_sieving_cost=False):
     """ total expected cost in bits """
     cost = float(cost)
     prec = 0.0001
@@ -273,17 +303,15 @@ def get_parameters_A_log2B(cost, ell_wordsize, aut, deg_h, cst_filtering=cst_fil
     # time_linalg is a function
     # whose optional parameters are
     # cst_filtering, weight_per_row
-    # maybe pu these functions inside the class...
-    B_min, t1, B_max, t2 = find_log2_B(cost-1, start_log2B, prec, time_linalg, ell_wordsize, cst_filtering=cst_filtering, weight_per_row=weight_per_row)
+    # maybe put these functions inside the class...
+    B_min, t1, B_max, t2 = find_log2_B(cost-1, start_log2B, prec, time_linalg, ell_wordsize, weight_per_row, cst_filtering, barbulescu_duquesne, aut)
     log2_B = B_max
-    
-    start_A = floor(2**(RR((cost-1)/(2*deg_h))))
-    if count_sieving_cost:
-        A_min, s1, A_max, s2 = find_A(deg_h, float(cost-1.0+log(aut,2.0)), start_A, log2_B=log2_B)
-    else:
-        A_min, s1, A_max, s2 = find_A(deg_h, float(cost-1.0+log(aut,2.0)), start_A)
+
+    # starting rough estimate: |a_i|^(2*deg_h) ~= 2^(cost-1) => |a_i| ~= 2^((cost-1)/(2*deg_h))
+    start_A = floor(float(2.0)**(float((cost-1)/(2*deg_h))))
+    A_min, s1, A_max, s2 = find_A(deg_h, float(cost-1.0+log(aut,2.0)), start_A, count_sieving_cost, log2_B)
     # we should have s1 <= cost <= s2 and A_min + 1 = A_max
-    if abs(cost-1+log(aut,2)-s1) < 0.5 and abs(cost-1.0+log(aut,2.0)-s2) >= 0.5:
+    if abs(cost-1+log(aut,2.0)-s1) < 0.5 and abs(cost-1.0+log(aut,2.0)-s2) >= 0.5:
         A = A_min
     else:
         A = A_max
@@ -355,12 +383,12 @@ def statistics(h,f,g,Rxy,A,logB,samples,alpha_f=0.0,alpha_g=0.0,check_coprime=Tr
             # when there is an automorphism in Kh y->-y or y->1/y, then uses
             #a=random_vec_positive_lc(A,deg_h)
             if check_coprime:
-                coprime_ideal=(Kh.ideal(Kh(a))+Kh.ideal(Kh(b))==Kh.ideal(1))
+                coprime_ideal = (gcd(gcd(a), gcd(b)) == 1) and (Kh.ideal(Kh(a))+Kh.ideal(Kh(b))==Kh.ideal(1))
                 while not coprime_ideal:
                     a=random_vec(A,deg_h)
                     b=random_vec_positive_lc(A,deg_h)
                     duplicates += 1
-                    coprime_ideal=(Kh.ideal(Kh(a))+Kh.ideal(Kh(b))==Kh.ideal(1))
+                    coprime_ideal = (gcd(gcd(a), gcd(b)) == 1) and (Kh.ideal(Kh(a))+Kh.ideal(Kh(b))==Kh.ideal(1))
             #a=(1+2*y) + (3+4*y)*x
             #h.resultant(a.resultant(f0,x),y)
             log_Nf[i] = log(RR(abs(h.resultant(f.resultant(Rxy([a,b]))))))
@@ -417,7 +445,8 @@ class Simulation_TNFS():
                  A=None, log2_B=None,
                  alpha_f=None, alpha_g=None,
                  weight_per_row=None, cst_filtering=None,
-                 count_sieving=None, label=None):
+                 count_sieving=None, label=None,
+                 barbulescu_duquesne=False):
         """
         :param              p: prime integer defining a prime field
         :param            ell: prime integer, order of the mult. subgroup where to compute DL
@@ -438,6 +467,7 @@ class Simulation_TNFS():
         :param  cst_filtering: default 20, reduction factor of the matrix
         :param  count_sieving: True/False, whether count the cost of sieving or not
         :param          label: for printing (such as name of the pairing-friendly curve)
+        :param barbulescu_duquesne: if True, apply Conjecture 1 in [BD19] claiming a variable cst_filtering = (log B)^(1+o(1))
 
         put somewhere else:
         :param            A: bound on the coefficients a_ij in the relation collection -A <= a_ij <= A
@@ -450,7 +480,7 @@ class Simulation_TNFS():
         """
         self.p = p     #
         self.ell = ell
-        self.ell_wordsize = ceil(RR(log(ell,2))/64)
+        self.ell_wordsize = ceil(RR(log(ell,2))/machine_wordsize)
         self.Fp = Fp   # there 3 params are in classes PFCurve (BN, BLS12, KSS18, BLS24)
         self.Fpz = Fpz #
         self.h = h
@@ -473,20 +503,30 @@ class Simulation_TNFS():
 
         self.check_coprime_ideals=True
         #
-        if weight_per_row == None or weight_per_row <= 0:
-            self.weight_per_row = default_weight_per_row
+        if weight_per_row is None or weight_per_row <= 0:
+            if not barbulescu_duquesne:
+                self.weight_per_row = default_weight_per_row
+            else:
+                self.weight_per_row = barbulescu_duquesne_weight_per_row
         else:
             self.weight_per_row = weight_per_row
         
-        if cst_filtering == None or cst_filtering <= 0:
-            self.cst_filtering = default_cst_filtering
+        if cst_filtering is None or cst_filtering <= 0:
+            if not barbulescu_duquesne:
+                self.cst_filtering = default_cst_filtering
+            else:
+                self.cst_filtering = barbulescu_duquesne_cst_filtering
         else:
             self.cst_filtering = cst_filtering
-        if count_sieving == None:
-            self.count_sieving = True
+        if count_sieving is None:
+            if not barbulescu_duquesne:
+                self.count_sieving = True
+            else:
+                self.count_sieving = False
         else:
             self.count_sieving = count_sieving
         self.label = label
+        self.barbulescu_duquesne = barbulescu_duquesne
 
         self.set_cost(log2_cost, A=A, log2_B=log2_B)
 
@@ -496,7 +536,7 @@ class Simulation_TNFS():
             raise ValueError("Error cost: should be a strictly positive number but received {}".format(cost))        
         self.cost = cost
         if A is None or log2_B is None:
-            self.A, self.log2_B = get_parameters_A_log2B(self.cost, self.ell_wordsize, self.aut, self.deg_h, self.cst_filtering, self.weight_per_row,count_sieving_cost=self.count_sieving)
+            self.A, self.log2_B = get_parameters_A_log2B(self.cost, self.ell_wordsize, self.aut, self.deg_h, self.weight_per_row, self.cst_filtering, self.barbulescu_duquesne, self.count_sieving)
         else:
             self.A = A
             self.log2_B = log2_B
@@ -504,28 +544,47 @@ class Simulation_TNFS():
         self.core_vol_sieving_space, self.w, self.inv_zeta_Kh = core_volume_sieving_space(self.A, self.h, inv_zeta_Kh=self.inv_zeta_Kh)
         self.vol_sieving_space = volume_sieving_space(self.A, self.deg_h)
         self.log2_vol_sieving_space = lg2_volume_sieving_space(self.A, self.deg_h)
-        self.logB = float(self.log2_B*log(2))
+        self.logB = float(self.log2_B*log2)
         self.loglogB = float(log(self.logB)) # log(log(B)) where log=ln is natural logarithm (log(e) = 1)
-        self.log2_time_relcol = float(log(self.vol_sieving_space/self.aut,2))
-        self.log2_time_relcol_sieving_cost = self.log2_time_relcol + float(log(self.loglogB,2))
-
-        self.log2_time_linalg, self.log2_size_FB, self.log2_size_matrix, self.logB = time_linalg(self.log2_B, self.ell_wordsize)
-        self.log2_total_time = float(log(2**self.log2_time_relcol + 2**self.log2_time_linalg ,2))
-        self.log2_total_time_sieving_cost = float(log(2**self.log2_time_relcol_sieving_cost + 2**self.log2_time_linalg ,2))
-        
-        self.log2_size_FB_BD = float(1+self.log2_B - log(self.logB,2))
-        self.log2_size_matrix_BD = float(self.log2_size_FB_BD - log(self.log2_B,2))
-        self.log2_time_linalg_BD = float(5 + 2*self.log2_size_matrix_BD)
-        self.log2_total_time_BD = float(log(2**self.log2_time_relcol + 2**self.log2_time_linalg_BD ,2))
+        self.log2_time_relcol = float(log(self.vol_sieving_space/self.aut,2.0))
+        if self.count_sieving:
+            self.log2_time_relcol_sieving_cost = self.log2_time_relcol + float(log(self.loglogB,2.0))
+        else:
+            self.log2_time_relcol_sieving_cost = self.log2_time_relcol
+
+        self.log2_time_linalg, self.log2_size_FB, self.log2_size_matrix, self.logB = time_linalg(self.log2_B, self.ell_wordsize, self.weight_per_row, self.cst_filtering, self.barbulescu_duquesne, self.aut)
+        self.log2_total_time = float(log(2.0**self.log2_time_relcol + 2.0**self.log2_time_linalg, 2.0))
+        self.log2_total_time_sieving_cost = float(log(2.0**self.log2_time_relcol_sieving_cost + 2.0**self.log2_time_linalg, 2.0))
+
+        #in [BD19], the size of the factor basis uses the rough estimate 2*B/log B whose log_2 is 1 + log2_B - log_2(log B)
+        # Equation 2 of the model of cost in Barbulescu-Dquesne
+        # cost = {number of relations to collect} + {cost of linear algebra}
+        # cost = 2*B/(aut*log(B)) / (dickman_rho(log_2(Nf)/log_2(B)) * dickman_rho(log_2(Ng)/log_2(B)))
+        #        + 2^7 * (B/(aut*log(B)*log_2(B)))^2
+        # where aut = number of automorphisms, <= eta*kappa/gcd(eta*kappa)
+        #                                         where eta*kappa=k and eta=deg(h)
+        # size of the matrix is expected to be (size of factor basis)/(filtering_factor * aut)
+        # = (2*B/log(B))/(log_2(B) * aut)
+        # the cost of linear algebra is expected to be
+        # the size of the matrix squared, times the weight per row (the non-zero entries) times some factor c_linalg that can be modeled to be 1/4 (not 1 but 1/4), finally
+        # ((2*B/log(B))/(log_2(B) * aut))^2 * 128 / 4
+        # there is no factor (log2_r/wordsize) in Barbulescu-Duquesne cost model
+        self.log2_size_FB_BD = float(1.0 + self.log2_B - log(self.logB, 2.0))
+        # consider c_filter = log2_B + automorphisms
+        self.log2_size_matrix_BD = float(self.log2_size_FB_BD - log(self.log2_B, 2.0) - log(self.aut, 2.0))
+        self.log2_time_linalg_BD = float(log(barbulescu_duquesne_weight_per_row,2.0) + log(barbulescu_duquesne_cst_linalg,2.0) + 2.0*self.log2_size_matrix_BD)
+        # where for Barbulescu-Duquesne, weight_per_row = 128 and c_linalg = 1/4, hence 128/4 = 32 = 2^5, explaining why this is 2^5 and not 2^7.
+        self.log2_total_time_BD = float(log(2.0**self.log2_time_relcol + 2.0**self.log2_time_linalg_BD, 2.0))
             
     def print_params(self):
         #def run_simulation_tnfs(p,Fp, Fpz,h,f,g,Rxy,A,log2_B,ell_wordsize,samples=10**5,inv_zeta_Kh=1.0,aut=1,alpha_f=0.0,alpha_g=0.0,check_coprime=True):
-        
+        if self.barbulescu_duquesne:
+            print("barbulescu_duquesne=True, weight_per_row={}, count_sieving={}".format(self.weight_per_row, self.count_sieving))
         print("cost = {0} bits, deg_h = {1}, A = {2}, ell_wordsize = {3}, B = 2^{4:.3f} \n".format(self.cost, self.deg_h, self.A, self.ell_wordsize, float(self.log2_B)))
         print("vol sieving space = 2^{0:.2f}, core vol sieving = 2^{1:.2f}, w = {2}, 1/zeta_Kh(2) = {3:.6f}".format(float(log(self.vol_sieving_space,2)), float(log(self.core_vol_sieving_space,2)), self.w, float(self.inv_zeta_Kh)))
         if self.with_alpha:
             print("alpha_f ={0:7.4f} (basis e), alpha_g ={1:7.4f} (basis e), sum ={2:7.4f}".format(float(self.alpha_f), float(self.alpha_g), float(self.alpha_f+self.alpha_g)))
-            print("alpha_f ={0:7.4f} (basis 2), alpha_g ={1:7.4f} (basis 2), sum ={2:7.4f}".format(float(self.alpha_f/log(2)), float(self.alpha_g/log(2)), float(self.alpha_f/log(2)+self.alpha_g/log(2))))
+            print("alpha_f ={0:7.4f} (basis 2), alpha_g ={1:7.4f} (basis 2), sum ={2:7.4f}".format(float(self.alpha_f/log2), float(self.alpha_g/log2), float(self.alpha_f/log2+self.alpha_g/log2)))
 
     def simulation(self, samples):
         self.samples = samples
@@ -555,52 +614,64 @@ class Simulation_TNFS():
             # do not use the total space for the relation collection, stop when enough rels are obtained
             # time constraint: stop when the time consumed is cost/2 (log_2 cost - 1)
             # where to put the gain thanks to an automorphism?
-            log2_time_relcol_sieving_cost = RR(self.cost-1)
-            log2_cost_sieving = RR(log(self.loglogB,2))
-            log2_vol_sieving_space = log2_time_relcol_sieving_cost - log2_cost_sieving + RR(log(self.aut,2))
-            log2_rels = RR(log(self.inv_zeta_Kh,2) + log2_vol_sieving_space + log(self.avrg_proba,2))
-            log2_total_time_sieving_cost = float(log(2**log2_time_relcol_sieving_cost + 2**self.log2_time_linalg,2))
+            log2_target_relcol_cost = RR(self.cost-1)
+            if self.count_sieving:
+                log2_cost_sieving = RR(log(self.loglogB,2.0))
+            else:
+                log2_cost_sieving = float(0.0)
+            log2_vol_sieving_space = log2_target_relcol_cost - log2_cost_sieving + RR(log(self.aut,2.0))
+            log2_rels = RR(log(self.inv_zeta_Kh,2.0) + log2_vol_sieving_space + log(self.avrg_proba,2.0))
+            log2_total_time_relcol_cost = float(log(2.0**log2_target_relcol_cost + 2.0**self.log2_time_linalg,2.0))
             print("adjusting volume sieving space")
-            print("time relcol = 2^{0:.4f}".format(float(log2_time_relcol_sieving_cost)))
+            print("time relcol = 2^{0:.4f}".format(float(log2_target_relcol_cost)))
             print("time linalg = 2^{0:.4f}".format(float(self.log2_time_linalg)))
-            print("total time  = 2^{0:.4f} (counting sieving cost as ln(ln(B)))".format(log2_total_time_sieving_cost))
+            if self.count_sieving:
+                print("total time  = 2^{0:.4f} (counting sieving cost as ln(ln(B)))".format(log2_total_time_relcol_cost))
+            else:
+                print("total time  = 2^{0:.4f} (counting sieving cost as free (c_sieve=1))".format(log2_total_time_relcol_cost))
             print("relations obtained:                                  2^{:.4f}".format(float(log2_rels)))
             # now minimal time to ensure enough relations
             # minimal_vol_sieving_space*inv_zeta_Kh*avrg_proba = needed_relations
-            # minimal_vol_sieving_space * cost_sieving = minimal_time_relcol
-            log2_minimal_vol_sieving_space = RR(self.log2_size_FB - log(self.inv_zeta_Kh,2) - log(self.avrg_proba,2))
-            log2_minimal_time_relcol = RR(log2_minimal_vol_sieving_space - RR(log(self.aut,2)) + log2_cost_sieving)
+            # minimal_vol_sieving_space/aut * cost_sieving = minimal_time_relcol
+            log2_minimal_vol_sieving_space = RR(self.log2_size_FB - log(self.inv_zeta_Kh,2.0) - log(self.avrg_proba,2.0))
+            log2_minimal_time_relcol = RR(log2_minimal_vol_sieving_space - RR(log(self.aut,2.0)) + log2_cost_sieving)
             # rels_min = inv_zeta_Kh * time_rel_col_sieving_cost/time_sieving * avrg_proba
-            log2_rels_min = RR(log(self.inv_zeta_Kh,2) + log2_minimal_vol_sieving_space + log(self.avrg_proba,2))
+            log2_rels_min = RR(log(self.inv_zeta_Kh,2.0) + log2_minimal_vol_sieving_space + log(self.avrg_proba,2.0))
             print("time relcol = 2^{0:.4f}".format(float(log2_minimal_time_relcol)))
             print("time linalg = 2^{0:.4f}".format(float(self.log2_time_linalg)))
-            log2_total_time_sieving_cost = float(log(2**log2_minimal_time_relcol + 2**self.log2_time_linalg ,2))
-            print("total time  = 2^{0:.4f} (counting sieving cost as ln(ln(B)))".format(log2_total_time_sieving_cost))
+            log2_total_time_relcol_cost = float(log(2.0**log2_minimal_time_relcol + 2.0**self.log2_time_linalg,2.0))
+            if self.count_sieving:
+                print("total time  = 2^{0:.4f} (counting sieving cost as ln(ln(B)))".format(log2_total_time_relcol_cost))
+            else:
+                print("total time  = 2^{0:.4f} (counting sieving cost as free (c_sieve=1))".format(log2_total_time_relcol_cost))
             print("relations obtained:                                  2^{:.4f}".format(float(log2_rels_min)))
     
     def print_results(self):
         print("experimental 1/zeta_Kh(2) = {:.8f}".format(float(self.inv_zeta_Kh_exp)))
         print("theoretical  1/zeta_Kh(2) = {:.8f}".format(float(self.inv_zeta_Kh)))
-        #print("avrg_log2_Nf = {}\navrg_log2_Ng = {}\navrg_smooth_proba_f = {}\navrg_smooth_proba_g = {}\navrg_proba  = {}\n".format(float(self.avrg_log_Nf/log(2)), float(self.avrg_log_Ng/log(2)), self.avrg_smooth_proba_f, self.avrg_smooth_proba_g, self.avrg_proba))
+        #print("avrg_log2_Nf = {}\navrg_log2_Ng = {}\navrg_smooth_proba_f = {}\navrg_smooth_proba_g = {}\navrg_proba  = {}\n".format(float(self.avrg_log_Nf/log2), float(self.avrg_log_Ng/log2), self.avrg_smooth_proba_f, self.avrg_smooth_proba_g, self.avrg_proba))
         print("err_norm_f  = {}\nerr_norm_g  = {}\nerr_proba_f = {}\nerr_proba_g = {}\nerr_proba   = {}\n".format(self.err_norm_f, self.err_norm_g, self.err_proba_f, self.err_proba_g, self.err_proba))
         print("std_dev_norm_f = {}\nstd_dev_norm_g = {}\nstd_dev_proba_f = {}\nstd_dev_proba_g = {}\nstd_dev_proba   = {}\n".format(self.std_dev_norm_f, self.std_dev_norm_g, self.std_dev_proba_f, self.std_dev_proba_g, self.std_dev_proba))
         
         print("A = {0}, B = 2^{1:.3f}, deg(h) = {2}, #{{units}} = w = {3}, samples = {4}".format(self.A,float(self.log2_B),self.deg_h,self.w,self.samples))
         print("vol_sieving_space = (2*A+1)^(2*{0})/(2*{1}) = 2^{2:.3f}".format(self.deg_h, self.w, float(log(self.vol_sieving_space, 2))))
-        print("avrg_log2_Nf = {0:9.4f} (bits) avrg_log2_Ng = {1:9.4f} (bits) sum = {2:9.4f} (bits)".format(float(self.avrg_log_Nf/log(2)), float(self.avrg_log_Ng/log(2)), float((self.avrg_log_Nf+self.avrg_log_Ng)/log(2))))
-        print(" min_log2_Nf = {0:9.4f} (bits)  min_log2_Ng = {1:9.4f} (bits)".format(float(self.min_log_Nf/log(2)), float(self.min_log_Ng/log(2))))
-        print(" max_log2_Nf = {0:9.4f} (bits)  max_log2_Ng = {1:9.4f} (bits)".format(float(self.max_log_Nf/log(2)), float(self.max_log_Ng/log(2))))
+        print("avrg_log2_Nf = {0:9.4f} (bits) avrg_log2_Ng = {1:9.4f} (bits) sum = {2:9.4f} (bits)".format(float(self.avrg_log_Nf/log2), float(self.avrg_log_Ng/log2), float((self.avrg_log_Nf+self.avrg_log_Ng)/log2)))
+        print(" min_log2_Nf = {0:9.4f} (bits)  min_log2_Ng = {1:9.4f} (bits)".format(float(self.min_log_Nf/log2), float(self.min_log_Ng/log2)))
+        print(" max_log2_Nf = {0:9.4f} (bits)  max_log2_Ng = {1:9.4f} (bits)".format(float(self.max_log_Nf/log2), float(self.max_log_Ng/log2)))
         
-        print("smooth_proba(avrg_Nf,logB,alpha_f) = smooth_proba({0:.2f},2^{1:.2f},{2:.4f}) = {3} = 2^{4:.4f}".format(float(self.avrg_log_Nf*log(2)),float(self.log2_B),float(self.alpha_f), float(smoothness_proba(self.avrg_log_Nf,self.logB,self.alpha_f)),float(log(smoothness_proba(self.avrg_log_Nf,self.logB,self.alpha_f),2))))
-        print("smooth_proba(avrg_Ng,logB,alpha_g) = smooth_proba({0:.2f},2^{1:.2f},{2:.4f}) = {3} = 2^{4:.4f}".format(float(self.avrg_log_Ng*log(2)),float(self.log2_B),float(self.alpha_g), float(smoothness_proba(self.avrg_log_Ng,self.logB,self.alpha_g)),float(log(smoothness_proba(self.avrg_log_Ng,self.logB,self.alpha_g),2))))
-        print("prod of smooth proba of avrg norms: {0} = 2^{1:.4f}".format(float(self.proba_rough),float(log(self.proba_rough,2))))
-        print("avrg_proba = {0} = 2^{1:.4f}".format(float(self.avrg_proba),float(log(self.avrg_proba,2))))
+        print("smooth_proba(avrg_Nf,logB,alpha_f) = smooth_proba({0:.2f},2^{1:.2f},{2:.4f}) = {3} = 2^{4:.4f}".format(float(self.avrg_log_Nf*log2),float(self.log2_B),float(self.alpha_f), float(smoothness_proba(self.avrg_log_Nf,self.logB,self.alpha_f)),float(log(smoothness_proba(self.avrg_log_Nf,self.logB,self.alpha_f),2.0))))
+        print("smooth_proba(avrg_Ng,logB,alpha_g) = smooth_proba({0:.2f},2^{1:.2f},{2:.4f}) = {3} = 2^{4:.4f}".format(float(self.avrg_log_Ng*log2),float(self.log2_B),float(self.alpha_g), float(smoothness_proba(self.avrg_log_Ng,self.logB,self.alpha_g)),float(log(smoothness_proba(self.avrg_log_Ng,self.logB,self.alpha_g),2.0))))
+        print("prod of smooth proba of avrg norms: {0} = 2^{1:.4f}".format(float(self.proba_rough),float(log(self.proba_rough,2.0))))
+        print("avrg_proba = {0} = 2^{1:.4f}".format(float(self.avrg_proba),float(log(self.avrg_proba,2.0))))
 
-        print("relations obtained: core_vol_sieving_space * proba = 2^{0:.4f}".format(float(log(self.rels,2))))
-        print("relations obtained: core_vol_sieving_space * (rough proba) = 2^{0:.4f}".format(float(log(self.rels_rough,2))))
+        print("relations obtained: core_vol_sieving_space * proba = 2^{0:.4f}".format(float(log(self.rels,2.0))))
+        print("relations obtained: core_vol_sieving_space * (rough proba) = 2^{0:.4f}".format(float(log(self.rels_rough,2.0))))
         print("#{{factor base}} =~ 2LogIntegral(B) = 2^{0:.4f}".format(self.log2_size_FB))
         print("#{{factor base}} ~ 2*B/log(B)       = 2^{0:.4f}".format(self.log2_size_FB_BD))
-        print("size matrix = size factor base / {0} = 2^{1:.4f} (worst case)".format(self.cst_filtering,self.log2_size_matrix))
+        if not self.barbulescu_duquesne:
+            print("size matrix = size factor base / {0} = 2^{1:.4f} (worst case)".format(self.cst_filtering,self.log2_size_matrix))
+        else:
+            print("size matrix = size factor base / (log_2 B = {0}) = 2^{1:.4f} (worst case)".format(self.log2_B,self.log2_size_matrix))
         
         if not self.count_sieving:
             print("time relcol = 2^{0:.4f}".format(self.log2_time_relcol))
@@ -622,9 +693,12 @@ class Simulation_TNFS():
                 if self.log2_rels >= self.log2_size_FB - 0.8:
                     cost = self.cost+1
                 else:
-                    log2_cost_sieving = RR(log(self.loglogB,2))
-                    log2_minimal_vol_sieving_space = RR(self.log2_size_FB - log(self.inv_zeta_Kh,2) - log(self.avrg_proba,2))
-                    log2_minimal_time_relcol = RR(log2_minimal_vol_sieving_space - RR(log(self.aut,2)) + log2_cost_sieving)
+                    if self.count_sieving:
+                        log2_cost_sieving = RR(log(self.loglogB,2.0))
+                    else:
+                        log2_cost_sieving = float(0.0)
+                    log2_minimal_vol_sieving_space = RR(self.log2_size_FB - log(self.inv_zeta_Kh,2.0) - log(self.avrg_proba,2.0))
+                    log2_minimal_time_relcol = RR(log2_minimal_vol_sieving_space - RR(log(self.aut,2.0)) + log2_cost_sieving)
                     cost = ceil(log2_minimal_time_relcol+1)
                 if cost == cost_old:
                     cost = cost+1
@@ -635,9 +709,12 @@ class Simulation_TNFS():
                 if self.log2_rels <= self.log2_size_FB + 1.5:
                     cost = self.cost-1
                 else:
-                    log2_cost_sieving = RR(log(self.loglogB,2))
-                    log2_minimal_vol_sieving_space = RR(self.log2_size_FB - log(self.inv_zeta_Kh,2) - log(self.avrg_proba,2))
-                    log2_minimal_time_relcol = RR(log2_minimal_vol_sieving_space - RR(log(self.aut,2)) + log2_cost_sieving)
+                    if self.count_sieving:
+                        log2_cost_sieving = RR(log(self.loglogB,2.0))
+                    else:
+                        log2_cost_sieving = float(0.0)
+                    log2_minimal_vol_sieving_space = RR(self.log2_size_FB - log(self.inv_zeta_Kh,2.0) - log(self.avrg_proba,2.0))
+                    log2_minimal_time_relcol = RR(log2_minimal_vol_sieving_space - RR(log(self.aut,2.0)) + log2_cost_sieving)
                     cost = ceil(log2_minimal_time_relcol+1)
                 if cost == cost_old:
                     cost = cost-1
@@ -684,7 +761,7 @@ class Simulation_TNFS_BD18():
                 b=random_vec(A,deg_h)
                 # this test is not in Barbulescu--Duquesne paper:
             if check_coprime:
-                coprime_ideal=(Kh.ideal(Kh(a))+Kh.ideal(Kh(b))==Kh.ideal(1))
+                coprime_ideal=gcd(gcd(a), gcd(b)) == 1 and (Kh.ideal(Kh(a))+Kh.ideal(Kh(b))==Kh.ideal(1))
                 while not coprime_ideal:
                     if not corrected_randint:
                         a=random_vec_bad_randint(A,deg_h)
@@ -695,8 +772,8 @@ class Simulation_TNFS_BD18():
                     coprime_ideal=(Kh.ideal(Kh(a))+Kh.ideal(Kh(b))==Kh.ideal(1))
             #a=(1+2*y) + (3+4*y)*x
             #h.resultant(a.resultant(f0,x),y)
-            log_Nf[i] = log(2)*int(ceil(log(RR(abs(h.resultant(f.resultant(Rxy([a,b]))))),2)))
-            log_Ng[i] = log(2)*int(ceil(log(RR(abs(h.resultant(g.resultant(Rxy([a,b]))))),2)))
+            log_Nf[i] = log2*int(ceil(log(RR(abs(h.resultant(f.resultant(Rxy([a,b]))))),2)))
+            log_Ng[i] = log2*int(ceil(log(RR(abs(h.resultant(g.resultant(Rxy([a,b]))))),2)))
             # the average norm is computed over ceil values of log_2(Nf), log_2(Ng)
             proba_f[i] = smoothness_proba(RR(log_Nf[i]), logB, alpha_f)
             proba_g[i] = smoothness_proba(RR(log_Ng[i]), logB, alpha_g)
@@ -751,26 +828,26 @@ class Simulation_TNFS_BD18():
     
         avrg_log_Nf_BD, avrg_log_Ng_BD, avrg_smooth_proba_f_BD, avrg_smooth_proba_g_BD, avrg_proba_BD, err_norm_f_BD, err_norm_g_BD, err_proba_f_BD, err_proba_g_BD, err_proba_BD, std_dev_norm_f_BD, std_dev_norm_g_BD, std_dev_proba_f_BD, std_dev_proba_g_BD, std_dev_proba_BD, min_log_Nf, max_log_Nf, min_log_Ng, max_log_Ng = statistics_BD18(h,f,g,Rxy,A,logB,samples,corrected_randint=corrected_randint,check_coprime=check_coprime)
     
-        print("avrg_log2_Nf_BD = {}\navrg_log2_Ng_BD = {}\navrg_smooth_proba_f_BD = {}\navrg_smooth_proba_g_BD = {}\navrg_proba_BD  = {}\n".format(float(avrg_log_Nf_BD/log(2)), float(avrg_log_Ng_BD/log(2)), avrg_smooth_proba_f_BD, avrg_smooth_proba_g_BD, avrg_proba_BD))
+        print("avrg_log2_Nf_BD = {}\navrg_log2_Ng_BD = {}\navrg_smooth_proba_f_BD = {}\navrg_smooth_proba_g_BD = {}\navrg_proba_BD  = {}\n".format(float(avrg_log_Nf_BD/log2), float(avrg_log_Ng_BD/log2), avrg_smooth_proba_f_BD, avrg_smooth_proba_g_BD, avrg_proba_BD))
         print("err_norm_f_BD  = {}\nerr_norm_g_BD  = {}\nerr_proba_f_BD = {}\nerr_proba_g_BD = {}\nerr_proba_BD   = {}\n".format( err_norm_f_BD, err_norm_g_BD, err_proba_f_BD, err_proba_g_BD, err_proba_BD))
         print("std_dev_norm_f_BD = {}\nstd_dev_norm_g_BD = {}\nstd_dev_proba_f_BD = {}\nstd_dev_proba_g_BD = {}\nstd_dev_proba_BD   = {}\n".format(std_dev_norm_f_BD, std_dev_norm_g_BD, std_dev_proba_f_BD, std_dev_proba_g_BD, std_dev_proba_BD))
     
         print("A = {0}, B = 2^{1:.3f}, deg_h={2}, #{{roots of unity/+-1}} = w = {3}, samples = {4}".format(A,float(log2_B),deg_h,w,samples))
-        print("vol_sieving_space = (2*A+1)^(2*{0})/(2*{1}) = 2^{2:.3f}".format(deg_h, w, float(log(vol_sieving_space, 2))))
-        print("avrg_log2_Nf = {0:.4f} (bits) avrg_log2_Ng = {1:.4f} (bits)".format(float(avrg_log_Nf_BD/log(2)), float(avrg_log_Ng_BD/log(2))))
-        print(" min_log2_Nf = {0:.4f} (bits)  min_log2_Ng = {1:.4f} (bits)".format(float(min_log_Nf/log(2)), float(min_log_Ng/log(2))))
-        print(" max_log2_Nf = {0:.4f} (bits)  max_log2_Ng = {1:.4f} (bits)".format(float(max_log_Nf/log(2)), float(max_log_Ng/log(2))))
+        print("vol_sieving_space = (2*A+1)^(2*{0})/(2*{1}) = 2^{2:.3f}".format(deg_h, w, float(log(vol_sieving_space, 2.0))))
+        print("avrg_log2_Nf = {0:.4f} (bits) avrg_log2_Ng = {1:.4f} (bits)".format(float(avrg_log_Nf_BD/log2), float(avrg_log_Ng_BD/log2)))
+        print(" min_log2_Nf = {0:.4f} (bits)  min_log2_Ng = {1:.4f} (bits)".format(float(min_log_Nf/log2), float(min_log_Ng/log2)))
+        print(" max_log2_Nf = {0:.4f} (bits)  max_log2_Ng = {1:.4f} (bits)".format(float(max_log_Nf/log2), float(max_log_Ng/log2)))
 
-        print("smooth_proba(avrg_Nf,logB,alpha_f) = {0} = 2^{1:.4f}".format(float(smoothness_proba(avrg_log_Nf_BD,logB,alpha_f)),float(log(smoothness_proba(avrg_log_Nf_BD,logB,alpha_f),2))))
-        print("smooth_proba(avrg_Ng,logB,alpha_g) = {0} = 2^{1:.4f}".format(float(smoothness_proba(avrg_log_Ng_BD,logB,alpha_g)),float(log(smoothness_proba(avrg_log_Ng_BD,logB,alpha_g),2))))
+        print("smooth_proba(avrg_Nf,logB,alpha_f) = {0} = 2^{1:.4f}".format(float(smoothness_proba(avrg_log_Nf_BD,logB,alpha_f)),float(log(smoothness_proba(avrg_log_Nf_BD,logB,alpha_f),2.0))))
+        print("smooth_proba(avrg_Ng,logB,alpha_g) = {0} = 2^{1:.4f}".format(float(smoothness_proba(avrg_log_Ng_BD,logB,alpha_g)),float(log(smoothness_proba(avrg_log_Ng_BD,logB,alpha_g),2.0))))
         proba_rough_BD = float(smoothness_proba(avrg_log_Nf_BD,logB,alpha_f)*smoothness_proba(avrg_log_Ng_BD,logB,alpha_g))
-        print("prod of smooth proba of avrg norms: {0} = 2^{1:.4f}".format(proba_rough_BD, float(log(proba_rough_BD,2))))
-        print("avrg_proba = {0} = 2^{1:.4f}".format(float(avrg_proba_BD),float(log(avrg_proba_BD,2))))
+        print("prod of smooth proba of avrg norms: {0} = 2^{1:.4f}".format(proba_rough_BD, float(log(proba_rough_BD,2.0))))
+        print("avrg_proba = {0} = 2^{1:.4f}".format(float(avrg_proba_BD),float(log(avrg_proba_BD,2.0))))
         
         rels_BD = core_vol_sieving_space*avrg_proba_BD
         rels_rough_BD = core_vol_sieving_space*proba_rough_BD
-        print("relations obtained: vol_sieving_space * proba = 2^{0:.4f}".format(float(log(rels_BD,2))))
-        print("relations obtained: vol_sieving_space * (rough proba) = 2^{0:.4f}".format(float(log(rels_rough_BD,2))))
+        print("relations obtained: vol_sieving_space * proba = 2^{0:.4f}".format(float(log(rels_BD,2.0))))
+        print("relations obtained: vol_sieving_space * (rough proba) = 2^{0:.4f}".format(float(log(rels_rough_BD,2.0))))
         print("#{{factor base}} =~ 2LogIntegral(B) = 2^{0:.4f}".format(log2_size_FB))
         print("#{{factor base}} ~ 2*B/log(B)       = 2^{0:.4f}".format(log2_size_FB_BD))
         print("size matrix = size factor base / log2_B = 2^{0:.4f} ".format(log2_size_matrix_BD))
-- 
GitLab