text:
# Copyright 2018-2019, Carnegie Mellon University
# See LICENSE for details

Class(cos, AutoFoldExp, rec(
  computeType := (self) >> TReal,

));

Class(sin, AutoFoldExp, rec(
  computeType := (self) >> TReal,
));

Class(tan, AutoFoldExp, rec(
  computeType := (self) >> TReal,
));

Class(TIVReal, AtomicTyp, rec(
    hash := (_val, size) -> let(val := When(IsList(_val), _val, [_val, _val]),
        h := DoubleRep64(val[1]) + DoubleRep64(val[2]),
        1 + (h mod size)),
    check := v -> Cond(IsDouble(v) or IsList(v) and Length(v) = 2 and v[1] <= v[2], v, Error("<v> must be a Double or an ordered list of length 2")),
    realType    := self >> self,
    print := self >> Print(self.__name__, "(", self.t, ")"),
    base_t := self >> self.t,
    zero := self >> self.value(0.0),
    one := self >> self.value(1.0),
    __call__ := (self, t) >>
        WithBases(self, rec(
        t    := Checked(IsType(t), t),
        operations := TypOps))
));

Class(TIVBool, AtomicTyp, rec(
    hash := (val, size) -> 1 + (InternalHash(val) mod size),
    check := v -> Cond(IsInt(v) and v in [0, 1, -1], v, Error("<v> must be in [-1, 0, 1]")),
    zero := self >> self.value(0),
    one := self >> self.value(1)
));

Class(ivenv, chain);

Class(RulesHCOLIVArithType, RuleSet, rec(inType := "iCode", outType := "iCode"));
RewriteRules(RulesHCOLIVArithType, rec(
	ivenv_prop_types_bools := ARule(ivenv, [[@(1, assign), @(2,var).cond(e->e.t<>TInt),
		@(3).cond(e-> e.t<>@(2).val.t and e.t=TBool)], @(4)],
        e->let(
		nvar := var.fresh_t("n", TInt),
		#Print("New Type 2: ", @(1).val,"(",@(2).val, ",",@(3).val,")\n", @(2).val.t, " <> ", @(3).val.t, "\n New Var Type: ", nvar, ",",nvar.t," \n"),
	  [decl([nvar], SubstVars(chain(assign(nvar, @(3).val), @(4).val), rec((@(2).val.id) := nvar )))] ) ),

    ivenv_prop_types := ARule(ivenv, [[@(1, assign), @(2,var),
		@(3).cond(e-> e.t<>@(2).val.t and not (e.t in [TVect(T_Real(64), 2), TIVReal, TBool]))], @(4,...)],
        e->let(
		nvar := var.fresh_t("n", @(3).val.t),
		#Print("New Type: ", @(1).val,"(",@(2).val, ",",@(3).val,")", @(2).val.t, " <> ", @(3).val.t, " ", nvar, "\n", @(4).val, "\n"),
	  [decl([nvar], SubstVars(chain(assign(nvar, @(3).val), @(4).val), rec((@(2).val.id) := nvar )))] ) ),

));

Class(RulesHCOLIVArith, RuleSet, rec(inType := "iCode", outType := "iCode"));
RewriteRules(RulesHCOLIVArith, rec(
    ivenv_chain := Rule(@(1, ivenv, e->ForAny(e.cmds, f->ObjId(f) = chain)),
        e->ApplyFunc(ivenv, Flat(List(@(1).val.cmds, j->When(ObjId(j)=chain, j.cmds, j))))),
    ivenv_xyz_decl := ARule(ivenv, [@(1), @(2, decl)],
        e->[decl(@(2).val.vars, chain(@(1).val, @(2).val.cmd))]),
    ivenv_decl := Rule(@(1, ivenv, e->ObjId(e.cmds[1])=decl),
        e -> decl(@(1).val.cmds[1].vars, ivenv(Concat([@(1).val.cmds[1].cmd], Drop(@(1).val.cmds, 1))))),
    ivenv_drop_selfassign := ARule(ivenv, [@(1, assign, e->IsVar(e.loc) and IsVar(e.exp) and e.exp=e.loc), @(2)],
        e->[@(2).val]),
    ivenv_creturn := Rule(@(1, ivenv, e->ObjId(Last(e.cmds))=creturn),
        e->chain(ivenv(DropLast(@(1).val.cmds, 1)), Last(@(1).val.cmds))),
    loop_pull_const_array_tcast := Rule([@(1, loop), @(2),  @(3), [@(4, chain), [@(5, assign), @(6,var), [@(7, tcast), @(8), @(9,nth, e->IsValue(e.idx))]],...]],
        e->let(#Print(e, "\n"),
		chain(@(5).val, loop(@(1).val.var, @(1).val.range, chain(Drop(@(4).val.cmds,1)))))),
    ivenv_state := Rule([@(1, tcast), @(2), @@(3, nth, (e, cx)->IsBound(cx.opts.state) and e.loc.id=cx.opts.state.id)],
        (e, cx)->nth(cx.opts.ivstate, e.args[2].idx) )
));

RulesCodeUnrollHACMIVArth := CopyFields(MergedRuleSet( RulesHCOLIVArith, RulesUnrollHCOL, RulesStrengthReduce,  RulesCodeHCOL),
    rec(inType:="iCode", outType := "iCode"));

RulesCodeHACMIVArth := CopyFields(MergedRuleSet(RulesStrengthReduce, RulesCodeHCOL, RulesHCOLIVArith),
    rec(inType:="iCode", outType := "iCode"));

RulesCodeUnrollHACMIVArthType := CopyFields(MergedRuleSet(RulesCodeUnrollHACMIVArth, RulesHCOLIVArithType),
	rec(inType :="iCode", outType:="iCode"));

Class(IVArithMixin, rec(
    precision := "double",
    TRealCtype := "double",
    IVRealType := TIVReal(T_Real(64)),
    IVBoolType := TInt,
    IVValueType := TIVReal(T_Real(64)),
    codeRuleSet := RulesCodeUnrollHACMIVArth
));

Class(ivExp, errExp);

Class(RealEPS, errExp);

Class(RealMPS, errExp);

Class(IVArith, SumsBase, BaseContainer, rec(
    rng := meth(self) return self._children[1].rng(); end,
    dmn := meth(self) return self._children[1].dmn(); end,
    toOperator := self >> (vec -> self.child(1).toOperator()(vec))
));

Class(TIVArith, Tagged_tSPL_Container, rec(
    abbrevs :=  [ s -> Checked(IsSPL(s), [s]) ],
    transpose := self >> CopyFields(self, rec(transposed := not self.transposed)),
    dims := self >> [self.params[1].dims()[1], self.params[1].dims()[2]],
    isReal := True,
    toOperator := self >> (vec -> let(r := self.params[1].toOperator()(vec), List(r, i->ivExp(i))))
));

NewRulesFor(TIVArith, rec(
    TIVArith_Base := rec(
        applicable := True,
        forTransposition := false,
        children := nt -> [[nt.params[1]]],
        apply := (t, C, Nonterms) -> IVArith(C[1])
    )
));

HCOLSumsGen.IVArith := (self, o, opts) >> IVArith(self(o.child(1), opts));

Class(TypeUpdate, HierarchicalVisitor, rec(
	__call__ := meth(arg)
        local res;
        res := ApplyFunc(arg[1].visit, arg{[2..Length(arg)]});
        return res;
    end,
	func := (self, o, opts) >> let(func(o.ret, o.id, o.params, self(o.cmd, opts))),
	decl := (self, o, opts) >>  decl(o.vars, self(o.cmd, opts)),
	chain := (self, o, opts) >> chain(List(o.cmds, i-> self(i, opts))),
    ivenv := (self, o, opts) >> ivenv(List(o.cmds, i-> self(i, opts))),
    loop := (self, o, opts) >>  loop(o.var, o.range, self(o.cmd, opts)),

	assign := (self, o, opts) >> assign(o.loc, self(o.exp, opts)),
	add := (self, o, opts) >> add(self(o.args[1], opts), self(o.args[2], opts)),
	sub := (self, o, opts) >> sub(self(o.args[1], opts), self(o.args[2], opts)),
	mul := (self, o, opts) >> mul(self(o.args[1], opts), self(o.args[2], opts)),
	div := (self, o, opts) >> div(self(o.args[1], opts), self(o.args[2], opts)),
	cond := (self, o, opts) >> cond(o.args[1], self(o.args[2], opts), self(o.args[3], opts)),
	abs := (self, o, opts) >> abs(self(o.args[1], opts)),
	logic_and := (self, o, opts) >> logic_and(self(o.args[1], opts), self(o.args[2], opts)),
	logic_or := (self, o, opts) >> logic_or(self(o.args[1], opts), self(o.args[2], opts)),
	gt  := (self, o, opts) >> gt(self(o.args[1], opts), self(o.args[2], opts)),
	lt  := (self, o, opts) >> lt(self(o.args[1], opts), self(o.args[2], opts)),
	geq  := (self, o, opts) >> geq(self(o.args[1], opts), self(o.args[2], opts)),
	leq  := (self, o, opts) >> leq(self(o.args[1], opts), self(o.args[2], opts)),
	min := (self, o, opts) >> min(self(o.args[1], opts), self(o.args[2], opts)),
	max := (self, o, opts) >> max(self(o.args[1], opts), self(o.args[2], opts)),

	nth := (self, o, opts) >> o,
	Value := (self, o, opts) >> o,
	var := (self, o, opts) >> o,

	skip := (self, o, opts) >> o,
	tcast := (self, o, opts) >> o,

	creturn := (self, o, opts) >> o,
	creturnCond := (self, o, opts) >> o,

	eq := (self, o, opts) >> eq(self(o.args[1], opts), self(o.args[2], opts)),
));

HCOLCodegen.IVArith := (self, o, y, x, opts) >> let(
	#Generate non interval code
    cc := opts.codegen.OLCompose(o.child(1), y, x, opts),

	#Collect array of bools and make them into arrays of ints
	ba_vars := Filtered(FoldL(Collect(cc, var), (a,b)->When(not b.id in List(a, i->i.id), Concat([b],a), a), []),
    	    e->ObjId(e.t) = TArray and e.t.base_t()=TBool),
    ba_nvars := List(ba_vars, c-> [c, var.fresh_t("U", TArray(TInt, c.t.size))] ),
    ba_svars := FoldL(ba_nvars, (b, a)->CopyFields(rec((a[1].id):= a[2]), b), rec()),
    ba_cc1 := decl(List(ba_nvars, c->c[2]), SubstVars(cc, ba_svars)),

	#Cast bools variables into ints but make sure that they are not part of the input parameters
	b_vars :=Unique(Collect(ba_cc1, @(1,var).cond(e->e.t = TBool))),
	b_local_vars := Difference(b_vars, opts.params),
	b_global_vars := Difference(b_vars, b_local_vars),
	b_no_global := SubstBottomUp(ba_cc1, @(1,var).cond(e->e.t=TBool), f->let(Print(@(1).val.id, "\n"),
							tcast(TInt, @(1).val))),

	#Cast pointers to bools into pointers to ints
	b_cc2 := b_no_global,

	#Convert boolean constants into ints
	b_cc3 := SubstBottomUp(b_cc2, @(1,Value).cond(e->e.t=TBool),
        f->let(#Print(@(1).val, "\n"),
			Cond(@(1).val=V_true,  V(1),
			     @(1).val=V_false, V(0),
						   @(1).val))),

	#Collect array variables make them into arrays of intervals
    vars := Filtered(FoldL(Collect(b_cc3, var), (a,b)->When(not b.id in List(a, i->i.id), Concat([b],a), a), []),
    	    e->ObjId(e.t) = TArray and e.t.base_t() in [TReal, TDouble, T_Real(32), T_Real(64)]),
    nvars := List(vars, c-> [c, var.fresh_t("U", TArray(TIVReal(c.t.base_t()), c.t.size))] ),
    svars := FoldL(nvars, (b, a)->CopyFields(rec((a[1].id):= a[2]), b), rec()),
    cc1 := decl(List(nvars, c->c[2]), SubstVars( cc, svars)),

	#Cast doubles and reals variables into intervals
    cc11 := SubstBottomUp(cc1, @(1,var).cond(e->e.t in [TReal, TDouble, T_Real(32), T_Real(64)]), f->tcast(opts.IVRealType, @(1).val)),

	#Cast pointers to doubles and reals into intervals
    cc2 := SubstBottomUp(cc11, [@(1,nth), @(2, var, e->ObjId(e.t) = TPtr and e.t.base_t() in [TReal, TDouble, T_Real(32), T_Real(64)]), @(3)], f->tcast(opts.IVRealType, @(1).val)),

	#Cast constants into intervals
    cc3 := SubstBottomUp(cc2, @(1,Value).cond(e->e.t in [TReal, TDouble, T_Real(32), T_Real(64)]),
        f->let(tcast(opts.IVValueType, @(1).val))),

	result := ivenv(cc3),

	result2 := TypeUpdate(result, opts),

	result2
);

HCOLUnparser.RealEPS := (self,o,i,is) >> Print(Cond(o.args[1] = T_Real(64), "DBL_MIN", o.args[1] = T_Real(32), "FLT_MIN", "UNKNOWN_MIN"));
HCOLUnparser.RealMPS := (self,o,i,is) >> Print(Cond(o.args[1] = T_Real(64), "DBL_MAX", o.args[1] = T_Real(32), "FLT_MAX", "UNKNOWN_MAX"));
HCOLUnparser.ivenv := (self,o,i,is) >> self(chain(o.cmds), i, is);

HCOLUnparser.TIVReal := (self, t, vars, i, is) >>
		Print("interval_t ", self.infix(vars, ", ", i + is));
HCOLUnparser.TVect := (self, t, vars, i, is) >> Print("abc_interval_t ", self.infix(vars, ", ", i + is));

HCOLSSEUnparser.ivenv := (self,o,i,is) >> Print(Blanks(i), "{\n",
    Blanks(i+is), "unsigned _xm = _mm_getcsr();\n",
    Blanks(i+is), "_mm_setcsr(_xm & 0xffff0000 | 0x0000dfc0);\n",
    self(chain(o.cmds),i+is,is),
    Blanks(i+is),
    Cond(IsBound(self.opts.compiler) and self.opts.compiler = "IntelC", "__asm nop;\n",
        IsBound(self.opts.compiler) and self.opts.compiler = "GnuC", "asm volatile(\"\":::\"memory\");\n",
        "// BASIC BLOCK BARRIER\n"
    ),
    Blanks(i+is), "if (_mm_getcsr() & 0x0d) {\n",
    Blanks(i+2*is), "_mm_setcsr(_xm);\n",
	Blanks(i+2*is),
    Cond(self.opts.useCReturn, "return -1;\n", "Y[0] = -1;\n"),
    Blanks(i+is), "}\n",
    Blanks(i+is), "_mm_setcsr(_xm);\n",
    Blanks(i), "}\n");

HCOLSSEUnparser.RealEPS := HCOLUnparser.RealEPS;
HCOLSSEUnparser.RealMPS := HCOLUnparser.RealMPS;

Declare(ToIVArithBasic);

Class(ToIVArithBasic_Base, HierarchicalVisitor, rec(
    __call__ := meth(arg)
        local res;
        res := ApplyFunc(arg[1].visit, arg{[2..Length(arg)]});
        trace_log.addConversion(ObjId(arg[2]), arg[2], res, var);
       return res;
    end,
    func := (self, o, opts) >> let(Print("CONVERSION\n\n"),
		func(o.ret, o.id, FoldL(o.params, (a, b)->Concat(a, When(IsBound(opts.state) and b=opts.state, [opts.ivstate], [b])), []), self(o.cmd, opts))),
    creturn := (self, o, opts) >> o,
	creturnCond := (self, o, opts) >> o,
    decl := (self, o, opts) >> let(
        ovars := Filtered(o.vars, e->ObjId(e.t) = TIVReal),
        svars := FoldR(ovars, (b,a)->CopyFields(rec((a.id):=var.fresh_t("u", TVect(T_Real(64), 2))), b), rec()),
        vars := Filtered(o.vars, e->ObjId(e.t) <> TIVReal)::List(Filtered(RecFields(svars), i->not i in SystemRecFields), i->svars.(i)),
        cmd := SubstVars(o.cmd, svars),
        decl(vars, self(cmd, opts))),
    chain := (self, o, opts) >> chain(List(o.cmds, i-> self(i, opts))),
    ivenv := (self, o, opts) >> ivenv(List(o.cmds, i-> ToIVArithBasic(i, opts))),
    loop := (self, o, opts) >> loop(o.var, o.range, self(o.cmd, opts)),
    nth := (self, o, opts) >> o,

    Value := (self, o, opts) >> let(
        When(ObjId(o.t) = TIVReal,
          TVect(T_Real(64), 2).value([-o.v, o.v]),
          o)),

	skip := (self, o, opts) >> o,

	tcast := (self, o, opts) >> o,
	logic_and := (self, o, opts) >> o,
	logic_or := (self, o, opts) >> o,
	assign := (self, o, opts) >> o,
));

Class(ToIVArithBasic, HierarchicalVisitor, rec(
    __call__ := meth(arg)
        local res;
        res := ApplyFunc(arg[1].visit, arg{[2..Length(arg)]});
        trace_log.addConversion(ObjId(arg[2]), arg[2], res, var);
       return res;
    end,
    func := (self, o, opts) >> let(Print("CONVERSION\n\n"),
		func(o.ret, o.id, FoldL(o.params, (a, b)->Concat(a, When(IsBound(opts.state) and b=opts.state, [opts.ivstate], [b])), []), self(o.cmd, opts))),
    creturn := (self, o, opts) >> o,
	creturnCond := (self, o, opts) >> o,
    decl := (self, o, opts) >> let(
        ovars := Filtered(o.vars, e->ObjId(e.t) = TIVReal),
        svars := FoldR(ovars, (b,a)->CopyFields(rec((a.id):=var.fresh_t("u", TVect(T_Real(64), 2))), b), rec()),
        vars := Filtered(o.vars, e->ObjId(e.t) <> TIVReal)::List(Filtered(RecFields(svars), i->not i in SystemRecFields), i->svars.(i)),
        cmd := SubstVars(o.cmd, svars),
        decl(vars, self(cmd, opts))),
    chain := (self, o, opts) >> chain(List(o.cmds, i-> self(i, opts))),
    ivenv := (self, o, opts) >> ivenv(List(o.cmds, i-> self(i, opts))),
    loop := (self, o, opts) >> loop(o.var, o.range, self(o.cmd, opts)),
    nth := (self, o, opts) >> o,

    Value := (self, o, opts) >> let(
        When(ObjId(o.t) = TIVReal,
          TVect(T_Real(64), 2).value([-o.v, o.v]),
          o)),

	skip := (self, o, opts) >> o,

	tcast := (self, o, opts) >> o,

	logic_and := (self, o, opts) >> cond(eq(o.args[1],o.args[2]), o.args[1], mul(o.args[1], o.args[2])),
	logic_or := (self, o, opts) >> cond(logic_or(eq(o.args[1], 1), eq(o.args[2],1)), 1, min(o.args[1], o.args[2])),

    assign := (self, o, opts) >> let(
	  Cond(
	    ObjId(o.exp) = var,
			let(
				chain(
					assign(o.loc[1], tcast(o.loc[1].t,o.exp)),
					assign(o.loc[2], tcast(o.loc[2].t,o.exp))
				)
			),
	    ObjId(o.exp) = tcast and (o.exp.args[1] = TIVReal or o.exp.args[1] = TVect(T_Real(64), 2)),
				Cond(
					o.exp.args[2].t = T_Real(32),
						let(
							x0 := var.fresh_t("x", T_Real(64)),
							x1 := var.fresh_t("x", T_Real(64)),
							decl([x0, x1], chain(
								self(assign([x0,x1], o.exp.args[2]), opts),
								assign(o.loc[1], x0+RealEPS(T_Real(32))),
								assign(o.loc[2], x1-RealEPS(T_Real(32)))
							))
						),
					IsDouble(o.exp.args[2]),
						let(
							x0 := var.fresh_t("x", T_Real(64)),
							x1 := var.fresh_t("x", T_Real(64)),
							decl([x0, x1], chain(
								self(assign([x0,x1], o.exp.args[2]), opts),
								assign(o.loc[1], x0),
								assign(o.loc[2], -x1)
							))
						),
					Error("Don't know how to convert <o.args[2]> to an interval.")),
		ObjId(o.exp) = lt and ForAll(o.exp.args, i->i.t = TIVReal or i.t = TVect(T_Real(64), 2)),
            let(
				x10 := var.fresh_t("x", T_Real(64)),
				x11 := var.fresh_t("x", T_Real(64)),
				x20 := var.fresh_t("x", T_Real(64)),
				x21 := var.fresh_t("x", T_Real(64)),
                decl([x10, x11, x20, x21], chain(
					self(assign([x10, x11], o.exp.args[1]), opts),
					self(assign([x20, x21], o.exp.args[2]), opts),
					assign(o.loc, tcast(TInt, cond(lt(x10, x21), V(1), cond(lt(x20,x11), V(0),V(-1)))))
                ))
            ),
		ObjId(o.exp) = gt and ForAll(o.exp.args, i->i.t = TIVReal or i.t = TVect(T_Real(64), 2)),
            let(x10 := var.fresh_t("x", T_Real(64)),
				x11 := var.fresh_t("x", T_Real(64)),
				x20 := var.fresh_t("x", T_Real(64)),
				x21 := var.fresh_t("x", T_Real(64)),
                decl([x10, x11, x20, x21], chain(
					self(assign([x10, x11], o.exp.args[1]), opts),
					self(assign([x20, x21], o.exp.args[2]), opts),
					assign(o.loc, tcast(TInt, cond(gt(x11, x20), V(1), cond(gt(x21, x10), V(0), V(-1)))))
                ))
            ),

		ObjId(o.exp) = add  and ForAll(o.exp.args, i->i.t = TIVReal or i.t = TVect(T_Real(64), 2)),
			let(x10 := var.fresh_t("x", T_Real(64)),
				x11 := var.fresh_t("x", T_Real(64)),
				x20 := var.fresh_t("x", T_Real(64)),
				x21 := var.fresh_t("x", T_Real(64)),
				decl([x10, x11, x20, x21], chain(
					self(assign([x10, x11], o.exp.args[1]), opts),
					self(assign([x20, x21], o.exp.args[2]), opts),
					assign(o.loc[1], add(x11, x21)),
					assign(o.loc[2], add(x10, x20))
				))
		  ),
		ObjId(o.exp) = sub  and ForAll(o.exp.args, i->i.t = TIVReal or i.t = TVect(T_Real(64), 2)),
			let(x10 := var.fresh_t("x", T_Real(64)),
				x11 := var.fresh_t("x", T_Real(64)),
				x20 := var.fresh_t("x", T_Real(64)),
				x21 := var.fresh_t("x", T_Real(64)),
				decl([x10, x11, x20, x21], chain(
					self(assign([x10, x11], o.exp.args[1]), opts),
					self(assign([x20, x21], o.exp.args[2]), opts),
					assign(o.loc[1], sub(x11, x20)),
					assign(o.loc[2], sub(x10, x21))
				))
			),
		ObjId(o.exp) = abs  and ForAll(o.exp.args, i->i.t = TIVReal or i.t = TVect(T_Real(64), 2)),
            let(
				x10 := var.fresh_t("x", T_Real(64)),
				x11 := var.fresh_t("x", T_Real(64)),
				x2 := var.fresh_t("x", T_Real(64)),
				x3 := var.fresh_t("x", T_Real(64)),
                decl([x10, x11, x2, x3], chain(
					self(assign([x10, x11], o.exp.args[1]), opts),
					assign(x2, abs(x10)),
					assign(x3, abs(x11)),
					assign(o.loc[1], max(x2, x3)),
					assign(o.loc[2], min(x2, x3))
                ))
            ),
		assign(o.loc, self(o.exp, opts))
    )),
));

Class(ToIVArithSSE, HierarchicalVisitor, rec(
    __call__ := meth(arg)
        local res;
        res := ApplyFunc(arg[1].visit, arg{[2..Length(arg)]});
        trace_log.addConversion(ObjId(arg[2]), arg[2], res, var);
       return res;
    end,
    func := (self, o, opts) >> let(#Print("CONVERSION\n\n"),
		func(o.ret, o.id, FoldL(o.params, (a, b)->Concat(a, When(IsBound(opts.state) and b=opts.state, [opts.ivstate], [b])), []), self(o.cmd, opts))),
    creturn := (self, o, opts) >> o,
	creturnCond := (self, o, opts) >> o,
    decl := (self, o, opts) >> let(
        ovars := Filtered(o.vars, e->ObjId(e.t) = TIVReal),
        svars := FoldR(ovars, (b,a)->CopyFields(rec((a.id):=var.fresh_t("u", TVect(T_Real(64), 2))), b), rec()),
        vars := Filtered(o.vars, e->ObjId(e.t) <> TIVReal)::List(Filtered(RecFields(svars), i->not i in SystemRecFields), i->svars.(i)),
        cmd := SubstVars(o.cmd, svars),
        decl(vars, self(cmd, opts))),
    chain := (self, o, opts) >> chain(List(o.cmds, i-> self(i, opts))),
    ivenv := (self, o, opts) >> ivenv(List(o.cmds, i-> self(i, opts))),
    loop := (self, o, opts) >> loop(o.var, o.range, self(o.cmd, opts)),
    nth := (self, o, opts) >> o,
    Value := (self, o, opts) >> let(
        When(ObjId(o.t) = TIVReal,
          TVect(T_Real(64), 2).value([-o.v, o.v]),
          o)),
    var := (self, o, opts) >> o,
	skip := (self, o, opts) >> o,
    tcast := (self, o, opts) >> Cond(
        IsDouble(o.args[2]), vpack(-o.args[2], o.args[2]),
	IsValue(o.args[2]), TVect(T_Real(64), 2).value([-o.args[2].v, o.args[2].v]),
        o.args[2].t = T_Real(32), addsub_2x64f(vcvt_64f32f(vdup(RealEPS(T_Real(32)), 4)), vcvt_64f32f(vdup(o.args[2], 4))),
        o.args[2].t = T_Real(64), addsub_2x64f(vdup(RealEPS(T_Real(64))+RealEPS(T_Real(64)), 2), vdup(o.args[2], 2)),
	o.args[2].t = TIVReal(T_Real(64)),  TVect(T_Real(64),2).value(o.args[2]),
        Error("Don't know how to convert <o.args[2]> to an interval.")),

	logic_and := (self, o, opts) >> cond(eq(o.args[1],o.args[2]), o.args[1], mul(o.args[1], o.args[2])),
	logic_or := (self, o, opts) >> cond(logic_or(eq(o.args[1], 1), eq(o.args[2],1)), 1, min(o.args[1], o.args[2])),

    testc_4x32i:= (self, o, opts) >> o,

    assign := (self, o, opts) >> Cond(
	    ObjId(o.exp) = add and (o.loc.t = TIVReal or o.loc.t = TVect(T_Real(64), 2)),
		  let(
			x1 := var.fresh_t("x", TVect(T_Real(64), 2)), x2 := var.fresh_t("x", TVect(T_Real(64), 2)),
			decl([x1, x2], chain(
				comment("addition"),
				self(assign(x1, o.exp.args[1]), opts),
				self(assign(x2, o.exp.args[2]), opts),
				assign(o.loc, add(x1, x2)) ) )
		  ),
		ObjId(o.exp) = sub and (o.loc.t = TIVReal or o.loc.t = TVect(T_Real(64), 2)),
		  let(
			x1 := var.fresh_t("x", TVect(T_Real(64), 2)), x2 := var.fresh_t("x", TVect(T_Real(64), 2)),
			decl([x1, x2], chain(
				comment("sub"),
				self(assign(x1, o.exp.args[1]), opts),
				self(assign(x2, o.exp.args[2]), opts),
				assign(o.loc, add(x1, vushuffle_2x64f(x2, vparam([2,1])))) ) )
		  ),
        ObjId(o.exp) = mul and (o.loc.t = TIVReal or o.loc.t = TVect(T_Real(64), 2)),
            let(u := var.fresh_t("x", TVect(T_Real(64), 2)), a := var.fresh_t("x", TVect(T_Real(64), 2)),
                b := var.fresh_t("x", TVect(T_Real(64), 2)), c := var.fresh_t("x", TVect(T_Real(64), 2)),
                x1 := var.fresh_t("x", TVect(T_Real(64), 2)), x2 := var.fresh_t("x", TVect(T_Real(64), 2)),
                decl([u, a, b, c, x1, x2], chain(
					comment("mul"),
					self(assign(x1, o.exp.args[1]), opts),
					self(assign(x2, o.exp.args[2]), opts),
                    assign(u, addsub_2x64f(TVect(T_Real(64), 2).zero(), x1)),
                    assign(a, mul(u, x2)),
                    assign(b, mul(vushuffle_2x64f(u, vparam([2,1])), x2)),
                    assign(c, neg(min(a, b))),
                    assign(o.loc, add(max(max(a, b), vushuffle_2x64f(c, vparam([2,1]))), vdup(RealEPS(T_Real(64)), 2)))
                ))
            ),
        ObjId(o.exp) = abs and (o.loc.t = TIVReal or o.loc.t = TVect(T_Real(64), 2)),
            let(u := var.fresh_t("x", TVect(T_Real(64), 2)), x1 := var.fresh_t("x", TVect(T_Real(64), 2)),
                decl([u, x1], chain(
					comment("abs"),
                    self(assign(x1, o.exp.args[1]), opts),
                    assign(u, vushuffle_2x64f(x1, vparam([2,1]))),
                    assign(o.loc, vshuffle_2x64f(min(x1, u), max(x1, u), vparam([1,2])))
                ))
            ),
		ObjId(o.exp) = pow and (o.loc.t = TIVReal or o.loc.t = TVect(T_Real(64), 2)) and IsValue(o.exp.args[2]) and o.exp.args[2].v = 2,
			let( x1 := var.fresh_t("x", TVect(T_Real(64), 2)),
				decl([x1], chain(
					self(assign(x1, mul(o.exp.args[1],o.exp.args[1])), opts),
					assign(o.loc, x1)
				))
			),
		ObjId(o.exp) = pow and (o.loc.t = TIVReal or o.loc.t = TVect(T_Real(64), 2)) and IsValue(o.exp.args[2]) and IsInt(o.exp.args[2].v) and o.exp.args[2] > 2,
			let( x1 := var.fresh_t("x", TVect(T_Real(64), 2)),
				decl([x1], chain(
					self(assign(x1, mul(o.exp.args[1], pow(o.exp.args[1],o.exp.args[2].v-1))), opts),
					assign(o.loc, x1)
				))
			),

		ObjId(o.exp) = div and (o.loc.t = TIVReal or o.loc.t = TVect(T_Real(64), 2)) and not IsValue(o.exp.args[1]),
			let(
			Print("DIVISION - PART 1\n"),
			  x1 := var.fresh_t("x", TVect(T_Real(64), 2)),
			  decl([x1], chain(
				self(assign(x1, mul(o.exp.args[1], 1.0/o.exp.args[2])), opts),
				assign(o.loc, x1)
			  ))
			),

		ObjId(o.exp) = div and <#o.loc.t = TIVReal or o.loc.t = TVect(T_Real(64), 2)) and #>IsValue(o.exp.args[1]),
			let(
				x1 := var.fresh_t("x", TVect(T_Real(64), 2)),
				decl([x1], chain(
					assign(x1, vdiv_2x64f(TVect(T_Real(64), 2).value([-o.exp.args[1].v,o.exp.args[1].v]), o.exp.args[2])),
					assign(o.loc, vshuffle_2x64f(x1, x1, vparam([2,1])))
				))
			),
		ObjId(o.exp) = neg and (o.loc.t = TIVReal or o.loc.t = TVect(T_Real(64), 2)),
		    let(x1 := var.fresh_t("x", TVect(T_Real(64), 2)),x2 := var.fresh_t("x", TVect(T_Real(64), 2)),
				decl([x1], chain(
					self(assign(x1, TVect(T_Real(64), 2).value([0,0]), opts)),
					self(assign(x2, sub(x1, o.exp.args[1])), opts),
					assign(o.loc, x2)
				))
			),
        ObjId(o.exp) = max and (o.loc.t = TIVReal or o.loc.t = TVect(T_Real(64), 2)),
            let(x1 := var.fresh_t("x", TVect(T_Real(64), 2)), x2 := var.fresh_t("x", TVect(T_Real(64), 2)),
                decl([x1, x2], chain(
                    self(assign(x1, o.exp.args[1]), opts),
                    self(assign(x2, o.exp.args[2]), opts),
                    assign(o.loc, vshuffle_2x64f(min(x1, x2), max(x1, x2), vparam([1, 2])))
                ))
            ),
         ObjId(o.exp) = min and (o.loc.t = TIVReal or o.loc.t = TVect(T_Real(64), 2)),
            let(x1 := var.fresh_t("x", TVect(T_Real(64), 2)), x2 := var.fresh_t("x", TVect(T_Real(64), 2)),
                decl([x1, x2], chain(
                    self(assign(x1, o.exp.args[1]), opts),
                    self(assign(x2, o.exp.args[2]), opts),
                    assign(o.loc, vshuffle_2x64f(max(x1, x2), min(x1, x2), vparam([1, 2])))
                ))
            ),
        ObjId(o.exp) = geq <#and (o.loc.t = TBool or o.loc.t = TInt)#> and ForAll(o.exp.args, i->i.t = TIVReal or i.t = TVect(T_Real(64), 2)),
            let(x1 := var.fresh_t("x", TVect(T_Real(64), 2)), x2 := var.fresh_t("x", TVect(T_Real(64), 2)), u := var.fresh_t("x", TVect(T_Real(64), 2)),
                decl([x1, x2, u], chain(
					comment("geq"),
                    assign(x1, addsub_2x64f(TVect(T_Real(64), 2).zero(), self(o.exp.args[1], opts))),
                    assign(x2, addsub_2x64f(TVect(T_Real(64), 2).zero(), self(o.exp.args[2], opts))),
                    assign(u, cmpge_2x64f(x1, vushuffle_2x64f(x2, vparam([2, 1])))),
                    assign(o.loc, sub(
                        testc_4x32i(tcast(TVect(T_Int(32), 4), u), vhex(Replicate(4, "0xffffffff"))),
                        testnzc_4x32i(tcast(TVect(T_Int(32), 4), u), vhex(Replicate(4, "0xffffffff")))))
                ))
            ),
        ObjId(o.exp) = leq <#and o.loc.t = TBool#> and ForAll(o.exp.args, i->i.t = TIVReal or i.t = TVect(T_Real(64), 2)),
            let(x1 := var.fresh_t("x", TVect(T_Real(64), 2)), x2 := var.fresh_t("x", TVect(T_Real(64), 2)), u := var.fresh_t("x", TVect(T_Real(64), 2)),
                decl([x1, x2, u], chain(
                    assign(x1, addsub_2x64f(TVect(T_Real(64), 2).zero(), self(o.exp.args[1], opts))),
                    assign(x2, addsub_2x64f(TVect(T_Real(64), 2).zero(), self(o.exp.args[2], opts))),
                    assign(u, cmple_2x64f(x1, vushuffle_2x64f(x2, vparam([2, 1])))),
                    assign(o.loc, sub(
                        testc_4x32i(tcast(TVect(T_Int(32), 4), u), vhex(Replicate(4, "0xffffffff"))),
                        testnzc_4x32i(tcast(TVect(T_Int(32), 4), u), vhex(Replicate(4, "0xffffffff")))))
                ))
            ),
		ObjId(o.exp) = gt <#and (o.loc.t = TBool or o.loc.t = TInt)#> and ForAll(o.exp.args, i->i.t = TIVReal or i.t = TVect(T_Real(64), 2)),
            let(x1 := var.fresh_t("x", TVect(T_Real(64), 2)), x2 := var.fresh_t("x", TVect(T_Real(64), 2)), u := var.fresh_t("x", TVect(T_Real(64), 2)),
                decl([x1, x2, u], chain(
                    assign(x1, addsub_2x64f(TVect(T_Real(64), 2).zero(), self(o.exp.args[1], opts))),
                    assign(x2, addsub_2x64f(TVect(T_Real(64), 2).zero(), self(o.exp.args[2], opts))),
                    assign(u, cmpgt_2x64f(x1, vushuffle_2x64f(x2, vparam([2, 1])))),
                    assign(o.loc, sub(
                        testc_4x32i(tcast(TVect(T_Int(32), 4), u), vhex(Replicate(4, "0xffffffff"))),
                        testnzc_4x32i(tcast(TVect(T_Int(32), 4), u), vhex(Replicate(4, "0xffffffff")))))
                ))
            ),
        ObjId(o.exp) = lt <#and o.loc.t = TBool#> and ForAll(o.exp.args, i->i.t = TIVReal or i.t = TVect(T_Real(64), 2)),
            let(x1 := var.fresh_t("x", TVect(T_Real(64), 2)), x2 := var.fresh_t("x", TVect(T_Real(64), 2)), u := var.fresh_t("x", TVect(T_Real(64), 2)),
                decl([x1, x2, u], chain(
                    assign(x1, addsub_2x64f(TVect(T_Real(64), 2).zero(), self(o.exp.args[1], opts))),
                    assign(x2, addsub_2x64f(TVect(T_Real(64), 2).zero(), self(o.exp.args[2], opts))),
                    assign(u, cmplt_2x64f(x1, vushuffle_2x64f(x2, vparam([2, 1])))),
                    assign(o.loc, sub(
                        testc_4x32i(tcast(TVect(T_Int(32), 4), u), vhex(Replicate(4, "0xffffffff"))),
                        testnzc_4x32i(tcast(TVect(T_Int(32), 4), u), vhex(Replicate(4, "0xffffffff")))))
                ))
            ),
		ObjId(o.exp) = eq,
			let(Error("Eq not implemented for interval arithmetic.")),
		ObjId(o.exp) = neq,
			let(Error("Neq not implemented for interval arithmetic.")),

        assign(o.loc, self(o.exp, opts))
    ),
    mul := (self, o, opts) >> ApplyFunc(ObjId(o), List(o.args, i-> self(i, opts))),
    abs := (self, o, opts) >> ApplyFunc(ObjId(o), List(o.args, i-> self(i, opts))),
    max := (self, o, opts) >> ApplyFunc(ObjId(o), List(o.args, i-> self(i, opts))),
	add := (self, o, opts) >> ApplyFunc(ObjId(o), List(o.args, i-> self(i, opts))),
	sub := (self, o, opts) >> ApplyFunc(ObjId(o), List(o.args, i-> self(i, opts))),
    geq := (self, o, opts) >> ApplyFunc(ObjId(o), List(o.args, i-> self(i, opts))),
    leq := (self, o, opts) >> ApplyFunc(ObjId(o), List(o.args, i-> self(i, opts))),
	gt  := (self, o, opts) >> ApplyFunc(ObjId(o), List(o.args, i-> self(i, opts))),
    lt  := (self, o, opts) >> ApplyFunc(ObjId(o), List(o.args, i-> self(i, opts))),
	pow := (self, o, opts) >> ApplyFunc(ObjId(o), List(o.args, i-> self(i, opts))),
	div := (self, o, opts) >> ApplyFunc(ObjId(o), List(o.args, i-> self(i, opts))),
));

HCOLSSEUnparser.TIVReal := (self, t, vars, i, is) >> Print("__m128d ", self.infix(vars, ", ", i + is));

CoSynthesizeStrategies.IVArithSSE := [ (t, opts) -> RandomRuleTree(t, opts),
        (rt, opts) -> SPLRuleTree(rt),
        (s, opts) -> SumsSPL(s, opts),
        (s, opts) -> Rewrite(s, [RulesSumsHCOLv2a, RulesTerminateReductionHCOL, RulesSumsHCOLv2b], opts),
        (s, opts) -> HCOLProof_Codegen(s, opts),
        (c, opts) -> Rewrite(c, RulesCodeUnrollHACMIVArth, opts),
        (c, opts) -> HCOLProof_CodeConversion(c, ToIVArithSSE, opts),
        (c, opts) -> Rewrite(c, RulesCodeUnrollHACMIVArthType, opts) ];

CoSynthesizeStrategies.IVArithBasic := [ (t, opts) -> RandomRuleTree(t, opts),
        (rt, opts) -> SPLRuleTree(rt),
        (s, opts) -> SumsSPL(s, opts),
        (s, opts) -> Rewrite(s, [RulesSumsHCOLv2a, RulesTerminateReductionHCOL, RulesSumsHCOLv2b], opts),
        (s, opts) -> HCOLProof_Codegen(s, opts),
        (c, opts) -> Rewrite(c, RulesCodeUnrollHACMIVArth, opts),
        (c, opts) -> HCOLProof_CodeConversion(c, ToIVArithBasic_Base, opts),
        (c, opts) -> Rewrite(c, RulesCodeUnrollHACMIVArth, opts) ];

HCOLSSEUnparser.addsub_2x64f := (self, o, i, is) >> Checked(Length(o.args) = 2,
        CondPat(o,
           [addsub_2x64f, @TReal, @TVect], self(addsub_2x64f(vdup(o.args[1],o.t.size), o.args[2]), i, is),
           [addsub_2x64f, @TVect, @TReal], self(addsub_2x64f(o.args[1], vdup(o.args[2], o.t.size)), i, is),
           [addsub_2x64f, @TInt,  @TVect], self(addsub_2x64f(vdup(_toReal(o.args[1]), o.t.size), o.args[2]), i, is),
           [addsub_2x64f, @TVect, @TInt],  self(addsub_2x64f(o.args[1], vdup(_toReal(o.args[2]), o.t.size)), i, is),
           [addsub_2x64f, @TVect, @TVect], self.printf("_mm_addsub_pd($1, $2)", o.args),
           [addsub_2x64f, @TVect, @(1).cond(e->e.t=TIVReal)], self.printf("_mm_addsub_pd($1, $2)", o.args),
           [addsub_2x64f, @(1).cond(e->e.t=TIVReal), @TVect], self.printf("_mm_addsub_pd($1, $2)", o.args),
           Error("Don't know how to unparse <o>. Unrecognized type combination")
    ));

HCOLSSEUnparser.max := (self, o, i, is) >> let(n := Length(o.args), When(
	IsVecT(o.t) and n >2, self.printf("_mm_max_$1($2, $3)", [self.ctype_suffix(o.t, _isa(self)), o.args[1],
		ApplyFunc(max, Drop(o.args, 1))]),
	CondPat(o,
		[max, @TVect, @TVect], self.prefix("_mm_max_" :: self.ctype_suffix(o.t, _isa(self)),o.args),
	    [max, @TVect, @(1).cond(e->e.t=TIVReal)], self.printf("_mm_max_pd($1, $2)", o.args),
        [max, @(1).cond(e->e.t=TIVReal), @TVect], self.printf("_mm_max_pd($1, $2)", o.args),
		Inherited(o, i, is))
));

HCOLSSEUnparser.min := (self, o, i, is) >> let(
	CondPat(o,
		[min, @TVect, @TVect], self.prefix("_mm_min_" :: self.ctype_suffix(o.t, _isa(self)),o.args),
	    [min, @TVect, @(1).cond(e->e.t=TIVReal)], self.printf("_mm_min_pd($1, $2)", o.args),
        [min, @(1).cond(e->e.t=TIVReal), @TVect], self.printf("_mm_min_pd($1, $2)", o.args),
		Inherited(o, i, is))
);

geq.computeType := self >> TBool;
leq.computeType := self >> TBool;
eq.computeType := self >> TBool;
pow.computeType := (self) >> self.args[1].t;
neq.computeType := (self) >> TBool;

min.computeType := self >> let(
		types := List(self.args, e->e.t),
		When(Length(Collect(types, TIVReal))>0,
			TIVReal(UnifyTypes(List(Collect(types, TIVReal), i->i.t))),   <# TVect(T_Real(64),2),  #>
			UnifyTypes(List(self.args, x->x.t)))
	);
max.computeType := self >> let(
		types := List(self.args, e->e.t),
		When(Length(Collect(types, TIVReal))>0,
			TIVReal(UnifyTypes(List(Collect(types, TIVReal), i->i.t))),   <# TVect(T_Real(64),2),  #>
			UnifyTypes(List(self.args, x->x.t)))
	);

VecExp.computeType := self >> let(
        t       := self.args[1].t,
        deref_t := When(IsPtrT(t), t.t, t),
        el_t    := Cond(IsVecT(deref_t), deref_t.t, deref_t),
	test    := Cond(el_t = TIVReal, el_t.t, el_t),
        TVect(test, self.v));

add.computeType := meth(self)
	    local len, t, ptr_args, other_args, sum;
	    len := Length(self.args);
	    if   len=0  then return TInt;
	    elif len=1  then return self.args[1].t;
	    else
            [ptr_args, other_args] := SplitBy(self.args, x->IsPtrT(x.t) or IsArrayT(x.t));
            if Length(ptr_args)=0 then
	       if (self.args[1].t = TIVReal or self.args[2].t = TIVReal) then
	         return TIVReal(T_Real(64));
	       else
	         return UnifyTypesL(self.args);
               fi;
            elif Length(ptr_args)=1 then
                sum := Sum(other_args);
                if other_args<>[] and not IsIntT(sum.t) then Error("Can't add non-integer to a pointer"); fi;
		   return self._ptrPlusOfs(ptr_args[1].t, sum);
            elif Length(other_args)=0 then
                return self._addPtrT(ptr_args);
            else
                return Error("Addition of more than one pointer and integers is not defined");
            fi;
	    fi;
    end;

sub.computeType := meth(self)
	    local len, t, ptr_args, other_args, sum;
	    len := Length(self.args);
	    if   len=0  then return TInt;
	    elif len=1  then return self.args[1].t;
	    else
            [ptr_args, other_args] := SplitBy(self.args, x->IsPtrT(x.t) or IsArrayT(x.t));
            if Length(ptr_args)=0 then
	       if (self.args[1].t = TIVReal or self.args[2].t = TIVReal) then
	         return TIVReal(T_Real(64));
	       else
	         return UnifyTypesL(self.args);
               fi;
            elif Length(ptr_args)=1 then
                sum := Sum(other_args);
                if other_args<>[] and not IsIntT(sum.t) then Error("Can't add non-integer to a pointer"); fi;
		   return self._ptrPlusOfs(ptr_args[1].t, sum);
            elif Length(other_args)=0 then
                return self._addPtrT(ptr_args);
            else
                return Error("Addition of more than one pointer and integers is not defined");
            fi;
	    fi;
    end;

mul.computeType := meth(self)
        local len, t, ptr_t, ptr_args, other_args, prod, args;
	args := self.args;

	len := Length(args);
	if   len=0  then return TInt;
	elif len=1  then return args[1].t;
	else
	    [ptr_args, other_args] := SplitBy(args, x->IsPtrT(x.t));
	    if Length(ptr_args)=0 then
	        if (self.args[1].t = TIVReal or self.args[2].t = TIVReal) then
  		  return TIVReal(T_Real(64));
		else
		  return UnifyTypesL(args);
		fi;
	    elif Length(ptr_args) > 1 then Error("Can't multiply pointers");
	    else
		prod := Product(other_args);
		if other_args<>[] and not IsIntT(prod.t) then Error("Can't multiply a pointer by a non-integer"); fi;
		return  self._ptrMul(ptr_args[1].t, prod);
	    fi;
	fi;
    end;

div.computeType := self >> When(self.args[1].t = TIVReal or self.args[2].t = TIVReal, TIVReal(T_Real(64)), UnifyTypes(List(self.args, x->x.t)));

TempArrayType := (child, y, x, index) ->
let(When(IsBound(child.a.t_in), child.a.t_in[index],
		let( X := Flat([x])[1], Y:= Flat([y])[1],
			Cond(
				IsBound(X.t.t) and X.t.t = TComplex, 		TComplex,
				IsBound(Y.t.t) and Y.t.t = TComplex, 		TComplex,
				ObjId(child) = ISumReduction, 				When(IsRec(child.idval), child.idval.t, TReal),
				ObjId(child) = PointWise, 					child.op.expr.t,
				X.t.t))));
TempArray := (y, x, child) -> let(cols := Flat([ Cols(child) ]),
	newType := TempArrayType(child, y, x, 1),
      StripList(
        List([ 1 .. Length(cols) ], (i) -> TempVec(TArray(TempArrayType(child,
               y, x, i), cols[i])))) );

max_stars_repo_name: spiral-software/spiral-package-hcol

max_issues_repo_name: spiral-software/spiral-package-hcol

max_forks_repo_name: spiral-software/spiral-package-hcol


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

_md := function(lst)
    local i, str;
    str:="";
    if Length(lst) = 0  then
        return str;
    fi;
    str := Concat(str, StringInt(lst[1]));
    for i in [ 2 .. Length(lst) ]  do
        str := Concat(str, "x", StringInt(lst[i]));
    od;
    return str;
end;

#   doSimdDft argument parser
#
#   spiral> doSimdDft();
#   [ [ 4 ], SSE_2x64f, rec(
#      svct := true,
#      oddSizes := false ) ]
#   spiral> doSimdDft(1);
#   [ [ 4, 8 ], SSE_2x64f, rec(
#      svct := true,
#      oddSizes := false ) ]
#   spiral> doSimdDft([4,6]);
#   [ [ 4, 6 ], SSE_2x64f, rec(
#      svct := true,
#      oddSizes := true ) ]
#   spiral> doSimdDft(SSE_4x32f);
#   [ [ 16 ], SSE_4x32f, rec(
#      svct := true,
#      oddSizes := false ) ]
#   spiral> doSimdDft(1, SSE_4x32f);
#   [ [ 16, 32 ], SSE_4x32f, rec(
#      svct := true,
#      oddSizes := false ) ]
#   spiral> doSimdDft([4,6], SSE_4x32f);
#   [ [ 4, 6 ], SSE_4x32f, rec(
#      svct := true,
#      oddSizes := true ) ]
#   spiral> doSimdDft([4,6], SSE_4x32f, rec(svct:=false, oddSizes:=true));
#   [ [ 4, 6 ], SSE_4x32f, rec(
#      svct := false,
#      oddSizes := true ) ]

_parse := function(arglist, multiplier)
    local isa, n, sizes, optrec;

    isa := SIMD_ISA_DB.active()[1];
    optrec := rec(svct := true, oddSizes :=false);

    if Length(arglist) = 0 then
        n := 0;
    elif Length(arglist) = 1 then
        if IsInt(arglist[1]) then
            n := arglist[1];
        elif IsList(arglist[1]) then
            n := arglist[1];
            optrec := rec(svct := true, oddSizes :=true);
        else
            n := 0;
            isa := arglist[1];
        fi;
    elif Length(arglist) = 2 then
        n := arglist[1];
        isa := arglist[2];
        if IsList(arglist[1]) then
            optrec := rec(svct := true, oddSizes :=true);
        fi;
    else
        n := arglist[1];
        isa := arglist[2];
        optrec := arglist[3];
    fi;

    sizes := When(IsList(n), n, List([0..n], i->multiplier*isa.v^2*2^i));

    return [sizes, isa, optrec];
end;

codeSimdDft := function(arg)
    local sizes, opts, isa, res, t,rt, s, c, optrec;
    [sizes, isa, optrec] := _parse(arg, 1);
    opts := SIMDGlobals.getOpts(isa, optrec);

    t := ComplexT(DFT(sizes[1]), opts).withTags(opts.tags);
    rt := RandomRuleTree(t, opts);
    s := SumsRuleTreeOpts(rt, opts);
    c := CodeSumsOpts(s, opts);
    return rec(opts:=opts,t:=t,rt:=rt,s:=s,c:=c);
end;

_doSimd := function(sizes, isa, opts, dpr, name, trafo, mflops, stringint)
    local results, res, i;

    results := [];
    if not IsBound(dpr.hashTable) then dpr.hashTable := HashTableDP(); fi;

    Print(name, ": running ", isa, ": ", sizes, "...\n");
    for i in sizes do
      res := TimedAction(DP(trafo(i), dpr, opts));
      if Length(res[1]) > 0 then
        Add(results, res[1]);
        ExportCodeRuleTree(res[1][1].ruletree,
            Concat(name, "_", isa.name, "_", stringint(i)),
            CopyFields(opts, rec(verbosity := 0, fileinfo := rec(
                cycles := res[1][1].measured,
                flops := mflops(i),
                file := Concat(name, "_", isa.name, "_", stringint(i)),
                algorithm := res[1][1].ruletree
        ))));

        if opts.verify then
            _seqPerfStatsMflopsAcc(Concat(name, "_", isa.name, ".txt"),
                i, mflops(i), res[1][1].measured, res[2], IntDouble(d_log(1e-16+CVerifyRT(res[1][1].ruletree, opts))/d_log(10)));
        else
            _seqPerfStatsMflops(Concat(name, "_", isa.name, ".txt"),
                i, mflops(i), res[1][1].measured, res[2]);
        fi;
        HashSave(dpr.hashTable, Concat(name, "_", isa.name, ".hash"));
      fi;
    od;

    return results;
end;

#F doSimdDft(sizes, isa, optrec)
doSimdDft := function(arg)
    local sizes, opts, dpr, isa, optrec, trafos, dpbench, direction;

    [sizes, isa, optrec] := _parse(arg, 1);

	direction := When(IsBound(optrec.transInverse) and optrec.transInverse, 1, -1);

    opts := SIMDGlobals.getOpts(isa, optrec);
    opts.benchTransforms := List(sizes, k -> ComplexT(DFT(k, direction), opts).withTags(opts.tags));

    if not IsBound(optrec.globalUnrolling) and ForAny(sizes, i -> not IsInt(i/isa.v^2)) then
        optrec.globalUnrolling := 16384;
    fi;
    dpr  := rec(verbosity := 0, timeBaseCases:=true);
    dpbench := spiral.libgen.DPBench(rec((Concat(isa.name, When(opts.vector.conf.interleavedComplex, "_ic", "_sc"))) := opts), dpr);
    if IsBound(optrec.verify) then dpbench.matrixVerify := optrec.verify; fi;

	dpbench.sizes := sizes;

    return dpbench;
end;

#F doSimdWht(sizes, isa, optrec)
doSimdWht := function(arg)
    local sizes, isa, opts, optrec, defrec, dpr, dpbench;

    [sizes, isa, optrec] := _parse(arg, 1);

    defrec := rec(
        oddSizes := false,
        svct := true,
        stdTTensor := true,
        tsplPFA := false,
        splitL:=false,
        realVect := true,
        cplxVect := false,
        useConj := false,
    );
    opts := SIMDGlobals.getOpts(isa, CopyFields(optrec, defrec));
    opts.benchTransforms := List(sizes, k -> WHT(Log2Int(k)).withTags(opts.tags));
    #opts.benchTransforms := List(sizes, k -> GT(WHT(Log2Int(k)), GTPar, GTPar, [4]).withTags(opts.tags));
    #opts.benchTransforms := List(sizes, k -> GT(WHT(Log2Int(k)), GTPar, GTPar, [8]).withTags(opts.tags));
    dpr  := rec(verbosity := 0, timeBaseCases:=true);

    dpbench := spiral.libgen.DPBench(rec((isa.name) := opts), dpr);
    if IsBound(optrec.verify) then dpbench.matrixVerify := optrec.verify; fi;

	dpbench.sizes := sizes;

    return dpbench;
end;

doSimdMddft := function(arg)
    local sizes, opts, tags, dpr, isa, defrec, optrec, dpbench, direction;

    [sizes, isa, optrec] := _parse(arg, 1);

	direction := When(IsBound(optrec.transInverse) and optrec.transInverse, 1, -1);

    # for cannot do larger sizes at this point
    defrec := rec(
        realVect := true,
        cplxVect := false,
        useConj := false,
    );
    opts := SIMDGlobals.getOpts(isa, CopyFields(optrec, defrec));

    opts.benchTransforms := List(sizes, k -> ComplexT(MDDFT(k, direction), opts).withTags(opts.tags));
    dpr  := rec(verbosity := 0, timeBaseCases:=true);

    dpbench := spiral.libgen.DPBench(rec((Concat(isa.name, When(opts.vector.conf.interleavedComplex, "_ic", "_sc"))) := opts), dpr);
    if IsBound(optrec.verify) then dpbench.matrixVerify := optrec.verify; fi;

	dpbench.sizes := sizes;

    return dpbench;
end;

codeSimdSymDFT := function(arg)
    local sizes, opts, dpr, isa, defrec, optrec, t, r, s, c, ttype;

    ttype := arg[1];
    [sizes, isa, optrec] := _parse(Drop(arg, 1), 2);
    defrec := rec(oddSizes := false, svct := true, splitL:=false, realVect := true,
                  cplxVect := false, useConj := false, globalUnrolling := 1024);
    opts := SIMDGlobals.getOpts(isa, CopyFields(optrec, defrec));
    t := ttype(sizes[1]).withTags(opts.tags);
    r := RandomRuleTree(t, opts);
    s := SumsRuleTree(r, opts);
    c := CodeSums(s, opts);
    return rec(opts:=opts, t:=t, rt:=r, s:=s, c:=c);
end;

doSimdSymDFT := function(arg)
    local sizes, opts, dpr, isa, optrec, results, defrec, ur, ttype, dpbench;

    ttype := arg[1];
    [sizes, isa, optrec] := _parse(Drop(arg, 1), 2);

    defrec := rec(
        oddSizes := false,
        svct := true,
        splitL:=false,
        realVect := true,
        cplxVect := false,
        useConj := false,
        globalUnrolling := 1024
    );
    ur := rec(
        globalUnrollingMin := 1024,
        globalUnrollingMax := 1024
    );
    opts := SIMDGlobals.getOpts(isa, CopyFields(defrec, optrec));
    opts.benchTransforms := List(sizes, k -> ttype(k).withTags(Cond(IsBound(optrec.cplxVect) and optrec.cplxVect, opts.cxtags, opts.tags)));
    if not IsBound(optrec.globalUnrolling) and ForAny(sizes, i -> not IsInt(i/isa.v^2)) then
        optrec.globalUnrolling := 16384;
    fi;
    dpr  := rec(verbosity := 0, timeBaseCases:=true);

    dpbench := spiral.libgen.DPBench(rec((isa.name) := opts), dpr);
    if IsBound(optrec.verify) then dpbench.matrixVerify := optrec.verify; fi;

	dpbench.sizes := sizes;

    return dpbench;

#    results := _doSimd(sizes, isa, opts, dpr, arg[1], k -> ttype(k).withTags(opts.tags), i->2*i*d_log(i)/d_log(2) + i, StringInt);
#    return rec(opts:=opts, results:=results);
end;

doSimdSymMDDFT := function(arg)
    local sizes, isa, opts, optrec, defrec, dpr, dpbench, ttype;

    ttype := arg[1];
    [sizes, isa, optrec] := _parse(Drop(arg, 1), 2);

    defrec := rec(svct := true, splitL := false, stdTTensor := true, tsplRader:=false, tsplBluestein:=false,
        tsplPFA:=false, oddSizes:=false, interleavedComplex := false, pushTag := false, flipIxA := true);
    opts := SIMDGlobals.getOpts(isa, CopyFields(optrec, defrec));

    opts.benchTransforms := List(sizes, k -> TTensor(ttype(k), ttype(k)).withTags(opts.tags));
    dpr  := rec(verbosity := 0, timeBaseCases:=true);

    dpbench := spiral.libgen.DPBench(rec((isa.name) := opts), dpr);
    if IsBound(optrec.verify) then dpbench.matrixVerify := optrec.verify; fi;

	dpbench.sizes := sizes;

    return dpbench;
end;

doSimdMdConv := function(arg)
    local sizes, isa, opts, optrec, defrec, dpbench, transforms;

    [sizes, isa, optrec] := _parse(arg, 1);

    defrec := rec();

    opts := SIMDGlobals.getOpts(isa, CopyFields(optrec, defrec));

	transforms := Flat(List(sizes, n -> let(tt := TRConv2D(ImageVar([n, n])).withTags(opts.tags), [tt, tt.forwardTransform()])));

	dpbench := spiral.libgen.DPBench.build(transforms, opts);

	dpbench.sizes := sizes;

    return dpbench;
end;

ext: gi


text:
#############################################################################
####
##
#W  anupqid.gi              ANUPQ package                       Werner Nickel
#W                                                                Greg Gamble
##
##  This file installs functions to do with evaluating identities.
##
#Y  Copyright (C) 2001  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
##

#############################################################################
##
#F  PqEvalSingleRelation( <proc>, <r>, <instances> )
##
InstallGlobalFunction( PqEvalSingleRelation, function( proc, r, instances )
    local   w, datarec;

#    Print( instances, "\n" );

    datarec := ANUPQDataRecord(proc);
    datarec.nwords := datarec.nwords + 1;
    w := CallFuncList( r, instances );
    if w <> w^0 then
#        Print( w );
#        Print( "\n" );
#        PqSetOutputLevel( proc, 3 );
        PqCollect( proc, String(w) );
        PqEchelonise( proc );
#        PqSetOutputLevel( proc, 0 );
    fi;
end );

#############################################################################
##
#F  PqEnumerateWords( <proc>, <data>, <r> )
##
##    The parameters of PqEnumerateWords() have the following meaning:
##
##    <r>          a relation involving identical generators.
##    <instances>  the list of instances to be built up corresponding the
##                 identical generators.
##    <n>          the index of the current word in <instances>
##    <g>          the next generator in the current word
##    <wt>         the weight that can be spent on the next generators.
##
InstallGlobalFunction( PqEnumerateWords, function( proc, data, r )
    local   n,  g,  wt,  u,  save_wt,  save_u,  save_g;

    n  := data.currentinst;
    g  := data.currentgen;
    wt := data.weightleft;
    u  := data.instances[ n ];

    save_wt := wt;
    save_u  := u;
    save_g  := g;

    if wt = 0 then
        PqEvalSingleRelation( proc, r, data.instances );
        Info( InfoANUPQ, 3, "Instance: ", data.instances[ n ] );
        return;
    fi;

    if g > Length( data.pcgens ) then return; fi;

    while g <= data.nrpcgens and data.pcweights[g] <= wt do
        while data.pcweights[g] <= wt do
            u := u * data.pcgens[g];
            wt := wt - data.pcweights[g];

            data.instances[ n ] := u;
            data.weightleft     := wt;
            data.currentgen     := g+1;
            PqEnumerateWords( proc, data, r );

            if n < Length(data.instances) then
                data.currentinst := n+1;
                data.currentgen  := 1;
                PqEnumerateWords( proc, data, r );
                data.currentinst := n;
            fi;
        od;
        u  := save_u; wt := save_wt; g := g+1;
    od;
    data.instances[ n ] := save_u;
    data.weightleft     := save_wt;
    data.currentgen     := save_g;
end );

#############################################################################
##
#F  PqEvaluateIdentity( <proc>, <r>, <arity> )
##
InstallGlobalFunction( PqEvaluateIdentity, function( proc, r, arity )
    local   n,  class,  gens,  data,  c;

    n     := PqNrPcGenerators( proc );
    class := PqWeight( proc, n );
    if class > 1 then
        while n > 0 and PqWeight( proc, n ) = class do n := n-1; od;
    fi;

    if n = 0 then return; fi;

    gens := GeneratorsOfGroup( FreeGroup( n, "x" ) );
    data := rec( instances   := List( [1..arity], i->gens[1]^0 ),
                 currentinst := 1,
                 currentgen  := 1,
                 weightleft  := 0,
                 pcgens      := gens,
                 nrpcgens    := n,
                 pcweights   := List( [1..n], i->PqWeight( proc, i ) )
                 );

    for c in [1..class] do
#        Print( "words of class ", c, "\n" );
        data.weightleft := c;
        PqEnumerateWords( proc, data, r );
    od;

end );

#############################################################################
##
#F  PqWithIdentity( <G>, <p>, <Cl>, <identity> )
##
##  constructs a <p>-quotient <Q> of the fp or pc group <G> of class at  most
##  <Cl> that satisfies $<identity>(<w1>, \dots, <wn>) =  1$  for  all  words
##  $<w1>, \dots, <wn>$ in the pc generators of <Q>. The  following  examples
##  demonstrate its usage.
##
##  \beginexample
###  gap> F := FreeGroup(2);
##  <free group on the generators [ f1, f2 ]>
##  gap> f := w -> w^4;
##  function( w ) ... end
##  gap> PqWithIdentity( F, 2, 20, f );
##  #I  Evaluated 5 instances.
##  #I  Class 2 with 5 generators.
##  #I  Evaluated 18 instances.
##  #I  Class 3 with 7 generators.
##  #I  Evaluated 44 instances.
##  #I  Class 4 with 10 generators.
##  #I  Evaluated 95 instances.
##  #I  Class 5 with 12 generators.
##  #I  Evaluated 192 instances.
##  #I  Class 6 with 12 generators.
##  <pc group of size 4096 with 12 generators>
##  gap>
##  gap> G := F/[ F.1^11, F.2^11 ];
##  <fp group on the generators [ f1, f2 ]>
##  gap> f := function(u, v) return PqLeftNormComm( [u, v, v, v] ); end;
##  function( u, v ) ... end
##  gap> H := PqWithIdentity( G, 11, 20, f );
##  #I  Evaluated 14 instances.
##  #I  Class 2 with 3 generators.
##  #I  Evaluated 44 instances.
##  #I  Class 3 with 5 generators.
##  #I  Evaluated 122 instances.
##  #I  Class 4 with 5 generators.
##  <pc group of size 161051 with 5 generators>
##  gap> f( Random(H), Random(H) );
##  <identity> of ...
##  gap> f( H.1, H.2 );
##  <identity> of ...
##  \endexample
##
##  Compare    the    above    examples    with    those    generated     by:
##  `PqExample("B2-4-Id");' and `PqExample("11gp-3-Engel-Id");' which do  the
##  same as above using the function `Pq' with the `Identities' option.
##
##  `PqWithIdentity' and the functions it calls,  with  minor  modifications,
##  constitute the prototype provided by Werner Nickel,  for  constructing  a
##  quotient that satisfies an identity. The prototype  functions  have  been
##  merged into the single function `PQ_EVALUATE_IDENTITY' which is called by
##  `PQ_EVALUATE_IDENTITIES'    which    in     turn     is     called     by
##  `PQ_FINISH_NEXT_CLASS'   which   is   called   by   `PQ_NEXT_CLASS'   and
##  `PQ_EPI_OR_PCOVER',  if   the   `Identities'   option   has   been   set.
##  `PQ_EPI_OR_PCOVER' is the function called  by  `Pq',  `PqEpimorphism'  or
##  `PqPCover'.
##
InstallGlobalFunction( PqWithIdentity, function( G, p, Cl, identity )
    local   proc,  datarec,  prev_n,  class,  grp, arity;

    arity := NumberArgumentsFunction( identity );

    proc := PqStart( G : Prime := p );
    datarec := ANUPQData.io[ proc ];
    datarec.nwords := 0;

    prev_n := 0;
    Pq( proc : ClassBound := 1 );    class := 1;

    PqEvaluateIdentity( proc, identity, arity );
    PqEliminateRedundantGenerators( proc );

    if PqNrPcGenerators( proc ) = 0 then
        return TrivialGroup( IsPcGroup );
    fi;

    while class < Cl and prev_n <> PqNrPcGenerators( proc ) do
        prev_n := PqNrPcGenerators( proc );

        PqSetupTablesForNextClass( proc );
        PqTails( proc, 0 );
        PqDoConsistencyChecks( proc, 0, 0 );
        PqCollectDefiningRelations( proc );

        PqDoExponentChecks( proc );

        datarec.nwords := 0;
        PqEvaluateIdentity( proc, identity, arity );
        Info(InfoANUPQ, 2, "Evaluated ", datarec.nwords, " instances." );

        PqEliminateRedundantGenerators( proc );

        class := class + 1;

        Info(InfoANUPQ, 1, "Class ", class, " with ",
                           PqNrPcGenerators(proc), " generators." );
    od;

    grp := PqCurrentGroup( proc );

    PqQuit( proc );

    return grp;
end );

#############################################################################
##
#F  PQ_EVALUATE_IDENTITY( <proc>, <identity> )
##
InstallGlobalFunction( PQ_EVALUATE_IDENTITY, function( proc, identity )
    local   EnumerateWords, data, datarec, nwords, arity, c;

    EnumerateWords := function()
      local   i,  g,  wt,  u,  w,  save_wt,  save_u,  save_g;

      i  := data.currentinst;
      g  := data.currentgen;
      wt := data.weightleft;
      u  := data.instances[ i ];

      save_wt := wt;
      save_u  := u;
      save_g  := g;

      if wt = 0 then
        #evaluate a single relation
        nwords := nwords + 1;
        w := CallFuncList( identity, data.instances );
        if w <> w^0 then
#         Print( w );
#         Print( "\n" );
#         PqSetOutputLevel( proc, 3 );
          PqCollect( proc, String(w) );
          PqEchelonise( proc );
#         PqSetOutputLevel( proc, 0 );
        fi;
        Info( InfoANUPQ, 3, "Instance: ", data.instances[ i ] );
        return;
      fi;

      while g <= data.nrpcgens and data.pcweights[g] <= wt do
        while data.pcweights[g] <= wt do
            u  := u * data.pcgens[g];
            wt := wt - data.pcweights[g];

            data.instances[ i ] := u;
            data.weightleft     := wt;
            data.currentgen     := g+1;
            EnumerateWords();

            if i < Length(data.instances) then
                data.currentinst := i+1;
                data.currentgen  := 1;
                EnumerateWords();
                data.currentinst := i;
            fi;
        od;
        u  := save_u; wt := save_wt; g := g+1;
      od;
      data.instances[ i ] := save_u;
      data.weightleft     := save_wt;
      data.currentgen     := save_g;
    end;

    data := rec( nrpcgens := PqNrPcGenerators( proc ) );
    if data.nrpcgens = 0 then
      return;
    fi;

    datarec  := ANUPQDataRecord( proc );
    if datarec.class > 1 then
      data.nrpcgens := datarec.ngens[ Length(datarec.ngens) - 1 ];
    fi;

    nwords := 0;
    arity  := NumberArgumentsFunction(identity);

    data.pcgens      := GeneratorsOfGroup( FreeGroup( data.nrpcgens, "x" ) );
    data.instances   := List( [1 .. arity], i -> data.pcgens[1]^0 );
    data.currentinst := 1;
    data.currentgen  := 1;
    data.pcweights   := List( [1 .. data.nrpcgens], i -> PqWeight( proc, i ) );

    for c in [1 .. datarec.class] do
#        Print( "words of class ", c, "\n" );
      data.weightleft := c;
      EnumerateWords();
    od;
    Info(InfoANUPQ, 2, "Evaluated ", nwords, " instances." );
end );

#E  anupqid.gi  . . . . . . . . . . . . . . . . . . . . . . . . . . ends here

max_stars_repo_name: gap-system/anupq

max_issues_repo_name: gap-system/anupq

max_forks_repo_name: gap-system/anupq


text:
CollatzSequence := function(n)
  local v;
  v := [ n ];
  while n > 1 do
    if IsEvenInt(n) then
      n := QuoInt(n, 2);
    else
      n := 3*n + 1;
    fi;
    Add(v, n);
  od;
  return v;
end;

CollatzLength := function(n)
  local m;
  m := 1;
  while n > 1 do
    if IsEvenInt(n) then
      n := QuoInt(n, 2);
    else
      n := 3*n + 1;
    fi;
    m := m + 1;
  od;
  return m;
end;

CollatzMax := function(a, b)
  local n, len, nmax, lmax;
  lmax := 0;
  for n in [a .. b] do
    len := CollatzLength(n);
    if len > lmax then
      nmax := n;
      lmax := len;
    fi;
  od;
  return [ nmax, lmax ];
end;

CollatzSequence(27);
# [ 27, 82, 41, 124, 62, 31, 94, 47, 142, 71, 214, 107, 322, 161, 484, 242, 121, 364, 182, 91, 274, 137, 412, 206,
#   103, 310, 155, 466, 233, 700, 350, 175, 526, 263, 790, 395, 1186, 593, 1780, 890, 445, 1336, 668, 334, 167, 502,
#   251, 754, 377, 1132, 566, 283, 850, 425, 1276, 638, 319, 958, 479, 1438, 719, 2158, 1079, 3238, 1619, 4858, 2429,
#   7288, 3644, 1822, 911, 2734, 1367, 4102, 2051, 6154, 3077, 9232, 4616, 2308, 1154, 577, 1732, 866, 433, 1300,
#   650, 325, 976, 488, 244, 122, 61, 184, 92, 46, 23, 70, 35, 106, 53, 160, 80, 40, 20, 10, 5, 16, 8, 4, 2, 1 ]
CollatzLength(27);
# 112

CollatzMax(1, 100);
# [ 97, 119 ]
CollatzMax(1, 1000);
# [ 871, 179 ]
CollatzMax(1, 10000);
# [ 6171, 262 ]
CollatzMax(1, 100000);
# [ 77031, 351 ]
CollatzMax(1, 1000000);
# [ 837799, 525 ]

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
# Using functions from 'RandomPerm.gi' and 'Orbit.gi'
Read("RandomPerm.gi"); Read("Orbit.gi");

#######################################################################
# Description: A function to calculate the pointwise image of a list. #
# Input: - list of integers B = [b_1,...,b_k]                         #
#        - permutation x as cycle                                     #
# Output: - list of integers [b_1^x,...,b_k^x]                        #
#######################################################################

BaseImage := function(B, x)
    local i, image;
    image := [];
    for i in [1..Length(B)] do
        Add(image, B[i]^x);
    od;
    return image;
end;

#######################################################################
#######################################################################

#######################################################################
# Description: A function to setup partial generating sets of         #
#              stabilizeres.                                          #
# Input: - list of integers B = [b_1,...,b_k]                         #
#        - generating system S as list of permutations                #
# Output: - list of generating system (lists of permutations),        #
#           one generating system for each stabilizer in chain        #
#######################################################################

GeneratingSets := function(B, S)
    local s, x, i, partialBase, sets;
    sets := [S];
    for i in [1..Length(B)-1] do
        Add(sets, []);
    od;
    for x in S do
        for i in [1..Length(B)-1] do
            partialBase := B{[1..i]};
            if BaseImage(partialBase, x) = partialBase then
                Add(sets[i+1], x);
            else
                break;
            fi;
        od;
    od;
    return sets;
end;

#######################################################################
#######################################################################

#######################################################################
# Description: A function to strip a permutation w.r.t. a base,       #
#              described in section 3.2.                              #
# Input: - permutation g as cycle                                     #
#        - list of integers B = [b_1,...,b_k]                         #
#        - generating system S as list of permutations                #
#        - data structure returned by 'GeneratingSets' function       #
# Output: - residue and dropout level of g                            #
#######################################################################

Strip := function(g, B, orbits, sgs)
    local m, i, u, b;
    m := g;
    for i in [1..Length(B)] do
        b := B[i]^m;
        if orbits[i][2][b] = 0 then
            return [m, i];
        fi;
        u := Transversal(b, orbits[i][2], sgs[i]);  m := m*(u^(-1));
    od;
    return [m, Length(B)+1];
end;

#######################################################################
#######################################################################

#######################################################################
# Description: The standard deterministic Schreier-Sims-Algorithm     #
#              described in section 3.3.                              #
# Input: - list B of base points (might be empty list [] first)       #
#        - list S of given generators of group G                      #
#        - integer n specifying the degree of symmetric group the     #
#          generated group G is embedded in                           #
# Output: - list 'BSGS' with three entries including                  #
#           -- the calculated base B as list in BSGS[1]               #
#           -- a strong generating set S as list of partial           #
#              generating sets for each stabilizer in BSGS[2]         #
#           -- a list of the sizes of the orbits belonging to B       #
#              and S in BSGS[3]                                       #
#######################################################################

SchreierSims := function(B, S, n)
    local x, u, y, h, b, i, j, k, t, orbits, sgs, sizes, stripped;

    # extending partial base so that no point of S
    # fixes all points
    for x in S do
        if BaseImage(B, x) = B then
            Add(B, MovedPoints(x)[1]);
        fi;
    od;
    # calculating generating sets for each stabilizer
    sgs := GeneratingSets(B, S);

    # calculate orbits and schreier vectors of base
    # points under action of generated stabilizers
    orbits := [];
    for i in [1..Length(B)] do
        Add(orbits, OrbitSV(B[i], sgs[i], n));
    od;
    i := Length(B);
    while i >= 1 do
        for b in orbits[i][1] do
            t := Transversal(b, orbits[i][2], sgs[i]);
            for x in sgs[i] do
                u := Transversal(b^x, orbits[i][2], sgs[i]);  y := true;
                if not t*x = u then
                    # check whether Schreier generator
                    # belongs to H^(i+1)
                    stripped := Strip(t*x*(u^(-1)), B, orbits, sgs);
                    j := stripped[2];  h := stripped[1];

                    # in case the Schreier generator does
                    # not belong to H^(i+1) do ...
                    if j <= Length(B) then
                        y := false;
                    elif not h = () then
                        Add(B, MovedPoints(h)[1]);
                        Add(sgs, []);  y := false;
                    fi;
                    if y = false then
                        # update generating sets, orbits
                        # and Schreier vectors
                        for k in [i+1..j] do
                            Add(sgs[k], h);
                            orbits[k] := OrbitSV(B[k], sgs[k], n);
                        od;
                        i := j;  break;
                    fi;
                fi;
            od;
            if y = false then
                break;
            fi;
        od;
        if y = false then
            continue;
        fi;
        i := i - 1;
    od;
    sizes := [];
    for i in [1..Length(B)] do
       Add(sizes, Length(orbits[i][1]));
    od;
    return [B, sgs, sizes];
end;

#######################################################################
#######################################################################

#######################################################################
# Description: The randomized Schreier-Sims-Algorithm,                #
#              described in section 4.2.                              #
# Input: - list B of base points (might be empty list [] first)       #
#        - list S of given generators of group G                      #
#        - an integer 'bound' capping the number of elements          #
#          getting stripped                                           #
#          with trivial residue in a row                              #
#        - integer n specifying the degree of symmetric group the     #
#          generated group G is embedded in                           #
# Output: - list 'BSGS' with three entries including                  #
#           -- the calculated base B as list in BSGS[1]               #
#           -- a strong generating set S as list of partial           #
#              generating sets for each stabilizer in BSGS[2]         #
#           -- a list of the sizes of the orbits belonging to B       #
#              and S in BSGS[3]                                       #
#######################################################################

RandomSchreierSims := function(B, S, bound, n)
    local x, g, i, h, X, c, j, k, y, sgs, orbits, sizes, SGS, stripped;

    for x in S do
        if BaseImage(B, x) = B then
            Add(B, MovedPoints(x)[1]);
        fi;
    od;

    sgs := GeneratingSets(B, S);  orbits := [];  SGS := ShallowCopy(S);

    for i in [1..Length(B)] do
        Add(orbits, OrbitSV(B[i], sgs[i], n));
    od;

   X := InitGenerator(S);  c := 0;
   while c < bound do
       g := RandomPerm(S, X);
       stripped := Strip(g, B, orbits, sgs);
       j := stripped[2];  h := stripped[1];  y := true;
       if j <= Length(B) then
           y := false;
       elif not h = () then
           y := false;  Add(B, MovedPoints(h)[1]);  Add(sgs, []);
       fi;
       if y = false then
           Add(SGS, h);
           for k in [2..j] do
               Add(sgs[k], h);
               orbits[k] := OrbitSV(B[k], sgs[k], n);
           od;
           c := 0;
       fi;
       c := c + 1;
   od;

   sizes := [];  sgs[1] := SGS;
   for i in [1..Length(B)] do
       Add(sizes, Length(orbits[i][1]));
   od;
   Add(sgs, [()]);
   return [B, sgs, sizes];
end;

#######################################################################
#######################################################################

#######################################################################
# Description: A function to remove redundant generators from         #
#              a strong generating set, described in section 4.2.     #
# Input: - BSGS data structure returned by RandomSchreierSims         #
#        - Generating set which defined the group. NOT the strong     #
#          generating set returned by RandomSchreierSims!             #
#        - integer n specifying the size of the set the group         #
#          operates on
#######################################################################

RemoveRedundantGens := function(BSGS, initialGenSet, n)
    local orbits, i, g, newGenSys, removeGen;
    orbits := [];

    removeGen := function(SGS, gen, level)
        local j;
        for j in [1..level] do
            RemoveSet(SGS[j], gen);
        od;
    end;
    for i in [1..Length(BSGS[1])] do
        BSGS[2][i] := Set(BSGS[2][i]);
        Add(orbits, Set(OrbitSimple(BSGS[1][i], BSGS[2][i])));
    od;
    for i in Reversed([1..Length(BSGS[1])-1]) do
        for g in BSGS[2][i] do
            if not g in BSGS[2][i+1]
               and not g in initialGenSet then
                newGenSys := Set(BSGS[2][i]);
                RemoveSet(newGenSys, g);
                if Set(OrbitSimple(BSGS[1][i], newGenSys)) = orbits[i]
                	then removeGen(BSGS[2], g, i);
                fi;
            fi;
        od;
    od;
    orbits := [];
    for i in [1..Length(BSGS[1])] do
        Add(orbits, OrbitSV(BSGS[1][i], BSGS[2][i], n));
    od;
    Add(BSGS, orbits);  MakeImmutable(BSGS);
end;

#######################################################################
#######################################################################

#######################################################################
# Description: The sequential deterministic verification,             #
#              described in section 4.3                               #
# Input: - BSGS data structure returned by RandomSchreierSims         #
#        - integer m specifying which layer H^{(m)}_{b_m} = H^{(m+1)} #
#          has to be verified                                         #
# Output: - boolean, true if H^{(m)}_{b_m} = H^{(m+1)} holds,         #
#           else false                                                #
#######################################################################

SimpleVerify := function(BSGS, m)
    local x, u, b, t, orbits, gen, stripped, sgs, B;

    B := BSGS[1]; sgs := BSGS[2]; orbits := BSGS[4];

    for b in orbits[m][1] do
         t := Transversal(b, orbits[m][2], sgs[m]);
         for x in sgs[m] do
            if not x in sgs[m+1] or orbits[m+1][2][b] = 0 then
                 u := Transversal(b^x, orbits[m][2], sgs[m]);
                 if not t*x = u then
                    gen := t*x*(u^(-1));
                    # check whether Schreier generator 'gen'
                    # belongs to H^(i+1)
                    stripped := Strip(gen, B, orbits, sgs);
                    # if not, return
                    if stripped[2] <= Length(B) then
                        return false;
                    elif not stripped[1] = () then
                        return false;
                    fi;
                 fi;
              fi;
         od;
    od;
    return true;
end;

#######################################################################
#######################################################################

#######################################################################
# Description: The deterministic verification, parallelized w.r.t.    #
#              the length of the basic orbits, described as method 2  #
#              in section 5.1.                                        #
# Input: - BSGS data structure returned by RandomSchreierSims         #
#        - integer m specifying which layer H^{(m)}_{b_m} = H^{(m+1)} #
#          has to be verified                                         #
# Output: - boolean, true if H^{(m)}_{b_m} = H^{(m+1)} holds,         #
#           else false                                                #
#######################################################################

ParallelVerify := function(BSGS, m)
    local b, i, orbits, sgs, B, Verify, tasks, results, numOfTasks,
          remaining, partition, pointsPerTask;

    # the well-known deterministic verify routine
    Verify := function(B, sgs, orbits, m, partition)
        local x, u, i, stripped, gen, t;
        for b in partition do
            t := Transversal(b, orbits[m][2], sgs[m]);
                 for x in sgs[m] do
                    if not x in sgs[m+1] or orbits[m+1][2][b] = 0 then
                         u := Transversal(b^x, orbits[m][2], sgs[m]);
                         if not t*x = u then
                            gen := t*x*(u^(-1));
                            stripped := Strip(gen, B, orbits, sgs);

                            if stripped[2] <= Length(B) then
                                return false;
                            elif not stripped[1] = () then
                                return false;
                            fi;
                         fi;
                      fi;
                 od;
        od;
        return true;
    end;

    B := BSGS[1]; sgs := BSGS[2];
    orbits := BSGS[4]; partition := [];

    pointsPerTask := Int(Length(orbits[m][1])*0.05) + 1;
    numOfTasks := Int(Length(orbits[m][1])/pointsPerTask) + 1;
    remaining := Length(orbits[m][1]) mod pointsPerTask;

    tasks := [1..numOfTasks]; results := [1..numOfTasks];

    # slice orbit into 'numOfTask' many nearly equally sized subsets
    for i in [1..numOfTasks-1] do
        Add(partition, orbits[m][1]{[(i-1)*pointsPerTask+1..
                                         i*pointsPerTask]});
    od;

    Add(partition, orbits[m][1]{[(numOfTasks-1)*pointsPerTask+1..
                    ((numOfTasks-1)*pointsPerTask+remaining)]});

    # start several verifications ...
    for i in [1..numOfTasks] do
        tasks[i] := RunTask(Verify, B, sgs, orbits, m, partition[i]);
    od;

    # ... and evaluate their results afterwards
    for i in [1..numOfTasks] do
        results[i] := TaskResult(tasks[i]);
    od;
    if false in results then
        return false;

    else
        return true;
    fi;
end;

max_stars_repo_name: AlexanderKlemps/parallel-schreier-sims

max_issues_repo_name: AlexanderKlemps/parallel-schreier-sims

max_forks_repo_name: AlexanderKlemps/parallel-schreier-sims


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

processCRP := meth(partA, partB)
    local size, switches, next, modA, modB, loc, setA, choice;
    size := Length(partA);

    modA := List([1..size], i-> Mod(partA[i], size));
    modB := List([1..size], i-> Mod(partB[i], size));

    setA := [];

    switches := List([1..size], i->-1);
    next := -1;

    while (-1 in switches) do

        # if next = -1, then we can make an arbitrary decision.
        if (next < 0) then
            # find the next switch that isn't set.
            choice := PositionProperty(modA, i-> i <> -1);

            switches[choice] := 0;
            Append(setA, [modA[choice]]);

            # We just assigned modA[choice] to partition A.
            # Thus, we precluded modB[choice] from being in parition A.
            # So, we need to check.  If modB[choice] is not already in
            # partition A, then we need to make it our next choice.

            if (modB[choice] in setA) then
                next := -1;
            else
                next := modB[choice];
            fi;
            modA[choice] := -1; modB[choice] := -1;

        else
            # if next >= 0, then it represents the 'mod' number we must
            # assign to partition A.

            # the next value must be in the set of unchosen values in one
            # of the two partitions.
            if (next in modA) then
                choice := Position(modA, next);
                switches[choice] := 0;
                Append(setA, [modA[choice]]);

                if (modB[choice] in setA) then
                    next := -1;
                else
                    next := modB[choice];
                fi;

            else
                choice := Position(modB, next);
                switches[choice] := 1;
                Append(setA, [modB[choice]]);

                if (modA[choice] in setA) then
                    next := -1;
                else
                    next := modA[choice];
                fi;
            fi;
            modA[choice] := -1; modB[choice] := -1;
        fi;

    od;

    return switches;
end;

# Given a permutation permIn, determine the switch settings for the
# un-simplified \Omega\Omega^{-1} network.
#
# See: K. Y. Lee.  On the rearrangeability of 2(\log_2n)-1 stage
# permutation networks.  IEEE Transactions on Computers,
# C-34(5):412-425, May 1985.

# We do not simplify the network, because my network generation code
# will simplify all unneeded switches.

routePerm := meth(permIn)
    local n, ss, ns, p, p2, col, numP, partSize, partA, partB, i, j, k, partStride, tmpSw, swCol, tmp;

#    Print("Routing perm: ", permIn, "\n");

    n := Length(permIn);

    # The perm comes in as row values.  I.e., PermRowVals[i] -> i.
    # We need to transpose this.
    p := List([1..n], i->0);
    for i in [0..n-1] do p[permIn[i+1]+1] := i; od;

    # This is suboptimal.  For now, I will simply let the unnecessary
    # switches be all-zero.
    ns := n/2 * (LogInt(n,2) + LogInt(n,2)-1);

    ss := [];

    for k in [0..LogInt(n,2)-2] do
        col := LogInt(n,2)-1-k;

        # Permute output values by L(n,n/2)
        p2 := List([1..n], i -> p[L(n,n/2).at(i-1).ev()+1]);

        # Now we need to set the switches for the E^{col}_{R} column.
        numP := 2^k;
        partSize := n/(2*numP);
        partStride := 2*numP;

        swCol := List([1..n/2], i -> -1);
        tmpSw := [];
        for i in [0..numP-1] do
            partA := [];
            partB := [];
            for j in [1..partSize] do
                Append(partA, [p2[partStride*(j-1)+1+(2*i)]]);
                Append(partB, [p2[partStride*(j-1)+2+(2*i)]]);
            od;

            Append(tmpSw, processCRP(partA, partB));

        od;

        swCol := List([1..n/2], i -> tmpSw[L(n/2,n/(2^(k+1))).at(i-1).ev()+1]);
        Append(ss, swCol);

        # Now we have to perform the permutation based on the
        # switches we set.
        for i in [0..n/2-1] do
            if (swCol[i+1] = 1) then
                p[i*2+1] := p2[i*2+2];
                p[i*2+2] := p2[i*2+1];
            else
                p[i*2+1] := p2[i*2+1];
                p[i*2+2] := p2[i*2+2];
            fi;
        od;
    od;

    # Permute output values by L(n,n/2)
    p2 := List([1..n], i -> p[L(n,n/2).at(i-1).ev()+1]);

    for k in [0..LogInt(n,2)-1] do

        swCol := [];
        # go through the n/2 switches.  at each switch, if the top input's value
        # has a 1 in bit k.
        for i in [0..n/2-1] do
            if (imod(idiv(p2[2*i+1], 2^k), 2).ev() = 1) then
                Append(swCol, [1]);
            else
                Append(swCol, [0]);
            fi;
        od;

        Append(ss, swCol);

        # Now we have to perform the permutation based on the
        # switches we set.
        for i in [0..n/2-1] do
            if (swCol[i+1] = 1) then
                p[i*2+1] := p2[i*2+2];
                p[i*2+2] := p2[i*2+1];
            else
                p[i*2+1] := p2[i*2+1];
                p[i*2+2] := p2[i*2+2];
            fi;
        od;

        # Permute output values by L(n,2)
        p2 := List([1..n], i -> p[L(n,2).at(i-1).ev()+1]);

    od;

    p := List([0..n-1], i->i);
    if (p <> p2) then
        Error("Failed at configuring switching network.");
    fi;

#    Print("Configuration: ", ss, "\n");

    return ss;
end;

# Given a permutation permIn, determine the switch settings for the
# \Omega^{-1} network.
#
# We will return -1 if the permutation cannot be mapped to that network.
# In this case, then, you should use routePerm to map it to the
# \Omega\Omega^{-1} network.

# For now, we will cheat and use the \Omega\Omega^{-1} routing algorithm.
# It appears that when we use this algorithm on a permutation that can be
# permformed on the \Omega^{-1}, then the configuration bits are simply
# all zero for the first \Omega portion of the network.  However, I haven't
# verified this and I haven't looked at the routePerm code in many years.

# I have not proven this will be the case, but it seems to work on the
# circular shift permutations I am currently using.

# All circular shifts will be Omega^{-1} passable by Theorems 3  and 13 of "Access
# and Alignment of Data in an Array Processor," Duncan H. Lawrie, IEEE Tr.
# on Computers, C-24 (12), 1975.
#
# Theorem 3: Duncan says that if you can represent your mapping as:
# P_N(a, b, c, d, e) = {(ax+b, cx+d) | 0 \leq x < e}
# where ax+b represents the source (mod N) and cx+d represents the
# destination (mod N),
# then your mapping is \Omega passable if gcd(a,N) \leq gcd(c,N) and
# e \leq N/(gcd(c,N).
# For a cyclic shift, a and c will always = 1, b=0, d equals the number
# of positions to shift, and e=N.
# Since a and c are 1, gcd(a,N) and gcd(c,N) are both 1 and hence equal.
# Further, N/gcd(c,N) = N, and therefore e = N/gcd(c,N).

# Lastly we note that Theorem 13 says that if a mapping is \Omega passable
# it must be \Omega^{-1} passable.

routePermOmegaInv := meth(permIn)
    local full, n, ns, numSwitchesZero;

    full := routePerm(permIn);
    n := Length(permIn);
    ns := n/2 * (LogInt(n,2) + LogInt(n,2)-1);

    # If this permutation is Omega^{-1} routable, then the first numSwitchesZero
    # must be configured to zero.
    numSwitchesZero := n/2 * (LogInt(n,2)-1);

    if (Sum(full{[1..numSwitchesZero]}) <> 0) then
	return -1;
    fi;

    return full{[numSwitchesZero+1 .. ns]};
end;

cyclicShiftList := meth(size, shift)
    local t, res, t2;
    t := [0..size-1];
    shift := Mod(shift, size);
    res := t{[shift+1..size]};
    if (shift = 1) then
	t2 := [0];
    else
	t2 := t{[1..shift]};
    fi;

    Append(res, t2);
    return res;
end;

Class(omegaInvNetwork, BaseMat, SumsBase, rec(

    abbrevs := [(n, p, it) -> [n, p, it]],

    new := (self, n, perms) >> SPL(WithBases(self, rec(
		    dimensions := [n, n],
		    _children := [n, perms]
	     ))),
    rChildren := self >> self._children,
    rSetChild := meth(self, n, what) self._children[n] := what; end,
    child     := (self, n) >> self._children[n],
    children  := self >> self._children,

    print      := (self,i,is) >> Print(self.name, "(", self.child(1), ", ", self.child(2), ")"),

    dims := self >> [self._children[1], self._children[1]],

));

overlapAddCyclicShifts := meth(num_points, step_size, num_steps)
    local res, i, it;

    res := [];
    for i in [0 .. num_steps-1] do
       Append(res, [[routePermOmegaInv(cyclicShiftList(num_points, step_size*i))]]);
    od;

    return omegaInvNetwork(num_points, res);

end;

# generage Verilog for Inverse Omega network for a family of cyclic shift permutation,
# given by C_n^(s*l), where l = 0 .. k.
# In code below, n is "num_points", s is "step_size", k is "num_steps".
#
# "bits" gives the number of bits per word.
# "name" gives the output filename, which will be created in /tmp/spiral/[PID] where
# [PID] is Spiral's process ID.  You can find the current PID in Spiral by typing
# getPid();
#
# If you want to change this diretory, it is easy to do
# so; just edit the file hacks.gi in spiral/paradigms/stream/hacks.gi
#
# To generate an 8 point network that will perform C_8^0, C_8^1, ..., C_8^7, and
# store it in file.v, run:
#    genOverlapAddShift(8, 1, 8, 16, "file");
genOverlapAddShift := function(num_points, step_size, num_steps, bits, name)
    local opts;
    opts := InitStreamUnrollHw();
    HDLGen(overlapAddCyclicShifts(num_points, step_size, num_steps), 1, 0, bits, 0, 0, name);
end;


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

nz := x->When(x=0,0,1);
ps := x -> PrintMat(MapMat(x,nz));
crdiag := rr -> List([1..DimensionsMat(rr)[1]/2],
    i->rr[2*i-1][2*i-1] + E(4)*rr[2*i][2*i-1]);

submat := (m, rbegin, colbegin, size) -> List(m{[rbegin..rbegin+size-1]}, x->x{[colbegin..colbegin+size-1]});

N := 8;
m := 4;
n := N/m;

BIGNT := PRDFT4;
RIGHTNT := PRDFT4;
LEFTNT := BIGNT;

# DTT -> K^N_m (dirsum SkewDTT(m)) (DST3(n) tensor I(m)) B^N_m
ck := MatSPL(K(N,m));
cm := ck^-1;
rk := RCMatCyc(ck);
rm := RCMatCyc(cm);

big := MatSPL(BIGNT(2*N));
lhs := Tensor(I(n), LEFTNT(2*m));
rhs := Tensor(RIGHTNT(n), I(2*m));
quo := MatSPL(lhs)^-1 * rm * big;  #pquo := quo ^ MatSPL(L(2*N, n)); squo := MapMat(pquo, nz);
rb := quo * MatSPL(rhs)^-1;

# N=2, m=2, n=4
#prb := (MatSPL(L(16,2)*Tensor(M(8,4),I(2)))*rb*MatSPL(L(16,8)));

#TensorProductMat(IdentityMat(2*m), MatSPL(IPRDFT2(n,-1)*Diag(2,2,1,1)));
#PRDFT3(n))^-1);
#IPRDFT2(n,-1)*Diag(2,2,1,1)));

#pquo := quo ^ MatSPL(L(2*N, 2*n));
#blocks := List([1..m], x -> pquo{[1+2*n*(x-1) .. 2*n*x]}{[1+2*n*(x-1) .. 2*n*x]});

#brfinv := MatSPL(PRDFT3(2*n))^-1;
#rdiag := List(blocks, b -> b * brfinv);
#cdiag := List(rdiag, crdiag);
#diag := MatSPL(L(2*N, m)) * quo * MatSPL(L(2*N, 2*n)) * TensorProductMat(IdentityMat(m), brfinv);


text:
MovingAverage := function(n)
  local sma, buffer, pos, sum, len;
  buffer := List([1 .. n], i -> 0);
  pos := 0;
  len := 0;
  sum := 0;
  sma := function(x)
    pos := RemInt(pos, n) + 1;
    sum := sum + x - buffer[pos];
    buffer[pos] := x;
    len := Minimum(len + 1, n);
    return sum/len;
  end;
  return sma;
end;

f := MovingAverage(3);
f(1);  #  1
f(2);  #  3/2
f(3);  #  2
f(4);  #  3
f(5);  #  4
f(5);  #  14/3
f(4);  #  14/3
f(3);  #  4
f(2);  #  3
f(1);  #  2

max_stars_repo_name: mullikine/RosettaCodeData