spread_classify2.cpp

Go to the documentation of this file.

// spread2.cpp
// 
// Anton Betten
// November 17, 2009
//
// 
//
// moved here from translation_plane.cpp: April 23, 2013
// moved to TOP_LEVEL: November 2, 2013
// renamed to spread2.cpp from translation_plane2.cpp: March 25, 2018
//
//
 
#include "orbiter.h"
 
using namespace std;
 
namespace orbiter {
namespace layer5_applications {
namespace spreads {
 
 
 
// spread_classify2.cpp
#if 0
static int spread_check_function_callback(int len, long int *S,
    void *data, int verbose_level);
static void spread_callback_report(isomorph *Iso, void *data, int verbose_level);
static void spread_callback_make_quotients(isomorph *Iso, void *data,
    int verbose_level);
static void callback_spread_print(std::ostream &ost, int len, long int *S, void *data);
#endif
 
void spread_classify::print_isomorphism_type(isomorph *Iso,
    int iso_cnt, groups::sims *Stab, groups::schreier &Orb,
    long int *data, int verbose_level)
// called from callback_print_isomorphism_type()
{
    int f_v = (verbose_level >= 1);
    string fname;
    string fname_klein;
    char str[1000];
 
 
 
 
    if (f_v) {
        cout << "spread_classify::print_isomorphism_type" << endl;
    }
 
    fname.assign(Iso->prefix);
    sprintf(str, "_%d.tex", iso_cnt);
    fname.append(str);
 
    //sprintf(fname, "%s_%d.tex", Iso->prefix, iso_cnt);
 
 
    fname_klein.assign(Iso->prefix);
    sprintf(str, "_%d_klein.tex", iso_cnt);
    fname_klein.append(str);
 
    //sprintf(fname_klein, "%s_%d_klein.tex", Iso->prefix, iso_cnt);
    
    int save_longinteger_f_print_scientific = orbiter_kernel_system::Orbiter->longinteger_f_print_scientific;
    orbiter_kernel_system::Orbiter->longinteger_f_print_scientific = FALSE;
    {
        ofstream file(fname);
        orbiter_kernel_system::latex_interface L;
        //ofstream file_klein(fname_klein);
        if (f_v) {
            cout << "opening file " << fname << " for writing" << endl;
            //cout << "opening file " << fname_klein << " for writing" << endl;
        }
    
        L.head(file,
            FALSE/* f_book*/, FALSE /* f_title */, NULL /*title*/,
            "Orbiter" /*char *author*/, FALSE/* f_toc*/,
            FALSE /* f_landscape*/,
            FALSE /* f_12pt */,
            TRUE /* f_enlarged_page */,
            TRUE /* f_pagenumbers */,
            NULL /* extra_praeamble */);
 
 
        print_isomorphism_type2(Iso, file, iso_cnt, Stab, Orb, data, verbose_level);
    
        L.foot(file);
 
        orbiter_kernel_system::Orbiter->longinteger_f_print_scientific = save_longinteger_f_print_scientific;
    }
 
    orbiter_kernel_system::file_io Fio;
 
    if (f_v) {
        cout << "spread_classify::print_isomorphism_type written file "
                << fname << " of size " << Fio.file_size(fname) << endl;
        cout << "spread_classify::print_isomorphism_type written file "
                << fname_klein << " of size "
                << Fio.file_size(fname_klein) << endl;
    }
 
 
 
}
 
 
void spread_classify::print_isomorphism_type2(isomorph *Iso,
        ostream &ost,
        int iso_cnt, groups::sims *Stab, groups::schreier &Orb,
        long int *data, int verbose_level)
{
    int f_v = (verbose_level >= 1);
    int i;
    ring_theory::longinteger_object so;
    long int *pt_list;
    int f, l, j, idx, pt;
    int p, h;
    int f_elements_exponential = FALSE;
    string symbol_for_print;
 
    number_theory::number_theory_domain NT;
    data_structures::sorting Sorting;
 
    if (f_v) {
        cout << "spread_classify::print_isomorphism_type2" << endl;
    }
 
 
    symbol_for_print.assign("\\alpha");
 
 
    NT.factor_prime_power(q, p, h);
    if (h > 1) {
        f_elements_exponential = TRUE;
    }
 
    ost << "\\parindent=0pt" << endl;
    ost << endl;
    ost << "\\section{Spread " << iso_cnt << "}" << endl;
    //int_vec_print(file, data, q + 1);
    //file << endl;
    
    Stab->group_order(so);
 
    data_structures::tally C;
 
 
    C.init(Orb.orbit_len, Orb.nb_orbits, FALSE, 0);
 
    ost << "Stabilizer has order " << so << "\\\\" << endl;
    
    ost << "\\bigskip" << endl;
 
    ost << "There are $" << Orb.nb_orbits
            << "$ orbits on the set.\\\\" << endl;
    ost << "The orbit type is $[";
    C.print_naked_tex(ost, FALSE /*f_backwards*/);
    ost << "]$\\\\" << endl;
    ost << "\\bigskip" << endl;
    
 
    if (k == 2) {
        klein(ost, /*file_klein,*/ Iso, iso_cnt,
                Stab, Orb, data, Iso->size,  verbose_level);
    }
 
    ost << "The set of " << k - 1 << "-Subspaces in "
            "PG($" << n - 1 << ", " << q << "$) is:\\\\" << endl;
    
    for (i = 0; i < Iso->size; i++) {
        ost << "$S_{" << i + 1 << "}$ has rank " << data[i]
            << " and is generated by\\\\" << endl;
        Grass->unrank_lint(data[i], 0);
        ost << "$$" << endl;
        ost << "\\left[" << endl;
        Mtx->GFq->latex_matrix(ost, f_elements_exponential, symbol_for_print,
            Grass->M, k, n);
        ost << "\\right]" << endl;
        ost << "$$" << endl << endl;
    }
 
    ost << "The stabilizer of order " << so
            << " is generated by:\\\\" << endl;
    for (i = 0; i < Stab->gens.len; i++) {
        
        int *fp, n;
        
        fp = NEW_int(A->degree);
        n = A->find_fixed_points(Stab->gens.ith(i), fp, 0);
        //cout << "with " << n << " fixed points" << endl;
        FREE_int(fp);
 
        ost << "$$ g_{" << i + 1 << "}=" << endl;
        A->element_print_latex(Stab->gens.ith(i), ost);
        ost << "$$" << endl << "with " << n << " fixed points" << endl;
    }
 
 
 
    if (f_v) {
        cout << "spread_classify::print_isomorphism_type "
                "calling induced_action_on_set_and_kernel" << endl;
    }
    Iso->induced_action_on_set_and_kernel(ost,
            A, Stab, Iso->size, data, verbose_level - 1);
 
    if (f_v) {
        cout << "spread_classify::print_isomorphism_type "
                "induced_action_on_set_and_kernel finished" << endl;
    }
 
 
#if 0
    for (i = 0; i < Stab->gens.len; i++) {
        file << "$g_{" << setw(2) << i + 1 << "} = ";
        AA->element_print_as_permutation_with_offset(Stab->gens.ith(i), file, 1);
        file << "$\\\\" << endl;
    }
#endif
    
    pt_list = NEW_lint(Iso->size);
    
    ost << "The orbits on the set are:\\\\" << endl;
    for (i = 0; i < Orb.nb_orbits; i++) {
        f = Orb.orbit_first[i];
        l = Orb.orbit_len[i];
        ost << "$O_{" << i << "}=\\{";
        for (j = 0; j < l; j++) {
            idx = f + j;
            pt = Orb.orbit[idx];
            pt_list[j] = pt;
        }
        
        Sorting.lint_vec_heapsort(pt_list, l);
        
        for (j = 0; j < l; j++) {       
            ost << pt_list[j] + 1;
            if (j < l - 1)
                ost << ",";
        }
        ost << "\\}$ (length $" << l << "$)\\\\" << endl;
    }
 
#if 0
    
    if (Orb.nb_orbits > 1) {
        file << "\\bigskip" << endl;
        file << "The actions induced on the orbits are:\\\\" << endl;
    
        for (i = 0; i < Orb.nb_orbits; i++) {
            induced_action_on_orbit(file, AA, Stab, Orb, i, verbose_level - 2);
        }
    }
#endif
 
 
 
    FREE_lint(pt_list);
 
 
    if (f_v) {
        cout << "spread_classify::print_isomorphism_type2 done" << endl;
    }
 
}
 
void spread_classify::save_klein_invariants(char *prefix,
    int iso_cnt, 
    long int *data, int data_size, int verbose_level)
{
    int f_v = (verbose_level >= 1);
    ring_theory::longinteger_object *R;
    long int **Pts_on_plane;
    int *nb_pts_on_plane;
    int nb_planes;
    int i, j;
 
    if (f_v) {
        cout << "spread_classify::klein_invariants" << endl;
    }
    
    if (Klein == NULL) {
        cout << "spread_classify::klein_invariants Klein == NULL" << endl;
        exit(1);
    }
    
    Klein->plane_intersections(data, data_size, 
        R,
        Pts_on_plane, 
        nb_pts_on_plane, 
        nb_planes, 
        verbose_level - 2);
 
 
    Vector v;
 
    v.m_l(3);
    v.m_ii(0, nb_planes);
    v.s_i(1).change_to_vector();
    v.s_i(2).change_to_vector();
 
    v.s_i(1).as_vector().m_l(nb_planes);
    v.s_i(2).as_vector().m_l(nb_planes);
    for (i = 0; i < nb_planes; i++) {
        v.s_i(1).as_vector().m_ii(i, R[i].as_int());
        //v.s_i(1).as_vector().s_i(i).change_to_longinteger();
        //v.s_i(1).as_vector().s_i(i).as_longinteger().allocate(1, R[i].rep());
        v.s_i(2).as_vector().s_i(i).change_to_vector();
        v.s_i(2).as_vector().s_i(i).as_vector().m_l(nb_pts_on_plane[i]);
        for (j = 0; j < nb_pts_on_plane[i]; j++) {
            v.s_i(2).as_vector().s_i(i).as_vector().m_ii(j, Pts_on_plane[i][j]);
        }
    }
 
    char fname[1000];
    
    sprintf(fname, "%s%d_klein_invariant.bin", prefix, iso_cnt);
    v.save_file(fname);
 
    delete [] R;
    for (i = 0; i < nb_planes; i++) {
        FREE_lint(Pts_on_plane[i]);
    }
    FREE_plint(Pts_on_plane);
    FREE_int(nb_pts_on_plane);
 
    if (f_v) {
        cout << "spread_classify::klein_invariants done" << endl;
    }
}
 
void spread_classify::klein(ostream &ost,
    isomorph *Iso, 
    int iso_cnt, groups::sims *Stab, groups::schreier &Orb,
    long int *data, int data_size, int verbose_level)
// Called from print_isomorphism_type if k == 2
{
    int f_v = (verbose_level >= 1);
    ring_theory::longinteger_object *R;
    long int **Pts_on_plane;
    int *nb_pts_on_plane;
    int nb_planes;
 
 
 
    int set_size = data_size;
    int a, i, j, h;
    
    Klein->plane_intersections(data, data_size, 
        R,
        Pts_on_plane, 
        nb_pts_on_plane, 
        nb_planes, 
        verbose_level);
 
    data_structures::tally C;
    int f_second = FALSE;
 
    C.init(nb_pts_on_plane, nb_planes, f_second, 0);
    if (f_v) {
        cout << "spread::klein: plane-intersection type: ";
        C.print(FALSE /*f_backwards*/);
    }
    
    ost << "Plane type of Klein-image is $( ";
    C.print_naked_tex(ost, FALSE /*f_backwards*/);
    ost << " )$" << endl << endl;
    ost << "\\bigskip" << endl << endl;
 
    int nb_blocks, f, l, m, u, uu, idx;
    int *Inc;
    
    m = 0;
    for (i = 0; i < C.nb_types; i++) {
        f = C.type_first[i];
        l = C.type_len[i];
        a = C.data_sorted[f];
        m = MAXIMUM(a, m);
    }
 
    nb_blocks = 0;
    for (i = 0; i < C.nb_types; i++) {
        f = C.type_first[i];
        l = C.type_len[i];
        a = C.data_sorted[f];
        if (a == m) {
            nb_blocks += l;
        }
    }
    if (f_v) {
        cout << "There are " << nb_blocks
                << " interesting planes" << endl;
    }
    Inc = NEW_int(set_size * nb_blocks);
    for (i = 0; i < set_size * nb_blocks; i++) {
        Inc[i] = 0;
    }
    j = 0;
    for (h = 0; h < C.nb_types; h++) {
        f = C.type_first[h];
        l = C.type_len[h];
        a = C.data_sorted[f];
        if (a == m) {
            for (u = 0; u < l; u++) {
                a = C.data_sorted[f + u];
                idx = C.sorting_perm_inv[f + u];
                for (uu = 0; uu < a; uu++) {
                    i = Pts_on_plane[idx][uu];
                    Inc[i * nb_blocks + j] = 1;
                }
                j++;
            } // next u
        } // if
    } // next h
 
    if (nb_blocks < 20) {
 
        cout << "we will draw an incidence picture" << endl;
        
        string fname_pic;
        char str[1000];
        geometry::incidence_structure *I;
        data_structures::partitionstack *Stack;
        
        I = NEW_OBJECT(geometry::incidence_structure);
        I->init_by_matrix(set_size, nb_blocks, Inc, 0 /* verbose_level */);
        Stack = NEW_OBJECT(data_structures::partitionstack);
        Stack->allocate(set_size + nb_blocks, 0 /* verbose_level */);
        Stack->subset_continguous(set_size, nb_blocks);
        Stack->split_cell(0 /* verbose_level */);
        Stack->sort_cells();
 
        fname_pic.assign(Iso->prefix);
        sprintf(str, "_%d_planes.tex", iso_cnt);
        fname_pic.append(str);
 
        {
            ofstream fp_pic(fname_pic);
 
            ost << "\\input " << fname_pic << endl;
            I->latex_it(fp_pic, *Stack);
            ost << "\\\\" << endl;
        }
        FREE_OBJECT(Stack);
        FREE_OBJECT(I);
    }
 
    // compute TDO:
    {
        geometry::incidence_structure *I;
        data_structures::partitionstack *Stack;
        int depth = INT_MAX;
        
        I = NEW_OBJECT(geometry::incidence_structure);
        I->init_by_matrix(set_size, nb_blocks, Inc, 0 /* verbose_level */);
        Stack = NEW_OBJECT(data_structures::partitionstack);
        Stack->allocate(set_size + nb_blocks, 0 /* verbose_level */);
        Stack->subset_continguous(set_size, nb_blocks);
        Stack->split_cell(0 /* verbose_level */);
        Stack->sort_cells();
 
        I->compute_TDO_safe(*Stack, depth, verbose_level + 2);
        
        char str[1000];
        string fname_row_scheme;
        string fname_col_scheme;
 
 
        fname_row_scheme.assign(Iso->prefix);
        sprintf(str, "_%d_planes_row_scheme.tex", iso_cnt);
        fname_row_scheme.append(str);
 
        fname_col_scheme.assign(Iso->prefix);
        sprintf(str, "_%d_planes_col_scheme.tex", iso_cnt);
        fname_col_scheme.append(str);
 
        {
            ofstream fp_row_scheme(fname_row_scheme);
            ofstream fp_col_scheme(fname_col_scheme);
            I->get_and_print_row_tactical_decomposition_scheme_tex(
                fp_row_scheme, FALSE /* f_enter_math */,
                TRUE /* f_print_subscripts */, *Stack);
            I->get_and_print_column_tactical_decomposition_scheme_tex(
                fp_col_scheme, FALSE /* f_enter_math */,
                TRUE /* f_print_subscripts */, *Stack);
        }
 
 
#if 0
        ost2 << "\\section*{Spread " << iso_cnt << "}" << endl << endl;
 
        longinteger_object go;
        Stab->group_order(go);
 
        ost2 << "Stabilizer has order " << go << "\\\\" << endl;
    
        ost2 << "\\bigskip" << endl;
 
 
        tally C0;
        C0.init(Orb.orbit_len, Orb.nb_orbits, FALSE, 0);
        
        ost2 << "There are $" << Orb.nb_orbits << "$ orbits on the set.\\\\" << endl;
        ost2 << "The orbit type is $[";
        C0.print_naked_tex(ost2);
        ost2 << "]$\\\\" << endl;
        ost2 << "\\bigskip" << endl;
 
        ost2 << "Plane type of Klein-image is $(";
        C.print_naked_tex(ost2);
        ost2 << ")$" << endl << endl;
#endif
 
 
        ost << "\\bigskip" << endl << endl;
        ost << "\\[" << endl;
        ost << "\\input " << fname_row_scheme << endl;
        ost << "\\quad" << endl;
        ost << "\\input " << fname_col_scheme << endl;
        ost << "\\]" << endl;
 
        FREE_OBJECT(Stack);
        FREE_OBJECT(I);
    }
 
    FREE_int(Inc);
    
    FREE_OBJECTS(R);
    for (i = 0; i < nb_planes; i++) {
        FREE_lint(Pts_on_plane[i]);
    }
    FREE_plint(Pts_on_plane);
    FREE_int(nb_pts_on_plane);
 
    if (f_v) {
        cout << "spread_classify::klein done" << endl;
    }
}
 
void spread_classify::plane_intersection_type_of_klein_image(
        geometry::projective_space *P3,
        geometry::projective_space *P5,
        geometry::grassmann *Gr,
    long int *data, int size,
    int *&intersection_type, int &highest_intersection_number, 
    int verbose_level)
{
    int f_v = (verbose_level >= 1);
    int f_v3 = FALSE; //(verbose_level >= 3);
    //finite_field *F;
    long int N;
    long int *the_set_out;
 
    if (f_v) {
        cout << "spread_classify::plane_intersection_type_of_klein_image" << endl;
    }
    the_set_out = NEW_lint(size);
    
    if (f_v3) {
        P3->print_line_set_numerical(data, size);
    }
    P3->klein_correspondence(P5, 
        data, size, the_set_out, 0/*verbose_level*/);
    if (f_v) {
        cout << "after Klein correspondence:" << endl;
        Lint_vec_print(cout, the_set_out, size);
        cout << endl;
    }
    if (f_v3) {
        P5->print_set_numerical(cout, the_set_out, size);
    }
 
 
    //F = P3->F;
 
    if (f_v) {
        cout << "spread_classify::plane_intersection_type_of_klein_image "
                "after P3->klein_correspondence" << endl;
    }
 
    
    N = P5->nb_rk_k_subspaces_as_lint(3);
    if (f_v) {
        cout << "spread_classify::klein N = " << N << endl;
        }
 
 
    P5->plane_intersection_type(Gr, 
        the_set_out, size, 
        intersection_type, highest_intersection_number, 
        verbose_level - 2);
 
    FREE_lint(the_set_out);
}
 
void spread_classify::czerwinski_oakden(int level, int verbose_level)
{
    int f_v = (verbose_level >= 1);
    int f_vv = (verbose_level >= 2);
    int sz = 26;
    long int data[26];
    int M[8];
    int h, u, i, a = 0, b = 0, c = 0, d = 0;
    orbiter_kernel_system::file_io Fio;
    int spreads[] =     
        {
            // S1:
            0,0,0,0,
            1,0,0,1,
            2,0,0,2,
            3,0,0,3,
            4,0,0,4,
            4,3,3,0,
            0,3,3,1,
            2,1,1,4,
            4,4,4,2,
            2,3,3,3,
            3,4,4,1,
            1,1,1,3,
            3,2,2,2,
            0,1,1,2,
            1,2,2,0,
            0,2,2,4,
            2,4,4,0,
            0,4,4,3,
            3,1,1,0,
            1,4,4,4,
            4,2,2,3,
            3,3,3,4,
            4,1,1,1,
            1,3,3,2,
            2,2,2,1,
 
            //S2:
            0,0,0,0,
            1,0,0,1,
            2,0,0,2,
            3,0,0,3,
            4,0,0,4,
            1,2,2,0,
            0,2,2,4,
            4,3,3,0,
            0,3,3,1,
            2,1,3,3,
            4,4,4,1,
            2,3,1,3,
            3,4,2,2,
            1,1,1,4,
            3,2,4,2,
            0,1,1,2,
            2,4,4,0,
            0,4,4,3,
            3,1,1,0,
            1,4,4,4,
            4,2,2,3,
            3,3,3,4,
            4,1,1,1,
            1,3,3,2,
            2,2,2,1,
 
            // S3:
            0,0,0,0,
            1,0,0,1,
            1,2,2,0,
            0,2,2,4,
            4,0,0,4,
            4,3,3,0,
            0,3,3,1,
            2,1,3,3,
            4,4,4,1,
            2,3,1,3,
            3,4,2,2,
            1,1,1,4,
            3,2,4,2,
            0,1,4,0,
            2,0,4,3,
            2,4,0,3,
            0,4,1,0,
            3,0,1,2,
            3,1,0,2,
            1,4,4,4,
            4,2,2,3,
            3,3,3,4,
            4,1,1,1,
            1,3,3,2,
            2,2,2,1,
 
            //S4:
            0,0,0,0,
            1,0,0,1,
            1,2,2,0,
            0,2,2,4,
            4,0,0,4,
            4,3,3,0,
            0,3,3,1,
            2,1,3,3,
            4,4,4,1,
            2,3,1,3,
            3,4,2,2,
            1,1,1,4,
            3,2,4,2, 
            0,1,1,2,
            2,0,0,2,
            2,4,4,0,
            0,4,4,3,
            3,0,0,3,
            3,1,1,0,
            1,4,3,4,
            4,2,1,1,
            3,3,3,2,
            4,1,2,1,
            1,3,4,4,
            2,2,2,3,
 
            // S5:
            0,0,0,0,
            1,0,0,1,
            2,0,0,2,
            3,0,0,3,
            4,0,0,4,
            0,2,1,0,
            0,3,4,0,
            1,1,3,1,
            1,4,2,1,
            4,1,3,4,
            4,4,2,4,
            1,2,1,4,
            1,3,4,4,
            2,1,2,0,
            2,2,4,3,
            2,3,3,3,
            2,4,1,3,
            3,1,4,2,
            3,2,2,2,
            3,3,1,2,
            3,4,3,0,
            4,2,1,1,
            4,3,4,1,
            0,1,2,3,
            0,4,3,2,
 
            // A1:
            0,0,0,0,
            1,0,0,1,
            2,0,0,2,
            3,0,0,3,
            4,0,0,4,
            1,1,3,1,
            4,4,2,4,
            1,4,2,1,
            4,1,3,4,
            3,2,1,3,
            2,2,1,2,
            3,3,4,3,
            2,3,4,2,
            4,2,4,4,
            1,3,1,1,
            3,4,3,0,
            2,1,2,0,
            0,4,3,3,
            0,1,2,2,
            1,2,4,0,
            4,3,1,0,
            2,4,3,2,
            3,1,2,3,
            0,2,4,1,
            0,3,1,4,
 
            //A2:
            0,0,0,0,
            1,0,0,1,
            2,0,0,2,
            3,0,0,3,
            4,0,0,4,
            2,2,1,2,
            0,3,4,0,
            0,1,3,0,
            4,1,3,4,
            2,4,2,2,
            3,4,2,3,
            3,2,3,1,
            4,3,2,1,
            3,1,4,2,
            2,1,4,1,
            4,4,1,0,
            0,4,1,1,
            1,1,1,3,
            1,4,4,4,
            2,3,3,3,
            1,3,3,2,
            0,2,2,4,
            1,2,2,0,
            4,2,4,3,
            3,3,1,4,
 
            //A3:
            0,0,0,0,
            1,0,0,1,
            2,0,0,2,
            3,0,0,3,
            4,0,0,4,
            1,3,4,1,
            2,3,4,2,
            4,2,1,4,
            2,4,4,2,
            1,4,2,1,
            2,2,3,0,
            3,2,3,1,
            3,3,2,0,
            0,1,4,4,
            4,1,4,3,
            4,4,2,3,
            0,4,2,4,
            1,1,3,2,
            1,2,1,3,
            0,2,1,2,
            2,1,1,0,
            3,1,1,1,
            3,4,4,0,
            0,3,3,4,
            4,3,3,3,
            
            // A4:
            0,0,0,0,
            1,0,0,1,
            2,0,0,2,
            3,0,0,3,
            4,0,0,4,
            2,2,1,2,
            0,3,4,0,
            4,1,3,4,
            3,4,2,3,
            3,2,3,1,
            3,1,4,2,
            2,1,4,1,
            4,4,1,0,
            1,4,4,4,
            2,3,3,3,
            0,2,2,4,
            1,2,2,0,
            4,2,4,3,
            3,3,1,4,
            2,4,2,1,
            0,1,1,3,
            4,3,2,2,
            1,1,3,0,
            1,3,1,1,
            0,4,3,2,
            
            // A5:
            0,0,0,0,
            1,0,0,1,
            2,0,0,2,
            3,0,0,3,
            4,0,0,4,
            0,1,3,0,
            1,2,1,4,
            2,3,3,3,
            3,4,1,2,
            1,3,4,4,
            4,3,4,2,
            0,4,2,4,
            3,1,3,4,
            0,2,4,0,
            3,2,4,3,
            2,1,2,2,
            1,4,3,1,
            4,1,2,1,
            1,1,2,3,
            2,2,4,1,
            0,3,1,1,
            2,4,1,3,
            4,4,1,0,
            3,3,2,0,
            4,2,3,2,
 
            //A6:
            0,0,0,0,
            1,0,0,1,
            2,0,0,2,
            3,0,0,3,
            4,0,0,4,
            4,2,1,1,
            1,3,4,4,
            3,4,2,2,
            3,2,1,0,
            0,2,1,2,
            2,2,2,3,
            3,3,3,2,
            1,1,1,4,
            2,3,3,1,
            4,3,3,3,
            0,1,3,0,
            0,4,2,0,
            0,3,4,0,
            4,4,2,4,
            1,4,2,1,
            2,4,1,3,
            3,1,4,2,
            1,2,3,4,
            2,1,4,1,
            4,1,4,3,
 
            //A7:
            0,0,0,0,
            1,0,0,1,
            2,0,0,2,
            3,0,0,3,
            4,0,0,4,
            2,2,2,3,
            1,3,3,0,
            1,4,4,3,
            2,1,1,0,
            4,2,1,1, 
            4,3,4,2,
            0,4,2,4,
            3,1,3,4,
            2,3,3,3,
            1,2,2,0,
            1,1,1,3,
            2,4,4,0,
            3,4,1,2,
            0,1,4,1,
            3,3,2,1,
            0,2,3,2,
            4,4,3,1,
            4,1,2,2,
            0,3,1,4,
            3,2,4,4,
 
            //A8:
            0,0,0,0,
            1,0,0,1,
            2,0,0,2,
            3,0,0,3,
            4,0,0,4,
            1,3,4,4,
            4,3,4,2,
            4,1,2,1,
            2,1,2,4,
            0,2,4,0,
            4,4,3,4,
            2,4,1,3,
            3,3,2,0,
            3,1,3,3,
            4,2,1,4,
            1,4,1,0,
            0,3,2,3,
            1,2,2,2,
            2,3,3,1,
            3,4,4,1,
            0,1,1,2,
            3,2,1,1,
            1,1,3,0,
            0,4,3,2,
            2,2,4,3,
 
            // B1:
            0,0,0,0,
            1,0,0,1,
            3,0,0,3,
            2,2,4,2,
            2,4,3,2,
            0,3,1,0,
            1,2,3,4,
            2,0,1,2,
            2,3,2,1,
            4,3,0,2,
            0,4,3,1,
            0,2,4,1,
            1,3,4,4,
            4,4,2,3,
            3,4,2,4,
            3,1,1,1,
            1,1,1,3,
            0,1,2,2,
            4,0,1,4,
            3,2,3,3,
            2,1,0,4,
            4,1,4,0,
            4,2,3,0,
            3,3,4,3,
            1,4,2,0,
 
            // B2:
            0,0,0,0,
            1,0,0,1,
            4,0,0,4,
            3,4,3,3,
            3,3,1,3,
            4,1,2,4,
            1,4,1,2,
            0,4,4,1,
            1,2,2,3,
            2,1,0,3,
            3,1,0,2,
            3,2,2,2,
            4,2,1,4,
            2,2,1,0,
            1,1,3,1,
            2,4,3,4,
            0,2,2,1,
            2,3,2,0,
            4,4,4,0,
            2,0,4,3,
            3,0,4,2,
            0,3,1,1,
            4,3,3,2,
            0,1,3,0,
            1,3,4,4,
            
            // B3:
            0,0,0,0,
            1,0,0,1,
            4,0,0,4,
            4,4,3,4,
            2,2,4,2,
            4,1,2,4,
            0,4,3,0,
            0,2,4,4,
            4,2,4,0,
            1,2,2,3,
            3,1,2,0,
            2,1,3,3,
            1,3,1,4,
            1,1,3,2,
            2,4,2,1,
            0,1,1,1,
            1,4,4,3,
            4,3,1,0,
            3,2,3,1,
            2,3,0,3,
            2,0,1,3,
            3,0,1,2,
            0,3,2,2,
            3,4,0,2,
            3,3,4,1,
 
            // B4:
            0,0,0,0,
            1,0,0,1,
            4,0,0,4,
            3,4,3,3,
            3,1,2,3,
            1,3,1,1,
            1,2,4,1,
            0,1,1,2,
            2,3,2,0,
            2,0,3,2,
            2,4,0,2,
            1,1,3,0,
            0,4,2,1,
            2,2,1,0,
            0,3,4,2,
            4,3,1,3,
            3,3,1,4,
            4,4,4,0,
            0,2,3,4,
            2,1,4,4,
            4,2,2,2,
            3,0,4,3,
            3,2,0,3,
            1,4,2,4,
            4,1,3,1,
            
            // B5:
            0,0,0,0,
            1,0,0,1,
            2,0,0,2,
            4,3,1,4,
            1,2,4,1,
            3,3,1,3,
            1,1,2,1,
            1,4,1,2,
            2,3,3,1,
            3,0,4,4,
            4,2,0,3,
            0,4,1,1,
            1,3,3,0,
            0,3,2,2,
            2,1,1,0,
            3,4,0,4,
            4,0,3,3,
            0,2,4,2,
            2,2,4,0,
            4,4,2,0,
            0,1,3,4,
            4,1,4,3,
            3,2,2,4,
            2,4,2,3,
            3,1,3,2,
 
            // B6:
            0,0,0,0,
            1,0,0,1,
            4,2,2,3,
            1,2,2,0,
            3,2,2,2,
            2,2,2,4,
            3,3,3,1,
            1,4,4,2,
            2,4,1,1,
            1,1,1,0,
            3,4,0,4,
            0,1,4,4,
            2,0,4,3,
            1,3,3,0,
            3,1,4,1,
            4,0,3,4,
            4,1,0,3,
            0,4,1,3,
            0,2,2,1,
            2,1,3,2,
            4,4,4,0,
            2,3,0,2,
            0,3,3,3,
            3,0,1,2,
            4,3,1,4,
 
            //B7:
            0,0,0,0,
            1,0,0,1,
            2,2,4,2,
            3,2,4,3,
            1,2,4,1,
            3,3,1,3,
            2,4,3,2,
            1,1,1,0,
            0,3,2,1,
            3,1,0,2,
            2,0,2,3,
            4,1,0,3,
            3,0,2,4,
            4,4,4,0,
            0,2,3,4,
            4,3,3,0,
            0,4,1,4,
            3,4,3,1,
            0,1,4,4,
            1,4,3,3,
            4,2,2,0,
            2,3,0,4,
            4,0,1,2,
            1,3,2,2,
            2,1,1,1,
 
            // B8:
            0,0,0,0,
            1,0,0,1,
            2,0,0,2,
            3,0,0,3,
            4,2,2,0,
            0,2,2,3,
            4,3,3,0,
            0,4,4,3,
            4,4,4,2,
            4,0,1,4,
            4,1,0,4,
            3,2,4,4,
            3,4,2,4,
            0,1,1,2,
            2,3,1,1,
            2,1,3,1,
            1,4,1,0,
            1,1,4,0,
            3,1,3,3,
            2,2,2,1,
            3,3,1,3,
            2,4,4,1,
            1,2,3,2,
            1,3,2,2,
            0,3,3,4,
 
        };
 
    if (f_v) {
        cout << "spread_classify::czerwinski_oakden" << endl;
        }
 
    const char *label[] = {
        "S1", 
        "S2", 
        "S3", 
        "S4", 
        "S5", 
        "A1", 
        "A2", 
        "A3", 
        "A4", 
        "A5", 
        "A6", 
        "A7", 
        "A8", 
        "B1", 
        "B2", 
        "B3", 
        "B4", 
        "B5", 
        "B6", 
        "B7", 
        "B8", 
        };
    char fname[] = "Czerwinski_Oakden.txt";
    string fname2;
    {
    ofstream fp(fname);
    
    for (h = 0; h < 21; h++) {
        for (u = 0; u < sz; u++) {
            for (i = 0; i < 8; i++) {
                M[i] = 0;
                }
            if (u == 0) {
                M[0 * 4 + 2] = 1;
                M[1 * 4 + 3] = 1;
                }
            else {
                M[0 * 4 + 0] = 1;
                M[1 * 4 + 1] = 1;
                a = spreads[h * 25 * 4 + (u - 1) * 4 + 0];
                b = spreads[h * 25 * 4 + (u - 1) * 4 + 1];
                c = spreads[h * 25 * 4 + (u - 1) * 4 + 2];
                d = spreads[h * 25 * 4 + (u - 1) * 4 + 3];
                M[0 * 4 + 2] = a;
                M[0 * 4 + 3] = b;
                M[1 * 4 + 2] = c;
                M[1 * 4 + 3] = d;
                }
            for (i = 0; i < 8; i++) {
                Grass->M[i] = M[i];
                }
            if (f_vv) {
                cout << "spread " << h << ", element " << u << ":" << endl;
                if (u) {
                    cout << "a=" << a << " b=" << b
                            << " c=" << c << " d=" << d << endl;
                    }
                }
            Int_matrix_print(Grass->M, 2, 4);
            data[u] = Grass->rank_lint(0);
 
            } // next u
 
        cout << "spread " << h << ":";
        Lint_vec_print(cout, data, sz);
        cout << endl;
 
        fp << "0 "; // a dummy
        for (i = 0; i < sz; i++) {
            fp << data[i] << " ";
            }
        fp << endl;
 
 
        fname2.assign("Czerwinski_Oakden_");
        fname2.append(label[h]);
        fname2.append(".txt");
 
        Fio.write_set_to_file(fname2, data, sz, 0/*verbose_level*/);
        cout << "Written file " << fname2 << " of size "
                << Fio.file_size(fname2) << endl;
        } // next h
    fp << -1 << endl;
    }
    cout << "Written file " << fname << " of size " << Fio.file_size(fname) << endl;
}
 
void spread_classify::write_spread_to_file(int type_of_spread, int verbose_level)
{
    int f_v = (verbose_level >= 1);
    long int *data;
    int sz = order + 1;
    char str[1000];
 
    if (f_v) {
        cout << "spread_classify::write_spread_to_file" << endl;
        }
    if (type_of_spread == SPREAD_OF_TYPE_FTWKB) {
        sprintf(str, "spread_q%d_FTW.txt", q);
        }
    else if (type_of_spread == SPREAD_OF_TYPE_KANTOR) {
        sprintf(str, "spread_q%d_Kantor.txt", q);
        }
    else if (type_of_spread == SPREAD_OF_TYPE_KANTOR2) {
        sprintf(str, "spread_q%d_Kantor2.txt", q);
        }
    else if (type_of_spread == SPREAD_OF_TYPE_GANLEY) {
        sprintf(str, "spread_q%d_Ganley.txt", q);
        }
    else if (type_of_spread == SPREAD_OF_TYPE_LAW_PENTTILA) {
        sprintf(str, "spread_q%d_Law_Penttila.txt", q);
        }
    else if (type_of_spread == SPREAD_OF_TYPE_DICKSON_KANTOR) {
        sprintf(str, "spread_q%d_DicksonKantor.txt", q);
        }
    else if (type_of_spread == SPREAD_OF_TYPE_HUDSON) {
        sprintf(str, "spread_q%d_Hudson.txt", q);
        }
 
    string fname;
 
    fname.assign(str);
 
    orbiter_kernel_system::file_io Fio;
 
    data = NEW_lint(sz);
    if (type_of_spread == SPREAD_OF_TYPE_DICKSON_KANTOR ||
        type_of_spread == SPREAD_OF_TYPE_HUDSON) {
        make_spread(data, type_of_spread, verbose_level);
    }
    else {
        make_spread_from_q_clan(data, type_of_spread, verbose_level);
    }
    Fio.write_set_to_file(fname, data, sz, 0/*verbose_level*/);
    cout << "Written file " << fname << " of size "
            << Fio.file_size(fname) << endl;
    FREE_lint(data);
}
 
void spread_classify::make_spread(long int *data,
        int type_of_spread, int verbose_level)
{
    int f_v = (verbose_level >= 1);
    int f_vv = (verbose_level >= 2);
    int sz = order + 1;
    int M[8];
    int h, i, h1, s, t, sq, tq, x, y, w, z, eta, exponent;
    number_theory::number_theory_domain NT;
    int q1 = NT.i_power_j(Mtx->GFq->p, (Mtx->GFq->e >> 1));
 
    if (f_v) {
        cout << "spread_classify::make_spread q=" << q << " q1=" << q1 << endl;
    }
    if (n != 4) {
        cout << "spread_classify::make_spread n != 4" << endl;
        exit(1);
    }
    if (EVEN(q)) {
        cout << "spread_classify::make_spread need q odd" << endl;
        exit(1);
    }
    if (k != 2) {
        cout << "spread_classify::make_spread k != 2" << endl;
        exit(1);
    }
    for (eta = q1; eta < q; eta++) {
        if (Mtx->GFq->negate(Mtx->GFq->power(eta, q1)) == eta) {
            if (f_v) {
                cout << "spread_classify::make_spread eta=" << eta << endl;
            }
            break;
        }
    }
    exponent = (Mtx->GFq->q + 2 * q1 - 1) >> 1;
    for (h = 0; h < sz; h++) {
        for (i = 0; i < 8; i++) {
            M[i] = 0;
        }
        if (h == 0) {
            M[0 * 4 + 2] = 1;
            M[1 * 4 + 3] = 1;
        }
        else {
            M[0 * 4 + 0] = 1;
            M[1 * 4 + 1] = 1;
            h1 = h - 1;
            s = h1 % q;
            t = (h1 - s) / q;
            x = s;
            y = t;
            sq = Mtx->GFq->power(s, q1);
            tq = Mtx->GFq->power(t, q1);
            if (type_of_spread == SPREAD_OF_TYPE_DICKSON_KANTOR) {
                w = Mtx->GFq->add(sq, tq);
            }
            else {
                w = Mtx->GFq->add(sq, Mtx->GFq->power(t, exponent));
            }
            z = Mtx->GFq->add(Mtx->GFq->mult(eta, sq), tq);
            M[0 * 4 + 2] = x;
            M[0 * 4 + 3] = y;
            M[1 * 4 + 2] = w;
            M[1 * 4 + 3] = z;
        }
        for (i = 0; i < 8; i++) {
            Grass->M[i] = M[i];
        }
        if (f_vv) {
            cout << "spread element " << h << ":" << endl;
            Int_matrix_print(Grass->M, 2, 4);
        }
        data[h] = Grass->rank_lint(0);
    } // next h
    if (check_function(sz, data, verbose_level - 2)) {
        if (f_v) {
            cout << "spread_classify::make_spread The set is a spread" << endl;
        }
    }
    else {
        cout << "spread_classify::make_spread The set is NOT a spread" << endl;
        exit(1);
    }
 
    if (f_v) {
        cout << "spread_classify::make_spread done" << endl;
    }
}
 
 
void spread_classify::make_spread_from_q_clan(long int *data,
        int type_of_spread, int verbose_level)
{
    int f_v = (verbose_level >= 1);
    int f_vv = (verbose_level >= 2);
    int sz = order + 1;
    int M[8];
    int h, h1, i, s, t, t2, t3, t4, t5, t7, t9;
    int a_t = 0, b_t = 0, c_t = 0, x, y, w, z, r;
    int three, five, nonsquare = 0, minus_nonsquare = 0, nonsquare_inv = 0;
    number_theory::number_theory_domain NT;
 
    if (f_v) {
        cout << "spread_classify::make_spread_from_q_clan" << endl;
    }
    if (n != 4) {
        cout << "spread_classify::make_spread_from_q_clan n != 4" << endl;
        exit(1);
    }
    if (k != 2) {
        cout << "spread_classify::make_spread_from_q_clan k != 2" << endl;
        exit(1);
    }
    three = Mtx->GFq->add3(1, 1, 1);
    five = Mtx->GFq->add3(three, 1, 1);
    if (type_of_spread == SPREAD_OF_TYPE_KANTOR || 
        type_of_spread == SPREAD_OF_TYPE_GANLEY || 
        type_of_spread == SPREAD_OF_TYPE_LAW_PENTTILA) {
        nonsquare = Mtx->GFq->alpha_power(1);
        nonsquare_inv = Mtx->GFq->inverse(nonsquare);
        minus_nonsquare = Mtx->GFq->negate(nonsquare);
    }
    for (h = 0; h < sz; h++) {
        for (i = 0; i < 8; i++) {
            M[i] = 0;
        }
        if (h == 0) {
            M[0 * 4 + 2] = 1;
            M[1 * 4 + 3] = 1;
        }
        else {
            M[0 * 4 + 0] = 1;
            M[1 * 4 + 1] = 1;
            h1 = h - 1;
            s = h1 % q;
            t = (h1 - s) / q;
 
            // create the q-clan:
 
            if (type_of_spread == SPREAD_OF_TYPE_FTWKB) {
                // Fisher Thas Walker Betten:
                // a_t = t, b_t = 3t^2, c_t = 3t^3
                // only when q \equiv 2 mod 3
                r = q % 3;
                if (r != 2) {
                    cout << "FTWKB needs q equiv 2 mod 3" << endl;
                    exit(1);
                }
                a_t = t;
                b_t = Mtx->GFq->product3(three, t, t);
                c_t = Mtx->GFq->mult(b_t, t);
            }
            else if (type_of_spread == SPREAD_OF_TYPE_KANTOR) {
                if (EVEN(q)) {
                    cout << "KANTOR needs q to be odd" << endl;
                    exit(1);
                }
                if (NT.is_prime(q)) {
                    cout << "KANTOR needs q to be a prime power" << endl;
                    exit(1);
                }
                a_t = t;
                b_t = 0;
                c_t = Mtx->GFq->mult(minus_nonsquare, Mtx->GFq->frobenius_power(t, 1));
            }
            else if (type_of_spread == SPREAD_OF_TYPE_KANTOR2) {
                if (EVEN(q)) {
                    cout << "KANTOR2 needs q to be odd" << endl;
                    exit(1);
                }
                if ((q % 5) != 2 && (q % 5) != 3) {
                    cout << "KANTOR needs q congruent 2 or 3 mod 5" << endl;
                    exit(1);
                }
                t3 = Mtx->GFq->product3(t, t, t);
                t5 = Mtx->GFq->product3(t3, t, t);
                a_t = t;
                b_t = Mtx->GFq->mult(five, t3);
                c_t = Mtx->GFq->mult(five, t5);
            }
            else if (type_of_spread == SPREAD_OF_TYPE_GANLEY) {
                int tmp1, tmp2;
                if (q % 3) {
                    cout << "GANLEY, q needs to be "
                            "a power of three" << endl;
                    exit(1);
                }
                t3 = Mtx->GFq->product3(t, t, t);
                t9 = Mtx->GFq->product3(t3, t3, t3);
                a_t = t;
                b_t = t3;
                tmp1 = Mtx->GFq->mult(nonsquare, t);
                tmp2 = Mtx->GFq->mult(nonsquare_inv, t9);
                c_t = Mtx->GFq->negate(Mtx->GFq->add(tmp1, tmp2));
            }
            else if (type_of_spread == SPREAD_OF_TYPE_LAW_PENTTILA) {
                int tmp1, tmp2;
                int n2, n3;
                if (q % 3) {
                    cout << "LAW_PENTTILA, q needs to "
                            "be a power of three" << endl;
                    exit(1);
                }
                t2 = Mtx->GFq->mult(t, t);
                t3 = Mtx->GFq->product3(t, t, t);
                t4 = Mtx->GFq->mult(t2, t2);
                t7 = Mtx->GFq->mult(t4, t3);
                t9 = Mtx->GFq->product3(t3, t3, t3);
                n2 = Mtx->GFq->mult(nonsquare, nonsquare);
                n3 =Mtx->GFq->product3(nonsquare, nonsquare, nonsquare);
                a_t = t;
                b_t = Mtx->GFq->negate(Mtx->GFq->add(t4, Mtx->GFq->mult(nonsquare, t2)));
                tmp1 = Mtx->GFq->add(t7, Mtx->GFq->mult(n2, t3));
                tmp2 = Mtx->GFq->negate(Mtx->GFq->add(Mtx->GFq->mult(nonsquare_inv, t9),
                        Mtx->GFq->mult(n3, t)));
                c_t = Mtx->GFq->add(tmp1, tmp2);
                if (h == 34) {
                    cout << "s=" << s << endl;
                    cout << "t=" << t << endl;
                    cout << "n=" << nonsquare << endl;
                    cout << "t2=" << t2 << endl;
                    cout << "t3=" << t3 << endl;
                    cout << "t4=" << t4 << endl;
                    cout << "t7=" << t7 << endl;
                    cout << "t9=" << t9 << endl;
                    cout << "a_t=" << a_t << endl;
                    cout << "b_t=" << b_t << endl;
                    cout << "c_t=" << c_t << endl;
                }
            }
            
 
            // create the spread element 
            // according to Gevaert-Johnson 1988:
 
            x = a_t;
            y = Mtx->GFq->add(b_t, s);
            w = Mtx->GFq->negate(s);
            z = c_t;
            M[0 * 4 + 2] = x;
            M[0 * 4 + 3] = y;
            M[1 * 4 + 2] = w;
            M[1 * 4 + 3] = z;
        }
        for (i = 0; i < 8; i++) {
            Grass->M[i] = M[i];
        }
        if (f_vv) {
            cout << "spread element " << h << ":" << endl;
            Int_matrix_print(Grass->M, 2, 4);
        }
        data[h] = Grass->rank_lint(0);
    }
    if (check_function(sz, data, verbose_level - 2)) {
        cout << "The set is a spread" << endl;
    }
    else {
        cout << "The set is NOT a spread" << endl;
        exit(1);
    }
}
 
void spread_classify::read_and_print_spread(std::string &fname, int verbose_level)
{
    long int *data;
    int sz;
    orbiter_kernel_system::file_io Fio;
 
    Fio.read_set_from_file(fname, data, sz, verbose_level);
    print_spread(cout, data, sz);
    FREE_lint(data);
}
 
void spread_classify::HMO(std::string &fname, int verbose_level)
{
    long int *data;
    int sz, i, h, h1;
    int *G, *H;
    int *Ge, *He;
    int *GG, *HH;
    int alpha, beta, omega, x, y, tmp1, tmp2, f, z;
    int M[8];
    orbiter_kernel_system::file_io Fio;
 
    if (order != q * q) {
        cout << "spread_classify::print_spread order != q * q" << endl;
        exit(1);
    }
    Fio.read_set_from_file(fname, data, sz, verbose_level);
    G = NEW_int(order);
    H = NEW_int(order);
    Ge = NEW_int(order);
    He = NEW_int(order);
    print_spread(cout, data, sz);
    get_spread_matrices(G, H, data, verbose_level);
 
 
    int q2;
    field_theory::finite_field *Fq2;
    field_theory::subfield_structure *Sub;
 
    q2 = q * q;
    Fq2 = NEW_OBJECT(field_theory::finite_field);
    Sub = NEW_OBJECT(field_theory::subfield_structure);
    Fq2->finite_field_init(q2, FALSE /* f_without_tables */, verbose_level);
    Sub->init(Fq2, Mtx->GFq, verbose_level);
    for (i = 0; i < q * q; i++) {
        Ge[i] = Sub->FQ_embedding[G[i]];
        He[i] = Sub->FQ_embedding[H[i]];
    }
    cout << "spread::HMO after embedding" << endl;
    cout << "Ge:" << endl;
    Int_matrix_print(Ge, q, q);
    cout << "He:" << endl;
    Int_matrix_print(He, q, q);
 
    GG = NEW_int(q2 * q2);
    HH = NEW_int(q2 * q2);
    omega = Sub->Basis[1];
    cout << "omega=" << omega << endl;
    for (alpha = 0; alpha < q2; alpha++) {
        for (beta = 0; beta < q2; beta++) {
            x = Sub->components[beta * 2 + 0];
            y = Sub->components[beta * 2 + 1];
            cout << "alpha=" << alpha << " beta=" << beta
                    << " x=" << x << " y=" << y << endl;
            tmp1 = Ge[x * q + y];
            tmp2 = Fq2->negate(Fq2->mult(He[x * q + y], omega));
            f = Fq2->add(tmp1, tmp2);
            z = Fq2->frobenius_power(alpha, Fq2->e / 2);
            GG[alpha * q2 + beta] = f;
            HH[alpha * q2 + beta] = z;
        }
    }
    cout << "GG:" << endl;
    Int_matrix_print(GG, q2, q2);
    cout << "HH:" << endl;
    Int_matrix_print(HH, q2, q2);
 
    geometry::grassmann *Gq2;
    long int *Data2;
    int Sz;
 
    Gq2 = NEW_OBJECT(geometry::grassmann);
    Gq2->init(n, k, Fq2, verbose_level);
 
    Sz = q2 * q2 + 1;
    Data2 = NEW_lint(Sz);
 
    for (h = 0; h < Sz; h++) {
        for (i = 0; i < 8; i++) {
            M[i] = 0;
        }
        if (h == 1) {
            M[0 * 4 + 2] = 1;
            M[1 * 4 + 3] = 1;
        }
        else {
            if (h > 1) {
                h1 = h - 1;
            }
            else {
                h1 = 0;
            }
            M[0 * 4 + 0] = 1;
            M[1 * 4 + 1] = 1;
            y = h1 % q2;
            x = (h1 - y) / q2;
            cout << "h=" << h << " x=" << x << " y=" << y << endl;
            M[0 * 4 + 2] = x;
            M[0 * 4 + 3] = y;
            M[1 * 4 + 2] = GG[x * q2 + y];
            M[1 * 4 + 3] = HH[x * q2 + y];
        }
        cout << "element " << h << ":" << endl;
        Int_matrix_print(M, 2, 4);
        for (i = 0; i < 8; i++) {
            Gq2->M[i] = M[i];
        }
        Data2[h] = Gq2->rank_lint(0);
        cout << "has rank " << Data2[h] << endl;
    }
 
    string fname2;
 
    fname2.assign("HMO_");
    fname2.append(fname);
 
    Fio.write_set_to_file(fname2, Data2, Sz, verbose_level);
    
    FREE_lint(Data2);
    FREE_OBJECT(Gq2);
    FREE_OBJECT(Fq2);
    FREE_OBJECT(Sub);
    FREE_lint(data);
    FREE_int(G);
    FREE_int(H);
    FREE_int(Ge);
    FREE_int(He);
    FREE_int(GG);
    FREE_int(HH);
}
 
void spread_classify::get_spread_matrices(int *G, int *H,
        long int *data, int verbose_level)
{
    int f_v = (verbose_level >= 1);
    int h, i, x, y, f, g;
    int sz = order + 1;
    int M[8];
 
    if (f_v) {
        cout << "spread_classify::get_spread_matrices" << endl;
    }
    for (i = 0; i < order; i++) {
        G[i] = -1;
        H[i] = -1;
    }
    for (h = 0; h < sz; h++) {
        Grass->unrank_lint(data[h], 0);
        for (i = 0; i < 8; i++) {
            M[i] = Grass->M[i];
        }
        //cout << "element " << h << ":" << endl;
        //int_matrix_print(M, 2, 4);
        if (M[0 * 4 + 0] == 0 && 
            M[0 * 4 + 1] == 0 && 
            M[1 * 4 + 0] == 0 && 
            M[1 * 4 + 1] == 0) {
            continue;
        }
        if (M[0 * 4 + 0] != 1) {
            cout << "generator matrix is not in standard form" << endl;
            exit(1);
        }
        if (M[1 * 4 + 1] != 1) {
            cout << "generator matrix is not in standard form" << endl;
            exit(1);
        }
        if (M[0 * 4 + 1] != 0) {
            cout << "generator matrix is not in standard form" << endl;
            exit(1);
        }
        if (M[1 * 4 + 0] != 0) {
            cout << "generator matrix is not in standard form" << endl;
            exit(1);
        }
        x = M[0 * 4 + 2];
        y = M[0 * 4 + 3];
        f = M[1 * 4 + 2];
        g = M[1 * 4 + 3];
        G[x * q + y] = f;
        H[x * q + y] = g;
    }
    if (f_v) {
        cout << "spread::get_FG_matrices" << endl;
        cout << "G:" << endl;
        Int_matrix_print(G, q, q);
        cout << "H:" << endl;
        Int_matrix_print(H, q, q);
    }
}
 
void spread_classify::print_spread(ostream &ost, long int *data, int sz)
{
    //int sz = order + 1;
    int h;
 
    for (h = 0; h < sz; h++) {
        Grass->unrank_lint(data[h], 0);
        ost << "Spread element " << h << ":" << endl;
        orbiter_kernel_system::Orbiter->Int_vec->matrix_print_ost(ost, Grass->M, k, n);
    }
}
 
void spread_classify::report2(isomorph &Iso, int verbose_level)
{
    int f_v = (verbose_level >= 1);
    char fname[1000];
 
    if (f_v) {
        cout << "spread_classify::report2" << endl;
    }
    sprintf(fname, "report_Spreads_q%d_k%d.tex", q, k);
 
    {
    ofstream f(fname);
    int f_book = TRUE;
    int f_title = TRUE;
    char title[1000];
    const char *author = "Orbiter";
    int f_toc = TRUE;
    int f_landscape = FALSE;
    int f_12pt = FALSE;
    int f_enlarged_page = TRUE;
    int f_pagenumbers = TRUE;
    orbiter_kernel_system::latex_interface L;
 
    sprintf(title, "$%d$-Spreads of PG($%d,%d$)", k - 1, 2 * k - 1, q);
    cout << "Writing file " << fname << " with "
            << Iso.Reps->count << " spreads:" << endl;
    L.head(f, f_book, f_title,
        title, author, 
        f_toc, f_landscape, f_12pt, f_enlarged_page, f_pagenumbers, 
        NULL /* extra_praeamble */);
 
    report3(Iso, f, verbose_level);
 
    L.foot(f);
 
    }
    orbiter_kernel_system::file_io Fio;
 
    cout << "Written file " << fname << " of size "
            << Fio.file_size(fname) << endl;
 
}
 
void spread_classify::report3(isomorph &Iso, ostream &ost, int verbose_level)
{
    int f_v = (verbose_level >= 1);
 
    if (f_v) {
        cout << "spread_classify::report3" << endl;
    }
 
    int target_size = order + 1;
    data_structures::sorting Sorting;
 
    ost << "\\chapter{Summary}" << endl << endl;
    ost << "There are " << Iso.Reps->count << " spreads." << endl << endl;
 
 
    //Iso.setup_and_open_solution_database(verbose_level - 1);
 
    int i, first, /*c,*/ id;
    int u, v, h, rep, tt;
    ring_theory::longinteger_object go;
    long int data[1000];
 
 
 
    ring_theory::longinteger_object *Ago, *Ago_induced;
    int *Ago_int;
 
    Ago = NEW_OBJECTS(ring_theory::longinteger_object, Iso.Reps->count);
    Ago_induced = NEW_OBJECTS(ring_theory::longinteger_object, Iso.Reps->count);
    Ago_int = NEW_int(Iso.Reps->count);
 
 
    for (h = 0; h < Iso.Reps->count; h++) {
        rep = Iso.Reps->rep[h];
        first = Iso.orbit_fst[rep];
        //c = Iso.starter_number[first];
        id = Iso.orbit_perm[first];     
        Iso.load_solution(id, data);
 
        groups::sims *Stab;
        
        Stab = Iso.Reps->stab[h];
 
        Iso.Reps->stab[h]->group_order(Ago[h]);
        Ago_int[h] = Ago[h].as_int();
        //f << "Stabilizer has order $";
        //go.print_not_scientific(f);
        if (f_v) {
            cout << "spread::print_isomorphism_types computing "
                    "induced action on the set (in data)" << endl;
        }
        Iso.induced_action_on_set(Stab, data, 0 /*verbose_level*/);
        
            
        Iso.AA->group_order(Ago_induced[h]);
    }
 
 
    data_structures::tally C_ago;
 
    C_ago.init(Ago_int, Iso.Reps->count, FALSE, 0);
    cout << "Classification by ago:" << endl;
    C_ago.print(FALSE /*f_backwards*/);
 
 
 
    ost << "\\chapter{Invariants}" << endl << endl;
 
    ost << "Classification by automorphism group order: $";
    C_ago.print_naked_tex(ost, FALSE /*f_backwards*/);
    ost << "$" << endl;
    ost << "\\\\" << endl;
 
    ost << "\\begin{center}" << endl;
    ost << "\\begin{tabular}{|c|l|}" << endl;
    ost << "\\hline" << endl;
    ost << "Ago & Isom. Types \\\\" << endl;
    ost << "\\hline" << endl;
    ost << "\\hline" << endl;
 
    int cnt, length, t, vv, *set;
 
    cnt = 0;
    for (u = C_ago.nb_types - 1; u >= 0; u--) {
        first = C_ago.type_first[u];
        length = C_ago.type_len[u];
        t = C_ago.data_sorted[first];
 
        ost << t << " & ";
 
        set = NEW_int(length);
        for (v = 0; v < length; v++, cnt++) {
            vv = first + v;
            i = C_ago.sorting_perm_inv[vv];
            set[v] = i;
        }
 
        Sorting.int_vec_heapsort(set, length);
 
        for (v = 0; v < length; v++, cnt++) {
 
            ost << set[v];
 
            if (v < length - 1) {
                ost << ",";
                if ((v + 1) % 10 == 0) {
                    ost << "\\\\" << endl;
                    ost << " & " << endl;
                }
            }
        }
        ost << "\\\\" << endl;
        if (u > 0) {
            ost << "\\hline" << endl;
        }
        FREE_int(set);
    }
    ost << "\\hline" << endl;
    ost << "\\end{tabular}" << endl;
    ost << "\\end{center}" << endl << endl;
 
 
    ost << "\\clearpage" << endl << endl;
 
    ost << "\\begin{center}" << endl;
    ost << "\\begin{tabular}{|r|r|r|r|r|r|}" << endl;
    ost << "\\hline" << endl;
    ost << "ID & Isom. Type & $|\\mbox{Aut}|$ & $|\\mbox{Aut}|$ "
            "(induced) & \\# Orbits & Orbit Type\\\\" << endl;
    ost << "\\hline" << endl;
    ost << "\\hline" << endl;
 
    cnt = 0;
    for (u = 0; u < C_ago.nb_types; u ++) {
        first = C_ago.type_first[u];
        length = C_ago.type_len[u];
        t = C_ago.data_sorted[first];
 
        set = NEW_int(length);
        for (v = 0; v < length; v++) {
            vv = first + v;
            i = C_ago.sorting_perm_inv[vv];
            set[v] = i;
        }
 
        Sorting.int_vec_heapsort(set, length);
 
 
        for (v = 0; v < length; v++) {
            vv = first + v;
            i = C_ago.sorting_perm_inv[vv];
            h = set[v];
            ost << setw(3) << first + v + 1 << " & ";
            ost << setw(3) << h << " & ";
            Ago[h].print_not_scientific(ost);
            ost << " & ";
            Ago_induced[h].print_not_scientific(ost);
 
            orbit_info_short(ost, Iso, h);
 
 
            ost << "\\\\" << endl;
            cnt++;
            if ((cnt % 30) == 0) {
                ost << "\\hline" << endl;
                ost << "\\end{tabular}" << endl;
                ost << "\\end{center}" << endl << endl;
                ost << "\\begin{center}" << endl;
                ost << "\\begin{tabular}{|r|r|r|r|r|r|}" << endl;
                ost << "\\hline" << endl;
                ost << "ID & Isom. Type & $|\\mbox{Aut}|$ & $|\\mbox{Aut}|$ "
                        "(induced) & \\# Orbits & Orbit Type\\\\" << endl;
                ost << "\\hline" << endl;
                ost << "\\hline" << endl;
            }
        }
        FREE_int(set);
    }
 
    ost << "\\hline" << endl;
    ost << "\\end{tabular}" << endl;
    ost << "\\end{center}" << endl << endl;
 
 
    ost << "\\chapter{The Spreads}" << endl << endl;
 
    ost << "\\clearpage" << endl << endl;
 
 
    for (h = 0; h < Iso.Reps->count; h++) {
        rep = Iso.Reps->rep[h];
        first = Iso.orbit_fst[rep];
        //c = Iso.starter_number[first];
        id = Iso.orbit_perm[first];     
        Iso.load_solution(id, data);
 
 
        ost << "\\section{Isomorphism type " << h << "}" << endl;
        ost << "\\bigskip" << endl;
 
        //f << "The ambient space has dimension "
        //<< Rk_of_span[h] << "\\\\" << endl;
 
        if (Iso.Reps->stab[h]) {
            Iso.Reps->stab[h]->group_order(go);
            ost << "Stabilizer has order $";
            go.print_not_scientific(ost);
            ost << "$\\\\" << endl;
        }
        else {
            //cout << endl;
        }
 
#if 0
        if (h == 35) {
            cout << "Isomorphism type 35:" << endl;
            cout << "rep=" << rep << endl;
            cout << "first=" << first << endl;
            cout << "id=" << id << endl;
            cout << "starter_number=" << c << endl;
            cout << "solution:" << endl;
            int_vec_print(cout, data, target_size);
            cout << endl;
        }
#endif
 
        groups::sims *Stab;
        
        Stab = Iso.Reps->stab[h];
 
        if (f_v) {
            cout << "spread::print_isomorphism_types computing "
                    "induced action on the set (in data)" << endl;
        }
        Iso.induced_action_on_set(Stab, data, 0 /*verbose_level*/);
        
        ring_theory::longinteger_object go1;
            
        Iso.AA->group_order(go1);
        cout << "action " << Iso.AA->label << " computed, "
                "group order is " << go1 << endl;
 
        ost << "Order of the group that is induced on the "
                "collection of subspaces is ";
        ost << "$";
        go1.print_not_scientific(ost);
        ost << "$\\\\" << endl;
        
 
        groups::schreier Orb;
        //longinteger_object go2;
        
        Iso.AA->compute_all_point_orbits(Orb,
                Stab->gens, verbose_level - 2);
        ost << "With " << Orb.nb_orbits
                << " orbits on the subspaces\\\\" << endl;
 
        data_structures::tally C_ol;
 
        C_ol.init(Orb.orbit_len, Orb.nb_orbits, FALSE, 0);
 
        ost << "Orbit lengths: $";
        //int_vec_print(f, Orb.orbit_len, Orb.nb_orbits);
        C_ol.print_naked_tex(ost, FALSE /*f_backwards*/);
        ost << "$ \\\\" << endl;
 
        ost << "Orbits: $";
        Orb.list_all_orbits_tex(ost);
        ost << "$ \\\\" << endl;
 
    
        tt = (target_size + 3) / 4;
 
        ost << "The subspaces by ranks:\\\\" << endl;
        ost << "\\begin{center}" << endl;
 
        for (u = 0; u < 4; u++) {
            ost << "\\begin{tabular}[t]{|c|c|}" << endl;
            ost << "\\hline" << endl;
            ost << "$i$ & Rank \\\\" << endl;
            ost << "\\hline" << endl;
            for (i = 0; i < tt; i++) {
                v = u * tt + i;
                if (v < target_size) {
                    ost << "$" << v << "$ & $" << data[v] << "$ \\\\" << endl;
                }
            }
            ost << "\\hline" << endl;
            ost << "\\end{tabular}" << endl;
        }
        ost << "\\end{center}" << endl;
 
        ost << "The subspaces:\\\\" << endl;
        for (i = 0; i < target_size; i++) {
            Grass->unrank_lint(data[i], 0/*verbose_level - 4*/);
            if ((i % 4) == 0) {
                if (i) {
                    ost << "$$" << endl;
                }
                ost << "$$" << endl;
            }
            ost << "\\left[" << endl;
            ost << "\\begin{array}{c}" << endl;
            for (u = 0; u < k; u++) {
                for (v = 0; v < n; v++) {
                    ost << Grass->M[u * n + v];
                }
                ost << "\\\\" << endl;
            }
            ost << "\\end{array}" << endl;
            ost << "\\right]" << endl;
        }
        ost << "$$" << endl;
 
 
        report_stabilizer(Iso, ost, h /* orbit */, 0 /* verbose_level */);
 
        if (k == 2) {
            klein(ost, &Iso, h, Stab, Orb, data, Iso.size,  verbose_level);
        }
 
        ost << "" << endl;
        ost << "\\bigskip" << endl;
        ost << "" << endl;
        ost << "Induced action on the orbits:\\\\" << endl;
        ost << "$$" << endl;
        ost << "\\begin{array}{|c|c|c|c|}" << endl;
        ost << "\\hline" << endl;
        ost << "\\mbox{Orbit} & \\mbox{Length} & \\mbox{Induced} "
                "& \\mbox{Kernel} \\\\" << endl;
        ost << "\\hline" << endl;
        for (i = 0; i < Orb.nb_orbits; i++) {
            int fst, len, j;
            fst = Orb.orbit_first[i];
            len = Orb.orbit_len[i];
 
            cout << "inducing action on orbit " << i
                    << " of length " << len << endl;
            ost << i << " & " << len << " & ";
            long int *set;
            actions::action *A1;
            ring_theory::longinteger_object go1, gok;
 
            set = NEW_lint(len);
            //A1 = NEW_OBJECT(action);
            for (j = 0; j < len; j++) {
                set[j] = data[Orb.orbit[fst + j]];
            }
            A1 = Iso.A->create_induced_action_by_restriction(
                    Stab,
                    len, set,
                    TRUE,
                    0/*verbose_level*/);
            if (f_v) {
                cout << "after induced_action_by_restriction" << endl;
            }
            A1->group_order(go1);
            A1->Kernel->group_order(gok);
            cout << "induced group order " << go1
                    << " kernel order " << gok << endl;
 
            ost << go1 << " & " << gok << "\\\\" << endl;
 
            FREE_lint(set);
            FREE_OBJECT(A1);
        }
        ost << "\\hline" << endl;
        ost << "\\end{array}" << endl;
        ost << "$$" << endl;
 
 
 
        for (i = 0; i < Orb.nb_orbits; i++) {
            int fst, len, j;
            fst = Orb.orbit_first[i];
            len = Orb.orbit_len[i];
 
            cout << "inducing action on orbit " << i
                    << " of length " << len << endl;
            ost << "Generators for the kernel in the action "
                    "on orbit " << i << " are:";
            long int *set;
            actions::action *A1;
            data_structures_groups::vector_ge *gens;
            int *tl;
            ring_theory::longinteger_object go1, gok;
 
            set = NEW_lint(len);
            //A1 = NEW_OBJECT(action);
            gens = NEW_OBJECT(data_structures_groups::vector_ge);
            tl = NEW_int(Iso.A_base->base_len());
            for (j = 0; j < len; j++) {
                set[j] = data[Orb.orbit[fst + j]];
            }
            A1 = Iso.A->create_induced_action_by_restriction(
                    Stab,
                    len, set,
                    TRUE,
                    0/*verbose_level*/);
            if (f_v) {
                cout << "after induced_action_by_restriction" << endl;
            }
            A1->group_order(go1);
            A1->Kernel->group_order(gok);
            cout << "induced group order " << go1
                    << " kernel order " << gok << endl;
            A1->Kernel->extract_strong_generators_in_order(*gens, tl, 0);
 
            if (gens->len == 0) {
                ost << " No generators\\\\" << endl;
            }
            for (j = 0; j < gens->len; j++) {
        
 
                ost << "$$ g_{" << j + 1 << "}=" << endl;
                A->element_print_latex(gens->ith(j), ost);
                ost << "$$" << endl;
            }
            
            FREE_lint(set);
            FREE_OBJECT(A1);
            FREE_OBJECT(gens);
            FREE_int(tl);
        }
 
 
    }
 
    //Iso.close_solution_database(verbose_level - 1);
 
    char prefix[1000];
    char label_of_structure_plural[1000];
 
    sprintf(prefix, "Spreads_%d_%d", q, k);
    sprintf(label_of_structure_plural, "Spreads");
 
 
    isomorph_global IG;
 
    IG.init(Iso.A_base, Iso.A, Iso.gen, verbose_level);
 
    IG.report_data_in_source_code_inside_tex(Iso,
        prefix, label_of_structure_plural, ost, verbose_level);
 
 
 
    FREE_int(Ago_int);
    FREE_OBJECTS(Ago);
    FREE_OBJECTS(Ago_induced);
 
    if (f_v) {
        cout << "spread_classify::report3 done" << endl;
    }
}
 
 
void spread_classify::all_cooperstein_thas_quotients(
        isomorph &Iso, int verbose_level)
{
    int f_v = (verbose_level >= 1);
    int h;
    char fname[1000];
    int cnt = 0;
 
    if (f_v) {
        cout << "spread_classify::all_cooperstein_thas_quotients" << endl;
    }
 
 
    //Iso.setup_and_open_solution_database(verbose_level - 1);
 
 
    sprintf(fname, "quotients_q%d.txt", order);
    {
        ofstream f(fname);
 
 
        for (h = 0; h < Iso.Reps->count; h++) {
 
            cooperstein_thas_quotients(Iso, f, h, cnt, verbose_level - 1);
 
        }
 
 
        f << -1 << endl;
    }
    orbiter_kernel_system::file_io Fio;
 
    cout << "Written file " << fname << " of size "
            << Fio.file_size(fname) << endl;
 
 
    //Iso.close_solution_database(verbose_level - 1);
 
    if (f_v) {
        cout << "spread::all_cooperstein_thas_quotients done" << endl;
    }
}
 
 
 
void spread_classify::cooperstein_thas_quotients(isomorph &Iso,
        ofstream &f, int h, int &cnt, int verbose_level)
{
    int f_v = (verbose_level >= 1);
    int f_vv = (verbose_level >= 2);
    int rep, first, /*c,*/ id;
    groups::sims *Stab;
    long int data[1000];
    long int data2[1000];
    int u, fst, orbit_length, i, idx, i0, j, the_point, pivot, rem;
    long int a, b;
    int *M;
    int *vec1;
    int *vec2;
    long int **Pts;
    int nb_points = 0;
    int *List;
    geometry::grassmann *Gr;
    ring_theory::longinteger_domain Dom;
    number_theory::number_theory_domain NT;
    data_structures::sorting Sorting;
    orbiter_kernel_system::file_io Fio;
 
 
 
    if (f_v ) {
        cout << "spread_classify::cooperstein_thas_quotients" << endl;
        cout << "spread_classify::cooperstein_thas_quotients h=" << h << endl;
    }
 
    rep = Iso.Reps->rep[h];
    first = Iso.orbit_fst[rep];
    //c = Iso.starter_number[first];
    id = Iso.orbit_perm[first];     
    Iso.load_solution(id, data);
 
 
    Stab = Iso.Reps->stab[h];
 
    groups::schreier Orb;
    ring_theory::longinteger_object go;
    ring_theory::longinteger_object stab_order;
        
    Stab->group_order(go);
 
    if (f_v) {
        cout << "spread_classify::cooperstein_thas_quotients "
                "Isomorphism type " << h << " with automorphism group "
                        "of order " << go << ":" << endl;
    }
 
    if (f_vv) {
        cout << "spread_classify::cooperstein_thas_quotients "
                "before compute_all_point_orbits" << endl;
    }
 
    Iso.A_base->compute_all_point_orbits(Orb,
            Stab->gens, 0 /*verbose_level - 2*/);
 
    if (f_vv) {
        cout << "spread_classify::cooperstein_thas_quotients There are "
                << Orb.nb_orbits << " orbits on points" << endl;
    }
 
 
    
    M = NEW_int(k * n);
    vec1 = NEW_int(n);
    vec2 = NEW_int(n);
    Pts = NEW_plint(order + 1);
    List = NEW_int(order);
    Gr = NEW_OBJECT(geometry::grassmann);
 
    Gr->init(n - 1, k, Mtx->GFq, 0 /* verbose_level */);
    for (i = 0; i < order + 1; i++) {
        Grass->unrank_lint_here(M, data[i], 0/*verbose_level - 4*/);
        Mtx->GFq->all_PG_elements_in_subspace(M, k, n,
                Pts[i], nb_points, 0 /* verbose_level */);
        Sorting.lint_vec_heapsort(Pts[i], nb_points);
    }
 
    for (u = 0; u < Orb.nb_orbits; u++) {
        fst = Orb.orbit_first[u];
        orbit_length = Orb.orbit_len[u];
        
        the_point = Orb.orbit[fst];
        Dom.integral_division_by_int(go, orbit_length, stab_order, rem);
 
        if (f_vv) {
            cout << "spread_classify::cooperstein_thas_quotients Orbit " << u
                    << " is represented by point " << the_point
                    << " orbit length = " << orbit_length
                    << " stabilizer order before " << go
                    << " after " << stab_order << endl;
        }
 
 
        Mtx->GFq->PG_element_unrank_modified(vec1, 1, n, the_point);
        Mtx->GFq->PG_element_normalize_from_front(vec1, 1, n);
        pivot = orbiter_kernel_system::Orbiter->Int_vec->find_first_nonzero_entry(vec1, n);
 
        for (i = 0; i < order + 1; i++) {
            if (Sorting.lint_vec_search(Pts[i], nb_points, the_point, idx, 0)) {
                break;
            }
        }
        if (i == order + 1) {
            cout << "spread_classify::cooperstein_thas_quotients "
                    "Did not find the point" << endl;
            exit(1);
        }
        i0 = i;
        j = 0;
        for (i = 0; i < order + 1; i++) {
            if (i != i0) {
                List[j++] = i;
            }
        }
        if (j != order) {
            cout << "spread_classify::cooperstein_thas_quotients j != order" << endl;
            exit(1);
        }
        for (i = 0; i < order; i++) {
            a = data[List[i]];
            Grass->unrank_lint_here(M, a, 0/*verbose_level - 4*/);
            Int_vec_copy(vec1, vec2, n);
            for (j = 0; j < k; j++) {
                Mtx->GFq->Linear_algebra->Gauss_step(vec2, M + j * n, n, pivot,
                        0 /* verbose_level*/);
                    // afterwards: v2[idx] = 0 and
                    // v1,v2 span the same space as before
                    // v1 is not changed if v1[idx] is nonzero
            }
 
#if 1
            // delete column 'pivot' in the k x n matrix M.
            // Afterwards, the matrix is k x (n - 1)
 
            orbiter_kernel_system::Orbiter->Int_vec->matrix_delete_column_in_place(M, k, n, pivot);
#endif
 
            if (f_vv) {
                cout << "spread_classify::cooperstein_thas_quotients the reduction "
                        "of the " << i << "-th matrix is:" << endl;
                Int_matrix_print(M, k, n - 1);
            }
 
            b = Gr->rank_lint_here(M, 0/*verbose_level - 4*/);
            data2[i] = b;
        }
 
        if (f_vv) {
            cout << "spread_classify::cooperstein_thas_quotients The quotient "
                    "system with respect to orbit " << u << " / "
                    << Orb.nb_orbits << " is:" << endl;
            Lint_vec_print(cout, data2, order);
            cout << endl;
        }
 
 
        f << h << " " << u << " " << the_point << " " << orbit_length << " ";
        go.print_not_scientific(f);
        f << " ";
        stab_order.print_not_scientific(f);
        f << " ";
        for (i = 0; i < order; i++) {
            f << " " << data2[i];
        }
        f << endl;
 
 
 
        char str[1000];
        sprintf(str, "quotient_q%d_iso%d_nb%d_orbit_length%d.txt",
                NT.i_power_j(q, k), h, u, orbit_length);
 
        string fname;
 
        fname.assign(str);
        Fio.write_set_to_file(fname, data2, order, 0 /* verbose_level*/);
        cnt++;
 
 
    } // next u
 
    FREE_int(M);
    FREE_int(vec1);
    FREE_int(vec2);
    for (i = 0; i < order + 1; i++) {
        FREE_lint(Pts[i]);
    }
    FREE_plint(Pts);
    FREE_int(List);
    FREE_OBJECT(Gr);
    if (f_v ) {
        cout << "spread_classify::cooperstein_thas_quotients done" << endl;
    }
}
 
void spread_classify::orbit_info_short(ostream &ost, isomorph &Iso, int h)
{
    int rep, first, /*c,*/ id;
 
    groups::sims *Stab;
        
    long int data[1000];
 
    rep = Iso.Reps->rep[h];
    first = Iso.orbit_fst[rep];
    //c = Iso.starter_number[first];
    id = Iso.orbit_perm[first];     
    Iso.load_solution(id, data);
 
 
    Stab = Iso.Reps->stab[h];
 
    Iso.induced_action_on_set(Stab, data, 0 /*verbose_level*/);
        
    ring_theory::longinteger_object go1;
            
    Iso.AA->group_order(go1);
    cout << "action " << Iso.AA->label << " computed, "
            "group order is " << go1 << endl;
 
 
#if 0
    f << "Order of the group that is induced on the "
            "collection of subspaces is ";
    f << "$";
    go1.print_not_scientific(f);
    f << "$\\\\" << endl;
#endif  
 
    groups::schreier Orb;
    //longinteger_object go2;
        
    Iso.AA->compute_all_point_orbits(Orb, Stab->gens, 0 /*verbose_level - 2*/);
    //f << "With " << Orb.nb_orbits << " orbits on the subspaces\\\\" << endl;
 
    ost << " & " << Orb.nb_orbits << " & ";
 
    data_structures::tally C_ol;
 
    C_ol.init(Orb.orbit_len, Orb.nb_orbits, FALSE, 0);
 
    //f << "Orbit lengths: ";
    //int_vec_print(f, Orb.orbit_len, Orb.nb_orbits);
    ost << "$";
    C_ol.print_naked_tex(ost, FALSE /*f_backwards*/);
    ost << "$";
    //f << " \\\\" << endl;
    
}
 
 
 
void spread_classify::report_stabilizer(isomorph &Iso,
        ostream &ost, int orbit, int verbose_level)
{
    groups::sims *Stab;
    ring_theory::longinteger_object go;
    int i;
 
    Stab = Iso.Reps->stab[orbit];
    Stab->group_order(go);
 
    ost << "The stabilizer of order $" << go << "$ is generated by:\\\\" << endl;
    for (i = 0; i < Stab->gens.len; i++) {
        
        int *fp, n, ord;
        
        fp = NEW_int(A->degree);
        n = A->find_fixed_points(Stab->gens.ith(i), fp, 0);
        //cout << "with " << n << " fixed points" << endl;
        FREE_int(fp);
 
        ord = A->element_order(Stab->gens.ith(i));
 
        ost << "$$ g_{" << i + 1 << "}=" << endl;
        A->element_print_latex(Stab->gens.ith(i), ost);
        ost << "$$" << endl << "of order $" << ord << "$ and with "
                << n << " fixed points" << endl;
    }
}
 
 
 
// #############################################################################
// global functions:
// #############################################################################
 
 
 
 
#if 0
static int spread_check_function_callback(int len,
        long int *S, void *data, int verbose_level)
{
    spread_classify *Spread = (spread_classify *) data;
    int f_OK;
    
    f_OK = Spread->check_function(len, S, verbose_level);
    return f_OK;
}
 
static void spread_callback_report(isomorph *Iso, void *data, int verbose_level)
{
    spread_classify *Spread = (spread_classify *) data;
    
    Spread->report2(*Iso, verbose_level);
}
 
static void spread_callback_make_quotients(isomorph *Iso,
        void *data, int verbose_level)
{
    spread_classify *Spread = (spread_classify *) data;
    
    Spread->all_cooperstein_thas_quotients(*Iso, verbose_level);
}
 
static void callback_spread_print(ostream &ost, int len, long int *S, void *data)
{
    spread_classify *Spread = (spread_classify *) data;
 
    Spread->print_spread(ost, S, len);
}
#endif
 
}}}