tensor_codes.cpp

Go to the documentation of this file.

/*
 * tensor_codes.cpp
 *
 *  Created on: Jun 24, 2021
 *      Author: betten
 */
 
 
 
 
#include "foundations.h"
 
using namespace std;
 
 
 
namespace orbiter {
namespace layer1_foundations {
namespace coding_theory {
 
 
void coding_theory_domain::twisted_tensor_product_codes(
    int *&H_subfield, int &m, int &n,
    field_theory::finite_field *F, field_theory::finite_field *f,
    int f_construction_A, int f_hyperoval,
    int f_construction_B, int verbose_level)
{
    int f_v = (verbose_level >= 1);
    int f_vv = (verbose_level >= 2);
    int index;
    int exponents[9];
    int *M;
    //int *H_subfield;
    int *C;
    int *C_inv;
 
    int q = f->q;
    int q2;
    int Q;
    //int m, n;
    int r;
    int beta, beta_q;
    int f_elements_exponential = TRUE;
    string symbol_for_print;
    string symbol_for_print_subfield;
 
 
    symbol_for_print.assign("\\alpha");
    symbol_for_print_subfield.assign("\\omega");
 
 
    if (f_v) {
        cout << "twisted_tensor_product_codes" << endl;
        cout << "f_construction_A=" << f_construction_A << endl;
        cout << "f_hyperoval=" << f_hyperoval << endl;
        cout << "f_construction_B=" << f_construction_B << endl;
        }
 
 
    q2 = q * q;
    Q = 0;
    if (f_construction_A) {
        Q = q2;
        }
    else if (f_construction_B) {
        Q = q2 * q;
        }
    index = (Q - 1) / (q - 1);
 
    if (Q != F->q) {
        cout << "twisted_tensor_product_codes Q != F->q" << endl;
        exit(1);
        }
 
 
    if (f_vv) {
        cout << "q = " << q << endl;
        cout << "Q = " << Q << endl;
        cout << "index = " << index << endl;
        }
 
#if 0
    F.init_override_polynomial(Q, override_poly_Q, verbose_level - 2);
 
    if (f_vv) {
        cout << "field of order " << Q << " initialized" << endl;
        }
 
    f.init_override_polynomial(q, override_poly_q, verbose_level - 2);
 
    if (f_vv) {
        cout << "field of order " << q << " initialized" << endl;
        cout << "index = " << index << endl;
        }
#endif
 
    F->compute_subfields(verbose_level - 2);
 
 
    create_matrix_M(
            M,
            F, f,
            m, n, beta, r, exponents,
            f_construction_A, f_hyperoval, f_construction_B,
            f_elements_exponential, symbol_for_print,
            verbose_level - 2);
 
    beta_q = F->power(beta, q);
 
    if (f_vv) {
        cout << "twisted_tensor_product_codes after create_matrix_M" << endl;
        cout << "m = " << m << endl;
        cout << "n = " << n << endl;
        cout << "Q = " << Q << endl;
        cout << "q2 = " << q2 << endl;
        cout << "beta = " << beta << endl;
        cout << "beta_q = " << beta_q << endl;
        cout << "Exponents: ";
        Int_vec_print(cout, exponents, m);
        cout << endl;
        }
 
    if (f_vv) {
        cout << "twisted_tensor_product_codes: M:" << endl;
        Int_vec_print_integer_matrix_width(cout, M, m, n, n, 2);
 
        F->latex_matrix(cout, f_elements_exponential, symbol_for_print, M, m, n);
        }
 
 
 
 
#if 0
    for (j = 0; j < n; j++) {
        PG_element_normalize(F, M + j, n, m);
        }
    cout << "column normalized M:" << endl;
    print_integer_matrix_width(cout, M, m, n, n, 2);
#endif
 
 
 
 
    C = NEW_int(m * m);
    C_inv = NEW_int(m * m);
    H_subfield = NEW_int(m * n);
 
 
    create_matrix_H_subfield(F, f,
        H_subfield, C, C_inv, M, m, n, beta, beta_q,
        f_elements_exponential, symbol_for_print, symbol_for_print_subfield,
        f_construction_A, f_hyperoval, f_construction_B,
        verbose_level - 2);
 
 
    if (f_v) {
        cout << "twisted_tensor_product_codes: after create_matrix_H_subfield" << endl;
        cout << "H_subfield:" << endl;
        Int_vec_print_integer_matrix_width(cout, H_subfield, m, n, n, 2);
        f->latex_matrix(cout, f_elements_exponential, symbol_for_print_subfield, H_subfield, m, n);
        }
 
    FREE_int(M);
    FREE_int(C);
    FREE_int(C_inv);
 
}
 
 
void coding_theory_domain::create_matrix_M(
    int *&M,
    field_theory::finite_field *F, field_theory::finite_field *f,
    int &m, int &n, int &beta, int &r, int *exponents,
    int f_construction_A, int f_hyperoval, int f_construction_B,
    int f_elements_exponential, std::string &symbol_for_print,
    int verbose_level)
// int exponents[9]
{
    int f_v = (verbose_level >= 1);
    //int f_vv = (verbose_level >= 2);
    int i, j, t, q, Q, q2;
 
    q = f->q;
    q2 = q * q;
 
 
    if (f_construction_A) {
 
#if 0
        //q = 2; override_poly_Q = ""; override_poly = "";
        //q = 3; override_poly_Q = ""; override_poly = "";
        q = 4; override_poly_Q = "19"; override_poly = "7";
            // F_256  generated by X^8 + X^4 + X^3 + X^2 + 1
            // F_16 generated by X^4+X+1 = 19
            // F_4 generated by X^2+X+1 = 7
        //q = 5; override_poly_Q = "47"; override_poly = "";
            // F_625  generated by X^4 + X^3 + 3X + 2
            // F_25  generated by X^2 + 4X + 2  = 47
        //q = 7; override_poly_Q = ""; override_poly = "";
        //q = 8; override_poly_Q = "97"; override_poly = "11";
            // F_4096 generated by x^12+x^6+x^4+x+1
            // F_64 generated by X^6+X^5+1 = 97
            // F_8 generated by X^3+X+1 = 11
        //q = 9; override_poly_Q = ""; override_poly = "17";
#endif
        Q = q2;
        beta = q;
        m = 9;
        if (f_hyperoval) {
            n = Q + 2;
            }
        else {
            n = Q + 1;
            }
        r = 5;
        if (q == 4)
            r = 7;
        if (q == 3)
            r = 9;
        // 3 orbits of length 1: 0, q+1, 2q+2
        exponents[0] = 0;
        exponents[1] = q + 1;
        exponents[2] = 2 * q + 2;
        //orbit (q,1)
        exponents[3] = q;
        exponents[4] = 1;
        // orbit (2q, 2)
        exponents[5] = 2 * q;
        exponents[6] = 2;
        // orbit (2q+1, q+2)
        exponents[7] = 2 * q + 1;
        exponents[8] = q + 2;
        //exponents[0] = 0;
        //exponents[1] = q;
        //exponents[2] = 2 * q;
        //exponents[3] = 1;
        //exponents[4] = q + 1;
        //exponents[5] = 2 * q + 1;
        //exponents[6] = 2;
        //exponents[7] = q + 2;
        //exponents[8] = 2 * q + 2;
        }
    else if (f_construction_B) {
 
#if 0
        //q = 2; override_poly_Q = ""; override_poly = ""; r = 9;
        //q = 3; override_poly_Q = ""; override_poly = ""; r = 5;
        q = 4; override_poly_Q = ""; override_poly = "7"; r = 4;
            // F_4096  generated by X^8 + X^4 + X^3 + X^2 + 1 = 4096
        //q = 5; override_poly_Q = ""; override_poly = ""; r = 4;
        //q = 7; override_poly_Q = ""; override_poly = ""; r = 4;
#endif
 
        beta = q;
        Q = q2 * q;
        m = 8;
        n = Q + 1;
 
        exponents[0] = 0;
        exponents[1] = q2 + q + 1;
        exponents[2] = 1;
        exponents[3] = q;
        exponents[4] = q2;
        exponents[5] = q + 1;
        exponents[6] = q2 + q;
        exponents[7] = q2 + 1;
        //exponents[0] = 0;
        //exponents[1] = q;
        //exponents[2] = 1;
        //exponents[3] = q + 1;
        //exponents[4] = q2;
        //exponents[5] = q2 + q;
        //exponents[6] = q2 + 1;
        //exponents[7] = q2 + q + 1;
        }
    else {
        cout << "coding_theory_domain::create_matrix_M please specify the construction using option -A or -B" << endl;
        exit(1);
        }
 
 
    // create matrix M:
    M = NEW_int(m * n);
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) {
            M[i * n + j] = 0;
            }
        }
    for (t = 0; t < Q; t++) {
        for (i = 0; i < m; i++) {
            M[i * n + t] = F->power(t, exponents[i]);
            }
        }
    if (f_construction_A) {
        M[2 * n + Q] = 1;
        if (f_hyperoval)
            M[1 * n + Q + 1] = 1;
        }
    else if (f_construction_B) {
        M[1 * n + Q] = 1;
        }
 
    if (f_v) {
        cout << "M:" << endl;
        Int_vec_print_integer_matrix_width(cout, M, m, n, n, 2);
 
        F->latex_matrix(cout, f_elements_exponential,
                symbol_for_print, M, m, n);
        }
 
 
    if (f_v) {
        int *all_one, *col_sum;
 
        all_one = NEW_int(n);
        col_sum = NEW_int(m);
        for (i = 0; i < n; i++)
            all_one[i] = 1;
        F->Linear_algebra->mult_matrix_matrix(M, all_one, col_sum, m, n, 1,
                0 /* verbose_level */);
        cout << "overall col_sum:" << endl;
        Int_vec_print_integer_matrix_width(cout, col_sum, m, 1, 1, 2);
        FREE_int(all_one);
        FREE_int(col_sum);
        }
 
 
}
 
 
 
void coding_theory_domain::create_matrix_H_subfield(
        field_theory::finite_field *F, field_theory::finite_field*f,
    int *H_subfield, int *C, int *C_inv, int *M,
    int m, int n, int beta, int beta_q,
    int f_elements_exponential, std::string &symbol_for_print,
    std::string &symbol_for_print_subfield,
    int f_construction_A, int f_hyperoval, int f_construction_B,
    int verbose_level)
{
    int f_v = (verbose_level >= 1);
    int f_vv = (verbose_level >= 2);
    int i, j;
    int q;
    //int *C;
    //int *C_inv;
    int *H;
    int *AA;
    linear_algebra::representation_theory_domain Rep;
 
    q = f->q;
 
 
    Rep.init(F, verbose_level);
    // matrix C is zero:
    H = NEW_int(m * n);
    AA = NEW_int(m * m);
    for (i = 0; i < m; i++) {
        for (j = 0; j < m; j++) {
            C[i * m + j] = 0;
            }
        }
 
 
    if (f_construction_A) {
        int nb_C_coeffs = 15;
        int k, aa;
        int C_coeffs[] = {
            0, 0, 1,
            1, 1, 1,
            2, 2, 1,
            3, 3, 1,
            3, 4, 1,
            4, 3, beta_q,
            4, 4, beta,
            5, 5, 1,
            5, 6, 1,
            6, 5, beta_q,
            6, 6, beta,
            7, 7, 1,
            7, 8, 1,
            8, 7, beta_q,
            8, 8, beta,
            };
        for (k = 0; k < nb_C_coeffs; k++) {
            i = C_coeffs[k * 3 + 0];
            j = C_coeffs[k * 3 + 1];
            aa = C_coeffs[k * 3 + 2];
            C[i * m + j] = aa;
            }
        }
    else if (f_construction_B) {
        int nb_C_coeffs = 20;
        int k, aa;
        int C_coeffs[] = {
            0, 0, 1,
            1, 1, 1,
            2, 2, 1,
            2, 3, 1,
            2, 4, 1,
            3, 2, beta,
            3, 3, beta_q,
            3, 4, Rep.beta_trinomial(q, beta, 1, 0, 0),
            4, 2, Rep.beta_trinomial(q, beta, 0, 0, 2),
            4, 3, Rep.beta_trinomial(q, beta, 0, 2, 0),
            4, 4, Rep.beta_trinomial(q, beta, 2, 0, 0),
            5, 5, 1,
            5, 6, 1,
            5, 7, 1,
            6, 5, Rep.beta_trinomial(q, beta, 0, 1, 1),
            6, 6, Rep.beta_trinomial(q, beta, 1, 1, 0),
            6, 7, Rep.beta_trinomial(q, beta, 1, 0, 1),
            7, 5, Rep.beta_trinomial(q, beta, 0, 2, 2),
            7, 6, Rep.beta_trinomial(q, beta, 2, 2, 0),
            7, 7, Rep.beta_trinomial(q, beta, 2, 0, 2),
            };
        for (k = 0; k < nb_C_coeffs; k++) {
            i = C_coeffs[k * 3 + 0];
            j = C_coeffs[k * 3 + 1];
            aa = C_coeffs[k * 3 + 2];
            C[i * m + j] = aa;
            }
        }
 
 
    if (f_v) {
        cout << "matrix C:" << endl;
        Int_vec_print_integer_matrix_width(cout, C, m, m, m, 2);
        F->latex_matrix(cout, f_elements_exponential,
                symbol_for_print, C, m, m);
        }
 
 
    F->Linear_algebra->invert_matrix(C, C_inv, m, 0 /* verbose_level */);
 
    if (f_vv) {
        cout << "C_inv:" << endl;
        Int_vec_print_integer_matrix_width(cout, C_inv, m, m, m, 2);
        }
 
    F->Linear_algebra->mult_matrix_matrix(C, C_inv, AA, m, m, m,
            0 /* verbose_level */);
 
    if (f_vv) {
        cout << "C * C_inv:" << endl;
        Int_vec_print_integer_matrix_width(cout, AA, m, m, m, 2);
        }
 
 
    F->Linear_algebra->mult_matrix_matrix(C, M, H, m, m, n,
            0 /* verbose_level */);
 
    if (f_v) {
        cout << "H = C * M:" << endl;
        Int_vec_print_integer_matrix_width(cout, H, m, n, n, 2);
        F->latex_matrix(cout, f_elements_exponential,
                symbol_for_print, H, m, n);
        }
 
 
#if 0
    rk = F.Gauss_int(M, FALSE /* f_special */, TRUE /* f_complete */, base_cols,
        FALSE /* f_P */, NULL, m /* m */, n /* n */, 0 /* Pn */,
        verbose_level - 2);
    cout << "has rank " << rk << endl;
#endif
 
    tt_field_reduction(*F, *f, m, n, H, H_subfield, verbose_level - 2);
 
    if (f_v) {
        cout << "H_subfield:" << endl;
        Int_vec_print_integer_matrix_width(cout, H_subfield, m, n, n, 2);
        f->latex_matrix(cout, f_elements_exponential,
                symbol_for_print_subfield, H_subfield, m, n);
        }
 
    FREE_int(H);
    FREE_int(AA);
}
 
 
 
void coding_theory_domain::tt_field_reduction(
        field_theory::finite_field &F, field_theory::finite_field &f,
        int m, int n, int *M, int *MM, int verbose_level)
{
    int f_v = (verbose_level >= 1);
    int f_vv = (verbose_level >= 2);
    int  i, j, a, z, b, c, Q, q;
    int index;
 
    Q = F.q;
    q = f.q;
    index = (Q - 1) / (q - 1);
    if (f_v) {
        cout << "field reduction, Q=" << Q
                << " q=" << q << " index=" << index << endl;
        }
    if (f_vv) {
        cout << "before:" << endl;
        Int_vec_print_integer_matrix_width(cout, M, m, n, n, 2);
        //print_integer_matrix(cout, M, m, n);
        cout << endl;
        F.print_integer_matrix_zech(cout, M, m, n);
        cout << endl;
        }
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) {
            a = M[i * n + j];
            if (a == 0) {
                c = 0;
                }
            else {
                if (f.e == 1) {
                    if (a >= q) {
                        cout << "field reduction: element does not "
                                "lie in the subfield: " << a << endl;
                        exit(1);
                        }
                    c = a;
                    }
                else {
                    z = F.log_alpha(a);
                    b = z / index;
                    if (b * index != z) {
                        cout << "b * index != z" << endl;
                        exit(1);
                        }
                    c = f.alpha_power(b);
                    }
                }
            MM[i * n + j] = c;
            }
        }
    if (f_vv) {
        cout << "after:" << endl;
        Int_vec_print_integer_matrix_width(cout, MM, m, n, n, 2);
        //print_integer_matrix(cout, MM, m, n);
        cout << endl;
        f.print_integer_matrix_zech(cout, MM, m, n);
        cout << endl;
        }
}
 
 
 
 
 
//############################################### old stuff:
 
void coding_theory_domain::make_tensor_code_9dimensional_as_point_set(
        field_theory::finite_field *F,
    int *&the_set, int &length,
    int verbose_level)
{
    int f_v = (verbose_level >= 1);
    int f_hyperoval = FALSE;
    const char *override_poly = "";
    const char *override_poly_Q = "";
    int i, t, q;
    int *code;
 
    if (f_v) {
        cout << "make_tensor_code_9dimensional_as_point_set" << endl;
        }
    q = F->q;
    if (q == 2) {
        override_poly_Q = ""; override_poly = ""; f_hyperoval = FALSE;
        }
    else if (q == 3) {
        override_poly_Q = ""; override_poly = ""; f_hyperoval = FALSE;
        }
    else if (q == 4) {
        override_poly_Q = "19"; override_poly = "7"; f_hyperoval = FALSE;
        // F_256  generated by X^8 + X^4 + X^3 + X^2 + 1
        // F_16 generated by X^4+X+1 = 19
        // F_4 generated by X^2+X+1 = 7
        }
    else if (q == 5) {
        override_poly_Q = "47"; override_poly = ""; f_hyperoval = FALSE;
        // F_625  generated by X^4 + X^3 + 3X + 2
        // F_25  generated by X^2 + 4X + 2  = 47
        }
    else if (q == 7) {
        override_poly_Q = ""; override_poly = ""; f_hyperoval = FALSE;
        }
    else if (q == 8) {
        override_poly_Q = "97"; override_poly = "11"; f_hyperoval = TRUE;
        // F_4096 generated by x^12+x^6+x^4+x+1
        // F_64 generated by X^6+X^5+1 = 97
        // F_8 generated by X^3+X+1 = 11
        }
    else if (q == 9) {
        override_poly_Q = ""; override_poly = "17"; f_hyperoval = FALSE;
    }
 
    string poly_Q, poly_q;
 
    poly_Q.assign(override_poly_Q);
    poly_q.assign(override_poly);
 
    make_tensor_code_9_dimensional(q, poly_Q, poly_q,
        f_hyperoval, code, length, verbose_level - 1);
 
    the_set = NEW_int(length);
 
    int pt[9], rk;
 
    for (t = 0; t < length; t++) {
        for (i = 0; i < 9; i++) {
            pt[i] = code[i * length + t];
            }
        F->PG_element_rank_modified(pt, 1, 9, rk);
        the_set[t] = rk;
        }
    FREE_int(code);
    if (f_v) {
        cout << "make_tensor_code_9dimensional_as_point_set done" << endl;
        cout << "created the set: ";
        Int_vec_print(cout, the_set, length);
        cout << endl;
        }
}
 
void coding_theory_domain::make_tensor_code_9_dimensional(int q,
    std::string &override_poly_Q, std::string &override_poly,
    int f_hyperoval,
    int *&code, int &length,
    int verbose_level)
{
    field_theory::finite_field F;
    field_theory::finite_field f;
    algebra::rank_checker rc;
    int exponents[9];
    int *M;
    int *C;
    int *C_inv;
    int *H;
    int *H_subfield;
    int index, Q, i, j, t, m, n, r, beta, beta_q;
    int f_v = (verbose_level >= 1);
    int f_vv = (verbose_level >= 2);
 
 
    if (f_v) {
        cout << "make_tensor_code_9_dimensional q=" << q << endl;
        }
 
    //q = 2; override_poly_Q = ""; override_poly = ""; int f_hyperoval = FALSE;
    //q = 3; override_poly_Q = ""; override_poly = ""; int f_hyperoval = FALSE;
    //q = 4; override_poly_Q = "19"; override_poly = "7"; int f_hyperoval = FALSE;
        // F_256  generated by X^8 + X^4 + X^3 + X^2 + 1
        // F_16 generated by X^4+X+1 = 19
        // F_4 generated by X^2+X+1 = 7
    //q = 5; override_poly_Q = "47"; override_poly = ""; int f_hyperoval = FALSE;
        // F_625  generated by X^4 + X^3 + 3X + 2
        // F_25  generated by X^2 + 4X + 2  = 47
    //q = 7; override_poly_Q = ""; override_poly = ""; int f_hyperoval = FALSE;
    //q = 8; override_poly_Q = "97"; override_poly = "11"; int f_hyperoval = TRUE;
        // F_4096 generated by x^12+x^6+x^4+x+1
        // F_64 generated by X^6+X^5+1 = 97
        // F_8 generated by X^3+X+1 = 11
    //q = 9; override_poly_Q = ""; override_poly = "17"; int f_hyperoval = FALSE;
    beta = q;
    Q = q * q;
    m = 9;
    if (f_hyperoval) {
        n = Q + 2;
        }
    else {
        n = Q + 1;
        }
    r = 5;
    if (q == 4)
        r = 7;
    if (q == 3)
        r = 9;
    exponents[0] = 0;
    exponents[1] = q + 1;
    exponents[2] = 2 * q + 2;
    exponents[3] = q;
    exponents[4] = 1;
    exponents[5] = 2 * q;
    exponents[6] = 2;
    exponents[7] = 2 * q + 1;
    exponents[8] = q + 2;
    //exponents[0] = 0;
    //exponents[1] = q;
    //exponents[2] = 2 * q;
    //exponents[3] = 1;
    //exponents[4] = q + 1;
    //exponents[5] = 2 * q + 1;
    //exponents[6] = 2;
    //exponents[7] = q + 2;
    //exponents[8] = 2 * q + 2;
 
    index = (Q - 1) / (q - 1);
 
 
    cout << "q = " << q << " override polynomial = " << override_poly << endl;
    cout << "Q = " << Q << endl;
 
    F.init_override_polynomial(Q, override_poly_Q, FALSE /* f_without_tables */, verbose_level);
    cout << "field of order " << Q << " initialized" << endl;
    beta_q = F.power(beta, q);
    f.init_override_polynomial(q, override_poly, FALSE /* f_without_tables */, verbose_level);
    cout << "field of order " << q << " initialized" << endl;
    cout << "n = " << n << endl;
    cout << "index = " << index << endl;
    cout << "beta = " << beta << endl;
    cout << "beta_q = " << beta_q << endl;
    F.compute_subfields(verbose_level - 3);
 
    M = NEW_int(m * n);
    C = NEW_int(m * m);
    C_inv = NEW_int(m * m);
    H = NEW_int(m * n);
    H_subfield = NEW_int(m * n);
 
    rc.init(&f, m, n, r + 1);
 
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) {
            M[i * n + j] = 0;
            }
        }
    for (i = 0; i < m; i++) {
        for (j = 0; j < m; j++) {
            C[i * m + j] = 0;
            }
        }
    for (t = 0; t < Q; t++) {
        for (i = 0; i < m; i++) {
            M[i * n + t] = F.power(t, exponents[i]);
            }
        }
    {
    M[2 * n + Q] = 1;
    if (f_hyperoval)
        M[1 * n + Q + 1] = 1;
    int nb_C_coeffs = 15;
    int k, aa;
    int C_coeffs[] = {
        0, 0, 1,
        1, 1, 1,
        2, 2, 1,
        3, 3, 1,
        3, 4, 1,
        4, 3, beta_q,
        4, 4, beta,
        5, 5, 1,
        5, 6, 1,
        6, 5, beta_q,
        6, 6, beta,
        7, 7, 1,
        7, 8, 1,
        8, 7, beta_q,
        8, 8, beta,
        };
    for (k = 0; k < nb_C_coeffs; k++) {
        i = C_coeffs[k * 3 + 0];
        j = C_coeffs[k * 3 + 1];
        aa = C_coeffs[k * 3 + 2];
        C[i * m + j] = aa;
        }
    }
 
    cout << "M:" << endl;
    Int_vec_print_integer_matrix_width(cout, M, m, n, n, 2);
 
    {
        int *all_one, *col_sum;
 
        all_one = NEW_int(n);
        col_sum = NEW_int(m);
        for (i = 0; i < n; i++) {
            all_one[i] = 1;
        }
        F.Linear_algebra->mult_matrix_matrix(M, all_one, col_sum, m, n, 1,
                0 /* verbose_level */);
        cout << "col_sum:" << endl;
        Int_vec_print_integer_matrix_width(cout, col_sum, m, 1, 1, 2);
        FREE_int(all_one);
        FREE_int(col_sum);
    }
 
#if 0
    for (j = 0; j < n; j++) {
        PG_element_normalize(F, M + j, n, m);
        }
    cout << "column normalized M:" << endl;
    print_integer_matrix_width(cout, M, m, n, n, 2);
#endif
 
 
    cout << "C:" << endl;
    Int_vec_print_integer_matrix_width(cout, C, m, m, m, 2);
 
    F.Linear_algebra->invert_matrix(C, C_inv, m, 0 /* verbose_level */);
 
    cout << "C_inv:" << endl;
    Int_vec_print_integer_matrix_width(cout, C_inv, m, m, m, 2);
 
    {
        int *AA;
        AA = NEW_int(m * m);
        F.Linear_algebra->mult_matrix_matrix(C, C_inv, AA, m, m, m,
                0 /* verbose_level */);
        cout << "C * C_inv:" << endl;
        Int_vec_print_integer_matrix_width(cout, AA, m, m, m, 2);
        FREE_int(AA);
    }
 
    F.Linear_algebra->mult_matrix_matrix(C, M, H, m, m, n,
            0 /* verbose_level */);
    cout << "H = C * M:" << endl;
    Int_vec_print_integer_matrix_width(cout, H, m, n, n, 2);
 
 
#if 0
    rk = F.Gauss_int(M, FALSE /* f_special */, TRUE /* f_complete */, base_cols,
        FALSE /* f_P */, NULL, m /* m */, n /* n */, 0 /* Pn */,
        FALSE, FALSE);
    cout << "has rank " << rk << endl;
#endif
 
    if (f_vv) {
        cout << "before field reduction:" << endl;
        Int_vec_print_integer_matrix_width(cout, H, m, n, n, 2);
        cout << endl;
        f.print_integer_matrix_zech(cout, H, m, n);
        cout << endl;
        }
    F.retract_int_vec(f, 2, H, H_subfield, m * n, 0 /* verbose_level */);
    //field_reduction(F, f, m, n, H, H_subfield, TRUE, TRUE);
    if (f_vv) {
        cout << "after field reduction:" << endl;
        Int_vec_print_integer_matrix_width(cout, H_subfield, m, n, n, 2);
        cout << endl;
        f.print_integer_matrix_zech(cout, H_subfield, m, n);
        cout << endl;
        }
    cout << "H_subfield:" << endl;
    Int_vec_print_integer_matrix_width(cout, H_subfield, m, n, n, 2);
 
    code = H_subfield;
    length = n;
 
    FREE_int(M);
    FREE_int(C);
    FREE_int(C_inv);
    FREE_int(H);
 
}
 
 
 
}}}