diff --git a/libff/algebra/curves/alt_bn128/README.md b/libff/algebra/curves/alt_bn128/README.md
new file mode 100644
index 00000000..7dd94cf2
--- /dev/null
+++ b/libff/algebra/curves/alt_bn128/README.md
@@ -0,0 +1,10 @@
+# Implementation of altbn128
+
+## Run the sage script to generate the curve parameters
+
+1. Make sure that you have [SageMath](https://www.sagemath.org/) installed
+
+2. Run:
+```bash
+sage alt_bn128.sage
+```
diff --git a/libff/algebra/curves/alt_bn128/alt_bn128.sage b/libff/algebra/curves/alt_bn128/alt_bn128.sage
new file mode 100644
index 00000000..0948ec63
--- /dev/null
+++ b/libff/algebra/curves/alt_bn128/alt_bn128.sage
@@ -0,0 +1,66 @@
+#!/usr/bin/env sage -python
+
+from sage.all import *
+import sys
+sys.path.append("../")
+import params_generator
+
+# Prime order of the subgroup we work in
+def r(x):
+    return 36*(x**4) + 36*(x**3) + 18*(x**2) + 6*x + 1
+
+# Prime used to generate the base finite field
+def q(x):
+    return 36*(x**4) + 36*(x**3) + 24*(x**2) + 6*x + 1
+
+# Compute G2 cofactor
+# See: Proposition 1, Section 3.3: https://eprint.iacr.org/2015/247.pdf
+def g2_h(x):
+    return 36*x^4+ 36*x^3+ 30*x^2+ 6*x + 1
+
+# Computes the order of G1, the safe subgroup of E/Fq
+def g1_order(curve_order):
+    decomposition = factor(curve_order)
+    # Factor returns the prime decomposition and orders prime
+    # factors from smaller to biggest
+    biggest_factor = decomposition[-1]
+    assert(biggest_factor[1] == 1)
+    return biggest_factor[0]
+
+def main():
+    print("Generating parameters for alt_bn128")
+    # Curve parameter
+    param = 0x44e992b44a6909f1
+
+    prime_r = r(param)
+    assert(prime_r == 21888242871839275222246405745257275088548364400416034343698204186575808495617)
+
+    prime_q = q(param)
+    assert(prime_q == 21888242871839275222246405745257275088696311157297823662689037894645226208583)
+    if (mod(prime_q, 6) != 1):
+        raise BaseException("Unexpected: q should be = 1 (mod 6). See: https://eprint.iacr.org/2007/390.pdf")
+
+    # Scalar field
+    print('prime_r = {}'.format(prime_r))
+    #params_generator.generate_libff_Fp_model_params(prime_r)
+    Fr = GF(prime_r)
+
+    # Base field
+    print('prime_q = {}'.format(prime_q))
+    #params_generator.generate_libff_Fp_model_params(prime_q)
+    Fq = GF(prime_q)
+
+    # E/Fq
+    curve = EllipticCurve(Fq, [0, 3])
+    curve_order = curve.order()
+
+    # Cofactors
+    h1 = curve_order // g1_order(curve_order)
+    # G1 cofactor should be 1
+    assert(h1 == 1)
+    print('h1 = {}'.format(h1))
+    h2 = g2_h(param)
+    print('h2 = {}'.format(h2))
+
+if __name__ == '__main__':
+    main()
diff --git a/libff/algebra/curves/alt_bn128/alt_bn128_g1.cpp b/libff/algebra/curves/alt_bn128/alt_bn128_g1.cpp
index 8be6340a..47309056 100755
--- a/libff/algebra/curves/alt_bn128/alt_bn128_g1.cpp
+++ b/libff/algebra/curves/alt_bn128/alt_bn128_g1.cpp
@@ -19,6 +19,7 @@ std::vector<size_t> alt_bn128_G1::fixed_base_exp_window_table;
 alt_bn128_G1 alt_bn128_G1::G1_zero = {};
 alt_bn128_G1 alt_bn128_G1::G1_one = {};
 bool alt_bn128_G1::initialized = false;
+bigint<alt_bn128_G1::h_limbs> alt_bn128_G1::h;
 
 alt_bn128_G1::alt_bn128_G1()
 {
@@ -361,6 +362,12 @@ alt_bn128_G1 alt_bn128_G1::dbl() const
     return alt_bn128_G1(X3, Y3, Z3);
 }
 
+alt_bn128_G1 alt_bn128_G1::mul_by_cofactor() const
+{
+    // Cofactor = 1
+    return *this;
+}
+
 bool alt_bn128_G1::is_well_formed() const
 {
     if (this->is_zero())
diff --git a/libff/algebra/curves/alt_bn128/alt_bn128_g1.hpp b/libff/algebra/curves/alt_bn128/alt_bn128_g1.hpp
index 232513ce..a49e702d 100755
--- a/libff/algebra/curves/alt_bn128/alt_bn128_g1.hpp
+++ b/libff/algebra/curves/alt_bn128/alt_bn128_g1.hpp
@@ -33,6 +33,11 @@ class alt_bn128_G1 {
     typedef alt_bn128_Fq base_field;
     typedef alt_bn128_Fr scalar_field;
 
+    // Cofactor
+    static const mp_size_t h_bitcount = 1;
+    static const mp_size_t h_limbs = (h_bitcount+GMP_NUMB_BITS-1)/GMP_NUMB_BITS;
+    static bigint<h_limbs> h;
+
     alt_bn128_Fq X, Y, Z;
 
     // using Jacobian coordinates
@@ -58,6 +63,7 @@ class alt_bn128_G1 {
     alt_bn128_G1 add(const alt_bn128_G1 &other) const;
     alt_bn128_G1 mixed_add(const alt_bn128_G1 &other) const;
     alt_bn128_G1 dbl() const;
+    alt_bn128_G1 mul_by_cofactor() const;
 
     bool is_well_formed() const;
 
diff --git a/libff/algebra/curves/alt_bn128/alt_bn128_g2.cpp b/libff/algebra/curves/alt_bn128/alt_bn128_g2.cpp
index 333c3421..e16c514b 100755
--- a/libff/algebra/curves/alt_bn128/alt_bn128_g2.cpp
+++ b/libff/algebra/curves/alt_bn128/alt_bn128_g2.cpp
@@ -19,6 +19,7 @@ std::vector<size_t> alt_bn128_G2::fixed_base_exp_window_table;
 alt_bn128_G2 alt_bn128_G2::G2_zero = {};
 alt_bn128_G2 alt_bn128_G2::G2_one = {};
 bool alt_bn128_G2::initialized = false;
+bigint<alt_bn128_G2::h_limbs> alt_bn128_G2::h;
 
 alt_bn128_G2::alt_bn128_G2()
 {
@@ -375,6 +376,11 @@ alt_bn128_G2 alt_bn128_G2::mul_by_q() const
                       (this->Z).Frobenius_map(1));
 }
 
+alt_bn128_G2 alt_bn128_G2::mul_by_cofactor() const
+{
+    return alt_bn128_G2::h * (*this);
+}
+
 bool alt_bn128_G2::is_well_formed() const
 {
     if (this->is_zero())
diff --git a/libff/algebra/curves/alt_bn128/alt_bn128_g2.hpp b/libff/algebra/curves/alt_bn128/alt_bn128_g2.hpp
index a22c06f3..498d30fb 100755
--- a/libff/algebra/curves/alt_bn128/alt_bn128_g2.hpp
+++ b/libff/algebra/curves/alt_bn128/alt_bn128_g2.hpp
@@ -34,6 +34,11 @@ class alt_bn128_G2 {
     typedef alt_bn128_Fq2 twist_field;
     typedef alt_bn128_Fr scalar_field;
 
+    // Cofactor
+    static const mp_size_t h_bitcount = 256;
+    static const mp_size_t h_limbs = (h_bitcount+GMP_NUMB_BITS-1)/GMP_NUMB_BITS;
+    static bigint<h_limbs> h;
+
     alt_bn128_Fq2 X, Y, Z;
 
     // using Jacobian coordinates
@@ -62,6 +67,7 @@ class alt_bn128_G2 {
     alt_bn128_G2 mixed_add(const alt_bn128_G2 &other) const;
     alt_bn128_G2 dbl() const;
     alt_bn128_G2 mul_by_q() const;
+    alt_bn128_G2 mul_by_cofactor() const;
 
     bool is_well_formed() const;
 
diff --git a/libff/algebra/curves/alt_bn128/alt_bn128_init.cpp b/libff/algebra/curves/alt_bn128/alt_bn128_init.cpp
index 32f6ba8b..3f5eff04 100755
--- a/libff/algebra/curves/alt_bn128/alt_bn128_init.cpp
+++ b/libff/algebra/curves/alt_bn128/alt_bn128_init.cpp
@@ -150,6 +150,9 @@ void init_alt_bn128_params()
                                     alt_bn128_Fq::one());
     alt_bn128_G1::initialized = true;
 
+    // Cofactor
+    alt_bn128_G1::h = bigint<alt_bn128_G1::h_limbs>("1");
+
     alt_bn128_G1::wnaf_window_table.resize(0);
     alt_bn128_G1::wnaf_window_table.push_back(11);
     alt_bn128_G1::wnaf_window_table.push_back(24);
@@ -215,6 +218,14 @@ void init_alt_bn128_params()
                                     alt_bn128_Fq2::one());
     alt_bn128_G2::initialized = true;
 
+    // Cofactor
+    // [Sage excerpt]
+    // See: https://eprint.iacr.org/2015/247.pdf
+    // u = 4965661367192848881
+    // h2 = (36 * u^4) + (36 * u^3) + (30 * u^2) + 6*u + 1; h2
+    // # 21888242871839275222246405745257275088844257914179612981679871602714643921549
+    alt_bn128_G2::h = bigint<alt_bn128_G2::h_limbs>("21888242871839275222246405745257275088844257914179612981679871602714643921549");
+
     alt_bn128_G2::wnaf_window_table.resize(0);
     alt_bn128_G2::wnaf_window_table.push_back(5);
     alt_bn128_G2::wnaf_window_table.push_back(15);
diff --git a/libff/algebra/curves/bn128/bn128_g1.cpp b/libff/algebra/curves/bn128/bn128_g1.cpp
index 6d1e8b2b..23eb1699 100755
--- a/libff/algebra/curves/bn128/bn128_g1.cpp
+++ b/libff/algebra/curves/bn128/bn128_g1.cpp
@@ -20,6 +20,7 @@ std::vector<size_t> bn128_G1::fixed_base_exp_window_table;
 bn128_G1 bn128_G1::G1_zero = {};
 bn128_G1 bn128_G1::G1_one = {};
 bool bn128_G1::initialized = false;
+bigint<bn128_G1::h_limbs> bn128_G1::h;
 
 bn::Fp bn128_G1::sqrt(const bn::Fp &el)
 {
@@ -337,6 +338,12 @@ bn128_G1 bn128_G1::dbl() const
     return result;
 }
 
+bn128_G1 bn128_G1::mul_by_cofactor() const
+{
+    // Cofactor = 1
+    return (*this);
+}
+
 bn128_G1 bn128_G1::zero()
 {
     return G1_zero;
diff --git a/libff/algebra/curves/bn128/bn128_g1.hpp b/libff/algebra/curves/bn128/bn128_g1.hpp
index 02a4ad0a..110dccce 100755
--- a/libff/algebra/curves/bn128/bn128_g1.hpp
+++ b/libff/algebra/curves/bn128/bn128_g1.hpp
@@ -37,6 +37,11 @@ class bn128_G1 {
     typedef bn128_Fq base_field;
     typedef bn128_Fr scalar_field;
 
+    // Cofactor
+    static const mp_size_t h_bitcount = 1;
+    static const mp_size_t h_limbs = (h_bitcount+GMP_NUMB_BITS-1)/GMP_NUMB_BITS;
+    static bigint<h_limbs> h;
+
     bn::Fp X, Y, Z;
     void fill_coord(bn::Fp coord[3]) const { coord[0] = this->X; coord[1] = this->Y; coord[2] = this->Z; return; };
 
@@ -62,6 +67,7 @@ class bn128_G1 {
     bn128_G1 add(const bn128_G1 &other) const;
     bn128_G1 mixed_add(const bn128_G1 &other) const;
     bn128_G1 dbl() const;
+    bn128_G1 mul_by_cofactor() const;
 
     bool is_well_formed() const;
 
diff --git a/libff/algebra/curves/bn128/bn128_g2.cpp b/libff/algebra/curves/bn128/bn128_g2.cpp
index fa59d3ae..662db110 100755
--- a/libff/algebra/curves/bn128/bn128_g2.cpp
+++ b/libff/algebra/curves/bn128/bn128_g2.cpp
@@ -20,6 +20,7 @@ std::vector<size_t> bn128_G2::fixed_base_exp_window_table;
 bn128_G2 bn128_G2::G2_zero = {};
 bn128_G2 bn128_G2::G2_one = {};
 bool bn128_G2::initialized = false;
+bigint<bn128_G2::h_limbs> bn128_G2::h;
 
 bn::Fp2 bn128_G2::sqrt(const bn::Fp2 &el)
 {
@@ -337,6 +338,11 @@ bn128_G2 bn128_G2::dbl() const
     return result;
 }
 
+bn128_G2 bn128_G2::mul_by_cofactor() const
+{
+    return bn128_G2::h * (*this);
+}
+
 bool bn128_G2::is_well_formed() const
 {
     if (this->is_zero())
diff --git a/libff/algebra/curves/bn128/bn128_g2.hpp b/libff/algebra/curves/bn128/bn128_g2.hpp
index df06bf6d..54d68296 100755
--- a/libff/algebra/curves/bn128/bn128_g2.hpp
+++ b/libff/algebra/curves/bn128/bn128_g2.hpp
@@ -38,6 +38,11 @@ class bn128_G2 {
     typedef bn128_Fq base_field;
     typedef bn128_Fr scalar_field;
 
+    // Cofactor
+    static const mp_size_t h_bitcount = 256;
+    static const mp_size_t h_limbs = (h_bitcount+GMP_NUMB_BITS-1)/GMP_NUMB_BITS;
+    static bigint<h_limbs> h;
+
     bn::Fp2 X, Y, Z;
     void fill_coord(bn::Fp2 coord[3]) const { coord[0] = this->X; coord[1] = this->Y; coord[2] = this->Z; };
 
@@ -63,6 +68,7 @@ class bn128_G2 {
     bn128_G2 add(const bn128_G2 &other) const;
     bn128_G2 mixed_add(const bn128_G2 &other) const;
     bn128_G2 dbl() const;
+    bn128_G2 mul_by_cofactor() const;
 
     bool is_well_formed() const;
 
diff --git a/libff/algebra/curves/bn128/bn128_init.cpp b/libff/algebra/curves/bn128/bn128_init.cpp
index 9b5f9eb7..caa3e92a 100755
--- a/libff/algebra/curves/bn128/bn128_init.cpp
+++ b/libff/algebra/curves/bn128/bn128_init.cpp
@@ -108,6 +108,9 @@ void init_bn128_params()
 
     bn128_G1::initialized = true;
 
+    // Cofactor
+    bn128_G1::h = bigint<bn128_G1::h_limbs>("1");
+
     bn128_G1::wnaf_window_table.resize(0);
     bn128_G1::wnaf_window_table.push_back(10);
     bn128_G1::wnaf_window_table.push_back(24);
@@ -173,6 +176,9 @@ void init_bn128_params()
 
     bn128_G2::initialized = true;
 
+    // Cofactor
+    bn128_G2::h = bigint<bn128_G2::h_limbs>("21888242871839275222246405745257275088844257914179612981679871602714643921549");
+
     bn128_G2::wnaf_window_table.resize(0);
     bn128_G2::wnaf_window_table.push_back(7);
     bn128_G2::wnaf_window_table.push_back(18);
diff --git a/libff/algebra/curves/mnt/README.md b/libff/algebra/curves/mnt/README.md
new file mode 100644
index 00000000..2ee942df
--- /dev/null
+++ b/libff/algebra/curves/mnt/README.md
@@ -0,0 +1,10 @@
+# Implementation of the MNT4/6 cycle
+
+## Run the sage script to generate the curve parameters
+
+1. Make sure that you have [SageMath](https://www.sagemath.org/) installed
+
+2. Run:
+```bash
+sage mnt.sage
+```
diff --git a/libff/algebra/curves/mnt/mnt.sage b/libff/algebra/curves/mnt/mnt.sage
new file mode 100644
index 00000000..dd15b393
--- /dev/null
+++ b/libff/algebra/curves/mnt/mnt.sage
@@ -0,0 +1,104 @@
+#!/usr/bin/env sage -python
+
+from sage.all import *
+import sys
+sys.path.append("../")
+import params_generator
+
+# Computes the order of G1, the safe subgroup of E/Fq
+def g1_order(curve_order):
+    decomposition = factor(curve_order)
+    # Factor returns the prime decomposition and orders prime
+    # factors from smaller to biggest
+    biggest_factor = decomposition[-1]
+    assert(biggest_factor[1] == 1)
+    return biggest_factor[0]
+
+def mnt4_298_params():
+    print("Generating parameters for MNT4-298")
+    # Scalar field
+    prime_r = 475922286169261325753349249653048451545124878552823515553267735739164647307408490559963137
+    params_generator.generate_libff_Fp_model_params(prime_r)
+    Fr = GF(prime_r)
+
+    # Base field characteristic
+    prime_q = 475922286169261325753349249653048451545124879242694725395555128576210262817955800483758081
+    params_generator.generate_libff_Fp_model_params(prime_q)
+    Fq = GF(prime_q)
+
+    # E/Fq
+    coeff_a = 2
+    coeff_b = 423894536526684178289416011533888240029318103673896002803341544124054745019340795360841685
+    curve = EllipticCurve(Fq, [coeff_a, coeff_b])
+
+    # E'/Fq2
+    non_residue = Fq(17)
+    Fqx.<j> = PolynomialRing(Fq, 'j')
+    assert(Fqx(j^2 + non_residue).is_irreducible())
+    Fq2.<u> = GF(prime_q^2, modulus=j^2 - non_residue)
+
+    twist_coeff_a = Fq2(coeff_a * non_residue)
+    twist_coeff_b = Fq2(coeff_b * non_residue * u)
+    twist = EllipticCurve(Fq2, [twist_coeff_a, twist_coeff_b])
+
+    # Cofactors
+    curve_order = curve.order()
+    print('curve_order = {}'.format(curve_order))
+    twist_order = twist.order()
+    print('twist_order = {}'.format(twist_order))
+
+    g1_h = curve_order // g1_order(curve_order)
+    assert(g1_h == 1)
+    assert(curve_order == prime_r)
+    print('g1_h = {}'.format(g1_h))
+    g2_h = twist_order // g1_order(curve_order)
+    print('g2_h = {}'.format(g2_h))
+
+def mnt6_298_params():
+    print("Generating parameters for MNT6-298")
+    # Scalar field
+    prime_r = 475922286169261325753349249653048451545124879242694725395555128576210262817955800483758081
+    params_generator.generate_libff_Fp_model_params(prime_r)
+    Fr = GF(prime_r)
+
+    # Base field characteristic
+    prime_q = 475922286169261325753349249653048451545124878552823515553267735739164647307408490559963137
+    params_generator.generate_libff_Fp_model_params(prime_q)
+    Fq = GF(prime_q)
+
+    # E/Fq
+    coeff_a = 11
+    coeff_b = 106700080510851735677967319632585352256454251201367587890185989362936000262606668469523074
+    curve = EllipticCurve(Fq, [coeff_a, coeff_b])
+
+    # E'/Fq3
+    non_residue = Fq(5)
+    print('non_residue = {}'.format(non_residue))
+    Fqx.<j> = PolynomialRing(Fq, 'j')
+    assert(Fqx(j^3 + non_residue).is_irreducible())
+    Fq3.<u> = GF(prime_q^3, modulus=j^3 - non_residue)
+
+    twist_coeff_a = Fq3(coeff_a * u^2)
+    twist_coeff_b = Fq3(coeff_b * non_residue)
+
+    twist = EllipticCurve(Fq3, [twist_coeff_a, twist_coeff_b])
+
+    # Cofactors
+    curve_order = curve.order()
+    print('curve_order = {}'.format(curve_order))
+    twist_order = twist.order()
+    print('twist_order = {}'.format(twist_order))
+
+    g1_h = curve_order // g1_order(curve_order)
+    assert(g1_h == 1)
+    assert(curve_order == prime_r)
+    print('g1_h = {}'.format(g1_h))
+    g2_h = twist_order // g1_order(curve_order)
+    print('g2_h = {}'.format(g2_h))
+
+def main():
+    mnt4_298_params()
+    mnt6_298_params()
+
+if __name__ == '__main__':
+    main()
diff --git a/libff/algebra/curves/mnt/mnt4/mnt4_g1.cpp b/libff/algebra/curves/mnt/mnt4/mnt4_g1.cpp
index 0863c778..4245beca 100755
--- a/libff/algebra/curves/mnt/mnt4/mnt4_g1.cpp
+++ b/libff/algebra/curves/mnt/mnt4/mnt4_g1.cpp
@@ -27,6 +27,7 @@ mnt4_G1 mnt4_G1::G1_one = {};
 bool mnt4_G1::initialized = false;
 mnt4_Fq mnt4_G1::coeff_a;
 mnt4_Fq mnt4_G1::coeff_b;
+bigint<mnt4_G1::h_limbs> mnt4_G1::h;
 
 mnt4_G1::mnt4_G1()
 {
@@ -349,6 +350,12 @@ mnt4_G1 mnt4_G1::dbl() const
     }
 }
 
+mnt4_G1 mnt4_G1::mul_by_cofactor() const
+{
+    // Cofactor = 1
+    return (*this);
+}
+
 bool mnt4_G1::is_well_formed() const
 {
     if (this->is_zero())
diff --git a/libff/algebra/curves/mnt/mnt4/mnt4_g1.hpp b/libff/algebra/curves/mnt/mnt4/mnt4_g1.hpp
index 4297cbf2..5555badd 100755
--- a/libff/algebra/curves/mnt/mnt4/mnt4_g1.hpp
+++ b/libff/algebra/curves/mnt/mnt4/mnt4_g1.hpp
@@ -40,6 +40,11 @@ class mnt4_G1 {
     typedef mnt4_Fq base_field;
     typedef mnt4_Fr scalar_field;
 
+    // Cofactor
+    static const mp_size_t h_bitcount = 1;
+    static const mp_size_t h_limbs = (h_bitcount+GMP_NUMB_BITS-1)/GMP_NUMB_BITS;
+    static bigint<h_limbs> h;
+
     mnt4_Fq X, Y, Z;
 
     // using projective coordinates
@@ -66,6 +71,7 @@ class mnt4_G1 {
     mnt4_G1 add(const mnt4_G1 &other) const;
     mnt4_G1 mixed_add(const mnt4_G1 &other) const;
     mnt4_G1 dbl() const;
+    mnt4_G1 mul_by_cofactor() const;
 
     bool is_well_formed() const;
 
diff --git a/libff/algebra/curves/mnt/mnt4/mnt4_g2.cpp b/libff/algebra/curves/mnt/mnt4/mnt4_g2.cpp
index 1f627c6e..114d21e0 100755
--- a/libff/algebra/curves/mnt/mnt4/mnt4_g2.cpp
+++ b/libff/algebra/curves/mnt/mnt4/mnt4_g2.cpp
@@ -28,6 +28,7 @@ mnt4_Fq2 mnt4_G2::coeff_b;
 mnt4_G2 mnt4_G2::G2_zero = {};
 mnt4_G2 mnt4_G2::G2_one = {};
 bool mnt4_G2::initialized;
+bigint<mnt4_G2::h_limbs> mnt4_G2::h;
 
 mnt4_Fq2 mnt4_G2::mul_by_a(const mnt4_Fq2 &elt)
 {
@@ -372,6 +373,11 @@ mnt4_G2 mnt4_G2::mul_by_q() const
                    (this->Z).Frobenius_map(1));
 }
 
+mnt4_G2 mnt4_G2::mul_by_cofactor() const
+{
+    return mnt4_G2::h * (*this);
+}
+
 bool mnt4_G2::is_well_formed() const
 {
     if (this->is_zero())
diff --git a/libff/algebra/curves/mnt/mnt4/mnt4_g2.hpp b/libff/algebra/curves/mnt/mnt4/mnt4_g2.hpp
index 24cb6a7a..7bf16988 100755
--- a/libff/algebra/curves/mnt/mnt4/mnt4_g2.hpp
+++ b/libff/algebra/curves/mnt/mnt4/mnt4_g2.hpp
@@ -42,6 +42,11 @@ class mnt4_G2 {
     typedef mnt4_Fq2 twist_field;
     typedef mnt4_Fr scalar_field;
 
+    // Cofactor
+    static const mp_size_t h_bitcount = 298;
+    static const mp_size_t h_limbs = (h_bitcount+GMP_NUMB_BITS-1)/GMP_NUMB_BITS;
+    static bigint<h_limbs> h;
+
     mnt4_Fq2 X, Y, Z;
 
     // using projective coordinates
@@ -71,6 +76,7 @@ class mnt4_G2 {
     mnt4_G2 mixed_add(const mnt4_G2 &other) const;
     mnt4_G2 dbl() const;
     mnt4_G2 mul_by_q() const;
+    mnt4_G2 mul_by_cofactor() const;
 
     bool is_well_formed() const;
 
diff --git a/libff/algebra/curves/mnt/mnt4/mnt4_init.cpp b/libff/algebra/curves/mnt/mnt4/mnt4_init.cpp
index dbc20d7a..5a0e3126 100755
--- a/libff/algebra/curves/mnt/mnt4/mnt4_init.cpp
+++ b/libff/algebra/curves/mnt/mnt4/mnt4_init.cpp
@@ -139,6 +139,9 @@ void init_mnt4_params()
                               mnt4_Fq::one());
     mnt4_G1::initialized = true;
 
+    // Cofactor
+    mnt4_G1::h = bigint<mnt4_G1::h_limbs>("1");
+
     mnt4_G1::wnaf_window_table.resize(0);
     mnt4_G1::wnaf_window_table.push_back(11);
     mnt4_G1::wnaf_window_table.push_back(24);
@@ -203,6 +206,9 @@ void init_mnt4_params()
                               mnt4_Fq2::one());
     mnt4_G2::initialized = true;
 
+    // Cofactor
+    mnt4_G2::h = bigint<mnt4_G2::h_limbs>("475922286169261325753349249653048451545124879932565935237842521413255878328503110407553025");
+
     mnt4_G2::wnaf_window_table.resize(0);
     mnt4_G2::wnaf_window_table.push_back(5);
     mnt4_G2::wnaf_window_table.push_back(15);
diff --git a/libff/algebra/curves/mnt/mnt6/mnt6_g1.cpp b/libff/algebra/curves/mnt/mnt6/mnt6_g1.cpp
index 9b527e47..7cb55459 100755
--- a/libff/algebra/curves/mnt/mnt6/mnt6_g1.cpp
+++ b/libff/algebra/curves/mnt/mnt6/mnt6_g1.cpp
@@ -27,6 +27,7 @@ mnt6_G1 mnt6_G1::G1_one = {};
 bool mnt6_G1::initialized = false;
 mnt6_Fq mnt6_G1::coeff_a;
 mnt6_Fq mnt6_G1::coeff_b;
+bigint<mnt6_G1::h_limbs> mnt6_G1::h;
 
 mnt6_G1::mnt6_G1()
 {
@@ -349,6 +350,12 @@ mnt6_G1 mnt6_G1::dbl() const
     }
 }
 
+mnt6_G1 mnt6_G1::mul_by_cofactor() const
+{
+    // Cofactor = 1
+    return (*this);
+}
+
 bool mnt6_G1::is_well_formed() const
 {
     if (this->is_zero())
diff --git a/libff/algebra/curves/mnt/mnt6/mnt6_g1.hpp b/libff/algebra/curves/mnt/mnt6/mnt6_g1.hpp
index 06e7b1ee..b9f281d3 100755
--- a/libff/algebra/curves/mnt/mnt6/mnt6_g1.hpp
+++ b/libff/algebra/curves/mnt/mnt6/mnt6_g1.hpp
@@ -40,6 +40,11 @@ class mnt6_G1 {
     typedef mnt6_Fq base_field;
     typedef mnt6_Fr scalar_field;
 
+    // Cofactor
+    static const mp_size_t h_bitcount = 1;
+    static const mp_size_t h_limbs = (h_bitcount+GMP_NUMB_BITS-1)/GMP_NUMB_BITS;
+    static bigint<h_limbs> h;
+
     mnt6_Fq X, Y, Z;
 
     // using projective coordinates
@@ -66,6 +71,7 @@ class mnt6_G1 {
     mnt6_G1 add(const mnt6_G1 &other) const;
     mnt6_G1 mixed_add(const mnt6_G1 &other) const;
     mnt6_G1 dbl() const;
+    mnt6_G1 mul_by_cofactor() const;
 
     bool is_well_formed() const;
 
diff --git a/libff/algebra/curves/mnt/mnt6/mnt6_g2.cpp b/libff/algebra/curves/mnt/mnt6/mnt6_g2.cpp
index e1412cac..a020c2a6 100755
--- a/libff/algebra/curves/mnt/mnt6/mnt6_g2.cpp
+++ b/libff/algebra/curves/mnt/mnt6/mnt6_g2.cpp
@@ -28,6 +28,7 @@ mnt6_Fq3 mnt6_G2::coeff_b;
 mnt6_G2 mnt6_G2::G2_zero = {};
 mnt6_G2 mnt6_G2::G2_one = {};
 bool mnt6_G2::initialized = false;
+bigint<mnt6_G2::h_limbs> mnt6_G2::h;
 
 mnt6_G2::mnt6_G2()
 {
@@ -378,6 +379,11 @@ mnt6_G2 mnt6_G2::mul_by_q() const
                    (this->Z).Frobenius_map(1));
 }
 
+mnt6_G2 mnt6_G2::mul_by_cofactor() const
+{
+    return mnt6_G2::h * (*this);
+}
+
 bool mnt6_G2::is_well_formed() const
 {
     if (this->is_zero())
diff --git a/libff/algebra/curves/mnt/mnt6/mnt6_g2.hpp b/libff/algebra/curves/mnt/mnt6/mnt6_g2.hpp
index 8c8d7dfa..48d276b0 100755
--- a/libff/algebra/curves/mnt/mnt6/mnt6_g2.hpp
+++ b/libff/algebra/curves/mnt/mnt6/mnt6_g2.hpp
@@ -42,6 +42,11 @@ class mnt6_G2 {
     typedef mnt6_Fq3 twist_field;
     typedef mnt6_Fr scalar_field;
 
+    // Cofactor
+    static const mp_size_t h_bitcount = 596;
+    static const mp_size_t h_limbs = (h_bitcount+GMP_NUMB_BITS-1)/GMP_NUMB_BITS;
+    static bigint<h_limbs> h;
+
     mnt6_Fq3 X, Y, Z;
 
     // using projective coordinates
@@ -71,6 +76,7 @@ class mnt6_G2 {
     mnt6_G2 mixed_add(const mnt6_G2 &other) const;
     mnt6_G2 dbl() const;
     mnt6_G2 mul_by_q() const;
+    mnt6_G2 mul_by_cofactor() const;
 
     bool is_well_formed() const;
 
diff --git a/libff/algebra/curves/mnt/mnt6/mnt6_init.cpp b/libff/algebra/curves/mnt/mnt6/mnt6_init.cpp
index 61599fa9..65555241 100755
--- a/libff/algebra/curves/mnt/mnt6/mnt6_init.cpp
+++ b/libff/algebra/curves/mnt/mnt6/mnt6_init.cpp
@@ -150,6 +150,9 @@ void init_mnt6_params()
                               mnt6_Fq::one());
     mnt6_G1::initialized = true;
 
+    // Cofactor
+    mnt6_G1::h = bigint<mnt6_G1::h_limbs>("1");
+
     mnt6_G1::wnaf_window_table.resize(0);
     mnt6_G1::wnaf_window_table.push_back(11);
     mnt6_G1::wnaf_window_table.push_back(24);
@@ -215,6 +218,9 @@ void init_mnt6_params()
                               mnt6_Fq3::one());
     mnt6_G2::initialized = true;
 
+    // Cofactor
+    mnt6_G2::h = bigint<mnt6_G2::h_limbs>("226502022472576270196498690498308461791828762732602586162207535351960270082712694977333372361549082214519252261735048131889018501404377856786623430385820659037970876666767495659520");
+
     mnt6_G2::wnaf_window_table.resize(0);
     mnt6_G2::wnaf_window_table.push_back(5);
     mnt6_G2::wnaf_window_table.push_back(15);
diff --git a/libff/algebra/curves/params_generator.sage b/libff/algebra/curves/params_generator.sage
new file mode 100755
index 00000000..5790285c
--- /dev/null
+++ b/libff/algebra/curves/params_generator.sage
@@ -0,0 +1,63 @@
+#!/usr/bin/env sage -python
+
+from sage.all import *
+
+def generate_libff_Fp_model_params(prime):
+    num_bits = ceil(log(prime, 2))
+    print('num_bits = {}'.format(num_bits))
+
+    euler = (prime-1)/2
+    print('euler = {}'.format(euler))
+
+    factorization = factor(prime-1)
+    t = 0
+    term_2 = factorization[0]
+    counter = 0
+    if term_2[0] != 2:
+        raise BaseException("The prime decomposition doesn't have any factor 2."
+        "The 'high 2-adicity' requirement isn't respected")
+    while not(is_odd(t)):
+        s = term_2[1] - counter
+        t = (prime-1)/(2**s)
+        counter = counter + 1
+    print('s = {}'.format(s))
+    is_odd(t); print('t = {}'.format(t))
+
+    t_minus_1_over_2 = (t-1)/2
+    print('t_minus_1_over_2 = {}'.format(t_minus_1_over_2))
+
+    multiplicative_generator = primitive_root(prime)
+    print('multiplicative_generator = {}'.format(multiplicative_generator))
+
+    root_of_unity = pow(multiplicative_generator, t, prime)
+    print('root_of_unity = {}'.format(root_of_unity))
+
+    nqr = least_quadratic_nonresidue(prime)
+    print('nqr = {}'.format(nqr))
+
+    nqr_to_t = pow(nqr, t, prime)
+    print('nqr_to_t = {}'.format(nqr_to_t))
+
+    word_len_64_bits = 64
+    W_64_bits = 2**(word_len_64_bits)
+    k_64_bits = ceil(num_bits/word_len_64_bits)
+    print('k_64_bits (nb limbs) = {}'.format(k_64_bits))
+    R_64_bits = mod(W_64_bits**k_64_bits, prime); R_64_bits
+    Rsquared_64_bits = R_64_bits**2
+    print('Rsquared_64_bits = {}'.format(Rsquared_64_bits))
+    Rcubed_64_bits = R_64_bits**3
+    print('Rcubed_64_bits = {}'.format(Rcubed_64_bits))
+    inv_64_bits = hex(int(mod((1/-prime), W_64_bits)))
+    print('inv_64_bits = {}'.format(inv_64_bits))
+
+    word_len_32_bits = 32
+    W_32_bits = 2**(32)
+    k_32_bits = ceil(num_bits/word_len_32_bits)
+    print('k_32_bits (nb limbs) = {}'.format(k_32_bits))
+    R_32_bits = mod(W_32_bits**k_32_bits, prime); R_32_bits
+    Rsquared_32_bits = R_32_bits**2
+    print('Rsquared_32_bits = {}'.format(Rsquared_32_bits))
+    Rcubed_32_bits = R_32_bits**3
+    print('Rcubed_32_bits = {}'.format(Rcubed_32_bits))
+    inv_32_bits = hex(int(mod(1/-prime, W_32_bits)))
+    print('inv_32_bits = {}'.format(inv_32_bits))
\ No newline at end of file