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

Declare(ISA_Bridge);

IsISABridge := (obj) -> IsRec(obj) and IsBound(obj._ISABridge) and obj._ISABridge=true;

#F VCvt(<n>, <ISA_Bridge>)
#F data type conversion from one ISA into another.
#F
# OBSOLETE, will be removed
Class(VCvt, TCast, rec(
    abbrevs := [],
    def := (n, brige) -> Perm((), n),

    dmn := self >> [ TArray(self.params[2].isa_from.gt(), self._spl().dims()[2]) ],
    rng := self >> [ TArray(self.params[2].isa_to.gt(),   self._spl().dims()[1]) ],

    toAMat := self >> self._spl().toAMat(),

    _spl := self >> self.params[2].splVCvt(self.params[1]),
));

Class(Cvt, SumsBase, Sym, rec(
    abbrevs := [  (isa_bridge) -> Checked(IsISABridge(isa_bridge), [isa_bridge]) ],
    def := (isa_bridge) -> Perm((), isa_bridge.granularity()),

    n         := self >> self.params[1].granularity(),
    isa_to    := self >> self.params[1].isa_to,
    isa_from  := self >> self.params[1].isa_from,

    dmn := self >> [ TArray(self.isa_from().t.base_t(), self.n()) ],
    rng := self >> [ TArray(self.isa_to().t.base_t(), self.n()) ],

    isPermutation := False,

    toAMat := self >> self.params[1].toAMat(),

    #doNotMeasure := true
));

Class(TCvt, Tagged_tSPL, rec(
    abbrevs := [
        (n, isa_to, isa_from) -> Checked(IsPosIntSym(n), IsISA(isa_to), IsISA(isa_from), [n, isa_to, isa_from, []]),
        (n, isa_to, isa_from, props) -> Checked(IsPosIntSym(n), IsISA(isa_to), IsISA(isa_from), IsList(props), [n, isa_to, isa_from, props]),
    ],

    #def := (n, isa_to, isa_from) -> Perm((), n),

    n         := self >> self.params[1],
    isa_to    := self >> self.params[2],
    isa_from  := self >> self.params[3],
    props     := self >> self.params[4],

    dmn := self >> [ TArray(self.isa_from().t.base_t(), self.n()) ],
    rng := self >> [ TArray(self.isa_to().t.base_t(), self.n()) ],

    toAMat := self >> I(self.n()).toAMat(),

    isReal := self >> true,

    isAuxNonTerminal := true,

    hashAs := self >> ObjId(self)(8*Lcm(self.isa_to().v, self.isa_from().v), self.isa_to(), self.isa_from(), self.props()).takeTA(self),
));

#
# ISA_Bridge: represents conversion from one SIMD_ISA to another.
# Each descendant implements particular conversion case and registered in
# the ISA_Bridge object by [source ISA, destination ISA] pair.
#

# <props> field
# "wraparound" | "saturation"
# "trunc" | "round" | "ceil" | "floor"

Class(ISA_Bridge, rec(
    _ISABridge := true,
    _table     := rec(),

    __call__ := meth(self)
        local inst;
        inst := WithBases(self);
        if not IsBound(self._table.(inst.isa_from.id())) then
            self._table.(inst.isa_from.id()) := rec();
        fi;
        if not IsBound(self._table.(inst.isa_from.id()).(inst.isa_to.id())) then
            self._table.(inst.isa_from.id()).(inst.isa_to.id()) := [];
        fi;
        Add(self._table.(inst.isa_from.id()).(inst.isa_to.id()), inst);
        return inst;
    end,

    add := (self, class) >> class(),

    printAvailableBridges := meth(self)
        local all;
        all := Flat(List(UserRecValues(self._table), r -> UserRecValues(r)));
        Sort(all, (a,b) -> let(a1 := StringPrint(a.isa_from), b1 := StringPrint(b.isa_from), When(a1=b1, StringPrint(a.isa_to)<StringPrint(b.isa_to), a1<b1)));
        Print(PrintPad("From", 27), PrintPad("To", 27), "Props\n", Replicate(60, '-'), "\n");
        DoForAll(all, e -> Print(PrintPad(StringPrint(e.isa_from), 26), " ", PrintPad(StringPrint(e.isa_to), 26), " ", e.props, "\n"));
    end,

    _applicable_cvt := (props, min_range, range_in, bridge) -> let(
        range_from  := bridge.isa_from.t.range(),
        range_to    := bridge.isa_to.t.range(),
        range_mid   := bridge.range(),

        range       := Cond( range_in=false, range_from, range_from * range_in),
        range_req   := Cond( range_in=false, min_range,  min_range  * range_in),
        range_mul   := range_to * range_mid * range_from,

        class_range := Cond( range.min >= range_to.min and range.max <= range_to.max and
                             range.min >= range_mid.min and range.max <= range_mid.max, [], ["saturation", "wraparound"]),
        class_prec  := Cond( range_from.eps = range_to.eps and range_from.eps = range_mid.eps,    [], ["round", "trunc", "ceil", "floor"]),
        bridge_class_range  := Intersection(class_range, props),
        bridge_class_prec   := Intersection(class_prec,  props),

        Intersection(bridge_class_range, bridge.props) = bridge_class_range
        and Intersection(bridge_class_prec, bridge.props) = bridge_class_prec
        and range_mul.min <= range_req.min and range_mul.max >= range_req.max
        and range_mul.eps <= range_req.eps
    ),

    # bridge: returns ISA_Bridge objects list which implement conversion from isa_from to isa_to

    find := (arg) >> let(
        self     := arg[1],
        isa_to   := arg[2],
        isa_from := arg[3],
        props    := arg[4],
        min_range:= Cond( Length(arg)>4, arg[5], isa_to.t.range() * isa_from.t.range() ),
        range_in := Cond( Length(arg)>5, arg[6], isa_from.t.range()),

        Cond( IsBound(self._table.(isa_from.id())) and IsBound(self._table.(isa_from.id()).(isa_to.id())),
            Filtered(self._table.(isa_from.id()).(isa_to.id()),
                e -> self._applicable_cvt(props, min_range, range_in, e)),
            [] )
    ),

    #######################
    # descendants should define the fields listed below.

    # isa_from: source SIMD_ISA
    isa_from := false,
    # isa_to: destination SIMD_ISA
    isa_to   := false,

    # minimal range
    range := (self) >> self.isa_to.t.range() * self.isa_to.t.range(),

    # NOTE: bridge conversion properties, explain
    # For example float-to-int conversion can be implemented with different rounding methods bla bla bla
    props       := [],

    # code(<y>, <x>, <opts>): called from codegen to get code that implements Cvt(...);
    code        := (self, y, x, opts) >> Error("not implemented"),

    # granularity: Cvt object height and width (which are the same).
    granularity := self >> Lcm(self.isa_to.v, self.isa_from.v),

    # toAMat: implements Cvt.toAMat()
    toAMat      := abstract(),

    # toSpl: conversion block formula NOTE: explain
    toSpl       := abstract(),

    ##############################################
    # helpers

    _x := (self, x, offs) >> vtref(self.isa_from.t, x, offs),
    _y := (self, y, offs) >> vtref(self.isa_to.t,   y, offs),

    _mkTL       := meth(self, isa, n, s)
        local t, hentry, sums;
        t := TL(n, s).withTags(isa.getTags());
        hentry := HashLookup(SIMD_ISA_DB.getHash(), t);
        if hentry=false then
            #return Error("SIMD_ISA_DB lookup failed, tried ", t);
            return Prm(L(n, s));
	else
            sums := _SPLRuleTree(hentry[1].ruletree).sums();
            sums := When(IsBound(sums.unroll), sums.unroll(), sums);
            return sums;
	fi;
    end,

));

# ISA_Bridge base class for bridges for which Cvt object is identity
Class(ISA_Bridge_I, ISA_Bridge, rec(
    toAMat      := (self) >> I(self.granularity()).toAMat(),
    toSpl       := (self) >> Cvt(self),
));

Class(ISA_Bridge_VvL, ISA_Bridge, rec(
    toAMat      := (self) >> L(self.granularity(), self.isa_from.v).toAMat(),
    toSpl       := (self) >> Cvt(self)*self._mkTL(self.isa_to, self.granularity(), div(self.granularity(), self.isa_from.v)),
));

Class(ISA_Bridge_VvR, ISA_Bridge, rec(
    toAMat      := (self) >> L(self.granularity(), self.isa_from.v).toAMat(),
    toSpl       := (self) >> self._mkTL(self.isa_to, self.granularity(), div(self.granularity(), self.isa_from.v)) * Cvt(self),
));

# plain conversion from one scalar data type to another,
# assuming isa_to type range includes isa_from type range
Class(ISA_Bridge_tcast, ISA_Bridge_I, rec(
    code  := (self, y, x, opts) >> assign( self._y(y,0), tcast(self.isa_to.t, self._x(x,0)) ),
));

# reinterpretation
Class(ISA_Bridge_wrap, ISA_Bridge_I, rec(
    range := (self) >> self.isa_to.t.range(),
    props := ["wraparound"],
    code  := (self, y, x, opts) >> assign( self._y(y,0), tcast(self.isa_to.t, self._x(x,0)) ),
));

_debug_print_cvt_chains := false;

DebugCVT := function(value)
    _debug_print_cvt_chains := value;
end;

# _cvt_build_chains( <TCvt non terminal> ) returns all data conversion chains available for given
#   TCvt object.
#

Declare(_cvt_build_chains_rec);
_cvt_build_chains_rec := function(chain, isa, min_range, range, props, available_isa)
    local ifrom, fr, r, b, dbg_print;
    ifrom := Last(chain).isa_to;
    if ifrom = isa then
        if _debug_print_cvt_chains then
            Print(GreenStr("CVT: "), ConcatSepList(chain, e -> e.__name__, " -> "), "\n");
        fi;
        return [chain];
    else
        r  := [];
        fr := ConcatList( available_isa, isa -> ISA_Bridge.find( isa, ifrom, props, min_range, range ));
        dbg_print := fr = [] and _debug_print_cvt_chains;
        for b in fr do
            Append(r, _cvt_build_chains_rec(chain :: [b], isa, min_range, b.range()*range, props, RemoveList(available_isa, b.isa_to)));
        od;
    fi;
    if dbg_print then
        Print(RedStr("Dead CVT: "), ConcatSepList(chain, e -> e.__name__, " -> "), "\n");
    fi;
    return r;
end;

_cvt_build_chains := function(t)
    local available_isa, r, range, fr, b, min_range;
    available_isa := Set(t.firstTag().params[1] :: [t.isa_from(), t.isa_to()]);#Set(RemoveList(t.firstTag().params[1], t.isa_from()) :: [t.isa_to()]);
    r         := [];
    range     := StripList(t.getA("r_in", [t.isa_from().t.range()]));
    min_range := t.isa_to().t.range() * t.isa_from().t.range();
    fr        := ConcatList( available_isa, isa -> ISA_Bridge.find( isa, t.isa_from(), t.props(), min_range, range ));
    for b in fr do
        Append(r, _cvt_build_chains_rec([b], t.isa_to(), min_range, b.range()*range, t.props(), RemoveList(available_isa, b.isa_to)));
    od;
    return Sort(r);
end;

NewRulesFor( TCvt, rec(
   TCvt_Term := rec(
       forTransposition := false,

       applicable := t -> t.firstTagIs(AMultiVec),
       freedoms   := t -> [ _cvt_build_chains(t) ],
       child      := (t, fr) -> [InfoNt(Reversed(fr[1]))],
       apply      := (self, t, C, nt) >> let(
           # merge adjasent loops when possible
           l := self._merge_loops(List(nt[1].params[1], b -> rec(
                   # NOTE: need ESReduce here to deal with symbolic expressions
                   #        this way can figure out divisibility at least some times
                   its := ESReduce(div(t.n(), b.granularity()), spiral.SpiralDefaults), bridge := b
                ))),

           Compose(List(l, e -> self._to_ISum(e)))
       ),

       _to_ISum := (self, r) >> let(
           ch   := Cond( IsBound(r.bridge), [r.bridge.toSpl()], List(r.children, e -> self._to_ISum(e))),
           cols := Cols(Last(ch)),
           rows := Rows(ch[1]),
           i    := Ind(r.its),
           ISum(i, Compose( [Scat(H(r.its*rows, rows, rows*i, 1))] :: ch :: [Gath(H(r.its*cols, cols, cols*i, 1))]))
       ),

       _merge_loops := meth(self, l)
           local c, e, new_ch, r;
           r := [];
           c := l[1];
           for e in Drop(l, 1) do
               if c.its <> 1 and _divides(c.its, e.its) then
                   new_ch := rec(its := ESReduce(div(e.its, c.its), spiral.SpiralDefaults), bridge := e.bridge);
                   if IsBound(c.bridge) then
                       # it's a single bridge, creating new node and making this one child
                       c := rec(
                           its      := c.its,
                           children := [rec(its := 1, bridge := c.bridge), new_ch]
                       );
                   else
                       Add(c.children, new_ch);
                   fi;
               elif e.its<>1 and _divides(e.its, c.its) then
                   new_ch := rec(its := ESReduce(div(c.its, e.its), spiral.SpiralDefaults));
                   if IsBound(c.bridge) then
                       new_ch.bridge := c.bridge;
                   else
                       new_ch.children := c.children;
                   fi;
                   c := rec(
                       its      := e.its,
                       children := [new_ch, rec(its := 1, bridge := e.bridge)]
                   );
               else
                   Add(r, c);
                   c := e;
               fi;
           od;

           Add(r, c);

           # recurse on children
           return List(r, e -> Cond(IsBound(e.children), CopyFields(e, rec(children := self._merge_loops(e.children))), e));
       end,

   ),
));


text:
# See section 7.5 of reference manual

# GAP has assertions levels. An assertion is tested if its level
# is less then the global level.

# Set global level
SetAssertionLevel(10);

a := 1;
Assert(20, a > 1, "a should be greater than one");
# nothing happens

a := 1;
Assert(4, a > 1, "a should be greater than one");
# error

# Show current global level
AssertionLevel();
# 10

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
Filtered([1 .. 10000], n -> Sum(DivisorsInt(n)) = 2*n);
# [ 6, 28, 496, 8128 ]


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

Declare(VBase);

Class(VBase, SumsBase, BaseContainer, rec(
    __call__ := (self, spl, vlen) >> Checked(IsSPL(spl), IsInt(vlen),
        SPL(WithBases(self, rec(_children := [spl],
                                vlen := vlen,
                                dimensions := spl.dimensions)))),
    rChildren := self >> [self._children[1], self.vlen],
    rSetChild := meth(self, n, newChild)
        if n=1 then self._children[1] := newChild; self.dimensions := newChild.dimensions;
        elif n=2 then self.vlen := newChild;
        else Error("<n> must be in [1..2]");
        fi;
    end,

    print := (self, i, is) >> Print(self.name, "(", self._children[1], ", ", self.vlen, ")"),

    transpose := self >> VBase(self.child(1).transpose(), self.vlen),

    sums := self >> self
));

ext: gi


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

#F RuleTreeOps.\=( <ruletree1>, <ruletree2> )
#F    returns true if <ruletree1> and <ruletree2> are equal, i.e., they
#F    represent the same breakdown strategy.  Otherwise, false is returned.
#F

RuleTreeOps := CantCopy(OperationsRecord("RuleTreeOps"));

IsRuleTree := T -> IsRec(T) and IsBound(T.isRuleTree) and T.isRuleTree;

#F Ruletrees
#F =========
#F
#F A ruletree is a record with the following mandatory fields:
#F
#F   isRuleTree  = "true"              # identifies ruletrees
#F   operations  = RuleTreeOps         # operations record
#F   node        = <spl>               # the spl expanded
#F   transposed  = false/true          # indicates whether <rule> or
#F                                       transpose - <rule> - transpose
#F   rule        = <rule>              # the rule chosen
#F   children    = <list>              # the children, a list of
#F                                       ruletrees/non-terminal spls
#F
#F Note that children can be a list of ruletrees as well as non-terminal spls.
#F The field children is empty iff the ruletree is a leaf.
#F In general, the .node is a non-terminal spl that is expanded by .rule.
#F
#F The following fields are optional:
#F   splOptions = [ <option1>, <option2>, .. ]
#F     valid <options> : "unrolled", "pcl"
#F

Class(RuleTreeClass, rec(
    spl  := self >> SPLRuleTree(self),
    measure := self >> SPLRuleTree(self).measure(),
    shortPrint := false,
    dims := self >> self.node.dims(),
    dmn  := self >> self.node.dmn(),
    rng  := self >> self.node.rng(),

    print := meth(self, i, is)
        local ch;

        if Length(self.children) = 0 then
           # print leaves in one line
            Print(self.rule, "( ");
            When (self.transposed, Print("\"T\","));
            self.node.print(i+is, is);
            When (IsBound(self.splOptions),
                Print(", ", self.splOptions, " )"),
                Print(" )"));
        else
            # print children on separate lines
            Print(self.rule, "( ",
                When(self.transposed, "\"T\", ", ""),
                When(not self.shortPrint, self.node.print(i+is, is), ""), ",\n", Blanks(i+is)
            );
            self.children[1].print(i+is, is);
            for ch in Drop(self.children, 1) do
                Print(",\n", Blanks(i+is));
                ch.print(i+is, is);
            od;

            When (IsBound(self.splOptions),
                Print("\n", Blanks(i), ", ", self.splOptions, " )"),
                Print(" )"));
        fi;
    end,

    # new with tSPL and opts.baseHash:
    # we cannot blindly transpose everything, as vector
    # permutations that are automatically generated
    # may not support transposition
    transpose := self >>
        When((not self.node.transposeSymmetric()) and self.node.isSymmetric(),
            self,
            RuleTree(
                self.rule,
                When(self.transposed, "", "T"),
                self.node.transpose(),
                List(self.children, x->x.transpose())
        )
    ),
    # RuleTreeRewriting enabled by default
    # one can disable it with "EnableRuleTreeRewriting(false)"
    rChildren := (self) >> [ self.node, self.children ],
    rSetChild := rSetChildFields("node", "children"),
    from_rChildren := (self, rch) >>
        RuleTree(self.rule, When(self.transposed, "T", ""), rch[1], rch[2])
));

RuleTreeOps.\= := function( T1, T2 )
    return
        IsRuleTree(T1) and IsRuleTree(T2) and
        T1.rule = T2.rule and
        T1.transposed = T2.transposed and
        IsIdenticalSPL(T1.node, T2.node) and
        T1.children = T2.children and
        (   (   IsBound(T1.splOptions) and
		        IsBound(T2.splOptions) and
                T1.splOptions = T2.splOptions
			) or
            (   not IsBound(T1.splOptions) and
			    not IsBound(T2.splOptions)
			)
        );
end;

#F RuleTreeNC( <rule>, <non-terminal spl>, <children> )
#F   returns the corresponding ruletree without any checking.

RuleTreeNC := function ( R, S, C )
    return WithBases(RuleTreeClass,
        rec(
            isRuleTree := true,
            operations := RuleTreeOps,
            node       := S,
            transposed := false,
            rule       := R,
            children   := C
		)
    );
end;

#F RuleTree(
#F   <rule>,
#F   [ "T" ,]
#F   <non-terminal spl>
#F   [, <children> ]
#F   [, <list-of-sploptions> ]
#F )
#F
#F   Creates a ruletree for <non-terminal spl>
#F   using <rule> with <children>. The <children> can be non-terminal spls
#F   or ruletrees. It is checked whether <children> can be derived
#F   from <non-terminal spl> using <rule>.
#F   If the string "T" is supplied, then instead of <rule>, the
#F   sequence transpose - <rule> - transpose is denoted. In this case
#F   <non-terminal spl> and <children> are also transposed. This construction
#F   avoids unnecessary recoding of rules that are just transposes of
#F   other rules. In this case, <rule> must have the field
#F   .forTransposition set to true.
#F
#F   The argument <list-of-sploptions> provides options used for spl compilation.
#F   See top of this file for valid options; examples are "unrolled" and "pcl".
#F   <children> might be empty (default) in which case the ruletree
#F   returned is a leaf.
#F
RuleTree := function ( arg )
    local a, R, t, S, C, u, S1, C1, RT;
    S := false;
    t := false;
    u := [ ];
    C := [ ];
    R := false;

    for a in arg do
        if IsString(a) and a = "T" then
	        t := true;
        elif IsList(a) and ForAll(a, IsString) then
	        u := a;
        elif IsSPL(a) then
			if S=false then
				S := a;
			else
				Add(C, a);
			fi;
		elif IsRuleTree(a) then
			Add(C, a);
		elif IsList(a) then
			Append(C, a);
		elif IsRule(a) and R=false then
			R := a;
		else Error(
			"usage: \n",
			"  RuleTree( \n",
			"    <rule>, [ \"T\" ,] <non-terminal spl> \n",
			"    [, <children> ] [, <list-of-sploptions> ] )");

		fi;
    od;

    Constraint(IsRule(R));
    Constraint(IsSPL(S));

    if not IsList(u) and ForAll(u, IsString) then
        Error("<u> must be a list of spl options");
    elif not (IsList(C) and ForAll(C, c -> IsSPL(c) or IsRuleTree(c))) then
        Error("<C> must be a list of ruletrees or non-terminal spls");
    fi;

    RT := RuleTreeNC(R, S, C);
    if u <> [ ]
	    then RT.splOptions := u;
	fi;
    return RT;
end;

#F ExtractChildrenSPL ( <spl> )
#F   returns the list of non-terminal spls contained in <spl>.
#F   The order is left first - depth first (non-terminals are
#F   always leaves.
#F
ExtractChildrenSPL := S -> Collect(S, @.cond(IsNonTerminal));

#F CopyRuleTree( <ruletree> )
#F   returns a copy of <ruletree>. All fields apart from .children
#F   are copied by reference.
#F
CopyRuleTree := function( orig )
    local new;
    new := ShallowCopy(orig);
    if IsRuleTree(new) then
	    new.children := List(orig.children, CopyRuleTree);
	fi;
    return new;
end;

_matchChildren := (C1, C2) ->
    Length(C1)=Length(C2) and
    ForAll([1..Length(C1)],
       j -> IsHashIdenticalSPL(HashAsSPL(C1[j]), HashAsSPL(C2[j])) or
            (ObjId(C1[j])=InfoNt and ObjId(C2[j])=InfoNt));

#F ApplyRuleTreeSPL( <ruletree>, <spl>, <opts> )
#F   creates a ruletree for <spl> using rules from the <ruletree>
#F   (if possible)
#F   Uses: opts.breakdownRules (for backwards compatibility of R.switch only)
ApplyRuleTreeSPL := function ( rtree, spl, opts )
    local i, R, C, C1, CC, children, L;
    Constraint(IsRuleTree(rtree));
    Constraint(IsSPL(spl));

    if rtree.transposed then
        return ApplyRuleTreeSPL(rtree.transpose(), spl.transpose(), opts).transpose();
    fi;

    R := rtree.rule;
    children:=[];

    # if <spl> was not compatible with <ruletree>.node to start with
    # then at some point there will be a split in <ruletree> that is not
    # applicable to the generated ruletree.
    if not IsApplicableRule( R, spl, opts.breakdownRules ) then
        Error("<rtree> can not be expanded from <spl>");
    fi;

    # determine the right set of children from allChildren to be used
    # in expansion
    C1 := _allChildren(R, rtree.node, opts);
    C := List([1..Length(rtree.children)], c -> rtree.children[c].node);
    i := PositionProperty([1..Length(C1)], j -> _matchChildren(C1[j], C));

    # DEBUG HACK: AppendTo("applytree", R, "\n", C, "\n", C1, "\n---------\n");
    if i = false then
        Error("<ruletree> can not be expanded from <spl>, could not find matching children");
	fi;

    # determine all children for <spl> and pick the right set
    L :=  _allChildren(R, spl, opts);
    if i > Length(L) then
        Error("<ruletree> can not be expanded from <spl>");
    else
        C1 := L[i];
    fi;

    # recurse and return
    children := List([1..Length(rtree.children)], i -> ApplyRuleTreeSPL(rtree.children[i], C1[i], opts));

    return RuleTreeNC(R, spl, children);
end;

#F Printing RuleTrees
#F ------------------
#F

#F RuleTreeOps.Print( <ruletree> [, <indent> , <indentStep> ] )
#F   prints <ruletree> with <indent>. Further indenting is done
#F   in steps of size <indentStep>. The default is
#F   indent = 0, indentStep = 2.
#F
RuleTreeOps.Print := function ( arg )
    local T, indent, indentStep;
	if Length(arg) = 1 then
        T          := arg[1];
        indent     := 0;
        indentStep := 2;
    elif Length(arg) = 3 then
        T          := arg[1];
        indent     := arg[2];
        indentStep := arg[3];
    else
        Error("usage: RuleTreeOps.Print( <ruletree> [, <indent> , <indentStep> ] )");
    fi;
    Constraint(IsRuleTree(T));
    Constraint(IsPosInt0(indent));
    Constraint(IsPosInt0(indentStep));
    T.print(indent, indentStep);
end;

PrintRuleTreeCustom := function (T, i, is, spl_print, ruletree_print)
    local ch;
    Constraint(IsRuleTree(T) or IsSPL(T));
    Constraint(IsPosInt0(i));
    Constraint(IsPosInt0(is));

    if IsSPL(T) then
        spl_print(T);
    else # T is a ruletree
        Print(Blanks(i), ruletree_print(T), "\n");
        for ch in T.children do
            PrintRuleTreeCustom(ch, i+is, is, spl_print, ruletree_print);
        od;
    fi;
end;

#F PrintNodesRuleTree(
#F   <ruletree/non-terminal spl> [, <indent> , <indentStep> ]
#F )
#F   pretty prints <ruletree> by displaying only the nodes (non-terminals)
#F   and its parameters. The tree structure is expressed by indentation.
#F   Non-terminals are being marked by (nt).
#F

PrintNodesRuleTree := function ( arg )
    local T, i, is;
    if Length(arg) = 1 then
        T := arg[1];
		i := 0;
		is := 2;
    elif Length(arg) = 3 then
        T := arg[1];
		i := arg[2];
		is := arg[3];
    else
        Error("usage: PrintNodesRuleTree( <ruletree> [, <indent> , <indentStep> ] )");
    fi;
    PrintRuleTreeCustom(T, i, is,
        T -> Print(T, " (nt)"),
        T -> Print(T.node));
end;

#F PrintRulesRuleTree(
#F   <ruletree/non-terminal spl> [, <indent> , <indentStep> ]
#F )
#F   pretty prints <ruletree> by displaying only the rules and
#F   the node sizes in parentheses. Transposed rules are marked
#F   by ^T. The tree structure is expressed by indentation.
#F   Non-terminals are denoted by nt.
#F
PrintRulesRuleTree := function(arg)
    local T, i, is;
    if Length(arg) = 1 then
        T := arg[1]; i := 0; is := 2;
    elif Length(arg) = 3 then
        T := arg[1]; i := arg[2]; is := arg[3];
    else
        Error("usage: PrintRulesRuleTree( <ruletree> [, <indent> , <indentStep> ] )");
    fi;
    PrintRuleTreeCustom(T, i, is,
    T -> Print("nt (", T.dimensions[1], ")"),
    T -> Print(T.rule, When(T.transposed," ^ T",""), "(", T.node.dimensions[1], ")"));
end;

#F PrettyPrintRuleTree(
#F   <ruletree/non-terminal spl> [, <whiteIndent> [, <linedIndent> ] ]
#F )
#F   pretty prints <ruletree> by displaying only the nodes (non-terminals),
#F   their parameters, and the rules. The tree structure is expressed by
#F   indentation.  Non-terminals are being marked by (nt).
#F

PrettyPrintRuleTree := function ( arg )
    local T, chIndent, whiteIndent, linedIndent, oldLinedIndent, newline, params, i;

    # new line plus indents
    newline := function ( )
    local i;

    Print("\n");
    for i in [1..whiteIndent] do
        Print(" ");
    od;
    if not linedIndent = "" then
        Print( linedIndent, "--" );
    fi;
    end;

    # decode arg
    if Length(arg) = 1 then
        T           := arg[1];
        whiteIndent := 0;
        linedIndent := "";
    elif Length(arg) = 2 then
        T           := arg[1];
        whiteIndent := arg[2];
        linedIndent := "";
    elif Length(arg) = 3 then
        T           := arg[1];
        whiteIndent := arg[2];
        linedIndent := arg[3];
    else
        Error("usage: PrintNodesRuleTree( <ruletree> [, <indent> , <indentStep> ] )");
    fi;

    # check arguments
    if not ( IsRuleTree(T) or IsSPL(T) and T.type = "nonTerminal" ) then
        Error("<T> must be a ruletree or a non-terminal spl");
    fi;
    if not IsInt(whiteIndent) and whiteIndent >= 0 then
        Error("<whiteIndent> must be pos-int");
    fi;
    if not IsString(linedIndent) then
        Error("<linedIndent> must be a string");
    fi;

    chIndent := whiteIndent+Length(linedIndent)+4;
    # spl case
    if IsSPL(T) then
        T.print(chIndent, 2); Print(" (nt)");
        return;
    fi;

    T.node.print(chIndent, 2);
    Print("     {" );
    Print(T.rule);
    if T.transposed then
        Print(" ^ T");
    fi;
    if IsBound(T.splOptions) then
        Print(", ", T.splOptions);
    fi;
    Print( "}" );

    if Length(T.children) > 0 then
        if linedIndent = "" then
            oldLinedIndent := linedIndent;
        else
            oldLinedIndent := Concatenation( linedIndent, "  " );
        fi;
        linedIndent := Concatenation( oldLinedIndent, " |" );
        for i in [1..Length(T.children)-1] do
            newline();
            PrettyPrintRuleTree(T.children[i], whiteIndent, linedIndent);
        od;
        linedIndent := Concatenation( oldLinedIndent, " \`" );
        newline();
        linedIndent := Concatenation( oldLinedIndent, "  " );
        PrettyPrintRuleTree(T.children[Length(T.children)],
                            whiteIndent, linedIndent);
        linedIndent := oldLinedIndent;
    fi;
    if whiteIndent = 0 and linedIndent = "" then
        newline();
    fi;
end;

#F Converting RuleTrees
#F --------------------
#F

#F AMatRuleTree( <ruletree> ) - returns an amat corresponding to ruletree.
#F
AMatRuleTree := T -> AMatSPL(SPLRuleTree(T));

#F MatRuleTree( <ruletree> ) -  returns the matrix corresponding to ruletree.
#F
MatRuleTree := T -> MatSPL(SPLRuleTree(T));

#F Creating Random RuleTrees
#F -------------------------
#F

Declare(SemiRandomRuleTree, _SemiRandomRuleTree);

#F RandomRuleTree( <spl> , <opts>)
#F    returns a random rule tree for <spl>
#F    or 'false' if no breakdown rule combination leads to a
#F    fully expanded ruletree.
#F
#F   Uses: opts.baseHashes
#F         opts.breakdownRules
#F   Optional:  (see Doc(_allChildren) ) for usage
#F         opts.restrictSplit
#F         opts.restrictSplitSize
#F
RandomRuleTree := (spl, opts) -> _SemiRandomRuleTree(spl, false, s->false, opts);

RandomRuleTreeCutoff := (spl, cutoff_func, opts) -> _SemiRandomRuleTree(spl, false, cutoff_func, opts);

RandomRuleTreeDP := (t, opts) -> let(res := spiral.search.DP(t, rec(measureFunction := (R,O)->Random([1..10000]), verbosity := 0), opts), When(res = [], false, res[1].ruletree));

#F SemiRandomRuleTree( <spl>, <top_rule>, opts )
#F    returns a random rule tree for <spl> with fixed top-level rule
#F    or 'false' if no breakdown rule combination leads to a
#F    fully expanded ruletree.
#F
#F   Uses: opts.baseHashes
#F         opts.breakdownRules
#F   Optional:  (see Doc(_allChildren) ) for usage
#F         opts.restrictSplit
#F         opts.restrictSplitSize
#F
SemiRandomRuleTree := (spl, top_rule, opts) -> _SemiRandomRuleTree(spl, top_rule, s->false, opts);

_SemiRandomRuleTree := function(spl, top_rule, cutoff_func, opts)
    local ch, i, R, rules, rule, rand, candidates, ch, children, ch_candidates;

    if not (IsSPL(spl)) then
        Error("usage: SemiRandomRuleTree( <spl>, <top_rule>, <cutoff>, <opts> )");
    elif cutoff_func(spl) then
	    return spl;
    elif MultiHashLookup(opts.baseHashes, spl) <> false then
        return MultiHashLookup(opts.baseHashes, spl)[1].ruletree;
    else
        if top_rule <> false then
            rules := [top_rule];
        else
            rules  := AllApplicableRulesDirect(spl, opts.breakdownRules);
            if rules = [ ] then
                return false;
            fi;
        fi;

        candidates := Set([1..(Length(rules))]);
        while candidates <> [] do
            rand := RandomList(candidates);
            RemoveSet(candidates, rand);

            if rand <= Length(rules) then
                rule := rules[rand];
                ch_candidates := _allChildren(rule, spl,opts);
                ch_candidates := Permuted(ch_candidates, Random(SymmetricGroup(Length(ch_candidates))));
                for children in ch_candidates do
                    children := List( children, s -> _SemiRandomRuleTree(s, false, cutoff_func, opts) );
                    if not ForAny(children, c -> IsBool(c) and c=false) then
                        return RuleTreeNC( rule, spl, children );
                    fi;
                od;
            fi;
        od;
        return false;
   fi;
end;

#F Creating all RuleTrees
#F ----------------------
#F
#ExpandSPLRules := function( S, ruleset )
ExpandSPLRules := function( arg )
    local i, L, R, S, ruleset, opts;
    S := arg[1];
    ruleset := arg[2];
    Constraint(IsSPL(S));

    if IsBound(arg[3]) then
	    opts := arg[3];
    else
	    opts := rec();
	fi;

    L  := [ ];
    for R in AllApplicableRulesDirect(S, ruleset) do
        Append(L, List(_allChildren(R,S,opts), c -> RuleTreeNC(R, S, c)));
    od;

    return L;
end;

#F ExpandSPL( <spl>, <opts> )
#F   performs one expansion of <spl> in all possible ways,
#F   and returns the list of the derived ruletrees.
#F   If there are no breakdown rules for <spl> (i.e., not a non-terminal)
#F   then the default base cases rule @_Base is used, which means that
#F   the recursion will continue on its children
#F
#F   Uses: opts.baseHashes
#F         opts.breakdownRules
#F
ExpandSPL := (S, opts) -> ExpandSPLRules(S, opts.breakdownRules, opts);

_AllRuleTrees := function ( S, cutoff, memohash, opts )
    local L, L1, T, Lc, i, cs, T1, lkup;

    # check cutoff, baseHashes, and our memoization hash
    if cutoff(S) then
		return [S];
	fi;
    lkup := MultiHashLookup(opts.baseHashes, S);
    if lkup <> false then
		return [lkup[1].ruletree];
	fi;

    lkup := HashLookup(memohash, S);
    if lkup <> false
		then return lkup;
	fi;

    L := ExpandSPL(S, opts);
    if ForAll(L, r -> r.children = [ ]) then
		HashAdd(memohash, S, L);
		return L;
    fi;

    L1 := [ ];
    for T in L do
        Lc := [ ];
		for i in [1..Length(T.children)] do
            Lc[i] := _AllRuleTrees(T.children[i], cutoff, memohash, opts);
		od;
		for cs in Cartesian(Lc) do
            T1 := ShallowCopy(T);
			T1.children := cs;
			Add(L1, T1);
		od;
    od;

    HashAdd(memohash, S, L1);
    return L1;
end;

#F AllRuleTrees( <spl>, <opts> )
#F   expands <spl> in all possible ways and returns the obtained list
#F   of fully expanded ruletrees (ruletrees containing no non-terminals).
#F   Note that subtrees are "by reference", i.e., modifying one subtree
#F   alters the same subtree in some other expansions.
#F
#F   Uses: opts.baseHashes
#F         opts.breakdownRules
#F   Optional:  (see Doc(_allChildren) ) for usage
#F         opts.restrictSplit
#F         opts.restrictSplitSize
#F
AllRuleTrees := (S, opts) -> _AllRuleTrees(S, x->false, HashTableSPL(), opts);

#F AllRuleTreesCutoff( <spl>, <cutoff_func>, <opts> )
#F   expands <spl> in all possible ways and returns the obtained list
#F   of fully expanded ruletrees (ruletrees containing no non-terminals).
#F   Note that subtrees are "by reference", i.e., modifying one subtree
#F   alters the same subtree in some other expansions.
#F   <cutoff_func> is of the form e -> e = <cutoff>
#F
#F   Uses: opts.baseHashes
#F         opts.breakdownRules
#F
AllRuleTreesCutoff := (S, cutoff_func, opts) -> _AllRuleTrees(S, cutoff_func, HashTableSPL(), opts);

_NofRuleTrees := function ( S, cutoff, memohash, opts, level, trace )
    local p, Cs, n, lkup, indentstr, i;
    Constraint(IsSPL(S));
    Constraint(IsRec(opts));

    if trace then
        indentstr := "";
        for i in [2..level] do
            indentstr := Concat(indentstr, ".  ");
        od;
        indentstr := Concat(indentstr, String(level), " ");
        Print(indentstr, ">> _NofRuleTrees(", S, ")\n");
    fi;

    # check cutoff, baseHashes, and our memoization hash
    if cutoff(S) then
        if trace then
            Print(indentstr, "<< (cutoff) 1\n");
        fi;
        return 1;
    fi;
    if MultiHashLookup(opts.baseHashes, S) <> false then
        if trace then
            Print(indentstr, "<< (MultiHashLookup) 1\n");
        fi;
        return 1;
    fi;
    lkup := HashLookup(memohash, S);
    if lkup <> false then
        if trace then
            Print(indentstr, "<< hashed: ", lkup, "\n");
        fi;
        return lkup;
    fi;

    # first apply rules as they are
    Cs := List(AllApplicableRulesDirect(S, opts.breakdownRules), r -> _allChildren(r,S, opts));
	if trace then
        Print(indentstr, "   Cs: ", Cs, "\n");
    fi;

    if Cs = [ [ ] ] then
        n := 1;
    else
        n := Sum(Cs, c -> Sum(c, l -> Product(l, x -> _NofRuleTrees(x, cutoff, memohash, opts, level+1, trace))));
    fi;

    HashAdd(memohash, S, n);

    if trace then
        Print(indentstr, "<< (", S, ") ", n, "\n");
    fi;

    return n;
end;

#F NofRuleTrees ( <spl>, <opts> )
#F   returns the number of ruletrees for <spl>.
#F
#F   Uses: opts.baseHashes
#F         opts.breakdownRules
#F   Optional:  (see Doc(_allChildren) ) for usage
#F         opts.restrictSplit
#F         opts.restrictSplitSize
#F
#F   if opts.verbosity > 3 prints trace info

NofRuleTrees := (spl, opts) -> _NofRuleTrees(spl, x->false, HashTableSPL(), opts, 1, IsBound(opts.verbosity) and (opts.verbosity > 3));

NofRuleTreesCutoff := (spl, cutoff_func, opts) -> _NofRuleTrees(spl, cutoff_func, HashTableSPL(), opts, 1, IsBound(opts.verbosity) and (opts.verbosity > 3));

#F RulesInRuleTree( <ruletree> )
#F   return the set of rules contained in <ruletree>.
#F
RulesInRuleTree := function ( T )
    local C;
    Constraint(IsRuleTree(T));

    # base case
    if T.children = [ ] then
	    return Set([T.rule]);
	fi;

    # recurse with ruletree children
    C := Filtered(T.children, IsRuleTree);
    C := Set(Concatenation(List(C, RulesInRuleTree)));
    AddSet(C, T.rule);
    return C;
end;

#F RulesInRuleTreeAll( <ruletree> )
#F   return the histogram of all rules contained in <ruletree> as set
#F   of pairs [rule, number].
#F
RulesInRuleTreeAll := function ( T )
    local C, C1, c, i;
    Constraint(IsRuleTree(T));

    # base case
    if T.children = [ ] then
	    return Set([ [T.rule, 1] ]);
	fi;

    # recurse with ruletree children
    C := Filtered(T.children, IsRuleTree);
    C := Concatenation(List(C, RulesInRuleTreeAll));

    # merge
    C1 := [ ];
    for c in C do
        i := PositionProperty(C1, p -> c[1] = p[1]);
        if i = false then
            AddSet(C1, c);
        else
            C1[i][2] := C1[i][2] + c[2];
        fi;
    od;

    i := PositionProperty(C1, p -> p[1] = T.rule);
    if i = false then
        AddSet(C1, [T.rule, 1]);
    else
        C1[i][2] := C1[i][2] + 1;
    fi;
        return C1;
end;

ApplyRuleTreeStep := function ( rt, recurse )
    local nt, C, Nonterms, S, c;
    if IsNonTerminal(rt) or IsSPL(rt) then
	    return rt;
	fi;

    Constraint(IsRuleTree(rt));

    C := List(rt.children, recurse);
    nt := rt.node;
    Nonterms := List(rt.children, c -> When(IsRuleTree(c), c.node, c));

    if rt.transposed then
        nt := TransposedSPL(nt);
        C := List(C, TransposedSPL);
        Nonterms := List(Nonterms, TransposedSPL);
    fi;

    S := _apply(rt.rule, nt, C, Nonterms);
    if rt.transposed then
	    S := TransposedSPL(S);
	fi;

    trace_log.addTreeExpansion(rt.rule,nt,S,rt.children,CopyFields(rt.node, rec(params := C)), var);

    S.root := rt.node;
    return S;
end;

SPLRuleTreeStep := rt -> ApplyRuleTreeStep(rt, c -> RecursStep(
    When(IsRuleTree(c), RTWrap(c), c)));

_SPLRuleTree := rt -> ApplyRuleTreeStep(rt, _SPLRuleTree);

SPLRuleTree := function(rt)
    local spl, t, tag;
    trace_log.beginStage("tSPL","SPL", rt.node);
    spl := ApplyRuleTreeStep(rt, _SPLRuleTree);

    # tags may inject container objects
    t := rt.node;
    if IsBound(t.getTags) then
        for tag in Reversed(t.getTags()) do
            if IsBound(tag.container)
			    then spl := tag.container(spl);
			fi;
        od;
    fi;
    trace_log.endStage("tSPL","SPL", spl);
    return spl;
end;

# forward declaration for function that does all the work.
Declare(_ruleTreeN);

#F RandomRuleTrees( <spl>, <num>, <opts>)
#F    returns <num> random rule trees for <spl>
#F    uniformly distributed amongst the potential trees
#F    or 'false' if no breakdown rule combination leads to a
#F    fully expanded ruletree.
#F
#F   Uses: opts.baseHashes
#F         opts.breakdownRules
#F   Optional:  (see Doc(_allChildren) ) for usage
#F
RandomRuleTrees := function(spl, num, opts)
    local n, rand, memohash;

    memohash := HashTableSPL();

    n := _NofRuleTrees(spl, s -> false, memohash, opts, 1, false);

    if num >= n then
        rand := [1..n];
    else
        rand := [];
        while Length(rand) <> num do
            rand := Set(Concat(rand, [RandomList([1..n])]));
        od;
    fi;

    return List(rand, e -> _ruleTreeN(spl, e-1, memohash, opts));
end;

# helper functions/arrays for _ruleTreeN and _ruleTree1

_getNumTrees := (s, memohash, opts) -> let(
    n := HashLookup(memohash, s),
    When(n = false,
        _NofRuleTrees(s, e -> false, memohash, opts, 1, false),
        n
    )
);
#F RuleTreeN(<spl>, <num>, <opts>)
#F   returns ruletree of index <num>
#F
#F   Uses: opts.basHashes, opts.breakdownRules.
#F
RuleTreeN := function(spl, num, opts)
    local memohash, n, r;

    memohash := HashTableSPL();
    n := _NofRuleTrees(spl, s -> false, memohash, opts, 1, false);

    if num in [1..n] then
        r := _ruleTreeN(spl, num-1, memohash, opts);

    else
        r := false;
    fi;

    return r;
end;

# fold list inside out, favoring the right hand side

_foldedMidFirst := function(origlist)
	local len, midpt, front, back, newlist;

	len := Length(origlist);
	if len < 2 then
		return origlist;
	fi;
	midpt := Int(len/2);
	front := Reversed(Sublist(origlist, [1..midpt]));
	back  := Sublist(origlist, [midpt+1..len]);

	newlist := [];

	if IsOddInt(len) then
		Append(newlist, [back[1]]);
		back := Drop(back,1);
	fi;

	while (front <> []) or (back <> []) do
		if back <> [] then
			Append(newlist, [back[1]]);
			back := Drop(back,1);
		fi;
		if front <> [] then
			Append(newlist, [front[1]]);
			front := Drop(front, 1);
		fi;
	od;

	return newlist;
end;

Declare(_ruleTreeMid);

#F RuleTreeMid(spl, opts)
#F
#F   Similar to RuleTree1, but tries to grab the middle of each range of choices
#F

RuleTreeMid := (spl, opts) -> _ruleTreeMid(spl, HashTableSPL(), opts);

_ruleTreeMid := function(spl, memohash, opts)
    local h, R, r, C, c, ch, npc, idx, idxlist, ridx, ridxlist;

    h := MultiHashLookup(opts.baseHashes, spl);

    if h <> false then
        return h[1].ruletree;
    fi;

    R := AllApplicableRulesDirect(spl, opts.breakdownRules);
	ridxlist := [1..Length(R)];
	if Length(ridxlist) > 1 then
		if Length(_allChildren(R[1], spl, opts)) < 2 then
			ridxlist := _foldedMidFirst(ridxlist);
		fi;
	fi;
    for ridx in ridxlist do
		r := R[ridx];
		C := _allChildren(r, spl, opts);
		idxlist := _foldedMidFirst([1..Length(C)]);
		for idx in idxlist do
			c := C[idx];
			npc := Product(c, e -> _getNumTrees(e, memohash, opts));
			if npc > 0 then
				ch := List(c, e -> _ruleTreeMid(e, memohash, opts));
				if not ForAny(ch, e -> IsBool(e) and e=false) then
					return RuleTreeNC(r, spl, ch);
				fi;
			fi;

        od;
    od;
	# no tree, return original spl
    return spl;
end;

Declare(_ruleTree1);

#F RuleTree1(spl, opts)
#F
#F   returns the first ruletree. this is often useful for
#F   testing whether a particular rule is firing, and you
#F   want to avoid the randomness of RandomRuleTree
#F
#F   NOTE: this function does not check if a tree can
#F   actually just be created, it just tries to return
#F   the first one. RuleTreeN(spl, 1, opts) does the
#F   same thing but first checks.
#F

RuleTree1 := (spl, opts) -> _ruleTree1(spl, HashTableSPL(), opts);

_ruleTree1 := function(spl, memohash, opts)
    local h, R, r, i, c, ch, npc;

    h := MultiHashLookup(opts.baseHashes, spl);

    if h <> false then
        return h[1].ruletree;
    fi;

    R := AllApplicableRulesDirect(spl, opts.breakdownRules);
    for r in R do
		for c in _allChildren(r, spl, opts) do
			npc := Product(c, e -> _getNumTrees(e, memohash, opts));
			if npc > 0 then
				ch := List(c, e -> _ruleTree1(e, memohash, opts));
				if not ForAny(ch, e -> IsBool(e) and e=false) then
					return RuleTreeNC(r, spl, ch);
				fi;
			fi;
        od;
    od;
	# no tree, return original spl
    return spl;
end;

#F _ruleTreeN(<spl>, <num>, <memohash>, <opts>)
#F   internal function which returns a ruletree for
#F   <spl> with index given by <num>. This is a recursive
#F   routine. Note that <num> is 0-based, unlike everything
#F   else in SPIRAL, which is 1-based.
#F
_ruleTreeN := function(spl, num, memohash, opts)
    local R, rtNC, h, r, c, npc, off, ch, i, n, p, offset;

    # lookup stuff like SSE base cases.
    h := MultiHashLookup(opts.baseHashes, spl);

    if h <> false then
       return h[1].ruletree;
    fi;

    # get list of normal rules
    R := AllApplicableRulesDirect(spl, opts.breakdownRules);
    offset := 0;

    for r in R do
        for c in _allChildren(r, spl, opts) do
            # figure out which group of children we should be looking at.
            npc := Product(c, e -> _getNumTrees(e, memohash, opts));

            # is this the correct group?
            if num < offset + npc then
                off := npc;
                ch := [];

                # step through the kids, and push down the ruletree number
                for i in [1..Length(c)] do
                    n := _getNumTrees(c[i], memohash, opts);

                    off := When(n > 0, off / n, off);

                    p := When(n <= 0 or num <= offset,
                        0,
                        RemInt(QuoInt(num - offset, off), n)
                    );

                    # all the kids get evaluated
                    Add(ch, _ruleTreeN(c[i], p, memohash, opts));
                od;

                # no kids should be false, but just in case.
                if not ForAny(ch, e -> IsBool(e) and e=false) then
                    return RuleTreeNC(r, spl, ch);
                fi;
            fi;

            # fixup offset
            offset := offset + npc;
        od;
    od;

	# no tree, return original spl
    return spl;
end;

Declare(_ruleTreeGetN);

#F RuleTreeGetN(<ruletree>, <opts>)
#F given a ruletree, determine the index
#F
#F uses:
#F      opts.breakdownRules
#F      opts.baseHashes
#F

RuleTreeGetN := function(ruletree, opts)
    local memohash, r;

    memohash := HashTableSPL();

    # prime the memohash.
    _NofRuleTrees(ruletree.node, s -> false, memohash, opts, 1, false);

    r := _ruleTreeGetN(ruletree, memohash, opts, 0);

    # ruleTreeGetN may return false if we couldn't match the tree in the passed
    # in breakdownRules/baseHashes.
    return When(IsInt(r), r + 1, false);
end;

_ruleTreeGetN := function(rt, memohash, opts, offset)
    local h, R, r, c, ch, off, n;

    # is the ruletree something from the baseHashes?
    h := MultiHashLookup(opts.baseHashes, rt.node);
    if IsList(h) and rt = h[1].ruletree then
        return offset;
    fi;

    R := AllApplicableRulesDirect(rt.node, opts.breakdownRules);

    for r in R do
        # does rule match?
        if r = rt.rule then
            for c in _allChildren(r, rt.node, opts) do
                # does child match?
                if c = List(rt.children, e -> e.node) then
                    off := Product(rt.children, e -> _getNumTrees(e.node, memohash, opts));
                    for ch in rt.children do
                        n := _getNumTrees(ch.node, memohash, opts);
                        off := When(n > 0, off / n, off);
                        offset := offset + (off * _ruleTreeGetN(ch, memohash, opts, 0));
                    od;
                    return offset;
                # child doesn't match, increase offset
                else
                    offset := offset + Product(c, e -> _getNumTrees(e, memohash, opts));
                fi;
            od;
        # rule doesn't match, so increase offset by all possible expansions we didn't use
        else
            offset := offset
                + Sum(_allChildren(r, rt.node, opts), C ->
                    Product(C, c -> _getNumTrees(c, memohash, opts))
                );
        fi;
    od;

    return offset;
end;


text:
# The tetrahedron with v=4, k=3, lm=2.
BindGlobal("TetrahedronGraph", CompleteGraph(SymmetricGroup(4)));

# The octahedron with v=6, k=4, lm=2, mu=4.
BindGlobal("OctahedronGraph", CocktailPartyGraph(3));

# The Petersen graph with v=10, k=3, lm=0, mu=1.
BindGlobal("PetersenGraph", OddGraph(2));

# The Clebsch graph with v=16, k=10, lm=6, mu=6.
BindGlobal("ClebschGraph", HalvedCubeGraph(5));

# The Schlaefli graph with v=27, k=16, lm=10, mu=8.
BindGlobal("SchlaefliGraph",
    Graph(DirectProduct(SymmetricGroup(6), SymmetricGroup(2)),
        [[-3,-3,1,1,1,1,1,1], [3,-1,-1,-1,-1,-1,-1,3]],
        OnRoots, RootAdjacency));

# The Hoffman-Singleton graph with v=50, k=7, lm=0, mu=1.
BindGlobal("HoffmanSingletonGraph", List([function()
        local G, dp, p1, p2, p3, p4;
        G := Group((1, 2, 3, 4, 5));
        dp := DirectProduct(G, G, G, Group((1, 2)));
        p1 := Projection(dp, 1);
        p2 := Projection(dp, 2);
        p3 := Projection(dp, 3);
        p4 := Projection(dp, 4);
        return Graph(dp, Cartesian(GF(5), GF(5), [1, 2]),
            function(x, g)
                local g4, j, s, u, v;
                u := 5^Image(p2, g) * Z(5)^0;
                v := 5^Image(p3, g) * Z(5)^0;
                g4 := Image(p4, g);
                s := x[3]^2;
                j := [1, 2*x[3]];
                return [1^g4*(x[1]+x[2]*(u-s*v)+s*(u^2-v^2)/2)+5^Image(p1, g),
                    j[1^g4]*(x[2]+s*u+v), x[3]^g4];
            end, function(x, y)
                return (x{[2,3]} = y{[2,3]} and y[1] in [x[1]-x[3], x[1]+x[3]])
                    or x[3] <> y[3] and y[1] = x[1] + x[2]*y[2]*x[3]^2;
            end, true);
    end])[1]());

# The Gewirtz graph with v=56, k=10, lm=0, mu=2.
BindGlobal("GewirtzGraph", Graph(MathieuGroup(21), [[1,2,3,7,10,20]],
                                    OnSets, DisjointSets));

# The strongly regular Witt graph with v=77, k=16, lm=0, mu=4.
BindGlobal("WittStronglyRegularGraph", Graph(MathieuGroup(22),
                                    [[1,2,3,7,10,20]], OnSets, DisjointSets));

# The graph with v=210, k=99, lm=48, mu=45 constructed by M. Klin
BindGlobal("KlinGraph", List([function()
        local G, H, N, V;
        H := Group([(1,2,3,4,5,6), (1,4)]);
        V := Set(List(SymmetricGroup(7), g -> H*g));
        N := Position(V, H*());
        G := EdgeOrbitsGraph(Action(SymmetricGroup(7), V, OnRight),
            List([(3,4,5,6,7), (2,4,6,3,5,7), (3,5,6,7), (2,4,5,6,7,3),
                        (2,3,5,6,7), (2,4,5,6), (2,3,5,6)],
                    g -> [N, Position(V, H*g)]));
        AssignVertexNames(G, V);
        return G;
    end])[1]());

# The cube with intersection array {3, 2, 1; 1, 2, 3}
BindGlobal("CubeGraph", HammingGraph(3, 2));

# The Heawood graph with intersection array {3, 2, 2; 1, 1, 3}.
BindGlobal("HeawoodGraph", DesarguesianPlaneIncidenceGraph(2));

# The icosahedron with intersection array {5, 2, 1; 1, 2, 5}.
BindGlobal("IcosahedronGraph", Graph(DirectProduct(AlternatingGroup(5),
                                                   Group((1, 2))),
                                     [(1,3,5,2,4)], function(p, g)
                                        return (p^g)^((-1)^(6^g));
                                     end, function(x, y)
                                        return Order(x*y) = 3;
                                     end));

# The Sylvester graph with intersection array {5, 4, 2; 1, 1, 4}
BindGlobal("SylvesterGraph", Graph(SymmetricGroup(6),
                    [[1, [[[1,2], [3,4], [5,6]],
                          [[1,3], [2,5], [4,6]],
                          [[1,4], [2,6], [3,5]],
                          [[1,5], [2,4], [3,6]],
                          [[1,6], [2,3], [4,5]]]]],
                    function(x, g)
                        return [x[1]^g,
                            Set(List(x[2], l -> OnSetsSets(l, g)))];
                    end, function(x, y)
                        return x[1] <> y[1] and x[2] <> y[2] and
                            Set([x[1], y[1]]) in Intersection(x[2], y[2])[1];
                    end));

# The Perkel graph with intersection array {6, 5, 2; 1, 1, 3}.
BindGlobal("PerkelGraph", Graph(PSL(2, 19),
                    Elements(Filtered(ConjugacyClassesSubgroups(PSL(2, 19)),
                        x -> Size(x) = 57 and Order(x[1]) = 60)[1]),
                    ConjugateGroup, function(x, y)
                        return Order(Intersection(x, y)) = 10;
                    end, true));

# The Gosset graph with intersection array {27, 10, 1; 1, 10, 27}.
BindGlobal("GossetGraph",
    Graph(DirectProduct(SymmetricGroup(8), SymmetricGroup(2)),
        [[-3,-3,1,1,1,1,1,1]], OnRoots, RootAdjacency));

# The truncated Witt graph with intersection array {15, 14, 12; 1, 1, 9}.
BindGlobal("Witt23Graph", Graph(MathieuGroup(23), [[1,4,8,12,13,19,21,23]],
                                OnSets, DisjointSets));

# The large Witt graph with intersection array {30, 28, 24; 1, 3, 15}.
BindGlobal("Witt24Graph", Graph(MathieuGroup(24), [[1,4,8,12,13,19,21,23]],
                                OnSets, DisjointSets));

# The bipartite graph associated to Higman's design
# with intersection array {50, 49, 36; 1, 14, 50}.
BindGlobal("HigmanGraph", Graph(Stabilizer(MathieuGroup(24), [23, 24], OnSets),
                            [[1,4,8,12,13,19,21,23]], OnSets,
                            function (x, y)
                                return Length(Difference(Intersection(x, y),
                                                         [23, 24])) in [0, 4];
                            end));

# The Coxeter graph with intersection array {3,2,2,1; 1,1,1,2}.
BindGlobal("CoxeterGraph", PolarGraphNOorth(1, 3, 7));

# The doubly truncated Witt graph with intersection array {7,6,4,4; 1,1,1,6}.
BindGlobal("Witt22Graph", Graph(MathieuGroup(22), [[1,2,3,4,5,10,18,21]],
                                OnSets, DisjointSets));

# The dodecahedron with intersection array {3,2,1,1,1; 1,1,1,2,3}.
BindGlobal("DodecahedronGraph", List([function()
        local G, O, act, dp, p1, p2, u, v;
        dp := DirectProduct(AlternatingGroup(5), Group((1, 2)));
        p1 := Projection(dp, 1);
        p2 := Projection(dp, 2);
        act := function(t, g)
            return OnTuples(t{Permuted([1,2], Image(p2, g))}, Image(p1, g));
        end;
        u := [1, 2];
        v := [3, 4];
        O := Arrangements([1..5], 2);
        G := EdgeOrbitsGraph(Action(dp, O, act),
                             [Position(O, u), Position(O, v)]);
        AssignVertexNames(G, O);
        return G;
    end])[1]());

# The Desargues graph with intersection array {3,2,2,1,1; 1,1,2,2,3}.
BindGlobal("DesarguesGraph", DoubledOddGraph(2));

# The Biggs-Smith graph with intersection array {3,2,2,2,1,1,1; 1,1,1,1,1,1,3}.
BindGlobal("BiggsSmithGraph", Graph(PSL(2, 17),
                    Elements(Filtered(ConjugacyClassesSubgroups(PSL(2, 17)),
                        x -> Size(x) = 102 and Order(x[1]) = 24)[1]),
                    ConjugateGroup, function(x, y)
                        return Order(Intersection(x, y)) = 8;
                    end, true));

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

Class(BarrierScratchUnparserProg, CScratchUnparserProg, rec(
    barrier_cmd := (self,o,i,is) >> Print(Blanks(i), self.opts.barrierCMD(self.opts), self.pinfix(o.args, ", "), ";\n"),
    nop_cmd := (self,o,i,is) >> Print(""),
	register := (self,o,i,is) >> Print(Blanks(i), "if(count == 0) ", self.opts.register(self.opts), self.pinfix(o.args, ", "), ";\n"),
    initialization := (self, o, i, is) >> Print(Blanks(i), self.opts.initialization(self.opts), self.pinfix(o.args, ", "), ";\n"),
	par_exec := (self,o,i,is) >> Print(Blanks(i), "parallel((void*)&sub_cpu) \n",
                              Blanks(i), "{\n", self(o.cmds[1],i+is,is),
                              Blanks(i),"}\n"),
));


text:
# According to section 21.18 of the reference manual, Sort is not stable (it's a Shell sort).
# However, SortingPerm is stable. We will see it on an example, showing indexes of elements after the sort.

n := 20;
L := List([1 .. n], i -> Random("AB"));
# "AABABBBABBABAABABBAB"

p := SortingPerm(L);
# (3,10,15,17,18,19,9,14,7,13,6,12,16,8,4)(5,11)

a := Permuted(L, p);;
b := Permuted([1 .. n], p);;

PrintArray(TransposedMat(List([1 .. n], i -> [a[i], b[i]])));
# [ [  'A',  'A',  'A',  'A',  'A',  'A',  'A',  'A',  'A',  'B',  'B',  'B',  'B',  'B',  'B',  'B',  'B',  'B',  'B',  'B' ],
#   [    1,    2,    4,    8,   11,   13,   14,   16,   19,    3,    5,    6,    7,    9,   10,   12,   15,   17,   18,   20 ] ]


text:
################################################################################
##
#W symmetric_data.gi              GroupTheoretical Package
##
#W Paul Bruillard, Cesar Galindo, Siu-Hung Ng, Julia Plavnik, Eric Rowell,
#W Zhenghan Wang
##
## Installation file for symmetric_data functions of the GroupTheoretical Package
##
#Y Copyright (C) 2016, Battelle Memorial Institute
##
################################################################################

################################################################################
##
#F symmetric_data(<group>,<group_element>). . . . compute the symmetric data for Rep(G,z)
##
InstallGlobalFunction(symmetric_data, function(G,z)
  local tbl, cc_sizes, order, chi, rank, n, m, k, chim_times_chin, N, FPdim, FPdimC, S, T;
  # first verify that z is central in G
  if not(z in Center(G)) then
    Error("z is not central in G\n");
  fi;

  # verify that z has order <= 2
  if Order(z) > 2 then
    Error("z has order > 2\n");
  fi;

  # first coompute the fusion matrices
  tbl:=CharacterTable(G);
  cc_sizes:=SizesConjugacyClasses(tbl);
  chi:=Irr(tbl);
  order:=Size(G);
  rank:=Size(cc_sizes);

  N:=List([1..rank], n->NullMat(rank,rank));

  for n in [1..rank] do
    for m in [n..rank] do
      chim_times_chin:=ComplexConjugate(chi[m]*chi[n]);
      for k in [1..rank] do
        N[n][m][k]:= Sum([1..rank],i->cc_sizes[i]*chi[k][i]*chim_times_chin[i])/order;
        N[m][n][k]:=N[n][m][k];
      od;
    od;
  od;

  # get the dimensions and form the S-matrix
  FPdim := List(chi,x->x[1]);
  FPdimC := Sum(List(FPdim,x->x^2));
  S:=List(FPdim,x->List(FPdim,y->x*y));

  if z = Identity(G) then
    T:=List(FPdim,x->1);
  else
    # VERIFY THIS
    Display("VERIFY THAT THIS IS CORRECT\n");
    n := Position(ConjugacyClasses(tbl),ConjugacyClass(G,z));
    T:=ShallowCopy(chi[n]);
  fi;

  return [chi,S,T,N,FPdim,FPdimC];
end);
#E symmetric_data.gi  . . . . . . . . . . . . . . . . . . . . . . . . ends here

max_stars_repo_name: pnnl/GroupTheoretical

max_issues_repo_name: pnnl/GroupTheoretical

max_forks_repo_name: pnnl/GroupTheoretical


text:
#############################################################################
##
#W  automaton.gi              automgrp package                 Yevgen Muntyan
#W                                                             Dmytro Savchuk
##
#Y  Copyright (C) 2003 - 2018 Yevgen Muntyan, Dmytro Savchuk
##

###############################################################################
##
#R  IsMealyAutomatonRep
##
DeclareRepresentation("IsMealyAutomatonRep",
                      IsComponentObjectRep and IsAttributeStoringRep,
                      ["table",     # transitions: table[i][j] is the state automaton goes to after processing letter j
                                    # if it was in state i
                       "perms",     # output: perms[i] is the output function at state i. It is either IsPerm or IsTransformation
                       "trans",     # images of elements of perms: i^perms[j] = trans[j][i]
                       "alphabet",  # by default it's [1..degree]; may be something different e.g. after taking "cross product"
                                    # it's used for pretty printing, internally only [1..degree] is used
                       "degree",    # Length(alphabet)
                       "states",    # names of the states
                       "n_states",  # Length(states), Length(table)
                       ]);
BindGlobal("AG_AutomatonFamily", NewFamily("AutomatonFamily", IsMealyAutomaton, IsMealyAutomaton, IsMealyAutomatonFamily));

BindGlobal("AG_IsInvertible",
function(t)
  return RankOfTransformation(t) = DegreeOfTransformation(t);
end);

BindGlobal("_AG_CreateAutomaton",
function(table, states, alphabet)
  local a, s, i, n_states, invertible, degree, perms;

  if not AG_IsCorrectAutomatonList(table, false) then
    Error("'", table, "' is not a correct table representing a finite automaton");
  fi;

  if states = fail then
    states := List([1..Length(table)], i -> Concatenation(AG_Globals.state_symbol, String(i)));
  else
    if Length(states) <> Length(table) then
      Error("number of state names is not equal to the number of states");
    fi;
    if Length(Set(states)) <> Length(states) then
      Error("duplicated states names");
    fi;
    states := List(states, s -> String(s));
  fi;

  if alphabet = fail then
    alphabet := [1 .. Length(table[1])-1];
  else
    if Length(alphabet) <> Length(table[1]) - 1 then
      Error("number of letters in alphabet does not agree with the automaton table");
    fi;
    alphabet := StructuralCopy(alphabet);
  fi;

  n_states := Length(table);
  degree := Length(alphabet);
  perms := List(table, s -> s[degree+1]);
  table := List(table, s -> s{[1..degree]});
  invertible := true;

  for i in [1..n_states] do
    if not IsPerm(perms[i]) and not AG_IsInvertible(perms[i]) then
      invertible := false;
      break;
    fi;
  od;

  if not invertible then
    for i in [1..n_states] do
      if IsPerm(perms[i]) then
        perms[i] := Transformation(List([1..degree], j -> j^perms[i]));
      fi;
    od;
  else
    for i in [1..n_states] do
      if not IsPerm(perms[i]) then
        perms[i] := PermList(ImageListOfTransformation(perms[i]));
      fi;
    od;
  fi;

  a := Objectify(NewType(AG_AutomatonFamily, IsMealyAutomaton and IsMealyAutomatonRep),
                 rec(table := table,
                     perms := perms,
                     trans := List([1..n_states], i -> List([1..degree], j -> j^perms[i])),
                     states := states,
                     n_states := n_states,
                     alphabet := alphabet,
                     degree := degree));

  SetIsInvertible(a, invertible);

  return a;
end);

BindGlobal("__AG_PrintPerm",
function(p)
  if IsPerm(p) then
    Print(p);
  else
    Print(ImageListOfTransformation(p));
  fi;
end);

BindGlobal("__AG_PermString",
function(p)
  if IsPerm(p) then
    return String(p);
  else
    return String(ImageListOfTransformation(p));
  fi;
end);

InstallMethod(MealyAutomaton, [IsList, IsList, IsList],
function(table, states, alphabet)
  return _AG_CreateAutomaton(table, states, alphabet);
end);

InstallMethod(MealyAutomaton, [IsList, IsList],
function(table, states)
  return _AG_CreateAutomaton(table, states, fail);
end);

InstallMethod(MealyAutomaton, [IsList],
function(table)
  return _AG_CreateAutomaton(table, fail, fail);
end);

InstallMethod(MealyAutomaton, [IsString],
function(string)
  local ret;
  ret := AG_ParseAutomatonString(string);
  return MealyAutomaton(ret[2], ret[1]);
end);

InstallMethod(MealyAutomaton, [IsTreeHomomorphism],
function(a)
  return MealyAutomaton([a]);
end);

InstallMethod(MealyAutomaton, [IsList and IsTreeHomomorphismCollection],
function(tree_hom_list)
  return MealyAutomaton(tree_hom_list, false);
end);

###############################################################################
##
##
##
InstallMethod(MealyAutomaton, [IsList, IsObject],
function(tree_hom_list, name_func)
  local a, states, names, MealyAutomatonLocal, aut_list;

  MealyAutomatonLocal := function(g)
    local cur_state;
    if g in states then return Position(states, g); fi;
    Add(states, g);
    if IsFunction(name_func) then
      Add(names, name_func(g));
    elif name_func = true then
      Add(names, Concatenation("<", String(g), ">"));
    fi;
    cur_state := Length(states);
    aut_list[cur_state] := List([1..g!.deg], x -> MealyAutomatonLocal(Section(g, x)));
    Add(aut_list[cur_state], g!.perm);
    return cur_state;
  end;

  names := [];
  states := [];
  aut_list := [];
  for a in tree_hom_list do
    MealyAutomatonLocal(a);
  od;

  if not IsEmpty(names) then
    return MealyAutomaton(aut_list, names);
  else
    return MealyAutomaton(aut_list);
  fi;
end);

InstallMethod(MealyAutomaton, [IsSelfSim],
function(a)
  local states, MealyAutomatonLocal, aut_list;

  MealyAutomatonLocal := function(g)
    local cur_state;
    if g!.word in states then return Position(states, g!.word); fi;
    Add(states, g!.word);
    cur_state := Length(states);
    aut_list[cur_state] := List([1..g!.deg], x -> MealyAutomatonLocal(Section(g, x)));
    Add(aut_list[cur_state], g!.perm);
    return cur_state;
  end;

  states := [];
  aut_list := [];
  MealyAutomatonLocal(a);
  return MealyAutomaton(aut_list);
end);

InstallMethod(SetStateName, [IsMealyAutomaton, IsInt, IsString],
function(aut, k, name)
  local new_names;
  if k < 1 or k > aut!.n_states then
    Error("wrong state number ", k);
  fi;
  new_names := List(aut!.states);
  new_names[k] := name;
  SetStateNames(aut, new_names);
end);

InstallMethod(SetStateNames, [IsMealyAutomaton, IsList],
function(aut, names)
  if not IsDenseList(names) or Length(names) <> aut!.n_states or
     not ForAll(names, IsString) or Length(AsSet(names)) <> aut!.n_states
  then
    Error("invalid state names list, it must be a dense list of strings of length ",
          aut!.n_states, " without duplicates");
  fi;
  aut!.states := MakeImmutable(List(names));
end);

InstallMethod(Display, [IsMealyAutomaton and IsMealyAutomatonRep],
function(a)
  local i, j;

  for i in [1..a!.n_states] do
    Print(a!.states[i], " = (");
    for j in [1..a!.degree] do
      Print(a!.states[a!.table[i][j]]);
      if j <> a!.degree then
        Print(", ");
      fi;
    od;
    Print(")");
    if not IsOne(a!.perms[i]) then
      __AG_PrintPerm(a!.perms[i]);
    fi;
    if i <> a!.n_states then
      Print(", ");
    fi;
  od;
end);

InstallMethod(String, [IsMealyAutomaton and IsMealyAutomatonRep],
function(a)
  local i, j, str;
  str := "";
  for i in [1..a!.n_states] do
    Append(str,Concatenation(String(a!.states[i]), " = ("));
    for j in [1..a!.degree] do
      Append(str,String(a!.states[a!.table[i][j]]));
      if j <> a!.degree then
        Append(str,", ");
      fi;
    od;
    Append(str, ")");
    if not IsOne(a!.perms[i]) then
      Append(str,__AG_PermString(a!.perms[i]));
    fi;
    if i <> a!.n_states then
      Append(str,", ");
    fi;
  od;
  return str;
end);

InstallMethod(ViewObj, [IsMealyAutomaton],
function(a)
  Print("<automaton>");
end);

InstallMethod(PrintObj, [IsMealyAutomaton],
function(a)
  Print("MealyAutomaton(\"",String(a),"\")");
end);

BindGlobal("AG_AutomatonTransform",
function(a, q, x)
  local i, j, nq, nx, t;

  nq := Length(q);
  nx := Length(x);

  for i in [1..nq] do
    for j in [1..nx] do
      t := x[j];
      x[j] := a!.trans[q[i]][t];
      q[i] := a!.table[q[i]][t];
    od;
  od;
end);

# InstallMethod(TransitionFunction, [IsMealyAutomaton and IsMealyAutomatonRep],
# function(a)
#   local func, table;
#
#   table := a!.table;
#
#   func := function(Q, X)
#
#   end;
#
#   return func;
# end);

InstallMethod(AutomatonList, [IsMealyAutomaton],
function(a)
  return List([1..a!.n_states], i -> Concatenation(a!.table[i], [a!.perms[i]]));
end);

InstallMethod(NumberOfStates, [IsMealyAutomaton],
function(A)
  return A!.n_states;
end);

InstallMethod(SizeOfAlphabet, [IsMealyAutomaton],
function(A)
  return A!.degree;
end);

InstallMethod(AG_MinimizedAutomatonList, "for [IsMealyAutomaton]", [IsMealyAutomaton],
function(A)
  return AG_AddInversesList(List(AutomatonList(A), x -> List(x)));
end);

InstallGlobalFunction(MinimizationOfAutomatonTrack,
function(a)
  local min_aut;
  min_aut := AG_MinimizationOfAutomatonListTrack(List(AutomatonList(a), x -> List(x)));
  return [MealyAutomaton(min_aut[1], a!.states{min_aut[2]}), min_aut[2], min_aut[3]];
end);

InstallGlobalFunction(MinimizationOfAutomaton,
function(a)
  local min_aut;
  min_aut := AG_MinimizationOfAutomatonListTrack(List(AutomatonList(a), x -> List(x)));
  return MealyAutomaton(min_aut[1], a!.states{min_aut[2]});
end);

# this is the method which uses PolynomialDegreeOfGrowthOfUnderlyingAutomaton
#InstallMethod(IsOfPolynomialGrowth, "for [IsMealyAutomaton]", true,
#              [IsMealyAutomaton],
#function(A)
#  local G, res;
#  G := AutomatonGroup(A);
#  res := IsGeneratedByAutomatonOfPolynomialGrowth(G);
#  if res then
#    SetPolynomialDegreeOfGrowth(A, PolynomialDegreeOfGrowthOfUnderlyingAutomaton(G));
#  fi;
#  SetIsBounded(A, IsGeneratedByBoundedAutomaton(G));
#  return res;
#end);

InstallMethod(IsOfPolynomialGrowth, "for [IsMealyAutomaton]", true,
              [IsMealyAutomaton],
function(A)
  local i, d, ver, nstates, cycles, cycle_of_vertex,
        IsNewCycle, known_vertices, aut_list, HasPolyGrowth,
        cycle_order, next_cycles, cur_cycles, cur_path, cycles_of_level,
        lev, ContainsTrivState, s;

  IsNewCycle := function(C)
    local i, l, cur_cycle, long_cycle;

    l := [2..Length(C)];
    Add(l, 1);
    long_cycle := PermList(l);

    for cur_cycle in cycles do
      if Intersection(cur_cycle, C) <> [] then
#        if Length(C)<>Length(cur_cycle) then return fail; fi;
#        for i in [0..Length(C)-1] do
#          if cur_cycle=Permuted(C, long_cycle^i) then return false; fi;
#        od;
        Info(InfoAutomGrp, 5, "cycle1=", cur_cycle, "cycle2=", C);
        return fail;
      fi;
    od;

    return true;
  end;

#  Example:
#  cycles = [[1, 2, 4], [3, 5, 6], [7]]
#  cur_cycles = [1, 3] (the first and the third cycles)
#  cycle_order = [[2, 3], [3], []] (means 1 -> 2 -> 3,  1 -> 3)

  HasPolyGrowth := function(v)
    local i, v_next, is_new, C, ver;
#    Print("v=", v, "\n");
    Add(cur_path, v);
    for i in [1..d] do
      v_next := aut_list[v][i];
      if not (v_next in known_vertices or v_next = nstates+1) then
        if v_next in cur_path then
          C := cur_path{[Position(cur_path, v_next)..Length(cur_path)]};
          is_new := IsNewCycle(C);
          if is_new = fail then
            return false;
          fi;

          Add(cycles, C);
          Add(cycle_order, []);
          for ver in C do
#              Print("next_cycles = ", next_cycles);
            UniteSet(cycle_order[Length(cycles)], next_cycles[ver]);
            cycle_of_vertex[ver] := Length(cycles);
            next_cycles[ver] := [Length(cycles)];
          od;
        else
          if not HasPolyGrowth(v_next) then
            return false;
          fi;
          if cycle_of_vertex[v] = 0 then
            UniteSet(next_cycles[v], next_cycles[v_next]);
          elif cycle_of_vertex[v] <> cycle_of_vertex[v_next] then
            UniteSet(cycle_order[cycle_of_vertex[v]], next_cycles[v_next]);
            Info(InfoAutomGrp, 5, "v=", v, "; v_next=", v_next);
            Info(InfoAutomGrp, 5, "cycle_order (local) = ", cycle_order);
          fi;
        fi;
      elif v_next in known_vertices then
        if cycle_of_vertex[v]=0 then
          UniteSet(next_cycles[v], next_cycles[v_next]);
        elif cycle_of_vertex[v] = cycle_of_vertex[v_next] then
          return false;
        else
          UniteSet(cycle_order[cycle_of_vertex[v]], next_cycles[v_next]);
        fi;

      fi;
    od;

    Remove(cur_path);
    Add(known_vertices, v);
    return true;
  end;

  aut_list := AG_MinimizationOfAutomatonList(List(AutomatonList(A), x -> List(x)));

# below we put the trivial state to the last position in the aut_list
  ContainsTrivState := false;

  for s in [1..Length(aut_list)] do
    if AG_IsTrivialStateInList(s, aut_list) then
      ContainsTrivState := true;
      if s < Length(aut_list) then
        aut_list := AG_PermuteStatesInList(aut_list, (s, Length(aut_list)));
      fi;
      break;
    fi;
  od;

  if not ContainsTrivState then
    SetIsBounded(A, false);
    return false;
  fi;

  nstates := Length(aut_list) - 1;
  d := A!.degree;
  cycles := [];
  cycle_of_vertex := List([1..nstates], x -> 0);  #if vertex i is in cycle j, then cycle_of_vertex[i]=j
  next_cycles := List([1..nstates], x -> []); #if vertex i is not in a cycle, next_cycles[i] stores the list of cycles, that can be reached immediately (with no cycles in between) from this vertex
  known_vertices := [];
  cur_path := [];
  cycle_order := [];

  while Length(known_vertices) < nstates do
    ver := Difference([1..nstates], known_vertices)[1];
    if not HasPolyGrowth(ver) then
      SetIsBounded(A, false);
      return false;
    fi;
  od;

# Now we find the longest chain in the poset of cycles
  cycles_of_level := [[]];
  for i in [1..Length(cycles)] do
    if cycle_order[i] = [] then
      Add(cycles_of_level[1], i);
    fi;
  od;

  lev := 1;

  while cycles_of_level[Length(cycles_of_level)] <> [] do
    Add(cycles_of_level, []);
    for i in [1..Length(cycles)] do
      if Intersection(cycles_of_level[lev], cycle_order[i]) <> [] then
        Add(cycles_of_level[lev+1], i);
      fi;
    od;
    lev := lev+1;
  od;

  if lev<=2 then
    SetIsBounded(A, true);
  else
    SetIsBounded(A, false);
  fi;
  SetPolynomialDegreeOfGrowth(A, lev-2);
  Info(InfoAutomGrp, 5, "Cycles = ", cycles);
  Info(InfoAutomGrp, 5, "cycle_order = ", cycle_order);
  Info(InfoAutomGrp, 5, "next_cycles = ", next_cycles);
  return true;
end);

InstallMethod(IsBounded, "for [IsMealyAutomaton]", true,
              [IsMealyAutomaton],
function(A)
  IsOfPolynomialGrowth(A);
  return IsBounded(A);
end);

InstallMethod(PolynomialDegreeOfGrowth, "for [IsMealyAutomaton]", true,
              [IsMealyAutomaton],
function(A)
  local res;

  res := IsOfPolynomialGrowth(A);

  if not res then
    Info(InfoAutomGrp, 1, "Error: the automaton <A> has exponential growth");
    return fail;
  fi;

  return PolynomialDegreeOfGrowth(A);
end);

InstallMethod(DualAutomaton, "for [IsMealyAutomaton]", true,
              [IsMealyAutomaton],
function(A)
  local list, dual_list, states, d, n;

  list := AutomatonList(A);
  d := Length(list[1]) - 1;
  n := Length(list);
  dual_list := List([1..d], i -> Concatenation(List([1..n], j -> i^list[j][d+1]),
                                               [Transformation(List([1..n], j -> list[j][i]))]));
  states := List([1..d], i -> Concatenation(AG_Globals.state_symbol_dual, String(i)));

  return MealyAutomaton(dual_list, states);
end);

InstallMethod(InverseAutomaton, "for [IsMealyAutomaton]", true,
              [IsMealyAutomaton],
function(A)
  local list, inv_list, states, n;

  if not IsInvertible(A) then
    Error("MealyAutomaton <A> is not invertible");
  fi;

  list := AutomatonList(A);
  n := Length(list);
  inv_list := AG_InverseAutomatonList(list);
  states := List([1..n], i -> Concatenation(AG_Globals.state_symbol, String(i)));

  return MealyAutomaton(inv_list, states);
end);

InstallMethod(Inverse, "for [IsMealyAutomaton]", true,
              [IsMealyAutomaton],
function(A)
  return InverseAutomaton(A);
end);

###############################################################################
##
#M  \^
##
##  Constructs the power <power> of an automaton, where <power> can be negative as well
##

InstallMethod(\^, "for [IsMealyAutomaton, IsCyclotomic]", [IsMealyAutomaton, IsCyclotomic],
function(A, power)
  if power>0 then
    TryNextMethod();
  elif power=0 then
    return MealyAutomaton([Concatenation(List([1..A!.degree],x->1),[()])]);
  else
    return InverseAutomaton(A)^-power;
  fi;
end);

###############################################################################
##
#M  \/
##
##  Constructs the power <power> of an automaton, where <power> can be negative as well
##

InstallMethod(\/, "for [IsMealyAutomaton, IsMealyAutomaton]", [IsMealyAutomaton, IsMealyAutomaton],
function(A1, A2)
  return A1*A2^-1;
end);

InstallMethod(IsBireversible, "for [IsMealyAutomaton]", true,
              [IsMealyAutomaton],
function(A)
  return IsInvertible(A) and IsInvertible(DualAutomaton(A)) and
         IsInvertible(DualAutomaton(InverseAutomaton(A)));
end);

InstallMethod(IsReversible, "for [IsMealyAutomaton]", true,
              [IsMealyAutomaton],
function(A)
  return IsInvertible(DualAutomaton(A));
end);

InstallMethod(IsIRAutomaton, "for [IsMealyAutomaton]", true,
              [IsMealyAutomaton],
function(A)
  return IsInvertible(A) and IsInvertible(DualAutomaton(A));
end);

InstallMethod(MDReduction, "for [IsMealyAutomaton]",
              [IsMealyAutomaton],
function(A)
  local B, DB, MDB;
  B:=MinimizationOfAutomaton(A);
  DB:=DualAutomaton(B);
  MDB:=MinimizationOfAutomaton(DB);
  while NumberOfStates(MDB)<NumberOfStates(DB) do
    B:=MDB;
    DB:=DualAutomaton(B);
    MDB:=MinimizationOfAutomaton(DB);
  od;
  return [B,MDB];
end);

InstallMethod(IsMDTrivial, "for [IsMealyAutomaton]", true,
              [IsMealyAutomaton],
function(A)
  local MDT;
  MDT:=MDReduction(A);
  if NumberOfStates(MDT[1])=1 then
    return true;
  fi;
  return false;
end);

InstallMethod(IsMDReduced, "for [IsMealyAutomaton]", true,
              [IsMealyAutomaton],
function(A)
  local B, DB, MDB;
  B:=MinimizationOfAutomaton(A);
  if NumberOfStates(B)<NumberOfStates(A) then
    return false;
  else
    DB:=DualAutomaton(B);
    MDB:=MinimizationOfAutomaton(DB);
    if NumberOfStates(MDB)<NumberOfStates(DB) then
      return false;
    fi;
  fi;
  return true;
end);

###############################################################################
##
#M  \*
##
##  Constructs a product of two noninitial automata
##
InstallMethod(\*, "for [IsMealyAutomaton, IsMealyAutomaton]", [IsMealyAutomaton, IsMealyAutomaton],
function(A1, A2)
  local n, m, d, i, j, aut_list, states;

  d := A1!.degree;
  if d <> A2!.degree then
    Error("The cardinalities of alphabets in <A1> and <A2> do not coincide");
  fi;

  n := A1!.n_states;
  m := A2!.n_states;
  aut_list := [];
  for i in [1..n] do
    for j in [1..m] do
#     (i, j)                                < ->     m*(i-1)+j
#     (Int((x-1)/m)+1, ((x-1) mod m) + 1 ) < ->     x
      Add(aut_list, Concatenation(List([1..d], x -> m*(A1!.table[i][x]-1) + A2!.table[j][x^A1!.perms[i]]), [A1!.perms[i]*A2!.perms[j]]));
    od;
  od;
  states := List([1..n*m], i -> Concatenation(AG_Globals.state_symbol, String(i)));
  return MealyAutomaton(aut_list, states);
end);

InstallMethod(IsTrivial, "for [IsMealyAutomaton]", [IsMealyAutomaton],
function(A)
  # XXX trivial transformation
  return AutomatonList(MinimizationOfAutomaton(A)) = [Concatenation(List([1..A!.degree], x -> 1), [()])];
end);

InstallMethod(DisjointUnion, "for [IsMealyAutomaton, IsMealyAutomaton]", [IsMealyAutomaton, IsMealyAutomaton],
function(A, B)
  local n, m, aut_list, states, perms, tableA, tableB;

  if A!.degree <> B!.degree then
    Error("The cardinalities of alphabets in <A> and <B> do not coincide");
  fi;

  n := A!.n_states;
  m := B!.n_states;
  tableA := List(A!.table, x -> List(x));
  tableB := List(B!.table, x -> List(x));
  Append(tableA, List(tableB, x -> List(x, y -> y+n)));
  perms := Concatenation(A!.perms, B!.perms);

  aut_list := List([1..n+m], i -> Concatenation(tableA[i], [perms[i]]));

  states := List([1..n+m], i -> Concatenation(AG_Globals.state_symbol, String(i)));

  return MealyAutomaton(aut_list, states);
end);

InstallMethod(AreEquivalentAutomata, "for [IsMealyAutomaton, IsMealyAutomaton]", [IsMealyAutomaton, IsMealyAutomaton],
function(A, B)
  local n, m, i, j, aut_list, found, Am, Bm, C, equiv_statesB;

  if A!.degree <> B!.degree then
    return false;
  fi;

  Am := MinimizationOfAutomaton(A);
  Bm := MinimizationOfAutomaton(B);
  n := Am!.n_states;
  m := Bm!.n_states;
  if m <> n then
    return false;
  fi;

  C := DisjointUnion(Am, Bm);
  aut_list := AutomatonList(C);

  equiv_statesB := [];

  for i in [1..n] do
    found := false;
    for j in [n+1..n+m] do
      if (not j in equiv_statesB) and AG_AreEquivalentStatesInList(i, j, aut_list) then
        found := true;
        Add(equiv_statesB, j);
        break;
      fi;
    od;
    if not found then
      return false;
    fi;
  od;

  return true;
end);

InstallMethod(SubautomatonWithStates, "for [IsMealyAutomaton, IsList]", [IsMealyAutomaton, IsList],
function(A, states)
  local G, d, i, g, perm, subautom_list, newsubautom_list;

  G := AutomatonList(A);
  d := A!.degree;

# first we add all states of <states> in the list
  for g in states do
    for i in [1..d] do
      if not (G[g][i] in states) then
        Add(states, G[g][i]);
      fi;
    od;
  od;

  Sort(states);
  subautom_list := G{states};
  newsubautom_list := [];

  for g in subautom_list do
    perm := g[d+1];
    g := List(g{[1..d]}, x -> Position(states, x));
    Add(g, perm);
    Add(newsubautom_list, g);
  od;

  return MealyAutomaton(newsubautom_list, A!.states{states});
end);

InstallMethod(AutomatonNucleus, "for [IsMealyAutomaton]", [IsMealyAutomaton],
function(A)
  local IsElemInNucleus, G, tmp, d, Nucl, i;

  IsElemInNucleus := function(g)
    local i, res;
    if g in tmp then
      for i in [Position(tmp, g)..Length(tmp)] do
        if not (tmp[i] in Nucl) then Add(Nucl, tmp[i]); fi;
      od;
      return g=tmp[1];
    fi;
    Add(tmp, g);
    res := false; i := 1;
    while (not res) and i<=d do
      res := IsElemInNucleus(G[g][i]);
      i := i+1;
    od;
    Remove(tmp);
    return res;
  end;

  G := AutomatonList(A);
  d := A!.degree;

  Nucl := [];
# first add elements of cycles
  for i in [1..Length(G)] do
    tmp := [];
    if not (i in Nucl) then IsElemInNucleus(i); fi;
  od;

  return SubautomatonWithStates(A, Nucl);
end);

InstallMethod(AdjacencyMatrix, "for [IsMealyAutomaton]", [IsMealyAutomaton],
function(A)
  local i, s, d, M, aut_list;
  aut_list := AutomatonList(A);
  M:=List([1..A!.n_states],x->List([1..A!.n_states],x->0));
  d := A!.degree;
  for s in [1..A!.n_states] do
    for i in [1..d] do
      M[s][aut_list[s][i]]:=M[s][aut_list[s][i]]+1;
    od;
  od;
  return M;
end);

InstallMethod(IsAcyclic, "for [IsMealyAutomaton]", [IsMealyAutomaton],
function(A)
  local i, j, M, N;
  M:=AdjacencyMatrix(A);
  N:=StructuralCopy(M);
  for i in [1..(A!.n_states-1)] do
    N:=N*M;
    for j in [1..A!.n_states] do
      if N[j][j]<>M[j][j]^(i+1) then return false; fi;
    od;
  od;
  return true;
end);

#E

max_stars_repo_name: gap-packages/automgrp

max_issues_repo_name: gap-packages/automgrp

max_forks_repo_name: gap-packages/automgrp


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

Import(code); # EvalScalar

_color := (c, str) -> When(LocalConfig.osinfo.useColor(), c(str), str);

_compute_cycles_from_mflops := (flop, flops) -> IntDouble(Double(EvalScalar(fdiv(flop,flops)))*LocalConfig.cpuinfo.freq);
_compute_mflops := (flop, cycles) -> Cond(cycles<>0, IntDouble(Double(EvalScalar(fdiv(flop,cycles)))*LocalConfig.cpuinfo.freq), "#NA");
_compute_gflops := (flop, cycles) -> Cond(cycles<>0, StringDouble("%.3f", (EvalScalar(fdiv(flop,cycles))/1000.0)*LocalConfig.cpuinfo.freq), "#NA");
_compute_gflops_freq := (flop, cycles, freq) -> Cond(cycles<>0, StringDouble("%.3f", (EvalScalar(fdiv(flop,cycles))/1000.0)*freq), "#NA");

_seqPerfStats := function(file, n, cycles, search_time)
    PrintLine("n ", n, "  seq ", BlueStr(cycles), " [cyc]  search ", search_time, " [s]");
   AppendTo(file,
    PrintLine("n ", n, "  seq ", cycles, " [cyc]  search ", search_time, " [s]"));
end;

_seqPerfStatsMflops := function(file, t, flops, cycles, search_time)
    local realFlops;
    if ObjId(flops)=ListClass then realFlops := flops[2]; flops := flops[1]; else realFlops := 0; fi;
    PrintLine(t, "  ", _color(BlueStr, cycles), " [cyc]  ", _color(DarkYellowStr, _compute_mflops(flops, cycles)), " [Mf/s]  ",
     _compute_mflops(realFlops, cycles), " [Mf/s, real]  ", search_time, " [search, s]");
    AppendTo(file,
    PrintLine(t, "  ", cycles, " [cyc]  ", _compute_mflops(flops, cycles), " [Mf/s]  ",  _compute_mflops(realFlops, cycles), " [Mf/s, real]  ", search_time, " [search, s]"));
end;

_seqPerfStatsGflops := function(file, t, freq, flops, cycles, search_time)
    local realFlops;

    if ObjId(flops)=ListClass then realFlops := flops[2]; flops := flops[1]; else realFlops := 0; fi;

    PrintLine(t, "  ", _color(BlueStr, cycles), " [cyc]  ", _color(DarkYellowStr, _compute_gflops_freq(flops, cycles, freq)), " [Gf/s]  ",
     _compute_gflops_freq(realFlops, cycles, freq), " [Gf/s, real]  ", search_time, " [search, s]");
    AppendTo(file,
    PrintLine(t, "  ", cycles, " [cyc]  ", _compute_gflops_freq(flops, cycles, freq), " [Gf/s]  ", _compute_gflops_freq(realFlops, cycles, freq), " [Gf/s, real]  ", search_time, " [search, s]"));
end;

_seqPerfStatsMflopsAcc := function(file, n, flops, cycles, search_time, acc)
    local realFlops;
    if ObjId(flops)=ListClass then realFlops := flops[2]; flops := flops[1]; else realFlops := 0; fi;
    PrintLine("n ", n, "  seq ", _color(BlueStr, cycles), " [cyc] ", _color(DarkYellowStr, _compute_mflops(flops, cycles)), " [Mflop/s] search ",
    _compute_mflops(realFlops, cycles), " [Mflop/s, real] search ", search_time, " [s] accuracy ", -acc, " [dig]");
    AppendTo(file,
    PrintLine("n ", n, "  seq ", cycles, " [cyc] ", _compute_mflops(flops, cycles), " [Mflop/s] search ", _compute_mflops(realFlops, cycles), " [Mflop/s, real] search ",
    search_time, " [s] accuracy ", -acc, " [dig]"));
end;

_seqPerfStatsGflopsAcc := function(file, t, flops, cycles, search_time, acc)
    local digits, realFlops;
    if ObjId(flops)=ListClass then realFlops := flops[2]; flops := flops[1]; else realFlops := 0; fi;
    digits := When(IsFloat(acc), -IntDouble(d_log(1e-16+AbsFloat(acc))/d_log(10)), acc);
    PrintLine(t, "  ", _color(BlueStr, cycles), " [cyc]  ", _color(DarkYellowStr,
    _compute_gflops(flops, cycles)), " [Gf/s]  ", _compute_gflops(realFlops, cycles), " [Gf/s, real]  ", search_time, " [search, s] ", digits, " [dig]");
    AppendTo(file,
    PrintLine(t, "  ", cycles, " [cyc]  ", _compute_gflops(flops, cycles), " [Gf/s]  ", _compute_gflops(realFlops, cycles), " [Gf/s, real]  ",
    search_time, " [search, s] ", digits, " [dig]"));
end;

_parPerfStats := function(file, n, threads, cycles, search_time)
    PrintLine("n ", n, "  threads ", threads, "  par ", _color(BlueStr, cycles), " [cyc]  search ", search_time, " [s]");
   AppendTo(file,
    PrintLine("n ", n, "  threads ", threads, "  par ", cycles, " [cyc]  search ", search_time, " [s]"));
end;

printcount := function(opcount, countrec)
    local n, i, retval;
    retval := "";
    n := 1;
    for i in opcount do
      retval := ConcatenationString(retval, " ", String(i), " ", countrec.printstrings[n]);
      n := n + 1;
    od;
    return(retval);
end;

_seqPerfStatsGflopsCount := function(file, t, flops, cycles, search_time, count, countrec)
    local realFlops;
    if ObjId(flops)=ListClass then realFlops := flops[2]; flops := flops[1]; else realFlops := 0; fi;

    PrintLine(t, "  ", _color(BlueStr, cycles), " [cyc]  ", _color(DarkYellowStr, _compute_gflops(flops, cycles)), " [Gf/s]  ",
     _compute_gflops(realFlops, cycles), " [Gf/s, real]  ", search_time, " [search, s] ", printcount(count, countrec));

    AppendTo(file,
    PrintLine(t, "  ", cycles, " [cyc]  ", _compute_gflops(flops, cycles), " [Gf/s]  ",
     _compute_gflops(realFlops, cycles), " [Gf/s, real]  ", search_time, " [search, s] ", printcount(count, countrec)) );
end;

_seqPerfStatsGflopsAccCount := function(file, t, flops, cycles, search_time, acc, count, countrec)
    local digits, realFlops;
    if ObjId(flops)=ListClass then realFlops := flops[2]; flops := flops[1]; else realFlops := 0; fi;
    digits := When(IsFloat(acc), -IntDouble(d_log(1e-16+AbsFloat(acc))/d_log(10)), acc);

    PrintLine(t, "  ", _color(BlueStr, cycles), " [cyc]  ", _color(DarkYellowStr, _compute_gflops(flops, cycles)), " [Gf/s]  ",
     _compute_gflops(realFlops, cycles), " [Gf/s, real]  ", search_time, " [search, s] ", digits, " [dig]", printcount(count, countrec));

    AppendTo(file,
    PrintLine(t, "  ", cycles, " [cyc]  ", _compute_gflops(flops, cycles), " [Gf/s]  ",
     _compute_gflops(realFlops, cycles), " [Gf/s, real]  ", search_time, " [search, s] ", digits, " [dig]", printcount(count, countrec)) );
end;


text:
# The Kneser graph on k-subsets of a set with n elements.
BindGlobal("KneserGraph", function(arg)
    local n, k, invt, vcs;
    if Length(arg) < 2 then
        Error("at least two arguments expected");
        return fail;
    else
        n := arg[1];
        k := arg[2];
        if Length(arg) > 2 then
            invt := arg[3];
        else
            invt := true;
        fi;
    fi;
    if invt then
        vcs := Combinations([1..n], k);
    else
        vcs := [[1..k]];
    fi;
    return Graph(SymmetricGroup(n), vcs, OnSets, DisjointSets);
end);

# The Odd graph of diameter d on 2*d+1 points.
BindGlobal("OddGraph", d -> KneserGraph(2*d+1, d, false));

# The doubled Odd graph on 2*d+1 points.
BindGlobal("DoubledOddGraph", function(d)
    local n, dp;
    n := 2*d+1;
    dp := DirectProduct(SymmetricGroup(n), SymmetricGroup(2));
    return Graph(dp, Union(Combinations([1..n], d), Combinations([1..n], d+1)),
                    OnDoubledOdd(n, dp), SymmetrizedInclusion, true);
end);

# The folded Johnson graph.
BindGlobal("FoldedJohnsonGraph",
    d -> AntipodalQuotientGraph(JohnsonGraph(2*d, d))
);

# The three Chang graphs with v=28, k=12, lm=6, mu=4
BindGlobal("ChangGraph", function(j)
    local J, S, Ss;
    Ss := [List([1..4], i -> [i, i+4]),
           Set(List([1..8], i -> Set([i, (i mod 8)+1]))),
           Union(List([1..3], i -> Set([i, (i mod 3)+1])),
                 List([1..5], i -> Set([i+3, (i mod 5)+4])))];
    J := JohnsonGraph(8, 2);
    S := List(Ss[j], x -> Position(J.names, x));
    return SwitchedGraph(J, S, Stabilizer(J.group, S, OnSets));
end);


text:
# Floating point infinity
inf := FLOAT_INT(1) / FLOAT_INT(0);

IS_FLOAT(inf);
#true;

# GAP has also a formal ''infinity'' value
infinity in Cyclotomics;
# true