text:
# GAP is case sensitive
ThreeDogs := function()
	local dog, Dog, DOG;
	dog := "Benjamin";
	Dog := "Samba";
	DOG := "Bernie";
	if dog = DOG then
		Print("There is just one dog named ", dog, "\n");
	else
		Print("The three dogs are named ", dog, ", ", Dog, " and ", DOG, "\n");
	fi;
end;

ThreeDogs();
# The three dogs are named Benjamin, Samba and Bernie

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

Class(DMACodegen, DefaultCodegen, rec(
    Formula := meth(self, o, y, x, opts)
        local icode, datas, prog, params, sub, initsub, io, t, initcode;
        [x, y] := self.initXY(x, y, opts);

        Add(x.t.qualifiers, opts.memModifier);
        Add(y.t.qualifiers, opts.memModifier);

        o := o.child(1);
        params := Set(Concatenation(Collect(o, param), Filtered(Collect(o, var), IsParallelLoopIndex)));

        datas := Collect(o, FDataOfs);
        [o,t] := UTimedAction(BlockSumsOpts(o, opts)); #PrintLine("BlockSums ", t);
        [icode,t] := UTimedAction(self(o, y, x, opts)); #PrintLine("codegen ", t);
        [icode,t] := UTimedAction(ESReduce(icode, opts)); #PrintLine("ESReduce ", t);
        icode := RemoveAssignAcc(icode);
        Unbind(Compile.times);
        [icode,t] := UTimedAction(BlockUnroll(icode, opts)); #PrintLine("BlockUnroll ", t);
        #PrintLine("---compile--");

        icode := DeclareHidden(icode);
        if IsBound(opts.isFixedPoint) and opts.isFixedPoint then
            icode := FixedPointCode(icode, opts.bits, opts.fracbits);
        fi;

        io := When(x=y, [x], [y, x]);
        sub := Cond(IsBound(opts.subNameDMA), opts.subNameDMA, "sub_dma");
        icode := func(TVoid, sub, Concatenation(io, params), icode);

        return icode;
    end,

    LSKernel := meth(self, o, y, x, opts)
        return chain(dma_signal(self.swp_var), cpu_wait(self.swp_var));
    end,

    DMAGath := meth(self, o, y, x, opts)
        local i, func, rfunc, ix, size;
        Add(self.loadbuffers, y);
        Add(self.membuffers, x);
        i := Ind();

        func := o.func;
        size := 1;
        if ObjId(func) = fTensor and ObjId(Last(func.children())) = fId then
            size := Last(func.children()).domain();
            func := fTensor(Concat(DropLast(func.children(), 1), [fBase(size, 0)]));
        fi;

        func := func.lambda();
        return loop(i, func.domain(), dma_load(y+(i*size), x+func.at(i), size));
    end,

    DMAScat := meth(self, o, y, x, opts)
        local i, func, rfunc, ix, size;

        Add(self.storebuffers, x);
        Add(self.membuffers, y);
        i := Ind();

        func := o.func;
        size := 1;
        if ObjId(func) = fTensor and ObjId(Last(func.children())) = fId then
            size := Last(func.children()).domain();
            func := fTensor(Concat(DropLast(func.children(), 1), [fBase(size, 0)]));
        fi;

        func := func.lambda();
        return loop(i, func.domain(), dma_store(y+func.at(i), x+(i*size), size));
    end,

    DMAFence := (self,o,y,x,opts) >> chain(self(o.child(1), y, x, opts), dma_fence()),

    swp_var := false,
    SWPSum := meth(self, o, y, x, opts)
        local old_swp, c;

        old_swp := self.swp_var;
        self.swp_var := o.var;
        c := swp_loop(o.var, o.domain, self(o.child(1), y, x, opts));
        self.swp_var := old_swp;
        return c;
    end

));

Class(CPUCodegen, DefaultCodegen, rec(
    Formula := meth(self, o, y, x, opts)
        local icode, datas, prog, params, sub, initsub, io, t, initcode;
        [x, y] := self.initXY(x, y, opts);

        o := o.child(1);
        params := Set(Concatenation(Collect(o, param), Filtered(Collect(o, var), IsParallelLoopIndex)));

        datas := Collect(o, FDataOfs);
        [o,t] := UTimedAction(BlockSumsOpts(o, opts)); #PrintLine("BlockSums ", t);
        [icode,t] := UTimedAction(self(o, y, x, opts)); #PrintLine("codegen ", t);
        [icode,t] := UTimedAction(ESReduce(icode, opts)); #PrintLine("ESReduce ", t);
        icode := RemoveAssignAcc(icode);
        Unbind(Compile.times);
        [icode,t] := UTimedAction(BlockUnroll(icode, opts)); #PrintLine("BlockUnroll ", t);
        #PrintLine("---compile--");

        icode := DeclareHidden(icode);
        if IsBound(opts.isFixedPoint) and opts.isFixedPoint then
            icode := FixedPointCode(icode, opts.bits, opts.fracbits);
        fi;

        io := When(x=y, [x], [y, x]);
        Add(params, Filtered(icode.free(), i->not IsBound(i.init)));
        sub := Cond(IsBound(opts.subNameCompute), opts.subNameCompute, "sub_cpu");

        icode := func(TVoid, sub, params, icode);

        return icode;

    end,

    DMAGath := meth(self, o, y, x, opts)
        local dma;
        if not IsBound(x.t.qualifiers) then x.t.qualifiers := []; fi;
        if not IsBound(y.t.qualifiers) then y.t.qualifiers := []; fi;
        Add(x.t.qualifiers, opts.memModifier);
        y.t.qualifiers := [ opts.scratchModifier ];
        Add(self.loadbuffers, y);
        return dma_wait(self.swp_var);
    end,

    DMAScat := meth(self, o, y, x, opts)
        local dma;
        if not IsBound(x.t.qualifiers) then x.t.qualifiers := []; fi;
        if not IsBound(y.t.qualifiers) then y.t.qualifiers := []; fi;
        x.t.qualifiers := [ opts.scratchModifier ];
        Add(y.t.qualifiers, opts.memModifier);
        Add(self.storebuffers, x);
        return cpu_signal(self.swp_var);
    end,

    LSKernel := (self, o, y, x, opts) >> self(o.child(1), y, x, opts),

    swp_var := false,
    SWPSum := meth(self, o, y, x, opts)
        local old_swp, c;

        old_swp := self.swp_var;
        self.swp_var := o.var;
        c := swp_loop(o.var, o.domain, self(o.child(1), y, x, opts));
        self.swp_var := old_swp;
        return c;
    end
));

_scratch_variables := function(vars, opts)
   local new_vars, i, nv;
   new_vars := [];
    for i in vars do
        nv := Copy(i);
        nv.t := TArray(nv.t, 2);
        nv.t.qualifiers := [opts.scratchModifier];
        Add(new_vars, nv);
    od;

    return new_vars;
end;

_double_buffer := function(cpu, dma, vars, opts)
    local new_cpu, new_dma, loops, new_loops, l, ld, st, s, w, c, i, s1, s2, ll, lst, new_vars, nv, svars, srec, v, vv, dvars, ass;

    svars := [];
    loops := Collect(dma, swp_loop);
    new_loops := [];
    new_dma := Copy(dma);
    new_cpu := cpu;

    for l in loops do
        i := l.var;

        ld := Collect(l, [loop, @(1), @(2), dma_load]);
        if Length(ld) = 0 then ld := Collect(l, dma_load)[1]; else ld := ld[1]; fi;
        st := Collect(l, [loop, @(1), @(2), dma_store]);
        if Length(st) = 0 then st := Collect(l, dma_store)[1]; else st := st[1]; fi;
        s := Collect(l, dma_signal)[1];
        w := Collect(l, cpu_wait)[1];
        s1 := Filtered(ld.free(), i->IsBound(i.t.qualifiers) and i.t.qualifiers=[opts.scratchModifier])[1];
        s2 := Filtered(st.free(), i->IsBound(i.t.qualifiers) and i.t.qualifiers=[opts.scratchModifier])[1];
        Add(svars, s1);
        Add(svars, s2);

        ll := Copy(l);
        lst := SubstVars(Copy(st), rec((s2.id) := nth(s2, bin_and(i, V(2)))));
        lst := SubstVars(lst, rec((i.id) := i-V(1)));
        lst := Collect(lst, dma_store)[1];

        ll.range := [Minimum(l.range)+1..Maximum(l.range)];
        ll := SubstTopDown(ll, @(1, dma_load), e->dma_load(nth(@(1).val.loc, bin_and(i, V(2))), @(1).val.exp, @(1).val.size));
        ll := SubstTopDown(ll, @(1, dma_store), e->lst);
        ll := SubstTopDown(ll, @(1, cpu_wait), e->cpu_wait(@(1).val.args[1]-1));

        c := chain([
            SubstVars(chain([Copy(ld), Copy(s)]), rec((i.id) := V(0), (s1.id) := nth(s1, 0))),
            ll,
            SubstVars(chain([Copy(w), Copy(st)]), rec((i.id) := i.range-V(1), (s2.id) := nth(s2, 1))),
        ]);
        c := RulesStrengthReduce(c);
        new_dma := SubstBottomUp(new_dma, @(1, swp_loop, e->e.var = i), e->c);
    od;

    svars := Set(svars);
    loops := Collect(cpu, swp_loop);
    new_cpu := Copy(cpu);
    new_loops := [];
    for l in loops do
        i := l.var;
        srec := rec();
        dvars := [];
        ass := [];
        for v in svars do
            vv := var.fresh_t("R", TPtr(v.t.t));
            Add(dvars, vv);
            Add(ass, assign(vv, nth(v, bin_and(i, V(2)))));
            srec.(v.id) := vv;
        od;
        c := SubstVars(l, srec);
        c := swp_loop(c.var, c.range, decl(dvars, chain(Concat(ass, [c.cmd]))));
        new_cpu := SubstBottomUp(new_cpu, @(1, swp_loop, e->e.var = i), e->c);
    od;

    new_vars := _scratch_variables(vars,opts);

    return [new_cpu, new_dma, new_vars];
end;

Class(ScratchMainCodegen, rec(
    genMain := (self, opts, sub, dmafunc, cpufunc, xy, params) >>
        func(TVoid, sub, Concatenation(xy, params), par_exec(
            ApplyFunc(call, Flat(Concatenation([dmafunc], xy))),
            call(cpufunc)
        ))
));

Class(ScratchCodegen, DefaultCodegen, rec(
    CPUCodegen := CPUCodegen,
    DMACodegen := DMACodegen,
    MainCodegen := ScratchMainCodegen,

    Formula := meth(self, o, y, x, opts)
        local tag, main, initfunc, cpufunc, dmafunc, prog, params, sub, initsub, memvars, scratchvars, v, io, loadvar, storevar, v, substrec, svars, initcode, datas, dvars, dv;

        datas := Collect(o, FDataOfs);
        params := Set(Concatenation(Collect(o, param), Filtered(Collect(o, var), IsParallelLoopIndex)));
        sub := Cond(IsBound(opts.subName), opts.subName, "transform");

        initsub := Cond(IsBound(opts.subName), Concat("init_", opts.subName), "init");
        if IsBound(opts.generateInitFunc) and opts.generateInitFunc then
            initcode := chain(List(datas, x -> SReduce(x.var.init, opts)));
            initfunc := func(TVoid, initsub, params :: Set(Collect(initcode, param)), initcode);
        else
            initfunc := func(TVoid, initsub, params, chain());
        fi;

        dvars := List(datas, x->x.var);
        for dv in dvars do
            dv.t.qualifiers := [opts.romModifier];
        od;

        self.DMACodegen.loadbuffers := Set([]);
        self.DMACodegen.storebuffers := Set([]);
        self.DMACodegen.membuffers := Set([]);
        dmafunc := self.DMACodegen.Formula(o, y, x, opts);
        self.DMACodegen.membuffers := Set(self.DMACodegen.membuffers);
        SubtractSet(self.DMACodegen.membuffers, Set([x, y]));
        for v in self.DMACodegen.membuffers do v.t.qualifiers := [opts.memModifier]; od;

        self.CPUCodegen.loadbuffers := Set([]);
        self.CPUCodegen.storebuffers := Set([]);
        cpufunc := self.CPUCodegen.Formula(o, y, x, opts);

        memvars := Set(Concat(
            Filtered(Flat(List(Collect(cpufunc, decl), i->i.vars)), j->IsBound(j.t.qualifiers) and opts.memModifier in j.t.qualifiers),
            self.DMACodegen.membuffers));
        scratchvars := Set(Filtered(Flat(List(Collect(cpufunc, decl), i->i.vars)), j->IsBound(j.t.qualifiers) and opts.scratchModifier in j.t.qualifiers));

        tag := opts.tags[1];

        loadvar := var.fresh_t("S", TArray(opts.XType.t, Cond(tag.isRegCx, 2 * tag.size, tag.size)));
        loadvar.t.qualifiers := [opts.scratchModifier];
        storevar := var.fresh_t("S", TArray(opts.XType.t, Cond(tag.isRegCx, 2 * tag.size, tag.size)));
        storevar.t.qualifiers := [opts.scratchModifier];

        substrec := rec();
        for v in self.CPUCodegen.loadbuffers do substrec.(v.id) := loadvar; od;
        for v in self.CPUCodegen.storebuffers do substrec.(v.id) := storevar; od;
        for v in self.DMACodegen.loadbuffers do substrec.(v.id) := loadvar; od;
        for v in self.DMACodegen.storebuffers do substrec.(v.id) := storevar; od;
        cpufunc :=  SubstVars(cpufunc, substrec);
        dmafunc :=  SubstVars(dmafunc, substrec);

        cpufunc := SubstTopDown(cpufunc, @(1, decl), e->decl(Filtered(@(1).val.vars, i->not i in Concat(memvars, scratchvars, [loadvar, storevar])), @(1).val.cmd));
        dmafunc := SubstTopDown(dmafunc, @(1, decl), e->decl(Filtered(@(1).val.vars, i->not i in Concat(memvars, scratchvars, [loadvar, storevar])), @(1).val.cmd));

        [cpufunc, dmafunc, svars] := When(opts.swp,_double_buffer(cpufunc, dmafunc, [loadvar, storevar], opts),[cpufunc,dmafunc,[loadvar,storevar]]);

        x.t.qualifiers := Set(x.t.qualifiers);
        y.t.qualifiers := Set(y.t.qualifiers);

        main := self.MainCodegen.genMain(opts, sub, dmafunc, cpufunc, [x,y], params);
        prog := program(
            decl(Concat(Set(memvars), svars, dvars), chain(
                initfunc,
                cpufunc,
		        dmafunc,
                main
        )));

        prog.dimensions := o.dims();
        return prog;
    end
));

ext: gi


text:
#############################################################################
##
##  ExternalPolymakePolytope.gi  JConvex package
##                               Martin Bies
##
##  Copyright 2021               University of Pennsylvania
##
##  A Gap package to do convex geometry by Polymake, Cdd and Normaliz
##
##  Chapter Polytopes in Polymake
##
#############################################################################

##############################################################################################
##
##  Section GAP category of PolymakePolytopes
##
##############################################################################################

DeclareRepresentation( "IsPolymakePolytopeRep", IsPolymakePolytope and IsAttributeStoringRep, [ ] );

BindGlobal( "TheFamilyOfPolymakePolytopes", NewFamily( "TheFamilyOfPolymakePolytopes" ) );

BindGlobal( "TheTypeOfPolymakePolytope", NewType( TheFamilyOfPolymakePolytopes, IsPolymakePolytopeRep ) );

BindGlobal( "MakePolymakePolytopeVRep", function( vertices, lineality )
    local poly, kwargs, new_vertices, new_lin, v;
    poly := Objectify( TheTypeOfPolymakePolytope,
            rec( vertices := Immutable( vertices ),
                 lineality := Immutable( lineality ),
                 number_type := "rational",
                 rep_type := "V-rep" ) );

    # add 1s at the beginning, since Polymake always considers polytopes as intersections with the hyperplane x0 = 1
    new_vertices := [];
    for v in vertices do
        v := ShallowCopy( v );
        Add( v, 1, 1 );
        Add( new_vertices, v );
    od;
    new_lin := [];
    for v in lineality do
        v := ShallowCopy( v );
        Add( v, 1, 1 );
        Add( new_lin, v );
    od;

    kwargs := rec();

    # TODO: verify that kwargs is set correctly, also if one or both lists are empty

    # check for degenerate case
    if Length( new_vertices ) > 0 then
        kwargs.POINTS := JuliaMatrixInt( new_vertices );

        # if we also have lineality, add them
        if Length( new_lin ) > 0 then
            kwargs.INPUT_LINEALITY := JuliaMatrixInt( new_lin );
        fi;
    elif Length( new_lin ) > 0 then
        kwargs.POINTS := JuliaMatrixInt( new_lin );
        kwargs.INPUT_LINEALITY := kwargs.POINTS;
    fi;

    poly!.pmobj := CallJuliaFunctionWithKeywordArguments( _Polymake_jl.polytope.Polytope, [], kwargs );
    return poly;
end);

BindGlobal( "MakePolymakePolytopeHRep", function( inequalities, equalities )
    local poly, kwargs;
    poly := Objectify( TheTypeOfPolymakePolytope,
            rec( inequalities := Immutable( inequalities ),
                 equalities := Immutable( equalities ),
                 number_type := "rational",
                 rep_type := "H-rep" ) );

    kwargs := rec();

    # TODO: verify that kwargs is set correctly, also if one or both lists are empty

    # check for degenerate case
    if Length( inequalities ) > 0 then
        kwargs.INEQUALITIES := JuliaMatrixInt( inequalities );

        # check if we also need equalities
        if Length( equalities ) > 0 then
            kwargs.EQUATIONS := JuliaMatrixInt( equalities );
        fi;
    elif Length( equalities ) > 0 then
        kwargs.INEQUALITIES := JuliaMatrixInt( equalities );
        kwargs.EQUATIONS := kwargs.INEQUALITIES;
    fi;

    poly!.pmobj := CallJuliaFunctionWithKeywordArguments( _Polymake_jl.polytope.Polytope, [], kwargs );
    return poly;
end);

##############################################################################################
##
##  Constructors for PolymakePolytopes
##
##############################################################################################

InstallGlobalFunction( Polymake_PolytopeByGenerators,
  function( arg )
    local poly, i, matrix, temp, dim;

    if Length( arg )= 0 or ForAll( arg, IsEmpty ) then

        Error( "Wrong input: Please provide some input!" );

    elif Length( arg ) = 1 and IsList( arg[1] ) then

        return Polymake_PolytopeByGenerators( arg[ 1 ], [ ] );

    elif Length( arg ) = 2 and IsList( arg[ 1 ] ) and IsList( arg[ 2 ] ) then

        if ( not IsEmpty( arg[ 1 ] ) ) and not ( IsMatrix( arg[ 1 ] ) ) then
            Error( "Wrong input: The first argument should be a Gap matrix!" );
        fi;

        if ( not IsEmpty( arg[ 2 ] ) ) and not ( IsMatrix( arg[ 2 ] ) ) then
            Error( "Wrong input: The second argument should be a Gap matrix!" );
        fi;

        poly := MakePolymakePolytopeVRep( arg[ 1 ], arg[ 2 ] );
        return Polymake_CanonicalPolytopeByGenerators( poly );

    fi;

end );

InstallGlobalFunction( Polymake_PolytopeFromInequalities,
  function( arg )
    local poly, i, temp, matrix, dim;

    if Length( arg ) = 0 or ForAll( arg, IsEmpty ) then

        Error( "Wrong input: Please provide some input!" );

    elif Length( arg ) = 1 and IsList( arg[ 1 ] ) then

        return Polymake_PolytopeFromInequalities( arg[ 1 ], [ ] );

    elif Length( arg ) = 2 and IsList( arg[ 1 ] ) and IsList( arg[ 2 ] ) then

        if ( not IsEmpty( arg[ 1 ] ) ) and not ( IsMatrix( arg[ 1 ] ) ) then
            Error( "Wrong input: The first argument should be a Gap matrix!" );
        fi;

        if ( not IsEmpty( arg[ 2 ] ) ) and not ( IsMatrix( arg[ 2 ] ) ) then
            Error( "Wrong input: The second argument should be a Gap matrix!" );
        fi;

        poly := MakePolymakePolytopeHRep( arg[ 1 ], arg[ 2 ] );
        return Polymake_CanonicalPolytopeFromInequalities( poly );

    fi;

end );

##############################################################################################
##
##  Canonicalize polytopes
##
##############################################################################################

InstallMethod( Polymake_CanonicalPolytopeByGenerators,
               [ IsPolymakePolytope ],
  function( poly )
    local vertices, v_copy, scaled_vertices, i, scale, lineality, scaled_lineality;

    if poly!.rep_type = "H-rep" then

        return fail;

    else

        # compute vertices
        vertices := PolymakeMatrixToGAP( poly!.pmobj.VERTICES );

        # sometimes, Polymake returns rational vertices - we turn them into integral vectors
        # also, Polymake requires x0 = 1 in affine coordinates - we remove this 1
        scaled_vertices := [];
        for i in [ 1 .. Length( vertices ) ] do
            v_copy := ShallowCopy( vertices[ i ] );
            Remove( v_copy, 1 );
            Add( scaled_vertices, v_copy );
        od;

        # extract lineality
        lineality := PolymakeMatrixToGAP( poly!.pmobj.LINEALITY_SPACE );

        # sometimes, Polymake returns rational lineality - we turn them into integral vectors
        scaled_lineality := [];
        for i in [ 1 .. Length( lineality ) ] do
            v_copy := ShallowCopy( lineality[ i ] );
            Remove( v_copy, 1 );
            Add( scaled_lineality, v_copy );
        od;

        # construct the new poly
        return MakePolymakePolytopeVRep( scaled_vertices, scaled_lineality );

    fi;

end );

InstallMethod( Polymake_CanonicalPolytopeFromInequalities,
               [ IsPolymakePolytope ],
  function( poly )
    local ineqs, eqs;

    if poly!.rep_type = "V-rep" then

        return fail;

    else

        # compute facets
        ineqs := PolymakeMatrixToGAP( poly!.pmobj.FACETS );

        # compute affine hull
        eqs := PolymakeMatrixToGAP( poly!.pmobj.AFFINE_HULL );

        # construct the new poly
        return MakePolymakePolytopeHRep( ineqs, eqs );

    fi;

end );

##############################################################################################
##
##  Conversion of polytopes
##
##############################################################################################

InstallMethod( Polymake_V_Rep,
               [ IsPolymakePolytope ],
  function( poly )
    local vertices, v_copy, scaled_vertices, i, lineality, scaled_lineality;

    if poly!.rep_type = "V-rep" then

        return poly;

    else

        # compute vertices
        vertices := PolymakeMatrixToGAP( poly!.pmobj.VERTICES );

        # sometimes, Polymake returns rational vertices - we turn them into integral vectors
        scaled_vertices := [];
        for i in [ 1 .. Length( vertices ) ] do
            v_copy := ShallowCopy( vertices[ i ] );
            Remove( v_copy, 1 );
            Add( scaled_vertices, v_copy );
        od;

        # compute lineality
        lineality := PolymakeMatrixToGAP( poly!.pmobj.LINEALITY_SPACE );

        # sometimes, Polymake returns rational lineality - we turn them into integral vectors
        scaled_lineality := [];
        for i in [ 1 .. Length( lineality ) ] do
            v_copy := ShallowCopy( lineality[ i ] );
            Remove( v_copy, 1 );
            Add( scaled_lineality, v_copy );
        od;

        # construct the new poly
        return MakePolymakePolytopeVRep( scaled_vertices, scaled_lineality );

    fi;

end );

InstallMethod( Polymake_H_Rep,
               [ IsPolymakePolytope ],
  function( poly )
    local ineqs, eqs;

    if poly!.rep_type = "H-rep" then

        return poly;

    else

        if poly!.rep_type = "V-rep" and poly!.vertices = [] then
            return Polymake_PolytopeFromInequalities( [ [ 0, 1 ], [ -1, -1 ] ] );
        fi;

        # compute inequalities
        ineqs := PolymakeMatrixToGAP( poly!.pmobj.FACETS );

        # compute equalities
        eqs := PolymakeMatrixToGAP( poly!.pmobj.AFFINE_HULL );

        # construct the new poly
        return MakePolymakePolytopeHRep( ineqs, eqs );

    fi;

end );

##############################################################################################
##
##  Attributes of PolymakeCones
##
##############################################################################################

InstallMethod( Polymake_AmbientSpaceDimension,
              "finding the dimension of the ambient space of the poly",
              [ IsPolymakePolytope ],
  function( poly )

    return Length( Polymake_V_Rep( poly )!.vertices[0] );

end );

InstallMethod( Polymake_Dimension,
              " returns the dimension of the poly",
            [ IsPolymakePolytope ],
  function( poly )

    if Polymake_IsEmpty( poly ) then
        return -1;
    fi;

    return poly!.pmobj.CONE_DIM - 1;

end );

InstallMethod( Polymake_Vertices,
              " return the list of generating vertices",
              [ IsPolymakePolytope ],
  function( poly )

    return Set( Polymake_V_Rep( poly )!.vertices );

end );

InstallMethod( Polymake_Linealities,
              " return the list of linealities",
              [ IsPolymakePolytope ],
  function( poly )

    return Set( Polymake_V_Rep( poly )!.linealities );

end );

InstallMethod( Polymake_Equalities,
              " return the list of equalities of a poly",
              [ IsPolymakePolytope ],
  function( poly )

    return Set( ( Polymake_H_Rep( poly ) )!.equalities );

end );

InstallMethod( Polymake_Inequalities,
              " return the list of inequalities of a poly",
              [ IsPolymakePolytope ],
  function( poly )

    return Set( ( Polymake_H_Rep( poly ) )!.inequalities );

end );

InstallMethod( Polymake_LatticePoints,
              " return the list of the lattice points of poly",
              [ IsPolymakePolytope ],
  function( poly )
    local point_list, lattice_points, i, copy;

    point_list := JuliaToGAP( IsList, JuliaMatrixInt( poly!.pmobj.LATTICE_POINTS_GENERATORS[1] ), true );

    # the point list is in affine coordinate, that is we have a 1 at position one, which should be removed
    lattice_points := [];
    for i in [ 1 .. Length( point_list ) ] do
        copy := ShallowCopy( point_list[ i ] );
        Remove( copy, 1 );
        Add( lattice_points, copy );
    od;

    # return result
    return lattice_points;

end );

InstallMethod( Polymake_Intersection,
               [ IsPolymakePolytope, IsPolymakePolytope ],
  function( poly1, poly2 )
    local poly1_h, poly2_h, new_ineqs, new_equ;

    # compute H-reps
    poly1_h := Polymake_H_Rep( poly1 );
    poly2_h := Polymake_H_Rep( poly2 );

    # add the inequalities and equalities
    new_ineqs := Concatenation( poly1_h!.inequalities, poly2_h!.inequalities );
    new_equ := Concatenation( poly1_h!.equalities, poly2_h!.equalities );

    # return result
    return Polymake_PolytopeFromInequalities( new_ineqs, new_equ );

end );

##############################################################################################
##
##  Properties of PolymakeCones
##
##############################################################################################

InstallMethod( Polymake_IsEmpty,
               "finding if the poly empty is or not",
               [ IsPolymakePolytope ],
  function( poly )

    return Length( Polymake_V_Rep( poly )!.vertices ) = 0;

end );

InstallMethod( Polymake_IsPointed,
               "finding if the poly is pointed or not",
               [ IsPolymakePolytope ],
  function( poly )
    return poly!.pmobj.POINTED;

end );

InstallMethod( Polymake_IsBounded,
              " returns if the polytope is bounded or not",
              [ IsPolymakePolytope ],
  function( poly )
    return poly!.pmobj.BOUNDED;

end );


text:
ZigZag := function(n)
  local a, i, j, k;
  a := NullMat(n, n);
  i := 1;
  j := 1;
  for k in [0 .. n*n - 1] do
    a[i][j] := k;
    if (i + j) mod 2 = 0 then
      if j < n then
        j := j + 1;
      else
        i := i + 2;
      fi;
      if i > 1 then
        i := i - 1;
      fi;
    else
      if i < n then
        i := i + 1;
      else
        j := j + 2;
      fi;
      if j > 1 then
        j := j - 1;
      fi;
    fi;
  od;
  return a;
end;

PrintArray(ZigZag(5));
# [ [   0,   1,   5,   6,  14 ],
#   [   2,   4,   7,  13,  15 ],
#   [   3,   8,  12,  16,  21 ],
#   [   9,  11,  17,  20,  22 ],
#   [  10,  18,  19,  23,  24 ] ]


text:
# A hyperconic in a Desarguesian projective plane of order 2^s.
BindGlobal("HyperconicProjectivePlane",
    s -> rec(points := Set(List(Union([[0,1,0], [1,0,0]]*Z(2^s)^0,
                                    List(GF(2^s), t -> [t, t^2, 1]*Z(2^s)^0)),
                                x -> Subspace(GF(2^s)^3, [x], "basis"))),
            group := HyperconicAutomorphismGroup(2^s), d := 2, q := 2^s)
);

# A classical conic in a Desarguesian projective plane of order q.
BindGlobal("ClassicalConicProjectivePlane",
    q -> rec(points := Set(List(Union([[0,1,0]]*Z(q)^0, List(GF(q),
            t -> [t, t^2, 1]*Z(q)^0)), x -> Subspace(GF(q)^3, [x], "basis"))),
            group := HyperconicAutomorphismGroup(q), d := 2, q := q)
);

# A nonclassical conic in a Desarguesian projective plane of order q.
BindGlobal("NonclassicalConicProjectivePlane",
    q -> rec(points := Set(List(Union([[1,0,0]]*Z(q)^0, List(GF(q),
            t -> [t, t^2, 1]*Z(q)^0)), x -> Subspace(GF(q)^3, [x], "basis"))),
            group := Group([DiagonalMat([Z(q), Z(q)^2, Z(q)^0]),
                            Z(q)*IdentityMat(3, GF(q))]),
            d := 2, q := q)
);

# An elliptic quadric in a d-dimensional projective space of order q.
BindGlobal("EllipticQuadricProjectiveSpace", function(q)
    local G, V;
    G := GO(-1, 4, q);
    V := GF(q)^4;
    return rec(points := IsotropicSpacesQuadraticForm(
                    InvariantQuadraticForm(G).matrix)(Subspaces(V, 1)),
                group := G, d := 3, q := q);
end);

max_stars_repo_name: jaanos/gap-graphs

max_issues_repo_name: jaanos/gap-graphs

max_forks_repo_name: jaanos/gap-graphs


text:
##  Copyright (c) 2018-2021, Carnegie Mellon University
##  See LICENSE for details

NewRulesFor(DFT, rec(

    DFT_PD_loop := rec(
        forTransposition := false,
        minSize := 7,
        maxSize := 13,

        TabPerm := r -> let(i:=Ind(r.domain()), f := FData(r.tolist()).at(i), Lambda(i, f).setRange(r.range())),

        applicable     := (self, nt) >> nt.params[1] > 2 and nt.params[1] in [self.minSize..self.maxSize] and IsPrime(nt.params[1]) and not nt.hasTags(),

        apply := (self, nt, C, cnt) >> let(
            N := nt.params[1],
            k := nt.params[2],
            root := PrimitiveRootMod(N),
            M:=(N-1)/2,
            i := Ind(M),
            j := Ind(M),
            k1 := Ind(M+1),
            m := MatSPL(DFT_PD.core(N, k, root, false) * DFT_PD.A(N)),
            m1 := Map(m{[2..(N+1)/2]}, r -> r{[2..(N+1)/2]}),
            m2 := Map(m{[(N+3)/2..N]}, r -> r{[(N+3)/2..N]}),
            d := Flat(m1)::Map(Flat(m2), c->im(ComplexAny(c)).v),
            fd := FData(d),

            s := Scat(self.TabPerm(RR(N, 1, root))),
            g := Gath(self.TabPerm(RR(N, 1, 1/root mod N))),

            gg := Gath(fCompose(fAdd(N, N-1, 1), fTensor(fId(2), fBase(M, k1-1)))),
            f2 := F(2),
            u := Ind(2),
            uf := Lambda(u, fd.at(u*M*M + M*i+(k1-1))),
            bb := DirectSum(Blk1(uf.at(0)), Scale(ImaginaryUnit(), Blk1(uf.at(1)))),

            krn0 := Mat([[1],[0]]) * Gath(fAdd(N, 1, 0)),
            krnn := bb * f2 * gg,
            krn := ScatAcc(fId(2)) * f2 * COND(eq(k1,0), krn0, krnn),

            q1 := L(2*M, 2) * IterVStack(i, ISumAcc(k1, krn)),
            #q2 := RowVec(fConst(13,V(1.0))),
            _i := Ind(N),
            q2 := ISumAcc(_i, ScatAcc(fAdd(1,1,0)) * Blk([[1]]) * COND(eq(_i, 0), Gath(fAdd(13,1,0)), Gath(fAdd(13,1,_i)))),
            q3 := VStack(q2, q1 * g),
            qq := s * q3,
            qq
        )
    )
));

RulesFor(PRDFT, rec(
   PRDFT_PD_loop := rec(
	forTransposition := false,
    minSize := 7,
	maxSize          := 13,
	isApplicable     := (self, P) >> P[1] > 2 and P[1] in [self.minSize..self.maxSize] and IsPrime(P[1]),

	rule := (self,P,C) >> let(N:=P[1], n:=N-1, k:=P[2], root:=PrimitiveRootMod(N),
        M:=n/2,
        i := Ind(M),
        j := Ind(M),
        k1 := Ind(M+1),
        m := MatSPL(DFT_PD.core(N, k, root, false) * DFT_PD.A(N)),
        m1 := Map(m{[2..(N+1)/2]}, r -> r{[2..(N+1)/2]}),
        m2 := Map(m{[(N+3)/2..N]}, r -> r{[(N+3)/2..N]}),
        d := Flat(m1)::Map(Flat(m2), c->im(ComplexAny(c)).v),
        fd := FData(d),
        #lfd1 := Lambda(j, fd.at(M*i+j)),
        #lfd1 := Lambda(k1, cond(eq(k1, V(0)), V(1.0), fd.at(M*i+(k1-V(1))))),
        #lfd2 := Lambda(j, fd.at(M*M+M*i+j)),

        gf := DFT_PD_loop.TabPerm(RR(N, 1, root)),
        g := Gath(gf),

        #kk1 := DirectSum(RowVec(lfd1), RowVec(lfd2)) * DirectSum(I(1), Tensor(F(2), I(M)) * OS(n, -1)),
        #q1 :=  IterVStack(i, BB(kk1)),

        #lbd := OS(12, -1).lambda(),
        _j := Ind(n),
        lbd := Lambda(_j, imod(V(n)-_j, V(n))),
        gg := Gath(fCompose(fAdd(N, N-1, 1), fCompose(lbd, fTensor(fId(2), fBase(M, k1-1))))),

        f2 := F(2),
        u := Ind(2),
        uf := Lambda(u, fd.at(u*M*M + M*i+(k1-1))),
        bb := DirectSum(Blk1(uf.at(0)), Blk1(uf.at(1))),

        krn0 := Mat([[1],[0]]) * Gath(fAdd(N, 1, 0)),
        krnn := bb * f2 * gg,
        krn := Grp(ScatAcc(fId(2))) * COND(eq(k1,0), krn0, krnn),

        q1 := IterVStack(i, ISumAcc(k1, krn)),
        q2a := RowVec(fConst(13,V(1.0))),
#        _i := Ind(N),
#        q2a := ISumAcc(_i, ScatAcc(fAdd(1,1,0)) * Blk([[1]]) * COND(eq(_i, 0), Gath(fAdd(13,1,0)), Gath(fAdd(13,1,_i)))),
        q2b := RowVec(fConst(13,V(0.0))),
        q3 := VStack(BB(VStack(q2a, q2b)), q1),
        sf := DFT_PD_loop.TabPerm(Refl((N+1)/2, N, (N+1)/2, RR(N,1,root))),
        s := Scat(fTensor(sf, fId(2))),

        u1 := Ind(N+1),
        df := Lambda(u1, cond(logic_and(eq(V(1), bin_and(u1, 1)), geq(gf.at(idiv(u1, 2)), (N+1)/2)), V(-1), V(1))),
        dfl := df.tolist(),
        m1s := Filtered([0..N], _i->dfl[_i+1] = V(-1)),
        m1eq := List(m1s, _i->eq(u1, V(_i))),
        cnd := ApplyFunc(logic_or, m1eq),
        lbdcnd := Lambda(u1, cond(cnd, V(-1), V(1))),
        diag := Diag(lbdcnd),
        sct := s * diag,
        qq := sct * q3 * g,
        qq
    )
)));

RulesFor(IPRDFT, rec(
   IPRDFT_PD_loop := rec(
	forTransposition := false,
    minSize := 7,
	maxSize          := 13,
	isApplicable     := (self, P) >> P[1] > 2 and P[1] in [self.minSize..self.maxSize] and IsPrime(P[1]),

	rule := (self,P,C) >> let(N:=P[1], n:=N-1, k:=P[2], root:=PrimitiveRootMod(N), M:=n/2,
        # loop variables
        ii := Ind(M),
        k1 := Ind(M+1),
        u := Ind(2),
        # scatter
        fstr := fDirsum(fId(1), L(n, n/2)),
        fos := fDirsum(fId(2), J(N-2)),
        frr := RR(N, 1, root),
        fsct := fCompose(frr, fos, fstr),
        fl := fsct.tolist(){[2..N]},
        flst := FList(N-1, fl),
        fst :=  DFT_PD_loop.TabPerm(flst),
        i := Ind(N),
        flbd := Lambda(i, cond(eq(i, V(0)), V(0), fst.at(i-V(1)))).setRange(N),
        sct := Scat(flbd),
        # diagonal sign flip and block matrix
        rr := RealRR_Out(N, root).transpose(),
        dd := rr.child(1),
        mc := TransposedSPL(DFT_PD.core(N, k, root, false) * DFT_PD.A(N)) * DirectSum(Mat([[1,0]]), 2*I(n)),
        dd2 := DirectSum(I(2), L(n, 2)) * dd * DirectSum(I(2), L(n, n/2)),
        m := MatSPL(mc * dd2),
        m1 := Map(m{[2..(N+1)/2]}, r -> r{[3..(N+1)/2+1]}),
        m2 := Map(m{[(N+3)/2..N]}, r -> Map(r{[(N+3)/2+1..N+1]}, c->-im(ComplexAny(c)).v)),
        d := Flat(m1)::Flat(m2),
        fd := FData(d),
        # gather
        gth := rr.child(2),
        gf := gth.func,
        gfl := gf.tolist(){2*[1..(N-1)/2]+1},
        gst := FList((N-1)/2, gfl),
        gftb := DFT_PD_loop.TabPerm(gst).setRange(N+1),
        glbd := Lambda(u, gftb.at(k1-1) + u).setRange(N+1),
        gath := Gath(glbd),
        # kernel
        uf := Lambda(u, fd.at(u*M*M + M*ii+(k1-1))),
        bb := DirectSum(Blk1(uf.at(0)), Blk1(uf.at(1))),
        f2 := F(2),
        krn0 := f2 * Mat([[1],[0]]) * Gath(fAdd(N+1, 1, 0)),
        krnn := f2 * bb * gath,
        krn := Grp(ScatAcc(fId(2))) * COND(eq(k1,0), krn0, krnn),
        # stack
        q1 := IterVStack(ii, ISumAcc(k1, krn)),
        r1 := RowVec(diagDirsum(fConst(1,V(1.0)), fConst(n/2,V(2.0)))) * Gath(fTensor(fId((N+1)/2), fBase(2, 0))),
        qq := VStack(r1, q1),
        qq1 := sct * qq,
        qq1
    ))
));

max_stars_repo_name: franzfranchetti/spiral-package-fftx

max_issues_repo_name: franzfranchetti/spiral-package-fftx

max_forks_repo_name: franzfranchetti/spiral-package-fftx


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

#F PrintActiveRules(<non-terminal>)
#F PrintActiveRules(<non-terminal.name>)
#F PrintActiveRules([<rulelist>])
#F
#F Prints out all applicable rules for for the transform <non-terminal>
#F together with their switches.
#F Accepts a nonterminal spl, its symbol, or a rule list as an argument.
#F

PrintActiveRules := function (nonterm)

local rules, switches, i;

# Check if rules are already found
if IsList(nonterm) and not IsString(nonterm) then
  if ForAll(nonterm, IsRule) then rules:=nonterm;
  else Error("List must be a list of rules");
  fi;
else rules := AllRules(nonterm);
fi;

switches :=
  List([1..Length(rules)],
    function(i)
      if rules[i].switch = false then return "OFF";
      else return "ON";
      fi;
    end
  );
Print("\n");
Print("No. \t\b\bSwitch \tRule\n");
Print("--------------------------------------------\n");
for i in [1..Length(rules)] do
  Print(i,".\t",switches[i],"\t",rules[i],"\n");
od;
Print("\n");
end;

#F SwitchRulesOn( <non-terminal>, <ind> )
#F SwitchRulesOn( <non-terminal.name>, <ind> )
#F SwitchRulesOn( [<rulelist>], <ind> )
#F
#F Switches on the rules with indices <ind> for a given <non-terminal>
#F   <non-terminal> can be either a nonterminal spl or its symbol.
#F   <ind> is either a single index or a list of indices
#F
#F Hint: Use PrintActiveRules(<non-terminal>) to find out more about
#F       rules and their numbering.
#F

SwitchRulesOn := function ( nonterm , index )
local i, L;

# Check if rules are already found
if IsList(nonterm) and not IsString(nonterm) then
  if ForAll(nonterm, IsRule) then L := nonterm;
  else Error("List must be a list of rules");
  fi;
else L := AllRules(nonterm);
fi;

if not IsList(index) then index := [index];fi;
if ForAny(index, k-> not (IsInt(k) and k <= Length(L) and k > 0)) then
     Error("Index to the rule list must be an integer in the proper range");
fi;
for i in index do
  L[i].switch:=true;
od;

PrintActiveRules(L);
end;

Declare(SwitchRulesQuiet, SwitchRulesByNameQuiet);

#F SwitchRules( <non-terminal>, <ind> )
#F SwitchRules( <non-terminal.name>, <ind> )
#F SwitchRules( [<rulelist>], <ind> )
#F
#F Switches on the rules with indices <ind> for a given <non-terminal>
#F   and off all the other rules for the same nonterm
#F   <non-terminal> can be either a nonterminal spl or its symbol.
#F   <ind> is either a single index or a list of indices
#F
#F Hint: Use PrintActiveRules(<non-terminal>) to find out more about
#F       rules and their numbering.
#F
SwitchRules := function ( nonterm , index )
    SwitchRulesQuiet(nonterm, index);
    PrintActiveRules(nonterm);
end;

SwitchRulesQuiet := function ( nonterm , index )
local i, L;

# Check if rules are already found
if IsList(nonterm) and not IsString(nonterm) then
  if ForAll(nonterm, IsRule) then L := nonterm;
  else Error("List must be a list of rules");
  fi;
else L := AllRules(nonterm);
fi;

if not IsList(index) then index := [index];fi;
if ForAny(index, k-> not (IsInt(k) and k <= Length(L) and k > 0)) then
     Error("Index to the rule list must be an integer in the proper range");
fi;

for i in [1..Length(L)] do
  L[i].switch:=false;
od;

for i in index do
  L[i].switch:=true;
od;
end;

SwitchRulesByName := function ( nonterm , rules )
    SwitchRulesByNameQuiet(nonterm, rules);
    PrintActiveRules(nonterm);
end;

SwitchRulesByNameQuiet := function ( nonterm , rules )
    local r, L;
    if not IsNonTerminal(nonterm) then Error("<nonterm> must be a non-terminal");
    elif not (IsList(rules) and ForAll(rules, IsRule)) then Error("<rules> must be a list of rules");
    fi;
    SwitchRulesQuiet(nonterm, []);
    for r in rules do
        r.switch:=true;
    od;
end;

#F SwitchRulesOff( <non-terminal>, <ind> )
#F SwitchRulesOff( <non-terminal.name>, <ind> )
#F SwitchRulesOff( [<rulelist>], <ind> )
#F
#F Switches off the rules with indices <ind> for a given <non-terminal>
#F   <non-terminal> can be either a nonterminal spl or its symbol.
#F   <ind> is either a single index or a list of indices
#F
#F Hint: Use PrintActiveRules(<non-terminal>) to find out more about
#F       rules and their numbering.
#F

SwitchRulesOff := function ( nonterm , index )
local i, L;

# Check if rules are already found
if IsList(nonterm) and not IsString(nonterm) then
  if ForAll(nonterm, IsRule) then L := nonterm;
  else Error("List must be a list of rules");
  fi;
else L := AllRules(nonterm);
fi;

if not IsList(index) then index := [index]; fi;
if ForAny(index, k-> not (IsInt(k) and k <= Length(L) and k > 0)) then
     Error("Index to the rule list must be an integer in the proper range");
fi;
for i in index do
  L[i].switch:=false;
od;

PrintActiveRules(L);
end;

# VIENNA ADDED:
#F SwitchRulesName( [<rulelist>], true|false )
#F
#F Switches the rules on and off by names
#F   <rulelist> is a list of names of rules
#F
SwitchRulesName := function(l,s)
local i;
for i in l do i.switch:=s; od;
end;


text:
# Copyright (c) 2018-2020, Carnegie Mellon University
# See LICENSE for details

# Radix N kernel for FFTE

buildScalarKernel_a := function(n, kk, conf, opts)
    local name, suffix, filesuffix, j, k, l, m, tmp1, tmp2, tmp3, i1, i2, gf, sf, gath, scat,
        gc, sc, dft, rts, rt, opcnts, dfts, dftc, i3, twt, twf, twl, tws, twc,
        cl, cc, c, cx, lp, mp;

    name := DFTnameStr(n, kk);
    suffix := "a";
    filesuffix := "c_";
#    suffix := "_j";
    Print("\n\n// == Generating \"", name, "\" ==================================================\n\n");

    # variables
    #t := var.fresh_t("t", TInt);
    j := var.fresh_t("j", TInt);
#    k := var.fresh_t("k", TInt);
    k := V(0);
    l := var.fresh_t("l", TInt);
 #   m := var.fresh_t("m", TInt);
   m := V(1);
    lp := var.fresh_t("lp", TPtr(TInt));
    mp := var.fresh_t("mp", TPtr(TInt));

    # temp arrays
    tmp1 := var.fresh_t("R", TArray(TReal, 2*n));
    tmp2 := var.fresh_t("S", TArray(TReal, 2*n));
    tmp3 := var.fresh_t("T", TArray(TReal, 2*n));

    # gather and scatter
    i1 := var.fresh_t("i", n);
    i2 := var.fresh_t("i", n);
    gf := fTensor(Lambda(i1, k + j*m + i1*l*m), fId(2));
    sf := fTensor(Lambda(i2, k+n*j*m + i2*m), fId(2));
    gath := Gath(gf);
    scat := Scat(sf);

    # generate gather and scatter code
    gc := BlockUnroll(unroll_cmd(opts.codegen(gath, tmp1, X, opts)), opts);
    sc := BlockUnroll(unroll_cmd(opts.codegen(scat, Y, tmp3, opts)), opts);

    # generate DFT kernel
    dft := RC(DFT(n, kk));
    rts := AllRuleTrees(dft, opts);
    opcnts := List(rts, r -> [ Length(Collect(CodeRuleTree(r, opts), @(1, [add, sub, mul], e->e.t=TReal))), r]);
    rt := Minimum(opcnts)[2];
    dfts := SumsRuleTree(rt, opts);
    dftc := BlockUnroll(opts.codegen(dfts, tmp2, tmp1, opts), opts);

    # twiddles as lookup table
    i3 := var.fresh_t("i", 2*n);
    twt := var.fresh_t("TW", TPtr(TReal));
    twf := Lambda(i3, cond(eq(i3, 0), 1, eq(i3, 1), 0, nth(twt, 2*((n-1)*j) + 2 * idiv(i3, 2) + (imod(i3, 2) - 2))));

    # debug twiddles
    twl := Map(twf.tolist(), RulesStrengthReduce);
    DoForAll([0..Length(twl)-1], _i->Unparse(assign(nth(Y, _i), twl[_i+1]), CUnparser, 0, 1));

    tws := RCDiag(twf);
    twc :=BlockUnroll(unroll_cmd(opts.codegen(tws, tmp3, tmp2, opts)), opts);

    # stitch code fragments together
    cl := func(TVoid, "transform", [Y, X, twt, j, k, l, m],
            decl([tmp1, tmp2, tmp3], chain(gc, dftc, twc, sc)));
    cc := Compile(cl, opts);

    # fixup to push decls inside (?)
    c := func(TVoid, "transform", [Y, X, twt, j, k, lp],
            decl(cc.vars::[l],
                chain(
                    assign(l, deref(lp)),
                    cc.cmd.cmds[1].cmd
                )
            )
         );

    # print Radix N FFTE kernel as function
#    PrintCode(name, c, opts);
#    PrintTo(name::".c", PrintCode(name, c, opts));

    # function with j and k loop in the function call
    cl := func(TVoid, "transform", [Y, X, twt, lp],
        decl(c.cmd.vars::c.free(),
            chain(
                assign(l, deref(lp)),
                loopn(j, l,
                    c.cmd.cmd)
            )
        )
    );

    # print Radix N FFTE kernel as function
    PrintCode(opts.FortranIze(name::suffix), cl, opts);
    PrintTo(name::filesuffix::suffix::".c", PrintCode(opts.FortranIze(name::suffix), cl, opts));
end;

max_stars_repo_name: spiral-software/spiral-package-ffte

max_issues_repo_name: spiral-software/spiral-package-ffte

max_forks_repo_name: spiral-software/spiral-package-ffte


text:
# The graph obtained from an adjacency function on the vertex set.
BindGlobal("AdjFunGraph", function(E, F)
    return Graph(Group(()), E, function(x, y) return x; end, F, true);
end);

# A generic product graph.
BindGlobal("ProductGraph", function(Gs, F)
    local G, GG, dp;
    dp := DirectProduct(List(Gs, H -> H.group));
    G := Graph(dp, Cartesian(List(Gs, H -> [1..H.order])),
        OnProduct(Length(Gs), dp), F, true);
    for GG in Gs do
        if not "names" in RecNames(GG) then
            GG.names := [1..GG.order];
        fi;
    od;
    AssignVertexNames(G, List(G.names,
        f -> List([1..Length(f)], i -> Gs[i].names[f[i]])));
    return G;
end);

# The box product of two or more graphs.
BindGlobal("BoxProductGraph", function(arg)
    local Gs;
    if Length(arg) = 1 then
        Gs := arg[1];
    else
        Gs := arg;
    fi;
    return ProductGraph(Gs, function(x, y)
        local l;
        l := List([1..Length(Gs)], i -> Distance(Gs[i], x[i], y[i]));
        return WeightVecFFE(l) = 1 and Sum(l) = 1;
    end);
end);

# The cross product of two or more graphs.
BindGlobal("CrossProductGraph", function(arg)
    local Gs;
    if Length(arg) = 1 then
        Gs := arg[1];
    else
        Gs := arg;
    fi;
    return ProductGraph(Gs, function(x, y)
        local l;
        l := List([1..Length(Gs)], i -> Distance(Gs[i], x[i], y[i]));
        return Minimum(l) = 1 and Maximum(l) = 1;
    end);
end);

# The strong product of two or more graphs.
BindGlobal("StrongProductGraph", function(arg)
    local Gs;
    if Length(arg) = 1 then
        Gs := arg[1];
    else
        Gs := arg;
    fi;
    return ProductGraph(Gs, function(x, y)
        return Maximum(List([1..Length(Gs)],
            i -> Distance(Gs[i], x[i], y[i]))) = 1;
    end);
end);

# The bipartite double of a graph.
BindGlobal("BipartiteDoubleGraph", function(G)
    local H;
    H := BipartiteDouble(G);
    CheckDualityFunctions(G);
    H.halfDuality := BipartiteDoubleDualityFunction(G.duality);
    H.halfPrimality := BipartiteDoubleDualityFunction(G.primality);
    return H;
end);

# The extended bipartite double of a graph.
BindGlobal("ExtendedBipartiteDoubleGraph", function(G)
    local dp, signs, H;
    signs := ["+", "-"];
    dp := DirectProduct(G.group, SymmetricGroup(2));
    H := Graph(dp, Cartesian(G.representatives, signs),
        OnSignedPoints(dp, signs), function(x, y)
            return x[2] <> y[2] and Distance(G, x[1], y[1]) <= 1;
        end);
    if "names" in RecNames(G) then
        AssignVertexNames(H, List(H.names, x -> [G.names[x[1]], x[2]]));
    fi;
    CheckDualityFunctions(G);
    H.halfDuality := BipartiteDoubleDualityFunction(G.duality);
    H.halfPrimality := BipartiteDoubleDualityFunction(G.primality);
    return H;
end);

# The halved graph of a bipartite graph. The optional second argument
# allows choosing between the first and second halves.
BindGlobal("HalvedGraph", function(arg)
    local n, G, G2, H, vs;
    if Length(arg) < 1 then
        Error("at least one argument expected");
        return fail;
    fi;
    G := arg[1];
    if Length(arg) > 1 then
        n := arg[2];
    else
        n := 1;
    fi;
    if not IsConnectedGraph(G) then
        Error("not a connected graph");
        return fail;
    fi;
    if not IsBipartite(G) then
        Error("not a bipartite graph");
        return fail;
    fi;
    G2 := DistanceGraph(G, 2);
    vs := ConnectedComponent(G2, 1);
    if n = 2 then
        vs := Difference([1..G.order], vs);
    fi;
    H := Graph(Stabilizer(G.group, vs, OnSets), vs, OnPoints,
        function(x, y) return Distance(G2, x, y) = 1; end, true);
    if "names" in RecNames(G) then
        AssignVertexNames(H, G.names{vs});
    fi;
    if "halfDuality" in RecNames(G) then
        if n = 2 then
            H.primality := G.halfDuality;
        else
            H.duality := G.halfDuality;
        fi;
    fi;
    if "halfPrimality" in RecNames(G) then
        if n = 2 then
            H.duality := G.halfPrimality;
        else
            H.primality := G.halfPrimality;
        fi;
    fi;
    return H;
end);

# The antipodal quotient of an antipodal cover.
BindGlobal("AntipodalQuotientGraph", function(G)
    local d, H;
    if not IsAntipodalCover(G) then
        Error("not an antipodal cover");
        return fail;
    fi;
    d := Diameter(G);
    H := Graph(G.group,
        Set(List(G.representatives, x -> DistanceSet(G, [0, d], x))),
        OnSets, function(x, y)
            return IsSubset(DistanceSet(G, 1, x), y);
        end);
    if "names" in RecNames(G) then
        AssignVertexNames(H, List(H.names, f -> G.names{f}));
    fi;
    return H;
end);

# A graph with the set of d-dimensional subspaces of V filtered by S
# as the vertex set, acted upon by the matrix group G,
# with two subspaces being adjacent iff their intersection has dimension d-1.
BindGlobal("SubspaceGraph", function(arg)
    local G, H, S, V, d, invt, vcs;
    if Length(arg) < 4 then
        Error("at least four arguments expected");
        return fail;
    fi;
    G := arg[1];
    S := arg[2];
    V := arg[3];
    d := arg[4];
    if Length(arg) > 4 then
        invt := arg[5];
    else
        invt := true;
    fi;
    if IsList(S) then
        vcs := S;
    else
        vcs := S(Subspaces(V, d));
    fi;
    H := Graph(G, vcs, OnSubspaces(V), function(x,y)
                    return Dimension(Intersection(x,y)) = d-1;
                end, invt);
    H.duality := Intersection;
    H.primality := Sum;
    return H;
end);

# The clique (dual geometry) graph of a collinearity graph. The optional second
# argument allows choosing a connected component of the resulting graph.
BindGlobal("CliqueGraph", function(arg)
    local C, G, H, n;
    if Length(arg) < 1 then
        Error("at least one argument expected");
        return fail;
    fi;
    G := arg[1];
    C := Cliques(G);
    if Length(arg) > 1 then
        n := arg[2];
        if not IsList(n) then
            n := [n];
        fi;
    else
        n := [1..Length(C)];
    fi;
    H := Graph(G.group, C{n}, OnSets,
                function(x, y)
                    return Size(Intersection(x,y)) = 1;
                end);
    if "names" in RecNames(G) then
        CheckDualityFunctions(G);
        H.duality := G.primality;
        H.primality := G.duality;
        AssignVertexNames(H, List(H.names, f -> G.duality(G.names{f})));
        if "halfDuality" in RecNames(G) then
            H.halfDuality := G.halfDuality;
        fi;
        if "halfPrimality" in RecNames(G) then
            H.halfPrimality := G.halfPrimality;
        fi;
    fi;
    return H;
end);

# The incidence graph of a collinearity graph.
BindGlobal("IncidenceGraph", function(arg)
    local C, G, H, n;
    if Length(arg) < 1 then
        Error("at least one argument expected");
        return fail;
    fi;
    G := arg[1];
    C := Cliques(G);
    if Length(arg) > 1 then
        n := arg[2];
        if not IsList(n) then
            n := [n];
        fi;
    else
        n := [1..Length(C)];
    fi;
    H := Graph(G.group, Union(G.representatives, C{n}),
                OnPointsOrLines(OnPoints, IsList),
                function(x, y)
                    return (IsList(y) and x in y) or (IsList(x) and y in x);
                end);
    if "names" in RecNames(G) and G.order > 0 then
        CheckDualityFunctions(G);
        if IsList(H.names[1]) then
            H.halfDuality := G.primality;
            H.halfPrimality := G.duality;
        else
            H.halfDuality := G.duality;
            H.halfPrimality := G.primality;
        fi;
        AssignVertexNames(H, List(H.names, function(f)
                                            if IsList(f) then
                                                return G.duality(G.names{f});
                                            else
                                                return G.names[f];
                                            fi;
        end));
    fi;
    return H;
end);


text:
#
#
#

Read("~/Workspace/Chevalley.gap/init.gi");

Read(Filename(home_dir,"lib/rsys.gd"));
Read(Filename(home_dir,"lib/rsys.gi"));

Read(Filename(home_dir,"lib/chvadj.gd"));
Read(Filename(home_dir,"lib/chvadj.gi"));

Read(Filename(home_dir,"lib/nilchv.gd"));
Read(Filename(home_dir,"lib/nilchv.gi"));


text:
LETTERS;
# "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz"
LETTERS{[5 .. 10]};
# "EFGHIJ"


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

memo := function(context, prefix, exp)
   local v, isum;
   exp := SReduce(toExpArg(exp), SpiralDefaults); ##Pass the real opts instead
   if IsValue(exp) then return exp.v;
   elif IsVar(exp) then return exp;
   else
       if IsBound(context.ISum) and context.ISum <> [] then
       isum := Last(context.ISum);
       v := var.fresh_t(prefix, TInt);
       if not IsBound(isum.memos) then isum.memos := []; fi;
       Add(isum.memos, [v, exp]);
       #Print(v.id , " -> ", exp, " [", isum.var, " in ", isum.domain, "]\n");
       v.mapping := exp;
       return v;
       else
       return exp;
       fi;
   fi;
end;

nomemo := (cx,pfx,exp) -> exp;

memo := nomemo; # memo feature is currently broken -- disable it

Class(ProcessMemos, RuleSet);
RewriteRules(ProcessMemos, rec(
   # pull out cmemos
   PullOutCmemo := Rule(cmemo,
     function(e,cx)
         local target;
	 target := cx.(e.target);
         if not IsBound(target.memos) then
	     target.memos := [ e.args[1], e.mapping ];
	 else Add(target.memos, [ e.args[1], e.mapping ]);
	 fi;
	 return e.args[1];
     end),

   # process memos at ISum
   ProcessISumMemos := Rule(@(1, [ISum, ISumAcc], e->IsBound(e.memos)),
     e -> ObjId(@(1).val)(e.var, e.domain,
	      FoldL(Reversed(e.memos),
		  (ex,m) -> Data(m[1], m[2], ex)), e.child(1)))
));


text:
#
# cls 66: A_2
#
handle66char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[66];
    info:=infos[66];

    rep:=Representative(o);
    #pr[19]=[0,0,1,1,1,0,0,0]
    #pr[17]=[0,1,1,1,0,0,0,0]
    #pr[37]=[1,1,1,2,1,0,0,0]
    tmp1:=ApplyRootsReflections(rep,[19,17,37]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,1]]);
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[3],[-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[19,17,37]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,1]]);
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[3],[-1]);

    return [tmp1,tmp2];
end;

#
# cls 64: A_2A_1
#
handle64char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[64];
    info:=infos[64];

    rep:=Representative(o);
    #pr[53]=[1,1,1,2,1,1,1,0]
    #pr[52]=[1,1,1,2,2,1,0,0]
    #pr[93]=[1,2,3,4,3,2,1,0]
    tmp1:=ApplyRootsReflections(rep,[53,52,93]);
    uu:=Unipotent(chevalleyAdj(rep),[[1,-1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[1],[-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[53,52,93]);
    uu:=Unipotent(chevalleyAdj(cent),[[1,-1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[1],[-1]);

    return [tmp1,tmp2];
end;

#
# cls 60: A_2^2
#
handle60char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[60];
    info:=infos[60];

    rep:=Representative(o);
    #pr[43]=[0,0,1,1,1,1,1,1]
    #pr[41]=[0,1,1,1,1,1,1,0]
    #pr[79]=[1,1,1,2,2,2,2,1]
    tmp1:=ApplyRootsReflections(rep,[43,41,79]);
    uu:=Unipotent(chevalleyAdj(rep),[[1,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[3],[-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[43,41,79]);
    uu:=Unipotent(chevalleyAdj(cent),[[1,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[3],[-1]);

    return [tmp1,tmp2];
end;

#
# cls 57: D_4(a_1)
#
handle57char5:=function(orbs,infos)
    local tmp1,tmp2,tmp3,tmp4,o,info,rep,cent,uu;
    o:=Classes(orbs)[57];
    info:=infos[57];

    rep:=Representative(o);
    #pr[4] =[0,0,0,1,0,0,0,0]
    #pr[51]=[1,1,2,2,1,1,0,0]
    #pr[52]=[1,1,1,2,2,1,0,0]
    tmp1:=ApplyRootsReflections(rep,[4,51,52]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,-1],[5,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[4,1,2,6],[-1,-1,-1,-1]);

    #pr[26]=[0,1,0,1,1,1,0,0]
    #pr[27]=[0,0,1,1,1,1,0,0]
    tmp3:=ApplyRootsReflections(rep,[26,27]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,1],[4,1],[5,-1],[11,1],[12,-1]],Ordering(rep));
    tmp3:=Conj(tmp3,uu);
    tmp3:=ConjugateByTori(tmp3,[2,4],[-1,-1,-1,-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[4,51,52]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,-1],[5,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[4,1,2,6],[-1,-1,-1,-1]);

    tmp4:=0;
    tmp4:=ApplyRootsReflections(cent,[26,27]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,1],[4,1],[5,-1],[11,1],[12,-1]],Ordering(cent));
    tmp4:=Conj(tmp4,uu);
    tmp4:=ConjugateByTori(tmp4,[2,4],[-1,-1,-1,-1]);

    return [tmp1,tmp3,tmp2,tmp4];
end;

#
# cls 53: D_4(a_1)A_1 Ross' rep
#
handle53char5:=function(orbs,infos)
    local tmp1,tmp2,tmp3,tmp4,o,info,rep,cent,uu;
    o:=Classes(orbs)[53];
    info:=infos[53];

    rep:=Representative(o);

    #pr[23]=[1,1,1,1,0,0,0,0]
    #pr[24]=[1,0,1,1,1,0,0,0]
    tmp1:=ApplyRootsReflections(rep,[23,24]);
    tmp1:=ConjNegUa(tmp1,4,1);
    uu:=Unipotent(chevalleyAdj(rep),[[2,-1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[2],[-1,-1,-1,-1]);

    #pr[65]=[1,1,2,2,1,1,1,1]
    #pr[67]=[1,1,1,2,2,1,1,1]
    tmp3:=ApplyRootsReflections(rep,[65,67,4]);
    uu:=Unipotent(chevalleyAdj(rep),[[12,-1],[11,-1]],Ordering(rep));
    tmp3:=Conj(tmp3,uu);
    tmp3:=ConjugateByTori(tmp3,[1,2,4,6,8],[-1,-1,-1,-1,-1,-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[23,24]);
    tmp2:=ConjNegUa(tmp2,4,1);
    uu:=Unipotent(chevalleyAdj(cent),[[2,-1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[2],[-1,-1,-1,-1]);

    tmp4:=0;
    tmp4:=ApplyRootsReflections(cent,[65,67,4]);
    uu:=Unipotent(chevalleyAdj(cent),[[12,-1],[11,-1]],Ordering(cent));
    tmp4:=Conj(tmp4,uu);
    tmp4:=ConjugateByTori(tmp4,[1,2,4,6,8],[-1,-1,-1,-1,-1]);

    return [tmp1,tmp3,tmp2,tmp4];
end;

#
# cls 52: A_3A_2 [3,4,2,6,7]
#
handle52char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[52];
    info:=infos[52];

    rep:=Representative(o);

    #pr[26]=[0,1,0,1,1,1,0,0]
    #pr[27]=[0,0,1,1,1,1,0,0]
    #pr[55]=[0,1,1,2,2,1,1,0]
    tmp1:=ApplyRootsReflections(rep,[26,27,55]);
    uu:=Unipotent(chevalleyAdj(rep),[[2,-1],[3,1],[6,-1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[4,6],[-1,-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[26,27,55]);
    uu:=Unipotent(chevalleyAdj(cent),[[2,-1],[3,1],[6,-1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[4,6],[-1,-1]);

    return [tmp1,tmp2];
end;

#
# cls 51: A_4 [1,2,3,4]
#
handle51char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[51];
    info:=infos[51];

    rep:=Representative(o);

    #pr[33]=[0,1,1,1,1,1,0,0]
    #pr[31]=[1,0,1,1,1,1,0,0]
    #pr[64]=[1,1,2,2,2,1,1,0]
    #pr[89]=[1,2,2,4,3,2,1,0]
    tmp1:=ApplyRootsReflections(rep,[33,31,64,89]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,1],[4,3],[2,2],[10,2],[11,-1],[17,2]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[1,3,4],[-1,-1,-1,-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[33,31,64,89]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,1],[4,3],[2,2],[10,2],[11,-1],[17,2]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[1,3,4],[-1,-1,-1,-1]);

    return [tmp1,tmp2];
end;
#vals[2]:=vals[1]+2*One(APR);
#vals[3]:=-(vals[1]-One(APR));
#vals[4]:=-(-vals[2]-One(APR));
#vals[10]:=vals[1]^2+2*vals[1]+vals[9]+3*One(APR);
#vals[11]:=-(-vals[1]-vals[9]+One(APR));
#vals[17]:=-(-vals[1]^2+3*vals[1]+vals[16]+3*One(APR));
# nu am terminat dar am gasit elementul

#
# cls 48: D_4(a_1)A_2
#
handle48char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[48];
    info:=infos[48];

    rep:=Representative(o);

    #pr[23]=[1,1,1,1,0,0,0,0]  asta merge si cu rep=[2,3,10,12,7,8] urmat de ConjNegUa(tmp1,4,1)
    #pr[24]=[1,0,1,1,1,0,0,0]
    tmp1:=ApplyRootsReflections(rep,[23,24]);
    tmp1:=ConjNegUa(tmp1,4,1);
    uu:=Unipotent(chevalleyAdj(rep),[[2,-1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[2],[-1,-1,-1,-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[23,24]);
    tmp2:=ConjNegUa(tmp2,4,1);
    uu:=Unipotent(chevalleyAdj(cent),[[2,-1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[2],[-1,-1,-1,-1]);

    return [tmp1,tmp2];
end;

#
# cls 47: A_4A_1
#
handle47char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[47];
    info:=infos[47];

    rep:=Representative(o);

    #pr[41]=[0,1,1,1,1,1,1,0]
    #pr[39]=[1,0,1,1,1,1,1,0]
    #pr[57]=[1,1,2,2,2,1,0,0]
    #pr[89]=[1,2,2,4,3,2,1,0]
    tmp1:=ApplyRootsReflections(rep,[41,39,57,89]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,-1],[4,2],[2,2],[10,3],[11,-1],[17,2]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[2,3,4],[-1,-1,-1,-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[41,39,57,89]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,-1],[4,2],[2,2],[10,3],[11,-1],[17,2]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[2,3,4],[-1,-1,-1,-1]);

    return [tmp1,tmp2];
end;
#vals[2]:=vals[1]+2*One(APR);
#vals[3]:=-(-vals[1]+One(APR));
#vals[4]:=-(vals[2]+One(APR));
#vals[10]:=-(vals[1]^2+2*vals[1]-vals[9]+3*One(APR));
#vals[11]:=-(-vals[1]+vals[9]+One(APR));
#vals[17]:=-(-vals[1]^2+3*vals[1]+vals[16]+3*One(APR));
#
# nu e gata dar am gasit elementul

#
# cls 45: D_5(a_1)
#
handle45char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[45];
    info:=infos[45];

    rep:=Representative(o);

    #pr[82]=[1,1,2,3,3,2,1,0]
    #pr[80]=[1,2,2,3,2,2,1,0]
    tmp1:=ApplyRootsReflections(rep,[82,80]);
    uu:=Unipotent(chevalleyAdj(rep),[[4,1],[12,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[2,4],[-1,-1,-1,-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[82,80]);
    uu:=Unipotent(chevalleyAdj(cent),[[4,1],[12,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[2,4],[-1,-1,-1,-1]);

    return [tmp1,tmp2];
end;

#
# cls 44: A_4A_1^2
#
handle44char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[44];
    info:=infos[44];

    rep:=Representative(o);

    #pr[71]=[1,1,2,2,2,2,1,0]
    #pr[72]=[1,1,2,2,2,1,1,1]
    #pr[112]=[1,3,3,5,4,3,2,1]
    #pr[111]=[2,2,3,5,4,3,2,1]
    tmp1:=ApplyRootsReflections(rep,[71,72,112,111]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,1],[4,3],[2,3],[10,3],[11,-1],[17,3]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[2,7],[-1,-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[71,72,112,111]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,1],[4,3],[2,3],[10,3],[11,-1],[17,3]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[2,7],[-1,-1]);

    return [tmp1,tmp2];
end;
#vals[2]:=-(-vals[1]+2*One(APR));
#vals[3]:=vals[1]+One(APR);
#vals[4]:=-vals[2]+One(APR);
#vals[10]:=-vals[1]^2+2*vals[1]+vals[9]+2*One(APR);
#vals[11]:=-vals[1]-vals[9]-One(APR);
#vals[17]:=-vals[1]^2+2*vals[1]-vals[16]+3*One(APR);
#
# nu e gata dar am gasit elementul

#
# cls 39: D_4A_2
#
handle39char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[39];
    info:=infos[39];

    rep:=Representative(o);

    #pr[34]=[0,1,0,1,1,1,1,0]
    #pr[35]=[0,0,1,1,1,1,1,0]
    #pr[68]=[0,1,1,2,2,2,1,1]
    tmp1:=ApplyRootsReflections(rep,[34,35,68]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,2],[4,-1],[5,-1],[8,-1],[11,-1],[12,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[4,6,7,8],[-1,-1,-1,-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[34,35,68]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,2],[4,-1],[5,-1],[8,-1],[11,-1],[12,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[4,6,7,8],[-1,-1,-1,-1]);

    return [tmp1,tmp2];
end;
#vals[3]:=-(-vals[2]+3*One(APR));
#vals[4]:=-(vals[2]+One(APR));
#vals[5]:=vals[4];
#vals[8]:=-vals[7]-One(APR);
#vals[11]:=-(vals[2]*vals[3]+vals[10]+One(APR));
#vals[12]:=-(-vals[2]+vals[10]-One(APR));
#
# nu e gata dar am gasit elementul

#
# cls 38: E_6(a_3) [1,4,6,17,18,19] this doesn't work in char2
#
handle38char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent;
    o:=Classes(orbs)[38];
    info:=infos[38];

    rep:=Representative(o);
    tmp1:=ConjugateByTorus(rep,4,-1);# h_1 --> h_4

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjugateByTorus(cent,4,-1);# h_1 --> h_4

    return [tmp1,tmp2];
end;

#
# cls 33: D_6(a_2)
#
handle33char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[33];
    info:=infos[33];

    rep:=Representative(o);

    #pr[42]=[0,1,0,1,1,1,1,0]
    #pr[43]=[0,0,1,1,1,1,1,0]
    tmp1:=ApplyRootsReflections(rep,[42,43]);
    uu:=Unipotent(chevalleyAdj(rep),[[4,1],[6,1],[10,1],[13,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[4,5,8],[-1,-1,-1,-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[42,43]);
    uu:=Unipotent(chevalleyAdj(cent),[[4,1],[6,1],[10,1],[13,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[4,5,8],[-1,-1,-1,-1]);

    return [tmp1,tmp2];
end;

#
# cls 32: E_6(a_3)A_1 [1,4,6,8,17,18,19] this doesn't work in char2
#
handle32char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent;
    o:=Classes(orbs)[32];
    info:=infos[32];

    rep:=Representative(o);
    tmp1:=ConjugateByTorus(rep,4,-1);# h_1 --> h_4

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjugateByTorus(cent,4,-1);# h_1 --> h_4

    return [tmp1,tmp2];
end;

#
# cls 31: E_7(a_5)
#
handle31char5:=function(orbs,infos)
    local tmp1,tmp2,tmp3,tmp4,o,info,rep,cent,uu;
    o:=Classes(orbs)[31];
    info:=infos[31];

    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[2,3,5]);
    uu:=Unipotent(chevalleyAdj(rep),[[4,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[3,4],[-1,-1]);

    tmp3:=ConjugateByTori(rep,[2,3,5],[bbb,bbb,bbb]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[2,3,5]);
    uu:=Unipotent(chevalleyAdj(cent),[[4,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[3,4],[-1,-1]);

    tmp4:=0;
    tmp4:=ConjugateByTori(cent,[2,3,5],[bbb,bbb,bbb]);

    return [tmp1,tmp3,tmp2,tmp4];
end;

#
# cls 29: E_8(a_7)                           # NOT FINISHED IN CHARACTERISTIC 5
#
handle29char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent;
    o:=Classes(orbs)[29];
    #info:=infos[29];

    rep:=Representative(o);
    tmp1:=ConjugateByTorus(rep,4,-1);# h_1 --> h_4

    tmp2:=0;
    #cent:=FromPositiveBorel(o,info[2]);
    #tmp2:=ConjugateByTorus(cent,4,-1);# h_1 --> h_4

    return [tmp1,tmp2];
end;

#
# cls 27: D_6(a_1)
#
handle27char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[27];
    info:=infos[27];

    rep:=Representative(o);

    #pr[42]=[0,1,0,1,1,1,1,0]
    #pr[43]=[0,0,1,1,1,1,1,0]
    tmp1:=ApplyRootsReflections(rep,[42,43]);
    uu:=Unipotent(chevalleyAdj(rep),[[4,1],[11,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[3,4],[-1,-1,-1,-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[42,43]);
    uu:=Unipotent(chevalleyAdj(cent),[[4,1],[11,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[3,4],[-1,-1,-1,-1]);

    return [tmp1,tmp2];
end;

#
# cls 25: E_7(a_4) [1,4,7,14,17,18,19,26] this doesn't work in char2
#
handle25char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent;
    o:=Classes(orbs)[25];
    info:=infos[25];

    rep:=Representative(o);
    tmp1:=ConjugateByTorus(rep,4,-1);# h_1 --> h_4

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjugateByTorus(cent,4,-1);# h_1 --> h_4

    return [tmp1,tmp2];
end;

#
# cls 24: E_6(a_1)
#
handle24char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[24];
    info:=infos[24];

    #pr[58]=[1,1,2,2,1,1,1,0]
    #pr[59]=[1,1,1,2,2,1,1,0]
    #pr[61]=[0,1,1,2,2,2,1,0]

    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[58,59,61]);
    tmp1:=ConjNegUa(tmp1,4,-1);
    uu:=Unipotent(chevalleyAdj(rep),[[2,1],[11,-1],[12,3],[10,-1],[5,2],[6,3],[13,3],[18,-1],[19,-1],[20,2],[26,2],[27,2],[32,1],[40,3],[48,-1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[2],[-1,-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[58,59,61]);
    tmp2:=ConjNegUa(tmp2,4,-1);
    uu:=Unipotent(chevalleyAdj(cent),[[2,1],[11,-1],[12,3],[10,-1],[5,2],[6,3],[13,3],[18,-1],[19,-1],[20,2],[26,2],[27,2],[32,1],[40,3],[48,-1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[2],[-1,-1]);

    return [tmp1,tmp2,uu];
end;
#vals[2]:=vals[1]+One(APR);
#vals[4]:=0;
#vals[3]:=0;
#vals[5]:=-(-vals[1]+3*One(APR));
#vals[6]:=vals[5]+One(APR);
#vals[10]:=-(-vals[1]+vals[5]-One(APR));
#vals[11]:=-vals[1]-One(APR);
#vals[12]:=-vals[2]+vals[10];
#vals[13]:=-(2*vals[1]^2+2*vals[9]+2*One(APR));
#vals[16]:=-(2*One(APR))^-1*(vals[1]^2+vals[1]-vals[9]);
#vals[17]:=-vals[1]^2-vals[1]-vals[16];
#vals[18]:=-(3*vals[1]^2+3*vals[1]+2*vals[9]-One(APR));
#vals[19]:=vals[1]-vals[16]+One(APR);
#vals[20]:=vals[1]^2+3*vals[1]+vals[18]+2*One(APR);
#
#vals[24]:=-(3*One(APR))^-1*(vals[1]^3+3*vals[1]*vals[9]-vals[1]+vals[23]);
#vals[25]:=vals[1]^3-vals[1]*vals[9]-vals[1]+2*vals[9]+vals[23]-vals[24];
#vals[26]:=-(2*vals[1]^2-vals[1]-vals[9]+vals[23]-One(APR));
#vals[27]:=-vals[1]^3+vals[1]*vals[9]+2*vals[1]+2*vals[9]+vals[24]+One(APR);
#
#vals[30]:=-vals[1]^4+3*vals[1]^3+2*vals[1]^2*vals[9]-vals[1]^2+2*vals[1]*vals[9]-vals[9]^2-vals[23];
#vals[32]:=vals[1]^3+2*vals[1]^2+vals[9]+vals[23]-One(APR);
#vals[33]:=2*vals[1]^4+vals[1]^3-vals[1]^2*vals[9]+2*vals[1]*vals[9]+2*vals[1]-vals[23]+vals[24]+vals[25]-vals[27]+vals[30]-vals[31]+One(APR);
#vals[37]:=-(vals[1]^4+3*vals[1]^3+2*vals[1]^2*vals[9]+vals[1]^2-vals[1]*vals[9]+vals[9]^2-vals[1]-vals[9]+2*vals[23]+vals[31]);
#
#vals[40]:=2*vals[1]^4-vals[1]^3+3*vals[1]^2*vals[9]+3*vals[1]*vals[9]+2*vals[1]-vals[9]+vals[37]-One(APR);
#vals[44]:=-(vals[1]^4+3*vals[1]^3*vals[9]+2*vals[1]^3+2*vals[1]^2*vals[9]+2*vals[1]*vals[9]^2+vals[1]^2-vals[1]*vals[9]+vals[23]-vals[38]);
#vals[48]:=3*vals[1]^5+3*vals[1]^3*vals[9]+3*vals[1]^3+vals[1]^2*vals[9]+vals[1]*vals[9]+vals[9]^2+2*vals[1]+2*vals[9]+vals[23]+vals[37]-vals[38]+vals[45]+3*One(APR);
#
#vals[51]:=-(3*vals[1]^3*vals[9]+2*vals[1]^3*vals[23]-vals[1]^2*vals[9]^2+3*vals[1]^2*vals[9]-vals[1]*vals[23]-vals[9]^2+2*vals[9]*vals[23]-vals[23]^2+vals[38]);
#
# nu e gata dar am gasit elementul

#
# cls 23: D_5A_2
#
handle23char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[23];
    info:=infos[23];

    rep:=Representative(o);

    #pr[87]=[1,1,2,3,3,2,1,1]
    #pr[86]=[1,2,2,3,2,2,1,1]
    #pr[119]=[2,3,4,6,5,4,3,1]
    tmp1:=ApplyRootsReflections(rep,[87,86,119]);
    uu:=Unipotent(chevalleyAdj(rep),[[8,-1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[1,4,7],[-1,-1,-1,-1,-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[87,86,119]);
    uu:=Unipotent(chevalleyAdj(cent),[[8,-1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[1,4,7],[-1,-1,-1,-1]);

    return [tmp1,tmp2];
end;

#
# cls 20: D_7(a_2)
#
handle20char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[20];
    info:=infos[20];

    rep:=Representative(o);

    #pr[113]=[2,3,3,5,4,3,2,1]
    #pr[114]=[2,2,4,5,4,3,2,1]
    tmp1:=ApplyRootsReflections(rep,[113,114]);
    uu:=Unipotent(chevalleyAdj(rep),[[4,1],[6,1],[11,-1],[13,-1],[20,-1],[27,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[5,6],[-1,-1,-1,-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[113,114]);
    uu:=Unipotent(chevalleyAdj(cent),[[4,1],[6,1],[11,-1],[13,-1],[20,-1],[27,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[5,6],[-1,-1,-1,-1]);

    return [tmp1,tmp2];
end;

#
# cls 18: E_6(a_1)A_1
#
handle18char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[18];
    info:=infos[18];

    #pr[106]=[1,2,2,4,4,3,2,1]
    #pr[105]=[1,2,3,4,3,3,2,1]
    #pr[104]=[2,2,3,4,3,2,2,1]

    # this is for E_6(a_1)A_1
    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[106,105,104]);
    tmp1:=ConjNegUa(tmp1,4,-1);
    uu:=Unipotent(chevalleyAdj(rep),[[2,1],[11,1],[12,2],[5,3],[10,-1],[6,2],[13,3],[18,-1],[19,-1],[20,2],[26,1],[27,3],[32,1],[40,2]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[1,2,3,5,6],[-1,-1,-1,-1,-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[106,105,104]);
    tmp2:=ConjNegUa(tmp2,4,-1);
    uu:=Unipotent(chevalleyAdj(cent),[[2,1],[11,1],[12,2],[5,3],[10,-1],[6,2],[13,3],[18,-1],[19,-1],[20,2],[26,1],[27,3],[32,1],[40,2]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[1,2,3,5,6],[-1,-1,-1,-1,-1]);

    return [tmp1,tmp2,uu];
end;
#vals[2]:=-(vals[1]-One(APR));
#vals[4]:=0;
#vals[3]:=0;
#vals[5]:=vals[1]+3*One(APR);
#vals[6]:=-(-vals[5]+One(APR));
#vals[10]:=-(vals[1]-vals[5]-One(APR));
#vals[11]:=-(vals[1]-One(APR));
#vals[12]:=vals[2]-vals[10];
#vals[13]:=-(2*One(APR))^-1*(-vals[1]^2-vals[9]-One(APR));
#
#vals[16]:=3*vals[1]+vals[9]-vals[13]+3*One(APR);
#vals[17]:=-(-vals[1]^2+vals[1]-vals[16]);
#vals[18]:=-vals[16]-One(APR);
#vals[19]:=-(-vals[1]^2+2*vals[1]+vals[9]-vals[16]-One(APR));
#vals[20]:=vals[1]^2+2*vals[1]-vals[18]+2*One(APR);
#vals[23]:=-(3*One(APR))^-1*(vals[1]^3+2*vals[1]^2+3*vals[1]*vals[9]-vals[1]+2*vals[9]+vals[13]+2*One(APR));
#vals[24]:=-(3*One(APR))^-1*(3*vals[1]^3-vals[1]*vals[9]+2*vals[1]);
#vals[25]:=-(vals[1]^3+vals[1]^2-vals[1]*vals[9]-vals[1]+vals[9]-vals[16]+vals[17]+vals[19]-vals[23]-vals[24]-One(APR));
#vals[26]:=-vals[1]+3*vals[9]+vals[13]-vals[23]+3*One(APR);
#vals[27]:=-(vals[1]^3-vals[1]*vals[9]+3*vals[1]+2*vals[9]-vals[24]+One(APR));
#
#vals[30]:=-(2*vals[1]^3-vals[1]^2*vals[9]+vals[1]^2*vals[13]+vals[1]^2+3*vals[1]*vals[9]+vals[9]*vals[13]-vals[13]^2+2*vals[9]+vals[13]-vals[23]+One(APR));
#vals[32]:=-(vals[1]^3+3*vals[1]^2-vals[9]-vals[23]+One(APR));
#vals[33]:=-(-vals[1]^4+vals[1]^3+vals[1]^2*vals[9]+vals[1]^2*vals[13]+2*vals[1]^2+2*vals[1]*vals[9]-vals[1]*vals[13]+2*vals[1]+2*vals[9]+vals[27]-vals[30]-vals[31]+One(APR));
#vals[37]:=vals[1]^4+3*vals[1]^3+2*vals[1]^2*vals[9]+vals[1]^2-vals[1]*vals[9]+vals[9]^2-vals[9]-vals[23]+vals[31];
#
#vals[40]:=-(2*vals[1]^3+vals[1]^2*vals[9]-vals[1]^2*vals[13]-vals[1]^2+3*vals[1]*vals[9]+3*vals[1]*vals[13]-vals[1]+3*vals[13]-vals[37]);
#
#vals[44]:=-vals[1]^4+3*vals[1]^3*vals[9]-vals[1]^3+3*vals[1]^2*vals[9]+2*vals[1]*vals[9]^2-vals[1]^2+3*vals[1]+3*vals[23]+vals[38];
#
#vals[48]:=-(3*vals[1]^5+3*vals[1]^3*vals[9]-vals[1]^2*vals[9]+2*vals[1]*vals[9]-vals[9]^2+3*vals[9]+3*vals[23]+vals[37]+vals[38]-vals[45]+2*One(APR));
#
#vals[51]:=-vals[1]^6+2*vals[1]^4*vals[9]-vals[1]^4+3*vals[1]^3*vals[9]+2*vals[1]^3*vals[23]+3*vals[1]^3*vals[24]-vals[1]^2*vals[9]^2+2*vals[1]^2*vals[9]+3*vals[1]*vals[9]^2+3*vals[1]*vals[9]*vals[23]-vals[1]*vals[9]*vals[24]+2*vals[1]^2-vals[1]*vals[9]+vals[1]*vals[23]+2*vals[1]*vals[24]-vals[9]^2-vals[9]*vals[23]+3*vals[23]*vals[24]-vals[24]^2+vals[38];
#
# nu e gata dar am gasit elementul

#
# cls 17: E_7(a_3) ([1,4,6,7,15,16,17] in E_7) [1,4,6,7,17,18,19] in E_8 this doesn't work in char2
#
handle17char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent;
    o:=Classes(orbs)[17];
    info:=infos[17];

    rep:=Representative(o);
    tmp1:=ConjugateByTorus(rep,4,-1);# h_1 --> h_4

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjugateByTorus(cent,4,-1);# h_1 --> h_4

    return [tmp1,tmp2];
end;

#
# cls 16: E_8(b_6)
#
handle16char5:=function(orbs,infos)
    local tmp1,tmp2,tmp3,tmp4,o,info,rep,cent,uu;
    o:=Classes(orbs)[16];
    info:=infos[16];

    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[1,2,5]);
    uu:=Unipotent(chevalleyAdj(rep),[[23,-1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[1,3,4,8],[-1,-1,-1,-1]);

    tmp3:=ConjugateByTori(rep,[1,2,5],[bbb,bbb,bbb^2]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[1,2,5]);
    uu:=Unipotent(chevalleyAdj(cent),[[23,-1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[1,3,4,8],[-1,-1,-1,-1]);

    tmp4:=0;
    tmp4:=ConjugateByTori(cent,[1,2,5],[bbb,bbb,bbb^2]);

    return [tmp1,tmp3,tmp2,tmp2];
end;

#
# cls 15: D_7(a_1)
#
handle15char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=Classes(orbs)[15];
    info:=infos[15];

    rep:=Representative(o);

    #pr[113]=[2,3,3,5,4,3,2,1]
    #pr[114]=[0,0,1,1,1,1,1,0]
    tmp1:=ApplyRootsReflections(rep,[113,114]);
    #uu:=Unipotent(chevalleyAdj(rep),[[4,1],[10,-1]],Ordering(rep)); si asta merge
    uu:=Unipotent(chevalleyAdj(rep),[[4,1],[11,-1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[2,3,4],[3*One(APR),3*One(APR),-1]);
    #tmp1:=ConjugateByTori(tmp1,[2,3,4],[ccc,ccc,-1]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[113,114]);
    uu:=Unipotent(chevalleyAdj(cent),[[4,1],[11,-1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[2,3,4],[3*One(APR),3*One(APR),-1]);
    #tmp2:=ConjugateByTori(tmp2,[2,3,4],[ccc,ccc,-1]);

    return [tmp1,tmp2];
end;

#
# cls 12: E_8(a_6)
#
handle12char5:=function(orbs,infos)
    local tmp1,tmp2,tmp3,tmp4,o,info,rep,cent,uu,vars;
    o:=Classes(orbs)[12];
    info:=infos[12];

    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[2,3,8,13]);
    uu:=Unipotent(chevalleyAdj(rep),[[18,-1],[29,-1],[31,1],[33,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[2,3,4,5,7,8],[-1,-1,-1,-1,-1,-1]);

    #tmp3:=ConjugateByTori(rep,[2,3,5,6,8],[aaa,aaa,aaa,aaa,aaa^2]);
    tmp3:=ConjugateByTori(rep,[2,3,5,6,8],[bbb,bbb,bbb,bbb,bbb^2]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[2,3,8,13]);
    uu:=Unipotent(chevalleyAdj(cent),[[18,-1],[29,-1],[31,1],[33,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[2,3,4,5,7,8],[-1,-1,-1,-1,-1,-1]);

    tmp4:=0;
    tmp4:=ConjugateByTori(cent,[2,3,5,6,8],[bbb,bbb,bbb,bbb,bbb^2]);

    return [tmp1,tmp3,tmp2,tmp4];
end;

#
# cls 10: E_8(b_5)
#
handle10char5:=function(orbs,infos)
    local tmp1,tmp2,tmp3,tmp4,o,info,rep,cent,uu;
    o:=Classes(orbs)[10];
    info:=infos[10];

    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[1,2,5]);
    uu:=Unipotent(chevalleyAdj(rep),[[19,-1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[1,3,4],[-1,-1,-1]);

    tmp3:=ConjugateByTori(rep,[1,2,5],[bbb,bbb,bbb^2]);

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[1,2,5]);
    uu:=Unipotent(chevalleyAdj(cent),[[19,-1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[1,3,4],[-1,-1,-1]);

    tmp4:=0;
    tmp4:=ConjugateByTori(cent,[1,2,5],[bbb,bbb,bbb^2]);

    return [tmp1,tmp3,tmp2,tmp4];
end;

# For representative deduced from Seitz' E_8(a_6) used in proof
#
# cls 8: E_8(a_5)
#
handle8char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent;
    o:=Classes(orbs)[8];
    info:=infos[8];

    rep:=Representative(o);
    tmp1:=ConjugateByTori(rep,[4,6,7],[-1,-1,-1,-1]);# should be Z(D_8)

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjugateByTori(cent,[4,6,7],[-1,-1,-1,-1]);# should be Z(D_8)

    return [tmp1,tmp2];
end;

#
# cls 7: E_8(b_4) [1,4,7,8,14,17,19,26]  this doesn't work in char2
#
handle7char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent;
    o:=Classes(orbs)[7];
    info:=infos[7];

    rep:=Representative(o);
    tmp1:=ConjugateByTorus(rep,4,-1);# h_1 --> h_4

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjugateByTorus(cent,4,-1);# h_1 --> h_4

    return [tmp1,tmp2];
end;

#
# cls 5: E_8(a_4) [1,4,6,8,17,18,19,21]  this doesn't work in char2
#
handle5char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent;
    o:=Classes(orbs)[5];
    info:=infos[5];

    rep:=Representative(o);
    tmp1:=ConjugateByTori(rep,[4,8],[-1,-1]);# h_4, h_8

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjugateByTori(cent,[4,8],[-1,-1]); # h_4, h_8

    return [tmp1,tmp2];
end;

#
# cls 4: E_8(a_3) [1,4,6,7,8,17,18,19]  this doesn't work in char2
#
handle4char5:=function(orbs,infos)
    local tmp1,tmp2,o,info,rep,cent;
    o:=Classes(orbs)[4];
    info:=infos[4];

    rep:=Representative(o);
    tmp1:=ConjugateByTorus(rep,4,-1);# h_1 --> h_4

    tmp2:=0;
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjugateByTorus(cent,4,-1);# h_1 --> h_4

    return [tmp1,tmp2];
end;


text:
# We will compute the sum exactly

# Computing an approximation of a rationnal (giving a string)
# Value is truncated toward zero
Approx := function(x, d)
	local neg, a, b, n, m, s;
	if x < 0 then
		x := -x;
		neg := true;
	else
		neg := false;
	fi;
	a := NumeratorRat(x);
	b := DenominatorRat(x);
	n := QuoInt(a, b);
	a := RemInt(a, b);
	m := 10^d;
	s := "";
	if neg then
		Append(s, "-");
	fi;
	Append(s, String(n));
	n := Size(s) + 1;
	Append(s, String(m + QuoInt(a*m, b)));
	s[n] := '.';
	return s;
end;

a := Sum([1 .. 1000], n -> 1/n^2);;
Approx(a, 10);
"1.6439345666"
# and pi^2/6 is 1.6449340668, truncated to ten digits