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

# ================================================================
# Integer Strength Reduction Rules
# ================================================================

Class(RulesStrengthReduce, RuleSet);

#F SReduceSimple(<c>) -- applies strength reduction rules to any object <c>
#F                       disables <opts> sensitive rules
SReduceSimple  := c -> ApplyStrategy(c, [RulesStrengthReduce], BUA, rec());

#F SReduce(<c>, <opts>) - applies strength reduction rules to any object <c>
#F
SReduce  := (c,opts) -> ApplyStrategy(c, [RulesStrengthReduce], BUA, opts);

#@1:=@(1);
#@2:=@(2);

@1int := @(1).cond(e -> IsValue(e) and IsInt(e.v));
@2int := @(2).cond(e -> IsValue(e) and IsInt(e.v));
@2ge_var := than@ -> @(2).cond(
    e -> IsValue(e) and IsBound(than@.val.range) and than@.val.range <= e.v);

isNegV := o -> o.v<0;
isNonZero := o -> o <> 0;

@negVal := @.target(Value).cond(isNegV);

@2Teq1 := @(2).cond( x -> x.t = @(1).val.t );

_0noneOrZero := (t) -> When(ObjId(@(0).val)=noneExp, noneExp(t), t.zero());
_sameTasOp := (opname, func) -> @@(0, [Value, noneExp]).cond((e, cx) ->
    let( op := Last(cx.(opname)),
        func(e) and ForAny(op.args, (a) -> not Same(a, e) and a.t=op.t)
    )
);

v0:=V(0);
v1:=V(1);

@vpos  := (n) -> @(n, Value, e -> IsPosInt(e.v));
@vpos0 := (n) -> @(n, Value, e -> IsPosInt0(e.v));

RewriteRules(RulesStrengthReduce, rec(

 # FF: SAR-specific rule. Not sure if they should be turned on by default.
 # YSV: below rule is invalid, because idiv(a+b, c) <> idiv(a) + idiv(b)
 #      example: idiv(2+3, 4) == 1
 #               idiv(2,4) + idiv(3,4) == 0
 # div_sar := Rule([@(1,idiv), @(2, add), @(3,Value) ], e -> ApplyFunc(add, List(@(2).val.args, i->idiv(i, @(3).val)))),

 ################
 # leq
 # ===========
 leq_single := Rule([leq, @(1)], e-> V_true),
 leq_eq := Rule([leq, @(1), @(2), @(3).cond(e->e=@(1).val)], e -> eq(@(1).val, @(2).val)),

 leq_val_mulindex     := Rule( [leq, @vpos0(1), [mul, @vpos0(2), @(3, [var, param], IsLoopIndex)]],
     e -> leq(idiv(@(1).val + @(2).val - 1, @(2).val), @(3).val)),

 leq_val_mulindex_val := Rule( [leq, @vpos0(1), [mul, @vpos(2), @(3, [var, param], IsLoopIndex)], @vpos0(4)],
     e -> leq(idiv(@(1).val + @(2).val - 1, @(2).val), @(3).val, idiv(@(4).val, @(2).val))),

 leq_mulindex_val     := Rule( [leq, [mul, @vpos0(2), @(3, [var, param], IsLoopIndex)], @vpos0(4)],
     e -> leq(@(3).val, idiv(@(4).val, @(2).val))),

 leq_index_neg := Rule([leq, @(2,var,IsLoopIndex), @(3, Value, x->x.v<0), ...], e -> V_false),
 leq_index_big := ARule(leq, [@(2,var,x->IsLoopIndex(x) and IsInt(x.range)), @(3, Value, x->x.v >= @(2).val.range-1)],
     e -> [@(2).val]),

 leq_neg_mul_index := Rule([leq, [mul, @(2, Value, x->x.v>0), @(3,var,IsLoopIndex)], @(4, Value, x->x.v<0), ...],
     e -> V_false),

 leq_index_0 := Rule([leq, @(2,var,IsLoopIndex), @(3, Value, x->x.v=0)], e -> eq(@(2).val, @(3).val)),
 leq_mul_index_0 := Rule([leq, [mul, @(2, Value, x->x.v<>0), @(3,var,IsLoopIndex)], @(4, Value, x->x.v=0)],
     e -> eq(@(3).val, @(4).val)),

 leq_neg_index_A := ARule(leq, [@(1, Value, x->x.v<=0), @(2,var,IsLoopIndex)], e -> [@(2).val]),
 leq_neg_mul_index_A := ARule(leq, [@(1, Value, x->x.v<=0), [@(0,mul), @(2, Value, x->x.v>0), @(3,var,IsLoopIndex)]],
     e -> [@(0).val]),

 leq_val_val_l := Rule([leq, @(1, Value), @(2, Value, x->@(1).val<=x), @, ...],
     e -> ApplyFunc(leq, Drop(e.args, 1))),
 leq_val_val_r := Rule([leq, ..., @, @(1, Value), @(2, Value, x->@(1).val<=x)],
     e -> ApplyFunc(leq, DropLast(e.args, 1))),

 ################
 # add/sub/mul 0
 # =============
 subzr := Rule([sub, @1, _sameTasOp("sub", _is0none)], e -> @1.val),
 subzl := Rule([sub, _sameTasOp("sub", _is0none), @1], e -> neg(@1.val)),

 div0 := Rule([@(1, [div, idiv, fdiv]), _0none, @], e->_0noneOrZero(e.t)),
 div1 := Rule([@(1, [div, idiv, fdiv]), @, @@(0).cond((e, cx) -> let( op := Last(cx.(ObjId(@(1).val).__name__)),
     Cond(IsValue(e), isValueOne(e), e=1) and ForAny(op.args, (a) -> not Same(a, e) and a.t=op.t)))],
     e -> @.val),

 add_assoc  := ARule(add, [ @(1,add) ], e -> @1.val.args),
 add_zero   := ARule(add, [_sameTasOp("add", _is0none)], e -> [ ]),
 add_single := Rule([add, @1],          e -> @1.val),

 add_cfold := Rule(@(1,add,e->Length(e.args)>2 and Length(Filtered(e.args, IsValue)) >= 2),
     e -> let(s := SplitBy(e.args, IsValue), values := s[1], exps := s[2],
          ApplyFunc(add, Concatenation([ApplyFunc(add,values)], exps)))),

 sub_sub := Rule( [@(1, sub), @(2), [sub, @(3), @(4)]], e -> sub(add(@(2).val, @(4).val), @(3).val)),

 ################
 # logic and/or
 # ============
 logic_single := Rule([@(1, [logic_and, logic_or]), @1], e->@1.val),
 logic_empty := Rule([@(1, [logic_and, logic_or])], e->V_true),

 logic_and_assoc := ARule(logic_and, [ @(1,logic_and) ], e -> @1.val.args),
 logic_and_false := Rule([logic_and, ..., _vfalse, ...],  e -> V_false),
 logic_and_true := ARule(logic_and, [_vtrue],  e -> [ ]),

 logic_or_assoc := ARule(logic_or, [ @(1,logic_or) ], e -> @1.val.args),
 logic_or_true := Rule([logic_or, ..., _vtrue, ...],  e -> V_true),
 logic_or_false := ARule(logic_or, [_vfalse],  e -> [ ]),

 logic_neg_neg := Rule([logic_neg, [logic_neg, @(1)]], e -> @(1).val),
 logic_neg_eq := Rule([logic_neg, [eq, @(1), @(2)]], e -> neq(@(1).val, @(2).val)),

 # 'maybe' is a "magic" boolean function that behaves as 'true' inside 'and'/'or' operators,
 # but also satisfies the uncertainty rule logic_not(maybe()) = maybe()
 #
 logic_neg_maybe := Rule([logic_neg, @(1, maybe)], e -> @(1).val),
 logic_and_maybe := ARule(logic_and, [maybe], e -> [ ]),
 logic_or_maybe := ARule(logic_or, [maybe], e -> [ ]),

 eq_logic_true := Rule( [eq, @(1), _vtrue], e -> @(1).val ),

 ################
 # mul
 # ===

 # mul_assoc is slow, and we enable it in RulesExpensiveStrength reduce
 #   mul_assoc := ARule(mul, [ @(1,mul) ], e -> @1.val.args),
 mul_one := ARule(mul, [ _sameTasOp("mul", _is1) ],  e -> [ ]),

 mul_single := Rule([mul, @1], e->@1.val),
 mul_zero := Rule([mul, ..., _v0none, ...], e -> _0noneOrZero(e.t)),  # must be @(0).val so that none*x -> none
 mul_mul := Rule([mul, @(1,Value), [mul, @(2,Value), @(3)]], e -> let(
     t := UnifyPair(@(1).val.t, @(2).val.t),
     val := t.value(t.product(@(1).val.v, @(2).val.v)),
     mul(val, @(3).val))),

 # apply distributivity a * (b+c) -> a*b + a*c, only if inside nth()
 mul_add_nth := Rule([@@(0,mul,(e,cx)->IsBound(cx.nth) and cx.nth<>[]), @(1), @(2,add)],
     e -> ApplyFunc(add, List(@(2).val.args, a->@(1).val * a))),

 # FF: NOTE: Need to guard this rule in the SAR case. Not sure how to properly do this.
 # Guard against nth indexing integer tables. This is at the moment only used in SAR.
 # NB: Without guard, there is an infinite loop with sar_sub rule
 mul_sub_nth := Rule([@@(0, mul, (e,cx) -> cx.isInside(nth) and cx.nth[1].loc.t.t <> TInt), @(1), @(2,sub)],
     e -> ApplyFunc(sub, List(@(2).val.args, a->@(1).val * a))),  # Above guard to prevent inf loop with sar_sub

 # NB: Below rule became obsolete due to change in segment length computation
 # sar_sub := Rule([sub, @(1,Value,e->e.t=TInt), [mul, @(2, Value, e->e.t=TInt and IsInt(@(1).val.v/e.v)), @(3)]],
 #    e->mul(@(2).val, sub(@(1).val/@(2).val, @(3).val))),

 # Value * (Value * a +/- Value * b) -> Value * a +/- Value * b
 mul_add_mul := Rule([mul, @(1, Value), [@(0,[add,sub]), [mul, @(2, Value), @(3)], [mul, @(4, Value), @(5)]]],
     e -> let(t1 := UnifyPair(@(1).val.t, @(2).val.t),
              t2 := UnifyPair(@(1).val.t, @(4).val.t),
              op := ObjId(@(0).val),
              op(t1.value(t1.product(@(1).val.v, @(2).val.v)) * @(3).val,
                 t2.value(t2.product(@(1).val.v, @(4).val.v)) * @(5).val))),

 #  mul_cfold := Rule(@(1,mul,e->Length(e.args)>2 and Length(Filtered(e.args, IsValue)) >= 2),
 #   e -> let(s := SplitBy(e.args, IsValue), values := s[1], vars := s[2],
 #            ApplyFunc(mul, Concatenation( [Product(values,v->v.v)], vars)))),

 fpmul_2power := Rule([fpmul, @(1, Value), @(2, Value, e->Is2Power(e.v)), @(3)], e ->
    let(sh := Log2Int(@(2).val.v),
        Cond(sh = @(1).val.v, @(3).val,
             sh > @(1).val.v, arith_shl(@(3).val, sh - @(1).val.v),
             sh < @(1).val.v, arith_shr(@(3).val, @(1).val.v - sh)))),

 # mul_add_expand := Rule([mul, @(1,Value,e->e.t=TInt), @(2,add)],
 #    e -> ApplyFunc(add, List(@2.val.args, a->mul(@1.val, a).eval()))),
 #
 # mul_sub_expand := Rule([mul, @(1,Value,e->e.t=TInt), @(2,sub)],
 #    e -> ApplyFunc(sub, List(@2.val.args, a->mul(@1.val, a).eval()))),

 ##############
 # neg
 # =====
 mul_negone := Rule([mul, _neg1, @2], e -> neg(@2.val)),
 neg_neg := Rule([neg, [neg, @1]], e -> @1.val), # -(-1) -> 1
 neg_sub := Rule([neg, [sub, @1, @2]], e -> sub(@2.val, @1.val)), # -(1-2) -> 2-1
 neg_zero := Rule([neg, _v0none], e -> e.args[1]),

 sub_neg1 := Rule([sub, @1, [neg, @2]], e -> add(@1.val, @2.val)),  # 1-(-2) -> 1+2
 sub_neg2 := Rule([sub, [neg, @1], @2], e -> neg(add(@1.val, @2.val))), # (-1)-2 -> -(1+2)
 add_neg1 := Rule([add, @1, [neg, @2]], e -> sub(@1.val, @2.val)), # 1+(-2) -> 1-2
 add_neg2 := Rule([add, [neg, @1], @2], e -> sub(@2.val, @1.val)), # (-1)+2 -> 2-1

 mul_neg1 := Rule([mul, @1, [neg, @2]], e -> neg(mul(@1.val, @2.val))), # 1*(-2) -> -(1*2)
 mul_neg2 := Rule([mul, [neg, @1], @2], e -> neg(mul(@1.val, @2.val))), # (-1)*2 -> -(1*2)
 mul_negc1 := Rule([mul, @1, @(2,Value,isNegV)], e -> neg(mul(Value.new(@2.val.t, -@2.val.v), @1.val))),
 mul_negc2 := Rule([mul, @(1,Value,isNegV), @2], e -> neg(mul(Value.new(@1.val.t, -@1.val.v), @2.val))),

 #################################
 # general constant folding
 # ==========================

 # The rules below can be used instead of generic cfold to speed up this ruleset
 # cfold_add  := Rule([add, Value,Value], e->Value.new(e.args[1].t, e.args[1].v + e.args[2].v)),
 # cfold_mul  := Rule([mul, Value,Value], e->Value.new(e.args[1].t, e.args[1].v * e.args[2].v)),
 # cfold_sub  := Rule([sub, Value,Value], e->Value.new(e.args[1].t, e.args[1].v - e.args[2].v)),
 # cfold_imod := Rule([imod, Value,Value], e->Value.new(e.args[1].t, e.args[1].v mod e.args[2].v)),
 # cfold_idiv := Rule([idiv, Value,Value], e->Value.new(e.args[1].t, QuoInt(e.args[1].v, e.args[2].v))),
 # cfold_nth  := Rule([nth, @(1,var,e->IsBound(e.value)), Value], e->e.loc.value.v[e.idx.v+1]),
 # cfold_nthv := Rule([nth, Value,Value],                         e->e.loc.v[e.idx.v+1]),
 # cfold_re   := Rule([re, Value], e->Value.new(TDouble, ReComplex(Complex(e.args[1].v)))),
 # cfold_im   := Rule([im, Value], e->Value.new(TDouble, ImComplex(Complex(e.args[1].v)))),
 # cfold_neg  := Rule([neg, Value], e->Value.new(e.args[1].t, -e.args[1].v)),

 # cfold_arith := Rule([@(1,[add,mul,sub]), Value,Value], e->e.eval()),
 # cfold_divmod := Rule([@(1,[imod,idiv]), Value,Value], e->e.eval()),
 # cfold_unary := Rule([@(1,[re,im,neg]), Value], e->e.eval()),
 # cfold_nth := Rule([nth, @(1,var,e->IsBound(e.value)), Value], e->e.eval()),
 # cfold_nthv := Rule([nth, Value, Value], e->e.eval()),

 cfold := Rule(@(1).cond(e -> IsExp(e) and e.can_fold()), e -> e.eval()),

 # In below rules @(1) is a type, and if t is a type,
 # t.value creates a value of this type, and thus performs the typecast
 tcast_varvalue := Rule([tcast, @(1), @(2,var,e->IsBound(e.value))],
     e -> Cond(IsPtrT(@(1).val), @(2).val, @(1).val.value(@(2).val.v))), # NOTE: IsPtrT branch is unsafe, and relies on correct
                                                                         #        scalarization/non-scalarization
# tcast_value    := Rule([tcast, @(1), @(2,Value)], e -> Cond(
#     # if an array value is typecast to a pointer of same type, just drop the typecast
#     IsPtrT(@(1).val) and IsArrayT(@(2).val.t) and @(1).val.t = @(2).val.t.t, @(2).val,
#     @(1).val.value(@(2).val.v))),

 # only applies to pointers, otherwise is not valid -- for example (double)((int)3.2) is not (double)3.2
 tcast_tcast := Rule([tcast, @(1).cond(t->IsPtrT(t)), [tcast, @(2).cond(t->IsPtrT(t)), @(3)]],
     e -> tcast(@(1).val, @(3).val)),

 # allow adding zero to ptrs.
 tptr_addzero := ARule(add, [@@(0, Value).cond( (e, cx) ->
    let( op := Last(cx.("add")),
        _is0none(e) and ObjId(op.t) = TPtr
    ))],
    ee -> []
 ),

 # cfold_vref := Rule(@(1).cond(e -> IsRec(e) and IsBound(e.isVref) and e.isVref and e.canScalarize() and not IsBound(e.didEval)),
 #     e -> e.eval()),

 # cfold_vref := Rule(@(1).cond(e -> IsRec(e) and IsBound(e.isVref) and e.isVref and e.canScalarize() and not IsBound(e.didEval)),
 #     e -> e.eval()),

 ##################
 # cond rules
 # ============
 cond_rid := Rule([cond, @(1), @(2), @(3).cond(e->e=@(2).val)], e -> @(2).val),

 cond_fold := Rule([cond, @(1,Value), @2, ...], e -> Cond(
         @1.val.v = false or @1.val.v = 0,
             When(Length(e.args) = 3, e.args[3], ApplyFunc(cond, Drop(e.args, 2))),
         @1.val.v = true or IsInt(@1.val.v),
             (@2.val),
         Error("non boolean condition"))),

 IF_fold := Rule([IF, @(1,Value), @(2), @(3)], e -> Cond(
         @(1).val.v = false or @1.val.v = 0,
             @(3).val,
         @(1).val.v = true or IsInt(@1.val.v),
             @(2).val,
         Error("non boolean condition"))),

 IF_skip := Rule([IF, @(1), skip, skip], e -> skip()),

 multi_if2 := Rule([multi_if, @(1), @(2)], e -> IF(@(1).val, @(2).val, skip())),
 multi_if3 := Rule([multi_if, @(1), @(2), @(3)], e -> IF(@(1).val, @(2).val, @(3).val)),
 multi_if_skip := Rule([multi_if, skip], e->skip()),
 multi_if_skip_skip := Rule([@(1, multi_if), @(2, e->ObjId(Last(@(1).val.rChildren()))=skip), skip, ...],
     e -> ApplyFunc(multi_if, Drop(@(1).val.rChildren(), 2))),
 #NOTE: code in autolib's _genPlan messing with 'args' directly, have to drop empty multi_if even though it could fold automatically
 multi_if_drop := ARule( chain, [[multi_if]], e -> [] ),

 ###################
 # chain rules
 # =============
 chain_empty := Rule([chain], e-> skip()),

 chain_skip := Rule([chain, skip], e-> skip()),

 ##########################
 # rules for functions
 # =====================
 frotate_rank1 := Rule([frotate, @(1).cond(e->e.rank()<=1)], e -> @(1).val),
 fcall_fexp  := Rule(
     [@(0,fcall),
	 @(1,[fcurry,flift,fsplit,frotate,Lambda], e -> e.rank()+2=Length(@(0).val.args)),
	 ...],
     e -> e.eval()),

 fexp_Lambda := Rule(
     [@(1,[fcurry,flift,fsplit,frotate]),
	 Lambda,
	 ...],
     e -> e.eval()),

 fcurry_void := Rule(
     [fcurry,
	 @(1),
	 @(2).cond(x -> IsValue(x) and ObjId(@(1).val.t)=TFunc and @(1).val.t.params[x.v] in [TVoid, TDummy]),
         @(3).cond(x -> x<>0)],
     x -> fcurry(@(1).val, @(2).val, 0)),

 fcall_void := Rule(
     [@(0,fcall),
	 @(1).cond(x -> ObjId(x.t)=TFunc and
	                let(pos := Position(x.t.params, TDummy),
			    pos <> false and
			    Length(@(0).val.args) >= 1+pos and
			    not IsValue(@(0).val.args[1+pos]))),
	 ...],
     e -> let(types := e.args[1].t.params,
	      len := Length(e.args),
	      ApplyFunc(fcall, [e.args[1]] :: List([2..len], i -> Cond(types[i-1]=TDummy, 0, e.args[i]))))),

 # if lambdaWrap(Lambda(..)) has Lambda of the form x -> exp, where exp does not depend on x, then eliminate the Lambda
 #lambdaWrap_dummyLambda := Rule([lambdaWrap, @(1, Lambda, x->x.t.params[1]=TDummy and Length(x.t.params)=2)], e -> @(1).val.at(0)),

 #lambdaWrap_rank0 := Rule([lambdaWrap, @(1).cond(x->x.rank()=0)], e -> fcall(@(1).val, 0)),

 ###############
 # omega
 # =======
 re_omega := Rule([re, [omega, @(1), @(2)]], e -> cospi(fdiv(2*@(2).val, @(1).val))),
 im_omega := Rule([im, [omega, @(1), @(2)]], e -> sinpi(fdiv(2*@(2).val, @(1).val))),

 re_omegapi := Rule([re, [omegapi, @(1)]], e -> cospi(@(1).val)),
 im_omegapi := Rule([im, [omegapi, @(1)]], e -> sinpi(@(1).val)),

 re_conj_omega := Rule([re, [conj, [omega, @(1), @(2)]]], e -> cospi(fdiv(2*@(2).val, @(1).val))),
 im_conj_omega := Rule([im, [conj, [omega, @(1), @(2)]]], e -> -sinpi(fdiv(2*@(2).val, @(1).val))),

 ######################
 # complex arith
 # ===============
 conj_conj := Rule([conj, [conj, @(1)]], e->@(1).val),
 re_conj := Rule([re, [conj, @(1)]], x->re(@(1).val)),
 im_conj := Rule([im, [conj, @(1)]], x->-im(@(1).val)),

 # the re_cond/im_cond rules prevent cond from having different types in different branches (ie TReal vs TComplex)
 #
 re_cond := Rule([re, @(1,cond)], x -> ApplyFunc(cond,
	 List(@(1).val.args, x->Cond(x.t=TBool, x, re(x))))),
 im_cond := Rule([im, @(1,cond)], x -> ApplyFunc(cond,
	 List(@(1).val.args, x->Cond(x.t=TBool, x, im(x))))),

 re_real := Rule([re, @.cond(e->e.t=TReal or ObjId(e.t)=T_Real)], e -> e.args[1]),
 im_real := Rule([im, @.cond(e->e.t=TReal or ObjId(e.t)=T_Real)], e -> e.t.zero()),

 re_mul := Rule([re, [@(1, mul), @(2).cond(e->e.t=TReal), ...]], e -> @(2).val * re(ApplyFunc(mul, Drop(@(1).val.args, 1)))),
 im_mul := Rule([im, [@(1, mul), @(2).cond(e->e.t=TReal), ...]], e -> @(2).val * im(ApplyFunc(mul, Drop(@(1).val.args, 1)))),

 re_cxpack := Rule([re, [cxpack, @(1), @]], e -> @(1).val),
 im_cxpack := Rule([im, [cxpack, @, @(1)]], e -> @(1).val),

 re_vdup := Rule([re, [vdup, @(1), @(2)]], e -> vdup(re(@(1).val), @(2).val)),
 im_vdup := Rule([im, [vdup, @(1), @(2)]], e -> vdup(im(@(1).val), @(2).val)),

 mul_cxpack := ARule( mul, [[cxpack, @(1), @(2)], [cxpack, @(3), @(4)]],
     e -> [cxpack( @(1).val*@(3).val - @(2).val*@(4).val,
                   @(1).val*@(4).val + @(2).val*@(3).val )]),

 add_cxpack := ARule( add, [[cxpack, @(1), @(2)], [cxpack, @(3), @(4)]],
     e -> [cxpack( add(@(1).val, @(3).val), add(@(2).val, @(4).val) )]),

 sub_cxpack := Rule( [sub, [cxpack, @(1), @(2)], [cxpack, @(3), @(4)]],
     e -> [cxpack( sub(@(1).val, @(3).val), sub(@(2).val, @(4).val) )]),

 # re( (a + i*b) * (a - i*b) ) -> a*a + b*b
 re_cxmulconj := Rule( [re, [mul, [cxpack, @(1), @(2)], [conj, [cxpack, @(3).cond(x->x=@(1).val), @(4).cond(x->x=@(2).val)]]]],
     e -> add(mul(@(1).val, @(1).val), mul(@(2).val, @(2).val))),
 # im( (a + i*b) * (a - i*b) ) -> 0
 im_cxmulconj := Rule( [im, [mul, [cxpack, @(1), @(2)], [conj, [cxpack, @(3).cond(x->x=@(1).val), @(4).cond(x->x=@(2).val)]]]],
     e -> e.t.zero() ),

 #############
 # imod
 # ======
 imod_by_one := Rule([imod, @1, _1], e -> v0),
 imod_zero := Rule([imod, _0none, @], e -> @(0).val),
 imod_one := Rule([imod, _1, @], e -> v1),

 imod_imod := Rule([imod, [imod, @(1,Value), @(2)], @(3,Value, e->@(2).val.v <= e.v )],
     e -> imod(@(1).val, @(2).val)),

 imod_small := Rule([imod, @(1,var,e->IsBound(e.range) and not IsSymbolic(e.range)),  @(2,Value,e->e.v >= @(1).val.range)],
     e -> @(1).val),

 # this is a less general version of imod_mul_reduce_const, commented out below
 # v1*exp mod v3 -> ((v1 mod v3)*exp) mod v3
 imod_mul := Rule([imod, [@(2, mul), @(1,Value), ...], @(3,Value, e->@(1).val.v >= e.v)],
     e -> imod(ApplyFunc(mul, [@(1).val.v mod @(3).val.v] :: Drop(@(2).val.args, 1)), @(3).val)),

 # (v1*exp + exp4) mod v3 -> ((v1 mod v3)*exp + exp4) mod v3
 imod_addmul_l := Rule([imod, [add, [@(2,mul), @(1,Value), ...], @(4)], @(3,Value, e->@(1).val.v >= e.v)],
     e -> imod(ApplyFunc(mul, [@(1).val.v mod @(3).val.v] :: Drop(@(2).val.args, 1)) + @(4).val, @(3).val)),
 imod_addmul_r := Rule([imod, [add, @(4), [@(2,mul), @(1,Value), ...]], @(3,Value, e->@(1).val.v >= e.v)],
     e -> imod(ApplyFunc(mul, [@(1).val.v mod @(3).val.v] :: Drop(@(2).val.args, 1)) + @(4).val, @(3).val)),

 no_mod_imod := Rule([no_mod, [imod, @1, @2]], e -> @1.val),
 no_mod_small_mod := Rule([no_mod, [small_mod, @1, @2]], e -> @1.val),

 # imod_mul_reduce_const := Rule( [imod, @(1,mul),
 #        @(2,Value, e -> Filtered(@1.val.args, a->IsValue(a) and a.v >= e.v) <>[]) ],
 #    e -> imod(
 #      ApplyFunc(mul, List(@1.val.args, a->Cond(IsValue(a), V(a.v mod @2.val.v), a))),
 #      @2.val)),

 # NOTE: These rules are very general but lead to an infinite loop.
 #
 # imod_mul := Rule([imod, @(1, mul, NOT_ALL_ARGS_ARE(imod)), @(2)],
 #     e -> imod(ApplyFunc(mul, List(@(1).val.args, a -> imod(a, @(2).val))), @(2).val)),
 #
 # imod_add := Rule([imod, @(1, [add,sub], NOT_ALL_ARGS_ARE(imod)), @(2)],
 #     e -> imod(ApplyFunc(ObjId(@(1).val), List(@(1).val.args, a -> imod(a, @(2).val))), @(2).val)),

 #############
 # idiv
 # ======
 idiv_small := Rule([idiv, @(1, var, e -> IsBound(e.range) and not IsSymbolic(e.range)), @(2, Value, e-> e.v >= @(1).val.range)],
     e -> V(0)),

 # exp / value ->  (1/value) * exp,  for floating-point types.  NOTE: use .t.isFloat()
 div_byV_tomul := Rule([@(0, [div, fdiv]), @(1), @(2, Value, x->x.t in [TReal, TComplex])],
     e -> div(@(2).val.t.one(), @(2).val) * @(1).val),

 # (v1*exp2 + exp3)/v4 -> v1/v4*exp2 + exp3/v4; v4 | v1;
 # NOTE: these rules are invalid if sign(exp2)<>sign(exp4) !!!

 #idiv_addmul_l := Rule([idiv, [add,  [mul, @(1, Value), @(2)], @(3)], @(4, Value, x -> _divides(x, @(1).val))],
 #    e -> add(mul(div(@(1).val, @(4).val), @(2).val), idiv(@(3).val, @(4).val))),

 #idiv_addmul_r := Rule([idiv, [add,  @(3), [mul, @(1, Value), @(2)]], @(4, Value, x -> _divides(x, @(1).val))],
 #    e -> add(mul(div(@(1).val, @(4).val), @(2).val), idiv(@(3).val, @(4).val))),

 # (v1*loop_idx + v2)/v3 -> v1/v3*loop_idx + v2/v3;  when  (v3 | v1) & (v2 >= 0) & (v1 >= 0);
 # assumption: loop index cannot be negative
 idiv_addmul_pos_l := Rule([idiv, [add,  [mul, @vpos0(1), @(2).cond(IsLoopIndex)], @vpos0(3)], @(4, Value, x -> _divides(x, @(1).val))],
     e -> add(mul(idiv(@(1).val, @(4).val), @(2).val), idiv(@(3).val, @(4).val))),

 idiv_addmul_pos_r := Rule([idiv, [add,  @vpos0(3), [mul, @vpos0(1), @(2).cond(IsLoopIndex)]], @(4, Value, x -> _divides(x, @(1).val))],
     e -> add(mul(idiv(@(1).val, @(4).val), @(2).val), idiv(@(3).val, @(4).val))),

 ############
 # pow, sqrt, rsqrt
 # =====
 pow_zero := Rule([pow, @1, _0none], e->v1),

 powmod_mul_const := Rule([powmod, @, Value, @(1,mul,e->ForAny(e.args,IsValue)), @],
    e -> let(split := SplitBy(@1.val.args, IsValue),
         value := Product(split[1], x->x.v),
         vars := split[2],
         powmod( e.args[1],
            (e.args[2].v ^ value) mod e.args[4].v,
                 ApplyFunc(mul, vars),
             e.args[4]))),

 powmod_cfold := Rule([powmod, @1, Value, Value, Value],
    e -> imod(((e.args[2].v ^ e.args[3].v) mod e.args[4].v)*e.args[1], e.args[4].v)),

 sqrt_pow2 := Rule( [sqrt, [pow, @(1), _2]], e -> abs(@(1).val) ),
 pow2_sqrt := Rule( [pow, [sqrt, @(1)], _2], e -> @(1).val ),

 rsqrt_pow2 := Rule( [rsqrt, [pow, @(1), _2]], e -> fdiv(e.t.one(), abs(@(1).val)) ),
 pow2_rsqrt := Rule( [pow, [rsqrt, @(1)], _2], e -> fdiv(e.t.one(), @(1).val) ),

 #############
 # bin
 # ====

 bin_xor_zero      := Rule([bin_xor, ..., _0none, ...], e -> When(@(0).val = e.args[2], e.args[1], e.args[2])),
 bin_and_zero      := Rule([bin_and, ..., _0none, ...], e -> @(0).val),

 bin_shr_zero      := Rule([bin_shr, _0none,      @], e -> e.args[1]),
 bin_shr_zeroshift := Rule([bin_shr,      @, _0none], e -> e.args[1]),
 bin_shl_zero      := Rule([bin_shl, _0none,      @], e -> e.args[1]),
 bin_shl_zeroshift := Rule([bin_shl,      @, _0none], e -> e.args[1]),

 absdiff2_to_abssub := Rule([@(0, absdiff2, x -> not IsOrdT(x.t.base_t())), @(1), @(2)],
     e -> Cond( IsValue(@(1).val) and isValueZero(@(1).val),
                    @(2).val,
                IsValue(@(1).val) and Length(Set(Flat([@(1).val.v])))=1,
                    sub(@(1).val, @(2).val),
                # else
                    abs(sub(@(1).val, @(2).val)))),

 # Memory allocation
 allocate_0 := Rule([allocate, @(1), [TArray, @(2), _0]], e -> assign(@(1).val, null())),
 deallocate_0 := Rule([deallocate, @(1), [TArray, @(2), _0]], e -> skip()),

 # funcExp rule

 expr_funcExp := Rule( [@(0, [mul, add, sub, div, idiv, imod]), ..., @(1, funcExp), ...],
     e -> @(1).val ),

 expr_cond_funcExp := Rule( [@(1, [mul, add]), ..., [@(2, cond), ..., funcExp, ...], ...],
     e -> let( args := RemoveList(e.args, @(2).val),
         _map_cond(@(2).val, p -> p, exp -> ApplyFunc(ObjId(e), args :: [exp])))),

 divmod_cond_funcExp := Rule( [@(1, [div, idiv, imod]), [@(2, cond), ..., funcExp, ...], @(3)],
     e -> _map_cond(@(2).val, p -> p, exp -> ObjId(e)(exp, @(3).val))),

 # always pull cond() from nth() for nth(funcExp) annihilation.
 nth_cond := Rule( [nth, @(1), @(2, cond)],
     e -> _map_cond(@(2).val, p->p, exp->nth(@(1).val, exp))),

 #############
 vdup_noneExp := Rule([vdup, [noneExp, @(1)], @(2)], e -> vdup(@(1).val.zero(), @(2).val) ),

 vdup_vdup := Rule([vdup, [vdup, @(1), @(2)], @(3)], e -> vdup(@(1).val, @(2).val * @(3).val)),

 vdup_vpack := Rule([vdup, @(1, vpack), @(3)], e -> ApplyFunc(vpack, ConcatList(@(1).val.args, x -> [x,x])))
));

Class(stickyNeg, neg, rec(doPeel := false, unparse := "neg"));

Class(RulesExpensiveStrengthReduce, RuleSet);

__groupSummands := function(lst, is_full)
    local e, grp, non, grpid, id, mult, i, one, ptrs;
    [grp, grpid, non] := [rec(), Set([]), []];
    [ptrs, lst] := SplitBy(lst, e -> IsPtrT(e.t) or IsArrayT(e.t));
    for e in lst do
        one := e.t.one();
        if ObjId(e) in [var, param] then
	    if not IsBound(grp.(e.id)) then  grp.(e.id) := [one, e];
	    else 		             grp.(e.id)[1] := grp.(e.id)[1] + one;
	    fi;
	    AddSet(grpid, e.id);
        elif ObjId(e) in [neg, stickyNeg] and ObjId(e.args[1]) in [var, param] then
	    id := e.args[1].id;
	    if not IsBound(grp.(id)) then  grp.(id) := [-one, e.args[1]];
	    else 		           grp.(id)[1] := grp.(id)[1] - one;
	    fi;
	    AddSet(grpid, id);
        elif ObjId(e)=mul and Length(e.args)=2 and ObjId(e.args[1])=Value and ObjId(e.args[2]) in [var, param] then
	    id := e.args[2].id;
	    if not IsBound(grp.(id)) then  grp.(id) := [e.args[1], e.args[2]];
	    else 		           grp.(id)[1] := grp.(id)[1] + e.args[1];
	    fi;
	    AddSet(grpid, id);
	elif ObjId(e) in [ neg, stickyNeg ] and ObjId(e.args[1]) = mul and Length(e.args[1].args) = 2
               and ObjId(e.args[1].args[1]) = Value and ObjId(e.args[1].args[2]) in [ var, param ]  then
            id := e.args[1].args[2].id;
            if not IsBound(grp.(id))  then grp.(id) := [ -e.args[1].args[1], e.args[1].args[2] ];
            else                           grp.(id)[1] := grp.(id)[1] - e.args[1].args[1];
            fi;
            AddSet(grpid, id);
        elif not is_full then
	    Add(non, [one, e]);
	else
	    mult := one;
	    if ObjId(e) in [neg, stickyNeg] then mult := -mult; e := e.args[1]; fi;
	    if ObjId(e) = mul and IsValue(e.args[1]) then mult := mult*e.args[1]; e := ApplyFunc(mul, Drop(e.args, 1)); fi;
	    if ObjId(e) in [neg, stickyNeg] then mult := -mult; e := e.args[1]; fi;

	    i := PositionProperty(non, x->x[2]=e);
	    if i=false then Add(non, [mult, e]);
	    else
		non[i][1] := non[i][1] + mult;
	    fi;
	fi;
    od;
    return ApplyFunc(add, ptrs :: List(grpid, id -> grp.(id)[1] * grp.(id)[2]) :: List(non, x->x[1]*x[2]));
end;

#F __groupSummandsVar(<lst>) is similar to __groupSummands but pulls out
#F variables. This helps to implement rolling pointers in Autolib
#F because grouping this way helps to get stride expression for a loop index variable.
#F

__groupSummandsVar := function(lst)
    local e, grp, non, grpid, id, mult, i, one, ptrs;
    [grp, grpid, non] := [rec(), Set([]), []];
    [ptrs, lst] := SplitBy(lst, e -> IsPtrT(e.t) or IsArrayT(e.t));
    for e in lst do
        one := e.t.one();
        if ObjId(e) in [var, param] then
	    if not IsBound(grp.(e.id)) then  grp.(e.id) := [one, e];
	    else 		             grp.(e.id)[1] := grp.(e.id)[1] + one;
	    fi;
	    AddSet(grpid, e.id);
        elif ObjId(e) in [neg, stickyNeg] and ObjId(e.args[1]) in [var, param] then
	    id := e.args[1].id;
	    if not IsBound(grp.(id)) then  grp.(id) := [-one, e.args[1]];
	    else 		           grp.(id)[1] := grp.(id)[1] - one;
	    fi;
	    AddSet(grpid, id);
        elif ObjId(e)=mul and ObjId(Last(e.args))=var then
	    id := Last(e.args).id;
	    if not IsBound(grp.(id)) then  grp.(id) := [ApplyFunc(mul, DropLast(e.args, 1)), Last(e.args)];
	    else 		           grp.(id)[1] := grp.(id)[1] + ApplyFunc(mul, DropLast(e.args, 1));
	    fi;
	    AddSet(grpid, id);
	elif ObjId(e) in [ neg, stickyNeg ] and ObjId(e.args[1]) = mul
               and ObjId(Last(e.args[1].args)) = var then
            id := Last(e.args[1].args).id;
            if not IsBound(grp.(id))  then grp.(id) := [ -ApplyFunc(mul, DropLast(e.args[1].args, 1)), Last(e.args[1].args) ];
            else                           grp.(id)[1] := grp.(id)[1] - ApplyFunc(mul, DropLast(e.args[1].args, 1));
            fi;
            AddSet(grpid, id);
        else
	    Add(non, [one, e]);
	fi;
    od;
    return ApplyFunc(add, ptrs :: List(grpid, id -> grp.(id)[1] * grp.(id)[2]) :: List(non, x->x[1]*x[2]));
end;

_groupSummands := lst -> __groupSummands(lst, true);
_groupSummandsSimple := lst -> __groupSummands(lst, false);

RewriteRules(RulesExpensiveStrengthReduce, rec(
  mul_add := Rule([mul, @(1,Value), [add, @(2, Value), @(3)]], e -> @(1).val*@(2).val + @(1).val*@(3).val),

  mul_assoc := ARule(mul, [ @(1,mul) ], e -> @1.val.args),

  mul_cfold := Rule(@(1,mul,e->Length(e.args)>2 and Length(Filtered(e.args, IsValue)) >= 2),
   e -> let(s := SplitBy(e.args, IsValue), values := s[1], exps := s[2],
            ApplyFunc(mul, Concatenation( [Product(values,v->v.v)], exps)))),

  add_assoc := Rule([add, ..., @(1,add), ... ], e -> let(
	  args := ConcatList(e.args, x->Cond(ObjId(x)=add, x.args, [x])),
	  #args)),
	  _groupSummands(args))),

  add_cfold := Rule(@(1,add,e->Length(e.args)>2 and Length(Filtered(e.args, IsValue)) >= 2),
        e -> let(s := SplitBy(e.args, IsValue), values := s[1], exps := _groupSummands(s[2]),
            ApplyFunc(add, Concatenation([Sum(values,v->v.v)], [exps])))),

  add_sub_val := Rule( [add, ..., [@(1, sub), @(2), @(3, Value)], ...],
      e -> ApplyFunc(add, [@(2).val, neg(@(3).val)] :: Filtered(e.args, x -> not Same(x, @(1).val)))),
  sub_add_val := Rule( [@(1, sub), @(2, add), @(3, Value)],
      e -> ApplyFunc(add, [neg(@(3).val)] :: @(2).val.args)),

  # super ugly:  ((a+b+c)-b) -> a+c
  # group summands could do this
  sub_add_same_exp := Rule( [sub, @(1, add), @(2).cond(x -> x in @(1).val.args)],
      e -> let( a := SplitBy(@(1).val.args, x -> x=@(2).val),
                ApplyFunc(add, a[2] :: Drop(a[1], 1) :: [a[1][1].t.zero()]))),
  add_sub_same_exp := Rule( [@(1, add), ..., [@(2, sub), @(3), @(4).cond(x -> x in @(1).val.args)], ...],
      e -> let( a := SplitBy(@(1).val.args, x -> x=@(4).val),
                b := SplitBy(a[2], x -> x=@(2).val),
                ApplyFunc(add, b[2] :: Drop(a[1], 1) :: Drop(b[1], 1) :: [@(3).val, a[1][1].t.zero()]))),

# fold the constants in an xor expression with >2 arguments of which >=2 arguments are constants
  xor_cfold := Rule(@(1,xor,e -> Length(e.args)>2 and Length(Filtered(e.args, IsValue))>=2),
    e -> let(
      s := SplitBy(e.args, IsValue),
      values := s[1],
      vars := s[2],

      ApplyFunc(xor, Concatenation( [Xor(values, v -> v.v)], vars ) )
    )
  ),

  xor_zero := ARule(xor, [_0none], e -> []),

  xor_one := Rule(@(1,xor, e -> Length(e.args) = 1),
    e -> e.args[1]
  ),

# NOTE: removed <idiv> since the rule is invalid for non-divisible numbers
#        can one create a new rule for idiv?

   mul_div := Rule([mul, ..., @(1, div), ...], e ->
       let(facs := e.args, divkind := ObjId(@(1).val), dd := Filtered(facs, f->ObjId(f)=divkind),
       divkind(Product(List(facs, f->When(ObjId(f)=divkind, f.args[1], f))),
           Product(List(dd, d->d.args[2]))))),

   imod_imod := Rule([imod, [imod, @(1), @(2, Value)], @(3, Value, e->@(2).val.v mod e.v = 0)],
       e -> imod(@(1).val, @(3).val)),

   imod_addmul := Rule([imod,[@(1,add),...,[@(2,mul),...,@(3,Value),...],...],@(4,Value,e->Mod(@(3).val.v,e.v)=0)],
       e->let(p1:=FirstPosition(@1.val.rChildren(),e->e=@(2).val),a:=ApplyFunc(add,ListWithout(@(1).val.rChildren(),p1)),imod(a,@(4).val))),

   imod_sub := Rule([imod,  [@(0,sub), @(1), @(2)],  @(3).cond(x -> _divides(x, @(1).val) or _divides(x, @(2).val))],
       e -> When(_divides(@(3).val, @(2).val), imod(@(1).val, @(3).val), imod(-@(2).val, @(3).val))),

   imod_add := Rule([@(1, imod), [ add, ..., @(2).cond(x -> _divides(@(1).val.args[2], x)), ... ], @],
       e -> imod(ApplyFunc(add, Filtered(e.args[1].args, x -> not _divides(e.args[2], x))), e.args[2])),

   imod_divides := Rule( [imod, @(1), @(2).cond( x -> _divides(x, @(1).val) )],
       e -> e.t.zero() ),

   # NOTE: idiv rule disabled, because its invalid if summands have different signs (a+2)/2 != a/2 + 1, if a<0
   #        can this be overcome?
   #idiv_addmul  := Rule([idiv,[@(1,add),...,[@(2,mul),...,@(3,Value),...],...],@(4,Value,e->Mod(@(3).val.v,e.v)=0)],
   #    e->let(p1:=FirstPosition(@1.val.rChildren(),e->e=@(2).val),a:=ApplyFunc(add,ListWithout(@(1).val.rChildren(),p1)),add(idiv(a,@(4).val),idiv(@(2).val,@(4).val)))),

   # NOTE: idiv rule disabled, because its invalid if summands have different signs (a+2)/2 != a/2 + 1, if a<0
   #        can this be overcome?
   divs_addmul2 := Rule([@(1, [div, fdiv]), [ add, ..., [ @(2, mul), ..., @(3).cond(x -> x = @(1).val.args[2]), ...], ... ], @],
       e -> add( ApplyFunc(mul, DropCond(@(2).val.args, x -> x=e.args[2], 1)),
                 ObjId(e)(ApplyFunc(add, DropCond(e.args[1].args, x -> x = @(2).val, 1)), e.args[2]))),

   # NOTE: idiv rule disabled, because its invalid if summands have different signs (a+2)/2 != a/2 + 1, if a<0
   #        can this be overcome?
   divs_addmul3 := Rule([@(1, [div, fdiv]), [ add, ..., @(2).cond(x -> x = @(1).val.args[2]), ...], @],
       e -> add( e.t.one(), ObjId(e)(ApplyFunc(add, DropCond(e.args[1].args, x -> x = @(2).val, 1)), e.args[2]))),

   div_div1 := Rule([div, @(1), [div, @(2), @(3)]], e -> div(@(1).val*@(3).val, @(2).val)),
   div_div2 := Rule([div, [div, @(1), @(2)], @(3)], e -> div(@(1).val, @(2).val*@(3).val)),

# NOTE: removed <idiv> since the rule is invalid for non-divisible numbers
#        can one create a new rule for idiv?
#   idiv_idiv1 := Rule([idiv, @(1), [idiv, @(2), @(3)]], e -> idiv(@(1).val*@(3).val, @(2).val)),
#   idiv_idiv2 := Rule([idiv, [idiv, @(1), @(2)], @(3)], e -> idiv(@(1).val, @(2).val*@(3).val)),

   fdiv_fdiv1 := Rule([fdiv, @(1), [fdiv, @(2), @(3)]], e -> fdiv(@(1).val*@(3).val, @(2).val)),
   fdiv_fdiv2 := Rule([fdiv, [fdiv, @(1), @(2)], @(3)], e -> fdiv(@(1).val, @(2).val*@(3).val)),

   mul_middleValue := Rule(@(1,mul,e -> PositionProperty(e.args, IsValue) > 1 and
                                    Length(Filtered(e.args, IsValue))=1),
       e -> let(i := PositionProperty(e.args, IsValue),
       Product(Permuted(e.args, (1,i))))),

   add_Value_1st := Rule(@(1,add,e->ForAny(e.args, IsValue) and not IsValue(e.args[1])),
    e -> let(s := SplitBy(e.args, IsValue), values := s[1], vars := s[2],
        ApplyFunc(add, [ApplyFunc(add, values)] :: vars))),

# NOTE: removed <idiv> since the rule is invalid for non-divisible numbers
#        can one create a new rule for idiv?

   div_add_Values := Rule([@(0,[div,fdiv]), [@(1,add), @(2, Value), ...],
                                             @(3, Value, e->IsInt(@(2).val.v) and (@(2).val.v mod e.v = 0))],
       e -> add(@(2).val.v / @(3).val.v,
            ObjId(@(0).val)(ApplyFunc(add, Drop(@(1).val.args, 1)), @(3).val))),

   # NOTE: this rule ignores the 3rd parameter of bin_shl (!!), which is 'bits'
   div_add_shl := Rule([@(0,[div,fdiv]), [@(1,add), ..., @(2, arith_shl)],
                                             @(3, Value, e->IsInt(@(2).val.args[2].v)
                                                 and Is2Power(e.v)
                                                 and Log2Int(e.v)<=@(2).val.args[2].v)],
       e -> add(
              arith_shl(@(2).val.args[1],@(2).val.args[2].v - Log2Int(@(3).val.v)),
            ObjId(@(0).val)(ApplyFunc(add, DropLast(@(1).val.args, 1)), @(3).val))),

   div_mul_mul_cfold1 := Rule([@(0,[idiv,div,fdiv]),
       [@(1,mul), @(2,Value), ...],
       [@(3,mul), @(4,Value,
           e -> (@(2).val.t = TInt and e.t=TInt and (@(2).val.v mod e.v)=0)), ...]],
   e -> ObjId(@(0).val)(ApplyFunc(mul, Concatenation([@(2).val.v/@(4).val.v], Drop(@(1).val.args, 1))),
                        ApplyFunc(mul, Drop(@(3).val.args, 1)))),

   div_mul_mul_cfold2 := Rule([@(0,[idiv,div,fdiv]),
       [@(1,mul), @(2,Value), ...],
       [@(3,mul), @(4,Value,
           e -> (@(2).val.t = TInt and e.t=TInt and (e.v mod @(2).val.v)=0)), ...]],
   e -> ObjId(@(0).val)(ApplyFunc(mul, Drop(@(1).val.args, 1)),
                        ApplyFunc(mul, Concatenation([@(4).val.v/@(2).val.v], Drop(@(3).val.args, 1))))),

   div_mul_Value_cfold1 := Rule([@(0,[div, idiv, fdiv]),
       [@(1,mul), @(2,Value), ...],
       @(4,Value, e->(@(2).val.t = TInt and e.t=TInt and (@(2).val.v mod e.v)=0))],
   e -> ApplyFunc(mul, Concatenation([@(2).val.v/@(4).val.v], Drop(@(1).val.args, 1)))),

   div_mul_Value_cfold2 := Rule([@(0,[div,idiv, fdiv]),
       [@(1,mul), @(2,Value), ...],
       @(4,Value, e->(@(2).val.t = TInt and e.t=TInt and (e.v mod @(2).val.v)=0))],
   e -> ObjId(@(0).val)(Product(Drop(@(1).val.args, 1)),
                        @(4).val.v/@(2).val.v)),

   div_mulmul_cfold_vars1 := Rule([@(0,[div,idiv,fdiv]),
       @(1, mul, e->Length(Filtered(e.args, IsLoc)) >= 1),
       @(2, mul, e->let(vv:=Filtered(e.args, IsLoc), Length(vv) >= 1 and Intersection(vv, @(1).val.args)<>[]))],
   e -> let(s1 := SplitBy(@(1).val.args, IsLoc), vars1 := s1[1], nonvars1 := s1[2],
            s2 := SplitBy(@(2).val.args, IsLoc), vars2 := s2[1], nonvars2 := s2[2],
            ObjId(@(0).val)(
                Product(nonvars1) * Product(ListDifference(vars1, vars2)),
                Product(nonvars2) * Product(ListDifference(vars2, vars1))))),

   div_mulmul_cfold_vars2 := Rule([@(0,[div,idiv,fdiv]),
       @(1, mul, e->Length(Filtered(e.args, IsLoc)) >= 1),
       @(2).cond(e->IsLoc(e) and e in @(1).val.args) ],
   e -> let(a := @(1).val.args, ApplyFunc(mul, ListWithout(a, Position(a, @(2).val))))),

   div_mulmul_cfold_vars3 := Rule([@(0,[div,idiv,fdiv]),
       @(1).cond(IsLoc),
       @(2, mul, e->@(1).val in e.args) ],
   e -> let(a := @(2).val.args, ObjId(@(0).val)(1, ApplyFunc(mul, ListWithout(a, Position(a, @(1).val)))))),

   div_cfold_vars4 := Rule([@(0,[div, idiv, fdiv]),
       @(1).cond(IsLoc),
       @(2).cond(x->IsLoc(x) and x=@(1).val)],
   e -> e.t.one()),

   pdiv_param := Rule([pdiv, @(0, param), @(1, Value)], e -> div(@(0).val, @(1).val)),

   pdiv_add := Rule([pdiv, @(0, add), @(1, Value)], e -> ApplyFunc(add, List(@(0).val.rChildren(), x->pdiv(x, @(1).val)))),

   #idiv is a trick to avoid inversion of this rule by mul_div
   pdiv_mul := Rule([pdiv, [mul, @(1, param), @(2)], @(3, Value)], e -> mul(idiv(@(1).val, @(3).val), @(2).val)),

   pdiv_const := Rule([pdiv, [mul, @(1, Value), @(2)], @(3, Value, e->Mod(@(1).val.v, e.v)=0)],
       e -> mul(div(@(1).val, @(3).val), @(2).val)),

   merge_logic_and := Rule([ logic_and, ..., logic_and, ... ],
       x -> ApplyFunc(logic_and, ConcatList(x.args, e->When(ObjId(e)=logic_and, e.args, [e]) )))
));

RulesMergedStrengthReduce := MergedRuleSet(RulesStrengthReduce, RulesExpensiveStrengthReduce);

ESReduce := (c, opts) -> ApplyStrategy(c, [RulesMergedStrengthReduce], BUA, opts);

Class(NormalizeLinearExp, RuleSet);
RewriteRules(NormalizeLinearExp, rec(
    # Normalization
    mul_zero := Rule([mul, _0none, @], e -> @(0).val),
    mul_one := Rule([mul, _1, @], e -> @.val),
    mul_assoc := ARule(mul, [ @(1,mul) ], e -> @1.val.args),

    mul_Value_Value := Rule([mul, @(1,Value), @(2,Value)], # UnifyTypes([@(1).val.t, @(2).val.t]).value
    e -> V(@(1).val.v * @(2).val.v)),

    mul_add_rt := Rule([mul, ..., @(2,add)],
        e -> ApplyFunc(add, List(@(2).val.args,
                a -> ApplyFunc(mul, Concatenation(Copy(DropLast(e.args, 1)),[a]))))), #Copy prevents aliasing
    mul_add_lft := Rule([mul, @(2,add), ...],
        e -> ApplyFunc(add, List(@(2).val.args,
                a -> ApplyFunc(mul, Concatenation([a], Copy(Drop(e.args, 1))))))), #Copy prevents aliasing

    mul_mul := Rule([mul, @(1,Value), [mul, @(2,Value), @(3)]], e -> mul(@(1).val.v * @(2).val.v, @(3).val)),

    sub_to_add := Rule([sub, @(1), @(2)], e -> add(@(1).val, stickyNeg(@(2).val))),

    add_assoc := Rule([add, ..., @(1,add), ... ], e -> let(
	  args := ConcatList(e.args, x->Cond(ObjId(x)=add, x.args, [x])),
	  #args)),
          _groupSummandsSimple(args))),

    add_zero := ARule(add, [_0none],  e -> [ ]),
    add_single := Rule([add, @(1)], e -> @(1).val),
    add_none := Rule([add], e -> V(0)),

    # Constant folding
    cfold := RulesStrengthReduce.rules.cfold,

    mul_cfold := RulesExpensiveStrengthReduce.rules.mul_cfold,
    add_cfold := RulesExpensiveStrengthReduce.rules.add_cfold,
    # put constants first
    mul_middleValue := RulesExpensiveStrengthReduce.rules.mul_middleValue,

    # negation
    neg        := Rule(neg, e -> stickyNeg(e.args[1])),
    mul_negone := Rule([mul, _neg1, @2], e -> stickyNeg(@2.val)),
    neg_neg := Rule([stickyNeg, [stickyNeg, @1]], e -> @1.val), # -(-1) -> 1

    mul_neg := Rule([mul, ..., @(1, [stickyNeg, Value], e->When(ObjId(e)=Value, isNegV(e), true)), ...],
        e -> let(cnt := Counter(0),
                 count := List(e.args, x -> When(ObjId(x)=stickyNeg or (ObjId(x)=Value and isNegV(x)), cnt.next(), 0)),
                 newargs := List(e.args, x -> Cond(
                         ObjId(x)=stickyNeg, x.args[1],
                         ObjId(x)=Value and isNegV(x), Value.new(x.t, -x.v),
                         x)),
                 res := ApplyFunc(mul, newargs),
                 When(IsOddInt(cnt.n), stickyNeg(res), res))),
    sink_stickyNeg_add := Rule([stickyNeg, @(1, add)], e -> ApplyFunc(add, List(@(1).val.args, stickyNeg))),

    # This rule is in conflict with mul_div rule from ESReduce but it simplifies index computations.
    # div implies divisibility but loop index takes different values from .range field -> other
    # arguments must be divisible by @(3)
    div_pull_out_idx := Rule([div, [@(1, mul), ..., @.cond(IsLoopIndex), ...], @(3, Value)],
        e -> let( m := SplitBy(@(1).val.args, IsLoopIndex),
                  ApplyFunc(mul, [div(ApplyFunc(mul, m[2]), @(3).val)] :: m[1]) )),
));

GroupSummandsExp := s -> NormalizeLinearExp(SubstBottomUp(NormalizeLinearExp(s), add, e->_groupSummands(e.args)));


text:
############################################################################
##
##  matmeths.gi                  IRREDSOL                   Burkhard Höfling
##
##  Copyright © 2003–2016 Burkhard Höfling
##

############################################################################
##
#M  Degree(<G>)
##
##  see IRREDSOL documentation
##
InstallMethod(Degree, "for matrix group", true, [IsMatrixGroup], 0,
    DegreeOfMatrixGroup);

############################################################################
##
#M  DegreeOfMatrixGroup(<G>)
##
##  see IRREDSOL documentation
##
InstallMethod(DegreeOfMatrixGroup, "for matrix group with dimension", true,
    [IsMatrixGroup and HasDimension], 0,
    Dimension);

############################################################################
##
#M  IsIrreducible(<G>)
##
##  see IRREDSOL documentation
##
InstallMethod(IsIrreducible, "for matrix group", true, [IsMatrixGroup], 0,
    function(G)
        return IsIrreducibleMatrixGroup(G, FieldOfMatrixGroup(G));
    end);

############################################################################
##
#M  IsIrreducible(<G>, <F>)
##
##  see IRREDSOL documentation
##
InstallMethod(IsIrreducible, "for matrix group and field", IsMatGroupOverFieldFam,
    [IsMatrixGroup, IsField], 0,
    function(G, F)
        return IsIrreducibleMatrixGroup(G, F);
    end);

############################################################################
##
#M  IsIrreducibleMatrixGroup(<G>)
##
##
InstallOtherMethod(IsIrreducibleMatrixGroup, "for matrix group", true, [IsMatrixGroup], 0,
    function(G)
        return IsIrreducibleMatrixGroup(G, FieldOfMatrixGroup(G));
    end);

############################################################################
##
#M  IsIrreducibleMatrixGroup(<G>, <F>)
##
InstallMethod(IsIrreducibleMatrixGroupOp, "for matrix group and finite field - use MeatAxe",
    IsMatGroupOverFieldFam, [IsFFEMatrixGroup, IsField and IsFinite], 0,
    function(G, F)

        if not IsSubset(F, FieldOfMatrixGroup(G)) then
            Error("G must be a matrix group over F");
        fi;
        if DegreeOfMatrixGroup(G) = 1 then
            return true;
        elif IsTrivial(G) then
            return false;
        elif IsSubset(F, FieldOfMatrixGroup(G)) then
            return MTX.IsIrreducible (GModuleByMats (GeneratorsOfGroup(G), F));
        else
            Error("G must be a matrix group over F");
        fi;
    end);

############################################################################
##
#M  IsIrreducibleMatrixGroup(<G>, <F>)
##
InstallMethod(IsIrreducibleMatrixGroupOp, "for matrix group and finite field - test attr IsIrreducibleMatrixGroup",
    IsMatGroupOverFieldFam, [IsFFEMatrixGroup and HasIsIrreducibleMatrixGroup,
        IsField and IsFinite], 0,
    function(G, F)

        if not IsSubset(F, FieldOfMatrixGroup(G)) then
            Error("G must be a matrix group over F");
        fi;
        if IsIrreducibleMatrixGroup(G) then
            if F = FieldOfMatrixGroup(G) then
                return true;
            fi;
        elif IsSubset(F, FieldOfMatrixGroup(G)) then
            return false;
        fi;
        TryNextMethod();
    end);

############################################################################
##
#M  IsIrreducibleMatrixGroup(<G>, <F>)
##
InstallMethod(IsIrreducibleMatrixGroupOp, "for matrix group and finite field - for absolutely irreducible matrix group",
    IsMatGroupOverFieldFam, [IsFFEMatrixGroup and IsAbsolutelyIrreducibleMatrixGroup,
        IsField and IsFinite], RankFilter(HasIsIrreducibleMatrixGroup),
    function(G, F)

        if not IsSubset(F, FieldOfMatrixGroup(G)) then
            Error("G must be a matrix group over F");
        fi;
        return true;
    end);

############################################################################
##
#M  IsAbsolutelyIrreducible(<G>)
##
InstallMethod(IsAbsolutelyIrreducible, "for matrix group", true, [IsMatrixGroup], 0,
    function(G)
        return IsAbsolutelyIrreducibleMatrixGroup(G);
    end);

############################################################################
##
#M  IsAbsolutelyIrreducibleMatrixGroup(<G>)
##
##  we use `InstallOtherMethod' because otherwise there will be a warning
##  generated by KeyDependentOperation
##
InstallOtherMethod(IsAbsolutelyIrreducibleMatrixGroup, "for mat group over finite field", true,
    [IsFFEMatrixGroup], 0,

    function(G)

    local M;

    if DegreeOfMatrixGroup(G) = 1 then
        return true;
    elif IsTrivial(G) then
        return false;
    else
        M := GModuleByMats (GeneratorsOfGroup(G), DefaultFieldOfMatrixGroup(G));
        return MTX.IsIrreducible (M) and MTX.IsAbsolutelyIrreducible (M);
    fi;
end);

############################################################################
##
#M  IsPrimitive(<G>)
##
##  see IRREDSOL documentation
##
InstallMethod(IsPrimitive, "for matrix group", true, [IsMatrixGroup], 0,
    function(G)
        return IsPrimitiveMatrixGroup(G,  FieldOfMatrixGroup(G));
    end);

############################################################################
##
#M  IsPrimitiveMatrixGroup(<G>)
##
##  we need an OtherMethod to avoid the warning generated by
##  KeyDependentOperation
##
InstallOtherMethod(IsPrimitiveMatrixGroup, "for matrix group", true, [IsMatrixGroup], 0,
    function(G)
        return IsPrimitiveMatrixGroup(G, FieldOfMatrixGroup(G));
    end);

############################################################################
##
#M  IsPrimitive(<G>, <F>)
##
##  see IRREDSOL documentation
##
InstallMethod(IsPrimitive, "for matrix group over field",
    IsMatGroupOverFieldFam, [IsMatrixGroup, IsField], 0,
    function(G, F)

        return IsPrimitiveMatrixGroup(G, F);
    end);

############################################################################
##
#M  IsPrimitiveMatrixGroupOp(<G>, <F>)
##
##  see IRREDSOL documentation
##
InstallMethod(IsPrimitiveMatrixGroupOp, "for matrix group over finite field, construct IsomorphismPcGroup",
    IsMatGroupOverFieldFam,
    [IsFFEMatrixGroup and IsSolvableGroup, IsField and IsFinite], 0,

    function(G, F)
        local iso, inv;
        if not IsSubset(F, FieldOfMatrixGroup(G)) then
            Error("G must be a matrix group over F");
        fi;
        iso := IsomorphismPcGroup(G);
        inv := InverseGeneralMapping (iso);
        SetIsBijective (inv, true);
        SetIsGroupHomomorphism (inv, true);
        return SmallBlockDimensionOfRepresentation (ImagesSource (iso), inv, F, DegreeOfMatrixGroup(G)) = DegreeOfMatrixGroup(G);
    end);

############################################################################
##
#M  IsPrimitiveMatrixGroupOp(<G>, <F>)
##
##  see IRREDSOL documentation
##
InstallMethod(IsPrimitiveMatrixGroupOp, "for matrix group over finite field, use RepresentationIsomorphism",
    IsMatGroupOverFieldFam,
    [IsFFEMatrixGroup and HasRepresentationIsomorphism, IsField and IsFinite],
    RankFilter (IsHandledByNiceMonomorphism) + 1, # rank higher than the nice mono. method

    function(G, F)
        if not IsSubset(F, FieldOfMatrixGroup(G)) then
            Error("G must be a matrix group over F");
        fi;
        return SmallBlockDimensionOfRepresentation (
            Source (RepresentationIsomorphism (G)), RepresentationIsomorphism (G), F, DegreeOfMatrixGroup(G)) = DegreeOfMatrixGroup(G);
    end);

############################################################################
##
#M  IsPrimitiveMatrixGroupOp(<G>, <F>)
##
##  see IRREDSOL documentation
##
InstallMethod(IsPrimitiveMatrixGroupOp, "for matrix group over finite field, use nice monomorphism",
    IsMatGroupOverFieldFam,
    [IsFFEMatrixGroup and IsHandledByNiceMonomorphism, IsField and IsFinite],
    0,

    function(G, F)
        local iso, inv;
        if not IsSubset(F, FieldOfMatrixGroup(G)) then
            Error("G must be a matrix group over F");
        fi;
        iso := NiceMonomorphism (G);
        inv := GroupHomomorphismByFunction(NiceObject (G), G,
            h -> PreImagesRepresentative(iso, h),
            g -> ImageElm(iso, g));
        SetIsBijective (inv, true);
        return SmallBlockDimensionOfRepresentation (NiceObject (G), inv, F, DegreeOfMatrixGroup(G)) = DegreeOfMatrixGroup(G);
    end);

############################################################################
##
#M  IsPrimitiveMatrixGroupOp(<G>, <F>)
##
##  see IRREDSOL documentation
##
InstallMethod(IsPrimitiveMatrixGroupOp, "for matrix group over finite field, try if IsPrimitive is set",
    IsMatGroupOverFieldFam,
    [IsFFEMatrixGroup and HasIsPrimitive, IsField and IsFinite],
    RankFilter (IsHandledByNiceMonomorphism) + 3, # rank higher than the nice mono. method

    function(G, F)
        if not IsSubset(F, FieldOfMatrixGroup(G)) then
            Error("G must be a matrix group over F");
        fi;
        if IsPrimitive(G) then
            if FieldOfMatrixGroup(G) = F then
                return true;
            fi;
        else
            return false;
        fi;
        TryNextMethod();
    end);

############################################################################
##
#F  SmallBlockDimensionOfRepresentation(G, hom, F, limit)
##
##  hom must be a homomorphism G -> GL(n, F), where G is a group and F a finite
##  field such that Image(hom, G) is irreducible over F. limit is an integer
##  The function returns an integer k such that Im hom has a block system
##  of block dimension k, where k < limit, or k >= limit and G has no
##  block system of block dimesnion < limit
##
InstallGlobalFunction(SmallBlockDimensionOfRepresentation, function(G, hom, F, limit)

    # computes a block dimension smaller than limit, if it exists,
    # or the smallest block dimension otherwise
    local max, min, dim, M, m, cf, i;

    max := AttributeValueNotSet(MaximalSubgroupClassReps, G);
    min := DegreeOfMatrixGroup(Range (hom));
    for M in max do
        if not IsTrivial(M) then
            m := GModuleByMats (List(GeneratorsOfGroup(M), x -> ImageElm(hom, x)), F);
            if not MTX.IsIrreducible (m) then
                cf := First (MTX.CompositionFactors(m),
                    cf -> MTX.Dimension (cf) * IndexNC(G, M) = DegreeOfMatrixGroup(Range (hom))
                            and Length(MTX.Homomorphisms (cf, m)) > 0);
                if cf <> fail then
                    dim := SmallBlockDimensionOfRepresentation (M,
                        GroupHomomorphismByImagesNC(M, GL(MTX.Dimension (cf), Size(F)),
                        GeneratorsOfGroup(M), MTX.Generators (cf)),
                        F, limit);
                    if dim < min then
                        min := dim;
                        if min < limit then
                            return min;
                        fi;
                    fi;
                fi;
            fi;
        fi;
    od;
    return min;
end);

############################################################################
##
#M  MinimalBlockDimension(<G>)
##
##  see IRREDSOL documentation
##
InstallMethod(MinimalBlockDimension, "for matrix group", true, [IsMatrixGroup], 0,
    function(G)
        return MinimalBlockDimensionOfMatrixGroup(G);
    end);

############################################################################
##
#M  MinimalBlockDimension(<G>, <F>)
##
##  see IRREDSOL documentation
##
InstallMethod(MinimalBlockDimension, "for matrix group and field",
    IsMatGroupOverFieldFam, [IsMatrixGroup, IsField], 0,
    function(G, F)

        if not IsSubset(F, FieldOfMatrixGroup(G)) then
            Error("G must be a matrix group over F");
        fi;
        return MinimalBlockDimensionOfMatrixGroup(G, F);
    end);

############################################################################
##
#M  MinimalBlockDimensionOfMatrixGroup(<G>)
##
##  see IRREDSOL documentation
##
InstallOtherMethod(MinimalBlockDimensionOfMatrixGroup, "for matrix group", true,
    [IsMatrixGroup], 0,
    function(G)
        return MinimalBlockDimensionOfMatrixGroup(G, FieldOfMatrixGroup(G));
    end);

############################################################################
##
#M  MinimalBlockDimensionOfMatrixGroupOp(<G>, <F>)
##
##  see IRREDSOL documentation
##
InstallMethod(MinimalBlockDimensionOfMatrixGroupOp, "for matrix group over finite field",
    IsMatGroupOverFieldFam,
    [IsFFEMatrixGroup and IsSolvableGroup, IsField and IsFinite], 0,

    function(G, F)
        local iso, inv;
        if not IsSubset(F, FieldOfMatrixGroup(G)) then
            Error("G must be a matrix group over F");
        fi;

        if not IsIrreducibleMatrixGroup(G, F) then
            TryNextMethod();
        fi;

        iso := IsomorphismPcGroup(G);
        inv := InverseGeneralMapping (iso);
        SetIsBijective (inv, true);
        SetIsGroupHomomorphism (inv, true);
        return SmallBlockDimensionOfRepresentation (ImagesSource (iso), inv, F, 2);
    end);

############################################################################
##
#M  MinimalBlockDimensionOfMatrixGroupOp(<G>, <F>)
##
##  see IRREDSOL documentation
##
InstallMethod(MinimalBlockDimensionOfMatrixGroupOp,
    "for matrix group over finite field with representation homomorphism",
    IsMatGroupOverFieldFam,
    [IsFFEMatrixGroup and HasRepresentationIsomorphism, IsField and IsFinite],
        RankFilter (IsHandledByNiceMonomorphism) + 1, # rank higher than the nice mono. method

    function(G, F)
        if not IsSubset(F, FieldOfMatrixGroup(G)) then
            Error("G must be a matrix group over F");
        fi;
        if not IsIrreducibleMatrixGroup(G, F) then
            TryNextMethod();
        fi;
        return SmallBlockDimensionOfRepresentation (
            Source (RepresentationIsomorphism (G)), RepresentationIsomorphism (G), F, 2) ;
    end);

############################################################################
##
#M  MinimalBlockDimensionOfMatrixGroupOp(<G>, <F>)
##
##  see IRREDSOL documentation
##
InstallMethod(MinimalBlockDimensionOfMatrixGroupOp, "for matrix group over finite field, use NiceMonomorphism",
    IsMatGroupOverFieldFam,
    [IsFFEMatrixGroup and IsHandledByNiceMonomorphism, IsField and IsFinite],
    0,

    function(G, F)
        local iso, inv;
        if not IsSubset(F, FieldOfMatrixGroup(G)) then
            Error("G must be a matrix group over F");
        fi;
        if not IsIrreducibleMatrixGroup(G, F) then
            TryNextMethod();
        fi;
        iso := NiceMonomorphism (G);
        inv := GroupHomomorphismByFunction(NiceObject (G), G,
            h -> PreImagesRepresentative(iso, h),
            g -> ImageElm(iso, g));
        SetIsBijective (inv, true);
        return SmallBlockDimensionOfRepresentation (NiceObject (G), inv, F, 2);
    end);

############################################################################
##
#M  MinimalBlockDimensionOfMatrixGroupOp(<G>, <F>)
##
##  see IRREDSOL documentation
##
InstallMethod(MinimalBlockDimensionOfMatrixGroupOp,
    "for matrix group over finite field which has MinimalBlockDimension",
    IsMatGroupOverFieldFam,
    [IsFFEMatrixGroup and HasMinimalBlockDimension, IsField and IsFinite],
        RankFilter (IsHandledByNiceMonomorphism) + 1, # rank higher than the nice mono. method

    function(G, F)
        if not IsSubset(F, FieldOfMatrixGroup(G)) then
            Error("G must be a matrix group over F");
        fi;
        if F = FieldOfMatrixGroup(G) then
            return MinimalBlockDimension (G);
        elif MinimalBlockDimension (G) = 1 then
            return 1;
        fi;
        TryNextMethod();
    end);

############################################################################
##
#M  CharacteristicOfField(<G>)
##
##  see IRREDSOL documentation
##
InstallMethod(CharacteristicOfField, "for matrix group", true, [IsMatrixGroup], 0,
    M -> Characteristic (DefaultFieldOfMatrixGroup(M)));

############################################################################
##
#M  Characteristic(<G>)
##
##  see IRREDSOL documentation
##
InstallMethod(Characteristic, "for matrix group", true, [IsMatrixGroup], 0,
    CharacteristicOfField);

############################################################################
##
#M  RepresentationIsomorphism(<G>)
##
##  see IRREDSOL documentation
##
InstallMethod(RepresentationIsomorphism, "for mat group handled by nice mono.", true,
    [IsMatrixGroup and IsHandledByNiceMonomorphism], 0,
    function(G)

        local nice, H;

        nice := NiceMonomorphism (G);
        H := NiceObject(G);

        if IsSolvableGroup(H) then
            nice := IsomorphismPcGroup(G);
            H := Range (nice);
        fi;

        return GroupHomomorphismByFunction(H, G,
            x -> PreImagesRepresentative(nice, x),
            x -> ImageElm(nice, x));
    end);

############################################################################
##
#M  RepresentationIsomorphism(<G>)
##
##  see IRREDSOL documentation
##
InstallMethod(RepresentationIsomorphism, "soluble group: inverse of IsomorphismPcGroup", true,
    [IsMatrixGroup], 0,
    function(G)

        local nice;

        nice := IsomorphismPcGroup(G);
        return GroupHomomorphismByFunction(Range(nice), G,
            x -> PreImagesRepresentative(nice, x),
            x -> ImageElm(nice, x));
    end);

############################################################################
##
#F  ImprimitivitySystemsForRepresentation(G, rep, F, limit)
##
##  G is a group, F a finite field, rep: G -> GL(n, F)
##
##  If G has no block system with block dimension <= limit, the function
##  computes a list of all imprimitivity systems of Im rep as a
##  subgroup of GL(n, F). Otherwise, the function computes systems of imprimitivity,
##  one of which will have block dimension <= limit.
##
##  Each imprimitivity system is represented by a record with the following entries:
##  bases: a list of lists of vectors, each list of vectors being a basis of a block
##            in the imprimitivity system
##  stab1: the stabilizer in G of the first block (i. e., the block with basis bases[1])
##  min:    true if the block system is a minimal block system amongst the systems returned
##
InstallGlobalFunction(ImprimitivitySystemsForRepresentation, function(G, rep, F, limit)

    local systems, max, M, gens, m, c, cf, hom, subsys, sys, newsys, homBasis, homSpace, bas, bas2, newbasis, pos, orb;

    if DegreeOfMatrixGroup(Range (rep))< limit then
        return     [rec(bases := [IdentityMat (DegreeOfMatrixGroup(Range (rep)), F)], stab1 := G, min := true)];
    fi;
    systems := [];
    max := AttributeValueNotSet(MaximalSubgroupClassReps, G);
    for M in max do
        if not IsTrivial(M) then # inducing up from the trivial rep gives a reducible representation
            gens := List(GeneratorsOfGroup(M), x -> ImageElm(rep, x));
            m := GModuleByMats (gens, F);
            if not MTX.IsIrreducible (m) then
                for c in MTX.CollectedFactors(m) do
                    cf := c[1];
                    if MTX.Dimension (cf) * IndexNC(G, M) = Degree (Range (rep)) then
                        homBasis := MTX.Homomorphisms (cf, m);
                        if Length(homBasis) > 0 then # submodule isomorphic with cf

                            # get imprimitivity systems for cf

                            hom := GroupHomomorphismByImagesNC(M, GL(MTX.Dimension (cf), Size(F)),
                                GeneratorsOfGroup(M), MTX.Generators (cf));
                            subsys := ImprimitivitySystemsForRepresentation (M, hom, F, limit);
                            Add(subsys, rec(bases := [IdentityMat (MTX.Dimension (cf), F)], stab1 := M,
                                min := Length(subsys) = 0));

                            # translate result back

                            homSpace := VectorSpace(F, homBasis{[2..Length(homBasis)]}, 0*homBasis[1], "basis");
                            for bas2 in Enumerator(homSpace) do
                                bas := homBasis[1] + bas2;
                                  for sys in subsys do
                                    newbasis := List(sys.bases[1]*bas, ShallowCopy);
                                    TriangulizeMat (newbasis);
                                    if ForAll (systems, sys -> not newbasis in sys.bases) then
                                        orb := Orbit (ImagesSet(rep, G), newbasis, OnSubspacesByCanonicalBasis);
                                        Assert(1, Length(orb) * Length(newbasis) = Degree (Range (rep)));
                                        Assert(1, Length(orb) = IndexNC(G, sys.stab1));
                                        if orb[1] <> newbasis then
                                            pos := Position (orb, newbasis);
                                            orb{[1, pos]} := orb{[pos, 1]};
                                        fi;
                                        Add(systems, rec(bases := orb, stab1 := sys.stab1, min := sys.min));
                                    fi;
                                od;
                            od;
                        fi;
                    fi;
                od;
            fi;
        fi;
    od;

    # add trivial system

    Add(systems, rec(bases := [IdentityMat (Degree (Range (rep)), F)], stab1 := G, min := Length(systems) = 0));
    return systems;
end);

############################################################################
##
#A  ImprimitivitySystemsOp(<G>, <F>)
##
##  see IRREDSOL documentation
##
InstallMethod(ImprimitivitySystemsOp, "for matrix group handled by nice mono. and finite field",
    IsMatGroupOverFieldFam,
    [IsFFEMatrixGroup and IsHandledByNiceMonomorphism, IsField and IsFinite], 0,
    function(G, F)
        local rep;
        rep := RepresentationIsomorphism (G);
        return ImprimitivitySystemsForRepresentation (Source (rep), rep, F, 0);
    end);

############################################################################
##
#M  ImprimitivitySystemsOp(<G>, <F>)
##
##  see IRREDSOL documentation
##
InstallMethod(ImprimitivitySystemsOp, "for matrix group handled by nice mono. and finite field",
    IsMatGroupOverFieldFam,
    [IsFFEMatrixGroup and IsHandledByNiceMonomorphism, IsField and IsFinite], 0,
    function(G, F)
        local rep;
        rep := RepresentationIsomorphism (G);
        return ImprimitivitySystemsForRepresentation (Source (rep), rep, F, 0);
    end);

############################################################################
##
#M  ImprimitivitySystems(<G>)
##
##  see IRREDSOL documentation
##
InstallOtherMethod(ImprimitivitySystems, "for matrix group: use FieldOfMatrixGroup", true,
    [IsFFEMatrixGroup], 0,
    function(G)
        return ImprimitivitySystems (G, FieldOfMatrixGroup(G));
    end);

############################################################################
##
#M  MinimalBlockDimensionOfMatrixGroupOp(<G>, <F>)
##
##  see IRREDSOL documentation
##
InstallMethod(MinimalBlockDimensionOfMatrixGroupOp, "for matrix group having imprimitivity systems",
    IsMatGroupOverFieldFam,
    [IsFFEMatrixGroup and HasComputedImprimitivitySystemss, IsField],
    0,
    function(G, F)

        local known, i, sys, d, l;

        known := ComputedImprimitivitySystemss (G);
        for i in [1,3..Length(known)-1] do
            if known[i] = F then
                d := DegreeOfMatrixGroup(G);
                for sys in known[i+1] do
                    l := Length(sys.bases[1]);
                    if l < d then
                        d := l;
                    fi;
                od;
                return d;
            fi;
        od;
        TryNextMethod();
    end);

############################################################################
##
#M  TraceField(<G>)
##
##  see IRREDSOL documentation
##
InstallMethod(TraceField, "for irreducible matrix group over finite field", true,
    [IsFFEMatrixGroup and IsFinite], 1,

    function(G)

    local ext, gens, q, q0, module, module2, gens2, c, i, p;

    q := Size(FieldOfMatrixGroup(G));
    Info(InfoIrredsol, 3, "TraceField: matrix group is over GF(",q,")");

    if IsTrivial(G) then
        Info(InfoIrredsol, 4, "TraceField: trivial group case");
        return FieldOfMatrixGroup(G);
    fi;

    # guess a field contained in the smallest field over which module can be realised

    q0 := Size( Field (List([1..LogInt(Size(G), 2) + 10], i -> TraceMat(PseudoRandom(G)))));

    if q = q0 then
        Info(InfoIrredsol, 3, "TraceField: trace field is GF(", q0, ")");
        return FieldOfMatrixGroup(G);
    fi;

    Info(InfoIrredsol, 3, "TraceField: trace field contains GF(",q0, ")");

    module := GModuleByMats (GeneratorsOfGroup(G), FieldOfMatrixGroup(G));

    if not MTX.IsIrreducible (module) then
        Info(InfoIrredsol, 3, "TraceField: trace field contains GF(",q0, ")");
        TryNextMethod();
    fi;

    for c in Collected (Factors(LogInt(q, q0))) do
        p := c[1];
        for i in [1..c[2]] do

            Info(InfoIrredsol, 1, "TraceField: trying if trace field is GF(",q0, ")");
            # compute Galois conjugate of module
            gens2 := List(MTX.Generators (module), g -> List(g, row -> List(row, a -> a^q0)));
            module2 := GModuleByMats (gens2, MTX.Field (module));

            # If module1 and module2 are conjugate, their traces are the same, and thus
            # invariant under the Galois automorphism. Therefore the traces must belong to GF(q0).
            # Thus by a theorem of Brauer, module1 (and module2) can be written over GF(q0).
            # Conversely, if module1^x is over GF(q0) for some matrix x,
            # then module1^x = (module1^x)^q0 = (module1^q0)^(x^q0) = module2^(x^q0)
            # which shows that module1 and module2 are conjugate.

            # see also S. P. Glasby, R. B. Howlett, Writing representations over mnimal fields,
            # Comm. Alg. 25 (1997), 1703--1711

             if MTX.Isomorphism (module, module2) <> fail then
                Info(InfoIrredsol, 1, "TraceField: trace field is GF(",q0, ")");
                return GF(q0);
            fi;

            q0 := q0^p; # size of minimal superfield of GF(q0)
        od;
    od;

    Info(InfoIrredsol, 1, "TraceField: not rewritable over subfield: trace field is GF(",q0, ")");

    return GF(q0);
end);

############################################################################
##
#M  TraceField(<G>)
##
InstallMethod(TraceField, "generic method for finite matrix groups via conjugacy classes", true,
    [IsMatrixGroup and IsFinite], 0,
    function(G)

        local F, rep;

        F := FieldOfMatrixGroup(G);
        if IsPrimeField (F) then
            return F;
        elif F = Field (List(GeneratorsOfGroup(G), TraceMat)) then
            return F;
        else
            rep := RepresentationIsomorphism (G);
            return Field (List(ConjugacyClasses (Source (rep)), cl -> TraceMat(ImageElm(rep, Representative(cl)))));
        fi;
    end);

RedispatchOnCondition (TraceField, true, [IsMatrixGroup], [IsFinite], 0);

############################################################################
##
#M  SplittingField(<G>)
##
InstallMethod(SplittingField, "use MeatAxe", true,
    [IsFFEMatrixGroup], 0,
    function(G)

        local F, module;

        F := FieldOfMatrixGroup(G);
        module := GModuleByMats (GeneratorsOfGroup(G), F);
        if MTX.IsIrreducible (module) then
            return GF(Characteristic (F)^MTX.DegreeSplittingField (module));
        else
            Error("G must be irreducible over FieldOfMatrixGroup(G)");
        fi;
    end);

############################################################################
##
#M  ConjugatingMatTraceField(<G>)
##
##  returns a matrix x such that the matrix entries of G^x lie in the
##  trace field of G.
##
##  The absolutely irreducible case is an impelemntation of an algorithm by
##  S. P. Glasby, R. B. Howlett, Writing representations over mnimal fields,
##  Comm. Alg. 25 (1997), 1703--1711
##
InstallMethod(ConjugatingMatTraceField, "for irreducible FFE matrix group",
    true,
    [IsFFEMatrixGroup], 0,

    function(G)

        local ext, q, q1, t, C, CC, D, Y, A, i, j, mu, nu,
            basis, moduleG, module, module2, module3, absirred;

        if Length(GeneratorsOfGroup(G)) = 0 or TraceField(G) = FieldOfMatrixGroup(G) then
            return One(G);
        fi;

        q := Size(TraceField (G));

        moduleG := GModuleByMats (GeneratorsOfGroup(G), FieldOfMatrixGroup(G));

        # reduce to the absolutely irreducible case

        if not MTX.IsIrreducible (moduleG) then
            TryNextMethod();
        fi;

        absirred := MTX.IsAbsolutelyIrreducible (moduleG);
        if absirred then
            module := moduleG;
            ext := FieldOfMatrixGroup(G);
        else
            ext := GF(CharacteristicOfField (G)^MTX.DegreeSplittingField (moduleG));
            module := GModuleByMats (GeneratorsOfGroup(G), ext);
            repeat
                basis := MTX.ProperSubmoduleBasis (module);
                module := MTX.InducedActionSubmodule (module, basis);
            until MTX.IsIrreducible (module);
            Assert(1, MTX.IsAbsolutelyIrreducible (module));
        fi;

        # moduleG can be rewritten over TraceField(G) but over no proper subfield
        # let GF(q1) be the trace field of module, then over GF(q1),
        # moduleG has block diagonal matrices which are conjugate via
        # a Galois automorphism of GF(q1) of order Deg(G)/Dim(module)
        # so GF(q1) has dimension Dim(moduleG)/Dim(module) over TraceField(G).

        # Thus, taking q1-th powers is a generator of Gal (MTX.Field(module)/GF(q1))

        q1 := q^(MTX.Dimension (moduleG)/MTX.Dimension(module));
        t := LogInt(Size(ext), q1);

        module2 := GModuleByMats (List(MTX.Generators (module),
            A -> List(A, row -> List(row, x -> x^q1))), ext);

        C := MTX.Isomorphism (module, module2);

        if C = fail then
            Error("panic: cannot rewrite matrix group over trace field!");
        fi;

        CC := C;
        D := C;
        for i in [1..t-1] do
            CC := List(CC, row -> List(row, x -> x^q1));
            D := D*CC;
        od;

        mu := D[1][1];

        repeat
            nu := Random(ext);
        until Norm(ext, GF(q1), nu) = mu;

        C := nu^-1 * C;

        repeat
            Y := RandomMat(MTX.Dimension(module), MTX.Dimension(module), ext);
            A := Y;
            for i in [2..t] do
                A := Y + C * List(A, row -> List(row, x -> x^q1));
            od;
        until Length(NullspaceMat(A)) = 0;

        Assert(1, ForAll(MTX.Generators(module), g -> IsSubset(GF(q1), Field(Flat(g^A)))));

        # now we have a solution for the absolutely irreducible case

        if absirred then
            MakeImmutable(A);
            ConvertToMatrixRep(A, ext);
            return A;
        fi;

        basis := CanonicalBasis(AsVectorSpace(GF(q), GF(q1)));

        module3 := GModuleByMats(List(MTX.Generators(module), x -> BlownUpMat(basis, x^A)), MTX.Field(moduleG));

        A := Immutable(MTX.Isomorphism(moduleG, module3));

        if A = fail then
            Error("could not find conjugating matrix");
        fi;
        Assert(1, ForAll(GeneratorsOfGroup(G), g -> IsSubset(GF(q), Field(Flat(g^A)))));
        ConvertToMatrixRep(A, FieldOfMatrixGroup(G));

        return A;

    end);

############################################################################
##
#E
##

ext: gi


text:
#
# Read("~/Workspace/groupsSB/epi/A2/testLieAlgebraModuleBasis.gi");
#

type:="A";
rank:=2;
nr_pos_roots:=3;

Read("~/Workspace/groupsSB/epi/group.gi");

highest_weight:=[1,1];# corresponding highest weight module is the Lie algebra

Read("~/Workspace/groupsSB/epi/ext.gi");

Vb:=Basis(V);

#
# action of elements in the chevalley basis cb0 on module
#
on_module:=List(cb0,b->List(Vb,v->Coefficients(Vb,b^v)));

#
# action of elements in the chevalley basis cb0 on module
#
on_algebra:=List(cb0,b->List(cb0,v->Coefficients(cb0,v*b)));

#
# Wights of torus cb0[7],cb0[8] on module
#
weights_module:=[];
Add(weights_module,[on_module[7][1][1],on_module[8][1][1]]);
Add(weights_module,[on_module[7][2][2],on_module[8][2][2]]);
Add(weights_module,[on_module[7][3][3],on_module[8][3][3]]);

Add(weights_module,[on_module[7][4][4],on_module[8][4][4]]);
Add(weights_module,[on_module[7][5][5],on_module[8][5][5]]);

Add(weights_module,[on_module[7][6][6],on_module[8][6][6]]);
Add(weights_module,[on_module[7][7][7],on_module[8][7][7]]);
Add(weights_module,[on_module[7][8][8],on_module[8][8][8]]);

#
# Let a=cb0[1] and b=cb0[2]
#

#
# Vb[1] corrsponds to a+b (heighest root) ---> cb0[3]
#

#
# Vb[2] and Vb[3] corrspond to a and b
#

#
# In the internal notation
#   y1 = cb0[4] corrseponding to -a corresponding to cb0[7]
#   y2 = cb0[5] corresponding to -b corresponding to cb0[8]
#   y3 = cb0[6] corresponding to -a-b

max_stars_repo_name: iuliansimion/groupsSB

max_issues_repo_name: iuliansimion/groupsSB

max_forks_repo_name: iuliansimion/groupsSB


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

ImportAll(paradigms.scratchpad);

Class(ScratchX86CMContext, rec (
	getOpts := meth(arg)
		local opts, elem, lssize, nrsgmts, vlen, size, swp, globalUnrolling, ttype;

		lssize := When (Length(arg) >= 2, arg[2], 2);
		nrsgmts := When (Length(arg) >= 3, arg[3], 1);
		vlen := When (Length(arg) >= 4, arg[4], 1);
		size:= When (Length(arg) >= 5, arg[5], 2);
         	ttype := When (Length(arg) >= 6, arg[6], 'R');
		swp := When (Length(arg) >= 7, arg[7], false);
		globalUnrolling := When(Length(arg) >=8, arg[8], 1);
		elem := ScratchpadGlobals.getOpts(lssize,nrsgmts,vlen,size,ttype,swp,globalUnrolling);

        	opts := CopyFields(elem);

		opts.unparser := CCMContextScratchUnparserProg;
		opts.codegen := CMContextScratchCodegen;

        	opts.profile.makeopts.CFLAGS := "-O2 -Wall -fomit-frame-pointer -msse4.1 -std=gnu99 -static";

        	opts.includes := [ "<include/omega64.h>" ];
        	Add(opts.includes, "\"scratchc.h\"");

		return opts;
	end,
));

Class(ScratchX86Globals, rec(
	getOpts := meth(arg)
		local opts, elem, lssize, nrsgmts, vlen, size, swp, globalUnrolling, ttype;

		lssize := When (Length(arg) >= 2, arg[2], 2);
		nrsgmts := When (Length(arg) >= 3, arg[3], 1);
		vlen := When (Length(arg) >= 4, arg[4], 1);
		size:= When (Length(arg) >= 5, arg[5], 2);
         	ttype := When (Length(arg) >= 6, arg[6], 'R');
		swp := When (Length(arg) >= 7, arg[7], false);
		globalUnrolling := When(Length(arg) >=8, arg[8], 1);
		elem := ScratchpadGlobals.getOpts(When(swp,lssize/2,lssize),nrsgmts,vlen,size,ttype,swp,globalUnrolling);

        opts := CopyFields(elem);

        opts.swp_var := swp;

        opts.profile.name := "linux-x86-thread";

		opts.unparser := BarrierScratchUnparserProg;
		opts.codegen := SWPBarrierScratchCodegen;

        opts.register := (self, opts) >> "REG_thread";

        opts.profile.makeopts.CFLAGS := "-O2 -Wall -fomit-frame-pointer -msse4.1 -std=gnu99 -static -lpthread";

        opts.barrierCMD := (self, opts) >> "BARRIER";
        opts.initialization := (self, opts) >> "INITIAL";

        opts.includes := [];
        Add(opts.includes, "\"scratch_barrier.h\"");
		Add(opts.includes, "<include/omega64.h>");

		return opts;
	end,
));


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

Declare(TDA, _ApplyAllRulesTopDown, limit_cx);

IsRewriteRule := r -> IsRec(r) and IsBound(r.isRewriteRule) and r.isRewriteRule;

Class(RewriteRuleBase, rec(
	isRewriteRule := true,

	regularize := self >> Error("Not implemented. This method needs to convert ",
	"a special kind of rewrite rule into a regular RewriteRule"),
));

#F RewriteRule(<rewrite-rule>)	-- returns <rewrite-rule> object unchanged
#F RewriteRule(<from>, <to_func>) -- creates an unnamed rewrite rule
#F RewriteRule(<from>, <to_func>, <name>) -- creates a rewrite rule
#F
Class(RewriteRule, RewriteRuleBase, rec(
	regularize := self >> self,

	__call__ := meth(arg)
		local op, from, to, name, self;

		self := arg[1];
		arg := Drop(arg, 1);

		if Length(arg)=1 and IsList(arg[1]) then
			arg := arg[1];
		fi;
		if Length(arg)=1 then
			return Checked(IsRewriteRule(arg[1]), arg[1]);
		elif Length(arg)=2 then
			from := arg[1];
			to   := Checked(IsCallableN(arg[2], 1) or IsCallableN(arg[2], 2), arg[2]);
			name := "unnamed";
		elif Length(arg)=3 then
			from := arg[1];
			to   := Checked(IsCallableN(arg[2], 1) or IsCallableN(arg[2], 2), arg[2]);
			name := Checked(IsString(arg[3]), arg[3]);
		else
			Error("Rule requires 1, 2 or 3 arguments. See Doc(ARule)");
		fi;

		if IsList(from) then
			if Length(from)=0 then
				Error("Rule must have non-empty left-side");
			fi;
			op := from[1];
		else
			op := from;
		fi;
		if Is@(op) or not IsRec(op) or not IsBound(op.name) then
			if IsBound(op._target) and Length(op._target)=1 then
				op := op._target[1];
			else
				op := @;
			fi;
		fi;
		return WithBases(self, rec(op := op, from := from, to := to, name := name, operations := PrintOps));
	end,

	print := self >> PrintEval("$1($2, $3, \"$4\")", self.__name__, self.from, self.to, self.name),

));

# Rule(..), alias for RewriteRule, see Doc(RewriteRule)
#
Rule := arg -> ApplyFunc(RewriteRule, arg);

#F ARule(<rewrite-rule>) - returns <rewrite-rule> object unchanged
#F ARule(<op>, <from>, <to_func>) - creates an unnamed rewrite rule for associative <op>
#F ARule(<op>, <from>, <to_func>, <name>) - creates a rewrite rule for associative <op>
#F
#F Create a rewrite rule for an associative operator
#F
Class(AssociativeRewriteRule, RewriteRuleBase, rec(
	__call__ := meth(arg)
		local op, from, to, name, lhs, rhs, self;

		self := arg[1];
		arg := Drop(arg, 1);

		if Length(arg)=1 and IsList(arg[1]) then
			arg := arg[1];
		fi;

		if Length(arg)=1 then
			return Checked(IsRewriteRule(arg[1]), arg[1]);
		elif Length(arg)=3 then
			op   := Checked(IsClass(arg[1]), arg[1]);
			from := Checked(IsList(arg[2]),  arg[2]);
			to   := Checked(IsCallableN(arg[3], 1) or IsCallableN(arg[3], 2), arg[3]);
			name := "unnamed";
		elif Length(arg)=4 then
			op   := Checked(IsClass(arg[1]), arg[1]);
			from := Checked(IsList(arg[2]),  arg[2]);
			to   := Checked(IsCallableN(arg[3], 1) or IsCallableN(arg[3], 2), arg[3]);
			name := Checked(IsString(arg[4]), arg[4]);
		else
			Error("ARule requires 1, 3 or 4 arguments. See Doc(ARule)");
		fi;

		return WithBases(self, rec(op := op, from := from, to := to, name := name, operations := PrintOps));
	end,

	regularize := meth(self)
		local op, from, to, lhs, rhs;

		[op, from, to] := [self.op, self.from, self.to];
		lhs := [op, ...] :: from :: [...];
		if NumArgs(to)=1 then
			rhs := DetachFunc(Subst(
				function(e)
				local rch;
				rch := _children(e);
					# ....left means O.left, where O is the '...' object, defined in rules.gi
				return _fromChildren(e, Concatenation(
				rch{[1 .. ....left]},
				$to(e),
				rch{[....right .. Length(rch)]}));
			end));
		else
			rhs := DetachFunc(Subst(
			function(e,cx)
				local rch;
				rch := _children(e);
					 # ....left means O.left, where O is the '...' object, defined in rules.gi
				return _fromChildren(e, Concatenation(
					rch{[1 .. ....left]},
					$to(e,cx),
					rch{[....right .. Length(rch)]}));
			end));
		fi;
		return RewriteRule(lhs, rhs, self.name);
	end,

	print := self >> PrintEval("$1($2, $3, $4, \"$5\")", self.__name__, self.op, self.from, self.to, self.name),
));

# ARule(..), alias for AssociativeRewriteRule, see Doc(AssociativeRewriteRule)
#
ARule := arg -> ApplyFunc(AssociativeRewriteRule, arg);

_ARule_Transparent := (op, transp_ops, from, to) -> Checked(
	IsClass(op), IsList(transp_ops), ForAll(transp_ops, IsClass),
	IsList(from), IsCallableN(to, 1) or IsCallableN(to, 2),
	List(transp_ops, tr_op ->
	ARule(op, List(from, f -> [tr_op, f]),
		When(IsCallableN(to, 1), Subst(x -> List($to(x), y->$tr_op(y))),
							 Subst((x,cx) -> List($to(x, cx), y->$tr_op(y))))))
);

IsRuleSet := x -> IsRec(x) and IsBound(x.isRuleSet) and x.isRuleSet;

Class(RuleSet, rec(
	isRuleSet	:= true,
	rules		:= rec(),
	# counts how many time this RuleSet changed
	_changes	 := 0,
	# change index of compiled ruleset, need to recompile if
	# '_comp_id' is not equal to '_changes'
	_comp_id	 := -1,

	__transparent__ := [],

	addRules := meth(self, rules)
		local f, rr;
		Constraint(IsRec(rules));
		self.rules := ShallowCopy(self.rules);
		rules := ShallowCopy(rules);

		for f in UserRecFields(rules) do
			rr := rules.(f);
			Constraint(IsRewriteRule(rr) or (IsList(rr) and ForAll(rr, IsRewriteRule)));
			rr := When(IsList(rr), rr, [rr]);
			rr := ConcatList(rr, r ->
			Cond(r _is AssociativeRewriteRule,
				[r] ::
				ConcatList(Filtered(self.__transparent__, t->t[1]=r.op),
				t -> _ARule_Transparent(r.op, t[2], r.from, r.to)),
				[r]));
			rules.(f) := rr;
		od;

		MergeIntoRecord(self.rules, rules);

		# ruleset changed
		self._changes := self._changes + 1;
		return self;
	end,

	compileRule := meth(self, name, rule)
		local op, from, to_func, owner, head;
		owner := rule.owner;
		rule := rule.regularize();
		[op, from, to_func] := [rule.op, rule.from, rule.to];
		op := op.__name__;
		head := When(IsList(from) and from<>[], from[1], from);

		if op = "@" and IsRec(head) and IsBound(head._target) then
			DoForAll(head._target, t ->
				self.compileRule(name, CopyFields(rule, rec(op := t))));
		else
			if not IsBound(self._compiled.(op)) then
				self._compiled.(op) := [ CopyFields(rule, rec(name := name)) ]; #[from, to_func, name, owner] ];
			else
				Add(self._compiled.(op), CopyFields(rule, rec(name := name))); #[from, to_func, name, owner]);
			fi;
		fi;
	end,

	compile := meth(self)
		self._compiled := tab();
		# save the owner ruleset, it will be used to produce warnings of duplicate
		# rule definitions, when dealing with merged rule sets
		DoForAll(self.rules, function(name, rule)
			local r;
			if IsRewriteRule(rule) then rule.owner := self;
			elif IsList(rule) and not IsString(rule) then  for r in rule do r.owner := self; od;
			fi;
		end);

		DoForAll(self.rules, (name, rule) ->
			Cond(IsRewriteRule(rule), self.compileRule(name, rule),
				IsList(rule) and not IsString(rule), DoForAll(rule, r -> self.compileRule(name, r)), 0)
		);
		# remember change number
		self._comp_id := self._changes;
		return self;
	end,

	__call__ := (self, s) >> TDA(s, self, rec()),

	apply1 := (self, s) >> _ApplyAllRulesTopDown(s, limit_cx(1), self),

	get_changes  := self >> self._changes,
	get_comp_ids := self >> self._comp_id,

	changed  := self >> self.get_comp_ids() <> self.get_changes(),

	compiled := self >> When(IsBound(self._locked) or not self.changed(), self, self.compile())._compiled,

	# locked ruleset doesn't check changes and doesn't recompile rules when asked for compiled() table
	locked   := self >> CopyFields(When(self.changed(), self.compile(), self), rec(_locked := true)),
));

Class(EmptyRuleSet, RuleSet);

AnonRuleSet := rules -> WithBases(RuleSet, rec()).addRules(rules);

Class(MergedRuleSet, RuleSet, rec(
	 _mergedConflictWarnings := [],
	 warnConflicts := false,
	 warnDuplicates := true,

	 __call__ := arg >> let(
		 self	 := arg[1],
		 rulesets := Drop(arg, 1),
		 comp_id  := Sum(rulesets, e -> e.get_changes())-1,
		 Checked(Length(rulesets) >= 1,
			 When( Length(rulesets)=1,
				 rulesets[1],
				 WithBases(self, rec(operations := PrintOps,
					 __call__ := RuleSet.__call__,
					 rulesets := rulesets,
					 _comp_id := comp_id))))),

	 get_changes  := self >> Sum(self.rulesets, e -> e.get_changes()),
	 get_comp_ids := self >> self._comp_id,

	 checkConflicts := meth(self, rules1, rules2, rulesets)
		 local conflicts, c;
		 conflicts := Intersection(UserRecFields(rules1), UserRecFields(rules2));
		 # make sure the origins of the rules are different
		 # origin == RuleSet where the rule first appears, i.e., not the MergedRuleSet
		 conflicts := Filtered(conflicts, c -> rules1.(c).owner.__name__ <> rules2.(c).owner.__name__);

		 # prevent duplicate warnings
		 if not self.warnDuplicates then
			 SubtractSet(conflicts, self._mergedConflictWarnings); fi;

		 if conflicts <> [] then
			 PrintErr("Warning: rewrite rule conflict when merging ", rulesets, "\n");
			 for c in conflicts do
				 Add(self._mergedConflictWarnings, c);
				 PrintErr(Blanks(9), c, "  ", rules1.(c).owner, "<->", rules2.(c).owner, "\n");
			 od;
		 fi;
	 end,

	 compile := meth(self)
		local rules, r;
		rules := rec();
		for r in self.rulesets do
			r.compile();
			if self.warnConflicts then
				self.checkConflicts(rules, r.rules, self.rulesets);
			fi;
			MergeIntoRecord(rules, r.rules);
		od;
		self.rules := rules;

		self._compiled := tab();

		DoForAll(self.rules, (name, rule) ->
			Cond(IsRewriteRule(rule), self.compileRule(name, rule),
				IsList(rule) and not IsString(rule),
				DoForAll(rule, r -> self.compileRule(name, r)),0)
		);

		self._comp_id := self.get_changes();
		return self;
	end,

	print := self >> Print(self.name, "(", PrintCS(self.rulesets), ")")
));

RewriteRules := (rule_set, rules) -> rule_set.addRules(rules);

_LookupRules := (expr, ruleset) -> let(
	name := ObjId(expr).__name__,
	R := When(IsBound(ruleset._locked), ruleset._compiled, ruleset.compiled()),
	When(IsBound(R.(name)), R.(name), []) :: When(IsBound(R.@), R.@, []));

#_AddRule2 := function(op, from, to_func)
#	 local rules;
#	 rules := _AddRule_old(op, from, to_func);
#	 if not IsBound(op._rules) then op._rules := BagAddr(rules); fi;
#end;
#_LookupRules2 := op -> When(IsBound(op._rules), BagFromAddr(op._rules), []);

ChkSPL := function(a,b, rule)
	local diff;
#	Error();
	if spiral.spl.IsSPL(a) and spiral.spl.IsSPL(b) then
		diff := spiral.code.InfinityNormMat(spiral.spl.MatSPL(a) -spiral.spl.MatSPL(b));
		if diff > 1E-4 then
			Print("Broken SPL rule: ", rule, "\nold:\n", a, "\nnew:\n", b);
			Error("bad SPL rule");
		fi;
	else
		return true;
	fi;
end;

RuleTrace := Ignore;
RuleStrategyTrace := Ignore;
RuleStrategyTiming := Ignore;
RuleStatus := Ignore;
RuleCheckSPL := false;

apply_rules := function(rules, expr, context)
	local rule, lhs, rhs, old;
	for rule in rules do
		lhs := rule.from;
		rhs := rule.to;
		while PatternMatch(expr, lhs, context) and context.rlimit <> 0 do
			context.rlimit := context.rlimit - 1;
			context.applied := context.applied + 1;
			old := Copy(expr);
			RuleTrace(rule);
			RuleStatus(rule, "OLD: ", [expr, "\n"]);
			if RuleCheckSPL then old := Copy(expr); fi;
			if NumArgs(rhs)=1 then expr := rhs(expr);
			else expr := rhs(expr, context);
			fi;
			RuleStatus(rule, "NEW: ", [expr, "\n"]);
			trace_log.addRewrite(rule.name,old,expr, []);
			if RuleCheckSPL then ChkSPL(old, expr, rule); fi;
		od;
	od;
	return expr;
end;

# non-iterative version (no more: while PatternMatch(...) do ...)
apply_rules_ni := function(rules, expr, context)
	local rule, lhs, rhs, old;
	for rule in rules do
		lhs := rule.from;
		rhs := rule.to;
		if PatternMatch(expr, lhs, context) and context.rlimit <> 0 then
			context.rlimit := context.rlimit - 1;
			context.applied := context.applied + 1;
			old := Copy(expr);
			RuleTrace(rule);
			RuleStatus(rule, "OLD: ", [expr, "\n"]);
			if NumArgs(rhs)=1 then expr := rhs(expr);
			else expr := rhs(expr, context);
			fi;
			RuleStatus(rule, "NEW: ", [expr, "\n"]);
			trace_log.addRewrite(rule.name,old,expr, []);
		fi;
	od;
	return expr;
end;

cx_enter := function(cx, expr)
	local opname;
	if IsRec(expr) and IsBound(expr.name) then
	opname := expr.name;
	if not IsBound(cx.(opname)) then cx.(opname) := [ expr ];
	else Add(cx.(opname), expr); fi;
	fi;
	Add(cx.parents, expr);
end;

cx_leave := function(cx, expr)
	# unupdate context back to original
	if IsRec(expr) and IsBound(expr.name) then
	RemoveLast(cx.(expr.name), 1);
	fi;
	RemoveLast(cx.parents, 1);
end;

map_children := function(expr, to_func)
	local ch, i;
	ch := _children(expr);
	for i in [1..Length(ch)] do
		_setChild(expr, i, to_func(ch[i]));
	od;
	return expr;
end;

map_children_safe := function(expr, to_func)
	local ch, i;
	ch := ShallowCopy(_children(expr));
	if ch=[] then return expr; fi;
	for i in [1..Length(ch)] do
		ch[i] := to_func(ch[i]);
	od;
	return _fromChildren(expr, ch);
end;

# Rule application context fields
#	parents - list of parents of current node,
#			  immediate parent is last, farthest parent is first,
#			  this field is updated by apply_rules()
#	rlimit  - maximum number of rules to apply, this is a parameter
#			  to apply_rules()
#	applied - number of rules applied so far, updated by apply_rules()
#
#	<objid> - for each <objid> list of parents with that id only,
#			  immediate parent is last. For instance Last(context.ISum)
#			  is the enclosing ISum.
#

_meth_cx_isInside := (self, oid) >> let(nam := Cond(IsString(oid), oid, oid.__name__),
	IsBound(self.(nam)) and self.(nam)<>[]
);

# construct an empty context, but set a limit on the number of applied rules
limit_cx := lim -> tab(
	isInside := _meth_cx_isInside,
	parents := [], rlimit := lim, applied := 0
);

# construct an empty initial context
empty_cx := () -> tab(
	isInside := _meth_cx_isInside,
	parents := [], rlimit := -1, applied := 0
);

_ApplyAllRulesTopDown := function(expr, context, ruleset)
	if (not IsRec(expr) or not IsBound(expr.name)) and (not IsList(expr) or BagType(expr) in [T_STRING, T_RANGE]) then
		return expr;
	fi;
	if IsBound(ruleset.__avoid__) and ObjId(expr) in ruleset.__avoid__ then
		return expr;
	fi;
	# apply rules
	expr := apply_rules(_LookupRules(expr, ruleset), expr, context);
	# recurse
	# NOTE: do not enter context if expr has no children!
	cx_enter(context, expr);
	expr := map_children(expr, c -> _ApplyAllRulesTopDown(c, context, ruleset));
	cx_leave(context, expr);
	return expr;
end;

ApplyAllRulesTopDown := (expr, ruleset) ->
	_ApplyAllRulesTopDown(expr, empty_cx(), ruleset.locked());

_ApplyAllRulesBottomUp := function(expr, context, ruleset)
	if (not IsRec(expr) or not IsBound(expr.name)) and (not IsList(expr) or BagType(expr) in [T_STRING, T_RANGE]) then
		return expr;
	fi;
	if IsBound(ruleset.__avoid__) and ObjId(expr) in ruleset.__avoid__ then
		return expr;
	fi;
	# recurse
	cx_enter(context, expr);
	expr := map_children(expr, c -> _ApplyAllRulesBottomUp(c, context, ruleset));
	cx_leave(context, expr);
	# apply rules
	expr := apply_rules(_LookupRules(expr, ruleset), expr, context);
	return expr;
end;

ApplyAllRulesBottomUp := (expr, ruleset) ->
	_ApplyAllRulesBottomUp(expr, empty_cx(), ruleset.locked());

_apply_strategy_step := function(expr, rset, apply_func, opts)
	if IsFunc(rset) then
		if NumArgs(rset)=2 then
			expr := rset(expr, opts);
		else
			expr := rset(expr);
		fi;
	else
		expr := apply_func(expr, rset.locked(), opts);
	fi;
	return expr;
end;

##
## ApplyStrategy(<expr>, <list-of-rulesets>, <apply_func>, <opts>)
##
## <expr>			 - expression to transform
## <list-of-rulesets> - list of rulesets to apply
## <apply_func>	   - of type (expr, rset) -> expr
##
ApplyStrategy := function(expr, strategy, apply_func, opts)
	local rset, l, r, i, t;
	# a hack to make this function reentrant
	# (e.g. to make it possible to call this function from within rhs of a rule)
	l := ....left; r := ....right; i := 1;
	for rset in strategy do
		RuleStrategyTrace(i, rset, expr);
	[expr, t] := UTimedAction(_apply_strategy_step(expr, rset, apply_func, opts));
		RuleStrategyTrace(i, "DONE", expr);
	RuleStrategyTiming(i, t);
		i := i + 1;
	od;
	....left := l;
	....right := r;
	return expr;
end;

TD := ApplyAllRulesTopDown;
BU := ApplyAllRulesBottomUp;

TDA := function(s, ruleset, opts)
	local cx;
	cx := empty_cx();
	cx.opts := opts;
	cx.applied := 1;
	while cx.applied > 0 do
		cx.applied := 0;
		s := _ApplyAllRulesTopDown(s, cx, ruleset);
	od;
	return s;
end;

BUA := function(s, ruleset, opts)
	local cx;

	cx := empty_cx();
	cx.opts := opts;
	cx.applied := 1;
	while cx.applied > 0 do
		cx.applied := 0;
		s := _ApplyAllRulesBottomUp(s, cx, ruleset);
	od;
	return s;
end;

UntilDone := TDA;

Rewrite := function(expr, rset, opts)
	local res;

	if IsList(rset) then
		return FoldL(rset, (a,b)->Rewrite(a, b, opts), expr);
	fi;

	trace_log.beginRuleset(rset, expr);
	res := UntilDone(expr, rset, opts);
	trace_log.endRuleset(rset, res);

	return res;
end;

#F TDA_1by1(<s>, <expr>, <after_each_rule>)
#F
#F This function applies rules in a top-down fashion to <s> until
#F it converges, just like TDA.  However after each application of the
#F rule it runs after_each_rule(s) on the newly transformed <s>.
#F This allows stepwise verification among other things.
#F
#F Example:
#F	r := RandomRuleTree(DFT(4));
#F	# following line (unlike SumsRuleTree) gives non-simplified Sigma-SPL
#F	s := SumsSPL(SPLRuleTree(r));
#F	TDA_1by1(s, MergedRuleSet(RulesSums, RulesDiag, RulesFuncSimp),
#F				x->PrintLine(MatSPL(x)-MatSPL(DFT(4))));
#F
TDA_1by1 := function(s, ruleset, after_each_rule)
	local cx;

	cx := limit_cx(1); cx.applied := 1;
	while cx.applied > 0 do
		cx := limit_cx(1); cx.applied := 0;
		s := _ApplyAllRulesTopDown(s, cx, ruleset);
		after_each_rule(s);
	od;
	return s;
end;