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

rSetChildFields := arg -> Checked(IsList(arg),
	Subst(
		meth(self, n, newChild)
			if not n in [1..numfields] then
				Error("<n> must be in ", [1..numfields]);
			else
				self.(RecName(fields[n])) := newChild;
			fi;
		end,
		numfields => Length(arg),
		fields => arg)
);

rSetChildFieldsCh := arg -> Checked(IsList(arg),
	Subst(
		meth(self, n, newChild)
			if not n in [1..chCount+numfields] then
				Error("<n> must be in ", [1..chCount+numfields]);
			elif n in [1..chCount] then
				self._children[n] := newChild;
			else
				self.(RecName(fields[n-chCount])) := newChild;
			fi;
		end,
		numfields => Length(arg)-1,
		chCount => arg[1],
		fields => Drop(arg, 1))
);

#  rChildrenFieldsCh: gives rChildren method where rChildren is self.children() :: [self.field1, self.field2, ... ],
#  arg[1] is "field1", arg[2] is "field2" etc.
#  ex: rChildren := rChildrenFieldsCh("stride", "v")
#
rChildrenFieldsCh := arg -> Checked(IsList(arg),
	Subst( (self) >> self._children :: List([1..numfields], i -> self.(RecName(fields[i]))),
		numfields => Length(arg),
		fields => arg )
);

rSetChildFieldsF := arg -> Checked(IsList(arg), Length(arg)>1,
	Subst(
		meth(self, n, newChild)
			if not n in [1..numfields] then
				Error("<n> must be in ", [1..numfields]);
			else
				self.(RecName(fields[n])) := func(newChild);
			fi;
		end,
		numfields => Length(arg)-1,
		fields => Drop(arg, 1),
		func => arg[1])
);

Local._print := function(arg)
	local a;
	for a in arg do
		if IsRec(a) and IsBound(a.__bases__) then
			Print(a.__bases__[1]);
		else
			Print(a);
		fi;
	od;
end;

Local._children := function(obj)
	if IsRec(obj) then
		if not IsBound(obj.rChildren) then
			return [];
		else
			return obj.rChildren();
		fi;
	elif BagType(obj) in [T_STRING, T_RANGE] then
		return [];
	elif IsList(obj) then
		return obj;
	else
		return [];
	fi;
end;

Local._setChild := function(obj, n, newchild)
	if IsRec(obj) then
		if not IsBound(obj.rSetChild) then
			Error("<obj> must have rSetChild(<n>, <newChild>) method");
		fi;
		obj.rSetChild(n, newchild);
	elif IsList(obj) then
		obj[n] := newchild;
	else
		Error("<obj> must be a record or a list");
	fi;
end;

Local._fromChildren := function(obj, newchildren)
	if IsRec(obj) then
		if not IsBound(obj.from_rChildren) then
			Error("<obj> must have from_rChildren(<newChildren>) method");
		fi;
		return obj.from_rChildren(newchildren);
	elif BagType(obj) in [T_STRING, T_RANGE] then
		return  obj;
	elif IsList(obj) then
		return newchildren;
	elif newchildren=[] then
		return obj;
	else
		return Error("<obj> must be a record or a list");
	fi;
end;

_PrintShape := function(obj, maxDepth, printCutoff, info, i, is)
	if maxDepth > 0 then
		Print(Blanks(i));
		_print(obj, " ", info(obj), "\n");
		DoForAll(_children(obj), c -> _PrintShape(c, maxDepth-1, printCutoff, info, i+is, is));
	elif printCutoff then
		Print(Blanks(i));
		Print("...\n");
	fi;
end;

PrintShape := (obj, maxDepth) -> Checked(IsInt(maxDepth) and maxDepth >= 0,
	_PrintShape(obj, maxDepth, true, x->"", 0, 4));

# info is a function : obj -> what to print
PrintShapeInfo := (obj, maxDepth, info) -> Checked(IsInt(maxDepth) and maxDepth >= 0,
	_PrintShape(obj, maxDepth, true, info, 0, 4));

PrintShape2 := (obj, maxDepth) -> Checked(IsInt(maxDepth) and maxDepth >= 0,
	_PrintShape(obj, maxDepth, false, x->"", 0, 4));

Local.id := ObjId; # ObjId now handles lists

#F ShapeObject(<obj>)
#F   Convert an object to a list-based pattern,
#F   i.e. add(X, Y) becomes [add, var, var]
#F
ShapeObject := obj -> Cond(
	BagType(obj) in [T_STRING,T_RANGE], obj,
	IsRec(obj) or IsList(obj), let(
		res := Concatenation([ObjId(obj)], List(_children(obj), ShapeObject)),
		When(Length(res)=1, res[1], res)),
	obj
);

#F ShapeObject@(<obj>, <func>)
#F
#F Same as ShapeObject(obj) but takes a boolean function <func>, and objects
#F for which <func> returns true are replaced by @.
#F
ShapeObject@ := (obj, func_replace_by_wildcard) -> Cond(
	BagType(obj) in [T_STRING,T_RANGE], obj,
	func_replace_by_wildcard(obj), @,
	IsRec(obj) or IsList(obj), let(
		res := Concatenation([ObjId(obj)], List(_children(obj), x -> ShapeObject@(x, func_replace_by_wildcard))),
		When(Length(res)=1, res[1], res)),
	obj
);

Class(..., rec(
	left := false,
	right := false
));

Left := z -> z.left;
Right := z -> z.right;

# @ : basic pattern object
# .parent is the pointer to the unique object in @.table
# That table does not store any precondition functions
#
Class(@, rec(
	is@ := true,

	table := WeakRef(tab()),

	__call__ := meth(arg)
		local res, self, num, target, precond;
		if Length(arg) < 2 or Length(arg) > 4 then
			Error("Usage: @(<num>, [<target>], [<cond>])");
		fi;

		self := arg[1];
		num := arg[2];
		if Length(arg) >= 3 then
			target := arg[3];
		fi;
		if Length(arg) = 4 then
			precond := arg[4];
		fi;

		if not IsBound(self.table.(num)) then
			res := self.new(num);
			self.table.(num) := res;
			if IsBound(target) then
				res := res.target(target);
			fi;
			if IsBound(precond) then
				res := res.cond(precond);
			fi;
			return res;
		else
			res := self.table.(num);
			Unbind(res._target);
			Unbind(res._precond);
			if IsBound(target) then
				res := res.target(target);
			fi;
			if IsBound(precond) then
				res := res.cond(precond);
			fi;
			return res;
		fi;
	end,

	new := (self, num) >> CantCopy(
		WithBases(self, rec(num := num, operations := PrintOps))),

	match := meth(self, obj, cx)
		if (not IsBound(self._target) or ObjId(obj) in self._target) and
		   (not IsBound(self._precond) or self._precond(obj))
		then
			if IsBound(self.parent) then
				self.parent.val := obj;
			else
				self.val := obj;
			fi;
			return true;
		else
			return false;
		fi;
	end,

	cond := (self, precond_func) >> Checked(IsCallableN(precond_func, 1), CantCopy(
		WithBases(self, rec(_precond := precond_func,
							 parent := When(IsBound(self.parent), self.parent, self))))),

	target := (self, target_id) >> CantCopy(
		WithBases(self, rec(#_taddr  := When(IsList(target_id), List(target_id,BagAddr), [BagAddr(target_id)]),
							_target := When(IsList(target_id), target_id,			   [target_id]),
							parent := When(IsBound(self.parent), self.parent, self)))),

	val := false,

	print := self >> Print("@(", self.num, ")",
						   When(IsBound(self._target),
								Print(".target(", self._target, ")"),
								""),
						   When(IsBound(self._precond),
								Print(".cond(", self._precond, ")"),
								"")),

	clear := meth(self)
		local e, l;
		l := TabToList(self.table);
		for e in l do
			if IsRec(e) and IsBound(e.val) then
				e.val := false;
			fi;
		od;
		if IsBound(self.val) then
			self.val := false;
		fi;
	end
));

# @@ : basic pattern object for context sensitive matching
# @@ has its own table.
# It still uses a .parent to point into @@.table where matches are stored
# That table does not store any precondition functions.
# Both lhs and rhs of rewrite rules must use either @ or @@. Otherwise wrong
# table will be used.

Class(@@, @, rec(
	table := WeakRef(tab()),

	match := meth(self, obj, cx)
		local base_match;
		base_match := @.match;
		if not IsBound(self._cxcond) then
			return base_match(self, obj, cx);
		else
			return base_match(self, obj, cx) and self._cxcond(obj, cx);
		fi;
	end,

	cond := (self, cxcond_func) >> CantCopy(Checked(IsCallableN(cxcond_func, 2),
		WithBases(self, rec(_cxcond := cxcond_func,
							 parent := When(IsBound(self.parent), self.parent, self))))),

	print := self >> Print("@@(", self.num, ")",
						   When(IsBound(self._target),
								Print(".target(", self._target, ")"),
								""),
						   When(IsBound(self._cxcond),
								Print(".cond(", self._cxcond, ")"),
								""))
));

Is@ := x -> IsRec(x) and IsBound(x.is@) and x.is@;

Local._midmatch := arg->false;
Local._normalmatch := arg->false;

Local._match_id := (obj, shape, cx) -> Cond(
	IsRec(shape) and IsBound(shape.is@),
		shape.match(obj, cx),
	BagType(shape) < T_FUNCTION or BagType(shape) in [T_STRING,T_RANGE],
		BagType(shape) = BagType(obj) and shape=obj,
	Same(obj, shape) or Same(ObjId(obj), shape));

##
Declare(cx_enter, cx_leave, empty_cx, apply_rules_ni, _SubstBottomUp, _SubstTopDown);
##

# PatternMatch(<obj>, <shape>, <cx>)
#	Returns true if <obj> matches the given <shape> in a given context <cx>.
#	For plain matches use empty context table empty_cx().
#
PatternMatch := function(obj, shape, cx)
	local ch, numch, shlen, res;

	if not IsList(shape) or BagType(shape) in [T_STRING,T_RANGE] then
		return _match_id(obj, shape, cx);
	elif shape = [] then
		return false;
	elif not _match_id(obj, shape[1], cx) then
		return false;
	else
		shlen := Length(shape);
		ch := _children(obj);
		numch := Length(ch);
		if shlen = 1
			then return numch = 0;
		else
			cx_enter(cx, obj);
			if shape[2] = ... then
				if shape[shlen] = ... then
					res := _midmatch(ch, shape{[3..shlen-1]}, cx);
				elif numch < shlen-2 then
					res := false;
				else
					....left := numch-shlen+2;
					....right := numch+1;
					res := _normalmatch(ch, numch-shlen+3, shlen-2,  shape, 3, shlen-2, cx);
				fi;
			elif Last(shape) = ... then
				if numch < shlen-2 then
					res := false;
				else
					....left := 0;
					....right := shlen-1;
					res := _normalmatch(ch, 1, shlen-2,  shape, 2, shlen-2, cx);
				fi;
			else
				....left := false;
				....right := false;
				res := _normalmatch(ch, 1, numch,  shape, 2, shlen-1, cx);
			fi;
			cx_leave(cx, obj);
			return res;
		fi;
	fi;
end;

_normalmatch := function(lst, lstart, llen, shape, sstart, slen, cx)
	local i, ch;
	if llen <> slen then
		return false;
	else
		for i in [0..slen-1] do
			if not PatternMatch(lst[lstart+i], shape[sstart+i], cx) then
				return false;
			fi;
		od;
		return true;
	fi;
end;

_midmatch := function(lst, shape, cx)
	local i, ch, shlen, res, llen;
	shlen := Length(shape);
	llen := Length(lst);
	if llen < shlen then
		return false;
	elif llen = shlen then
		res := _normalmatch(lst, 1, llen, shape, 1, shlen, cx);
		if res then
			....left := 0;
			....right := llen + 1;
		fi;
		return res;
	else
		for i in [1 .. llen - shlen + 1] do
			if PatternMatch(lst[i], shape[1], cx) then
				if shlen = 1 then
					....left := i-1;
					....right := i+shlen;
					return true;
				else
					if _normalmatch(lst, i+1, shlen-1,  shape, 2, shlen-1, cx) then
						....left := i-1;
						....right := i+shlen;
						return true;
					fi;
				fi;
			fi;
		od;
		return false;
	fi;
end;

#F AlternativesRewrite( <expr>, <from>, <to_func> )
#F
#F Creates all alternatives for substitution of an expression tree.
#F	<expr> - expression to substitute in
#F	<from> - shape to substitute
#F	<to_func> - a substitution function of the form e->f(e), where
#F			 as <e> will be passed the subtree matching <from>
#F

Declare(_AlternativesRewrite, _ConditionalAlternativesRewrite);

AlternativesRewrite := (expr, from, to_func) ->
	_AlternativesRewrite(expr, from, to_func, empty_cx(), []);

_AlternativesRewrite := function(expr, from, to_func, cx, list)
	local ch, i, clist;

	ch := _children(expr);
	if Length(ch) <> 0 then
		# not a leaf
		cx_enter(cx, expr);
		for i in [1..Length(ch)] do
			clist := _AlternativesRewrite(ch[i], from, to_func, cx, []);
			Append(list, List(clist, function(x) local a; a:=Copy(expr); _setChild(a, i, x); return a; end));
		od;
		cx_leave(cx, expr);
	fi;

	if PatternMatch(expr, from, cx) then
		Add(list, When(NumArgs(to_func)=2, to_func(cx, expr), to_func(expr)));
	fi;

	return list;
end;

ConditionalAlternativesRewrite := (expr, from, condition_to_func, rewrite_to_func) ->
	_ConditionalAlternativesRewrite(expr, from, condition_to_func, rewrite_to_func, empty_cx(), []);

_ConditionalAlternativesRewrite := function(expr, from, condition_to_func, rewrite_to_func, cx, list)
	local ch, i, clist;

	ch := _children(expr);
	if Length(ch) <> 0 then
		# not a leaf
		cx_enter(cx, expr);
		for i in [1..Length(ch)] do
			clist := _ConditionalAlternativesRewrite(ch[i], from, condition_to_func, rewrite_to_func, cx, []);
			Append(list, List(clist, function(x) local a; a:=Copy(expr); _setChild(a, i, x[2]); return [x[1], a]; end));
		od;
		cx_leave(cx, expr);
	fi;

	if PatternMatch(expr, from, cx) then
		Add(list, [When(NumArgs(condition_to_func)=2, condition_to_func(cx, expr), condition_to_func(expr)),
				When(NumArgs(rewrite_to_func)=2, rewrite_to_func(cx, expr), rewrite_to_func(expr))]);
	fi;

	return list;
end;

# Return a list of rewrite rule objects, from either a list of a record
# passing in a record has the advantage of clean rule naming
#
_parseRuleSet := rset -> Cond(
	IsList(rset), List(rset, Rule),
	IsRec(rset),  List(UserRecFields(rset),
	function(fld) local r; r := Rule(rset.(fld)); r.name := fld; return r; end),
	Error("<rset> must be a list of rules, or a record rec(rule1 := Rule(...), ...)")
);

#F SubstBottomUp( <expr>, <from>, <to_func> )
#F
#F Destructive bottom up substitution on an expression tree.
#F	<expr> - expression to substitute in
#F	<from> - shape to substitute
#F	<to_func> - a substitution function of the form e->f(e), where
#F			 as <e> will be passed the subtree matching <from>
#F
SubstBottomUp := (expr, from, to_func) ->
	_SubstBottomUp(expr, [ Rule(from, to_func, "unnamed(SubstBottomUp)") ], empty_cx());

SubstBottomUpRules := (expr, ruleset) ->
	_SubstBottomUp(expr, _parseRuleSet(ruleset), empty_cx());

_SubstBottomUp := function(expr, rules, cx)
	local ch, i;
	ch := _children(expr);
	if ch <> [] then
		cx_enter(cx, expr);
		for i in [1..Length(ch)] do
			_setChild(expr, i, _SubstBottomUp(ch[i], rules, cx));
		od;
		cx_leave(cx, expr);
	fi;
	return apply_rules_ni(rules, expr, cx);
end;

_SubstLeaves := function(expr, from, to_func, cx)
	local ch, i;
	ch := _children(expr);
	if Length(ch) <> 0 then
		# not a leaf
		cx_enter(cx, expr);
		for i in [1..Length(ch)] do
			_setChild(expr, i, _SubstLeaves(ch[i], from, to_func, cx));
		od;
		cx_leave(cx, expr);
		return expr;
	else
		# a leaf
		if PatternMatch(expr, from, cx) then
			return When(NumArgs(to_func)=2, to_func(cx, expr), to_func(expr));
		else
			return expr;
		fi;
	fi;
end;

SubstLeaves := (expr, from, to_func) -> _SubstLeaves(expr, from, to_func, empty_cx());

SubstChildren := function(expr, from, to_func)
	local ch, i;
	ch := _children(expr);
	#expr := map_children_safe(expr, c -> When(PatternMatch(c,from,empty_cx()), to_func(c), c));
	for i in [1..Length(ch)] do
		_setChild(expr, i,
			When(PatternMatch(ch[i],from,empty_cx()), to_func(ch[i]), ch[i]));
	od;
	return expr;
end;

SubstTopDown := (expr, from, to_func) ->
	_SubstTopDown(expr, [ Rule(from, to_func, "unnamed(SubstTopDown)") ], empty_cx());

SubstTopDown_named := (expr, from, to_func, name) ->
	_SubstTopDown(expr, [ Rule(from, to_func, name) ], empty_cx());

SubstTopDownRules := (expr, ruleset) ->
	_SubstTopDown(expr, _parseRuleSet(ruleset), empty_cx());

_SubstTopDown := function(expr, rules, cx)
	local ch, newch, res;
	expr := apply_rules_ni(rules, expr, cx);

	ch := _children(expr);
	if ch <> [] then
		cx_enter(cx, expr);
		newch :=  List(ch, c -> _SubstTopDown(c, rules, cx));
		res := _fromChildren(expr, newch);
		cx_leave(cx, expr);
	else
		res := expr;
	fi;
	return res;
end;

Declare(_SubstTopDownNR);

# same as SubstTopDown but doesn't recurse on the substitution
SubstTopDownNR := (expr, from, to_func) ->
	_SubstTopDownNR(expr, [ Rule(from, to_func, "unnamed(SubstTopDownNR)") ], empty_cx());

SubstTopDownNR_named := (expr, from, to_func, name) ->
	_SubstTopDownNR(expr, [ Rule(from, to_func, name) ], empty_cx());

SubstTopDownRulesNR := (expr, ruleset) ->
	_SubstTopDownNR(expr, _parseRuleSet(ruleset), empty_cx());

_SubstTopDownNR := function(expr, rules, cx)
	local ch, newch, res, n_applied;
	n_applied := cx.applied;
	cx.rlimit := 1;
	expr := apply_rules_ni(rules, expr, cx);
	if cx.applied > n_applied then
		return expr;
	fi;

	ch := _children(expr);
	if ch <> [] then
		cx_enter(cx, expr);
		newch :=  List(ch, c -> _SubstTopDownNR(c, rules, cx));
		res := _fromChildren(expr, newch);
		cx_leave(cx, expr);
	else
		res := expr;
	fi;
	return res;
end;

#F MatchSubst(<expr>, <rulelist>)
#F
#F Attempts matching 1 of the rules from the list to <expr> if it
#F succeeds the transformation is applied, otherwise original
#F expression is returned.
#F
#F MatchSubst is similar to SubstTopDown, but it does not recurse
#F on children.
#F
MatchSubst := function(expr, rules)
	local ch, i, r, cx;
	cx := empty_cx();
	for r in _parseRuleSet(rules) do
		if PatternMatch(expr, r.from, cx) then
			RuleTrace(r);
			return r.to(expr);
		fi;
	od;
	return expr;
end;

Declare(_Collect);

#F Collect(<expr>, <shape>)
#F
#F Returns a list of all subtrees of <expr> that match <shape>
#F
Collect := (expr, shape) -> _Collect(expr, shape, empty_cx(), true);

#F Contains(<expr>, <shape>)
#F
#F Returns 'true' if <expr> contains <shape>
#F NOTE: optimize this
Contains := (expr, shape) -> _Collect(expr, shape, empty_cx(), false)<>[];

#F CollectNR(<expr>, <shape>)
#F
#F Returns a list of all subtrees of <expr> that match <shape>.
#F Unlike 'Collect', once 'CollectNR' finds an object matching <shape>,
#F it does not search inside its children.
#F
CollectNR := (expr, shape) -> _Collect(expr, shape, empty_cx(), false);

_Collect := function(expr, shape, cx, do_recurse)
	local res, ch, c;
	if PatternMatch(expr, shape, cx) then
		if not do_recurse then
			return [expr];
		else
			cx_enter(cx, expr);
			res := Concatenation(List(_children(expr), c -> _Collect(c, shape, cx, do_recurse)));
			cx_leave(cx, expr);
			return Concatenation([expr], res);
		fi;
	else
		cx_enter(cx, expr);
		res := Concatenation(List(_children(expr), c -> _Collect(c, shape, cx, do_recurse)));
		cx_leave(cx, expr);
		return res;
	fi;
end;

Declare(_Pull);

#F Pull(<expr>, <shape>, <to_func>, <pull_func>)
#F
#F Pull is a hybrid of Collect and SubstBottomUp. First three parameters
#F specify substitution, and the last one collection:
#F
#F	<expr> - expression to substitute in
#F	<from> - shape to substitute
#F	<to_func> - a substitution function of the form e->f(e), where
#F			 as <e> will be passed the subtree matching <from>
#F	<pull_func> - applied to matching subtree <e> to collect data
#F
#F Pull retuns a tuple [<data>, <newtree>], where <data> is a list
#F of all matched subtrees (before substitutions) with <pull_func>
#F applied, and <newtree> is a tree obtained by substitution.
#F
## Note: This is Top-down recursive Pull, same as PullTD
Pull := (expr, shape, to_func, pull_func) -> _Pull(expr, shape, to_func, pull_func, empty_cx(),
	true, true);
# Top-down recursive Pull
PullTD := (expr, shape, to_func, pull_func) -> _Pull(expr, shape, to_func, pull_func, empty_cx(),
	true, true);

# Bottom-up recursive Pull (non-recursive bottom up does not exist)
PullBU := (expr, shape, to_func, pull_func) -> _Pull(expr, shape, to_func, pull_func, empty_cx(),
	true, false);

# Top-down non-recursive Pull
PullNR := (expr, shape, to_func, pull_func) -> _Pull(expr, shape, to_func, pull_func, empty_cx(),
	false, true);

_Pull := function(expr, shape, to_func, pull_func, cx, do_recurse, top_down)
	local ch, i, t, data, newdata, newexpr;
	data := [];

	if top_down then
		if PatternMatch(expr, shape, cx) then
			Add(data, When(NumArgs(pull_func)=2, pull_func(cx, expr), pull_func(expr)));
			expr := When(NumArgs(to_func)=2, to_func(cx, expr), to_func(expr));
			if not do_recurse then
				return [ data, expr ];
			fi;
		fi;
	fi;

	cx_enter(cx, expr);
	ch := _children(expr);
	for i in [1..Length(ch)] do
		t := _Pull(ch[i], shape, to_func, pull_func, cx, do_recurse, top_down);
		Append(data, t[1]);
		_setChild(expr, i, t[2]);
	od;
	cx_leave(cx, expr);

	if not top_down then
		if PatternMatch(expr, shape, cx) then
			Add(data, When(NumArgs(pull_func)=2, pull_func(cx, expr), pull_func(expr)));
			expr := When(NumArgs(to_func)=2, to_func(cx, expr), to_func(expr));
			if not do_recurse then
				return [ data, expr ];
			fi;
		fi;
	fi;
	return [data, expr];
end;

#F Harvest( <expr>, <shape-func-list> )
#F
#F Harvest is a generalization of Collect( <expr>, <shape> ).
#F It returns concatenated list of <func>(child) for all children
#F of <expr> that match <shape>, for each <shape>-<func> pair.
Harvest := function( expr, shape_func_list )
	local res, sf, c;
	res := [];
	for sf in shape_func_list do
		Append(res, List(Collect(expr, sf[1]), sf[2]));
	od;
	return res;
end;

#F SubstObj(<expr>, <obj>, <new_obj>
#F replacing all occurences of <obj> by <new_obj>
SubstObj := (expr, obj, new_obj) -> SubstTopDownNR_named(expr, @.cond(x -> x=obj), e -> new_obj, "SubstObj");


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

Class(swp_loop, loop);

Class(dma_signal, ExpCommand);
Class(dma_wait, ExpCommand);

Class(cpu_signal, ExpCommand);
Class(cpu_wait, ExpCommand);

Class(barrier_cmd, ExpCommand);
Class(nop_cmd, ExpCommand);

Class(dma_fence, skip);

Class(dma_transfer, Command, rec(
    isAssign := false, # transfers are not assignments, since they take >2 params
    __call__ := (self, loc, exp, size) >> WithBases(self, rec(
        operations := CmdOps,
        loc := toAssignTarget(loc),
                exp := toExpArg(exp),
        size := toExpArg(size)
   )),

   rChildren := self >> [self.loc, self.exp, self.size],
   rSetChild := rSetChildFields("loc", "exp", "size"),
   unroll := self >> self,

   print := (self,i,is) >> let(name := Cond(IsBound(self.isCompute) and self.isCompute,
                                            gap.colors.DarkYellow(self.__name__),
                                            IsBound(self.isLoad) and self.isLoad,
                                            gap.colors.DarkRed(self.__name__),
                                            IsBound(self.isStore) and self.isStore,
                                            gap.colors.DarkGreen(self.__name__),
                                            self.__name__),
                                Print(name, "(", self.loc, ", ", self.exp, ", ", self.size, ")"))
));

Class(dma_load, dma_transfer);
Class(dma_store, dma_transfer);

# wrapper for dma_xfer, includes the dma transfer size.
# dma_size( dma_xfer(a,b), c);
Class(dma_size, ExpCommand);
Class(par_exec, chain);

#The next two dma_jump and cpu_jump are used only in the CM stuff in ScratchX86
Class(dma_jump,ExpCommand);
Class(cpu_jump,ExpCommand);

#For thread registration using the Fast Barrier
Class(register,ExpCommand);
Class(initialization, ExpCommand);
Class(add_buffer,skip);

ext: gi


text:
InfiniteLoop := function()
	local n;
	n := 1;
	while true do
		Display(n);
		n := n + 1;
	od;
end;

# Prepare some coffee
InfiniteLoop();

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


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

_unTspl := (c, info) -> When(spiral.paradigms.common.IsTSPL(c), SumsSPL(c.toSpl(), info.opts), c);

SumsRuleTreeStep := function ( rtwrap, info )
  local S, rt, t, tag, tags, processed_tags;
  rt := rtwrap.rt;
  info.rsteps := info.rsteps + 1;

  tags := When(IsBound(rt.node.getTags), rt.node.getTags(), []);

  if IsBound(rt.node.noCodelet) then
      S := SumsSPL(_SPLRuleTree(rt), info.opts);

  else
      processed_tags := Union(info.processed_tags, List(tags, x->x.kind()));
      S := ApplyRuleTreeStep(rt, c -> Cond(
          not IsRuleTree(c) and IsBound(c.noCodelet) and c.noCodelet,
	      _unTspl(c, info),
	  not IsRuleTree(c),
	      RecursStep(_unTspl(c, info)).setA(processed_tags => processed_tags),
          IsBound(c.node.noCodelet) and c.node.noCodelet,
	      _SPLRuleTree(c),
	  # else
              RecursStep(RTWrap(c)).setA(processed_tags => processed_tags)
      ));
      S := SumsSPL(S, info.opts);
      S.root := rt.node;
  fi;

  # tags may inject container objects
  for tag in Reversed(tags) do
      if IsBound(tag.container) and not (tag.kind() in info.processed_tags) then
          S := tag.container(S);
      fi;
  od;
  return S;
end;

SumsRecursStep := function ( rstep, info )
  local rsteps;
  rsteps := info.rsteps;
  info.processed_tags := When(IsBound(rstep.a.processed_tags), rstep.a.processed_tags, []);
  rstep := SubstTopDownNR(rstep, RTWrap, e -> SumsRuleTreeStep(e, info));
  if rsteps = info.rsteps then # nothing has changed
      return rstep;
  else
      rstep := ApplyStrategy(rstep, info.strategy, UntilDone, info.opts);
      return rstep.child(1); # strip outer RecursStep container
  fi;
end;

Recurse := function(sums,info)
    local rt;
    if IsRuleTree(sums) then
	rt := sums;
	# keep track of tags that were converted into containers, via an attribute
	# this mechanism does not work well, and more VContainers will be created
	# than neeeded, however this is safe, because redundant ones are eliminated
	# by rewrite rules.
        return Recurse(RecursStep(RTWrap(rt)).setA(processed_tags=>[]), info);

    else
        return SubstTopDownNR(sums, RecursStep, x -> SumsRecursStep(x,info));
    fi;
end;

SumsRuleTreeStrategy := function ( rt, strategy, opts )
    local info;
    info := rec(rsteps := 0, strategy := strategy, opts := opts); #, cutoff := cutoff_func);
    rt := Recurse(rt, info);
    while info.rsteps > 0 do info.rsteps := 0; rt := Recurse(rt, info); od;
    return rt;
end;

#F SumsRuleTree(<rt>, <opts>)
#F
#F <opts> flags used:
#F   opts.formulaStrategies.sigmaSpl    Sigma-SPL rewriting strategy
#F   opts.formulaStrategies.rc          RC(.) rewriting strategy
#F   opts.generateComplexCode == bool   if set to false, then RC rewriting strategy is applied
#F
SumsRuleTree := function(rt, opts)
        local sums, rsums, t, tag;
        if IsNonTerminal(rt) then rt := RandomRuleTree(rt,opts); fi;
        sums := SumsRuleTreeStrategy(rt, [], opts);
        sums := ApplyStrategy(sums, opts.formulaStrategies.sigmaSpl, UntilDone, opts);
        sums := ApplyStrategy(sums, opts.formulaStrategies.preRC, UntilDone, opts);
        if (not opts.generateComplexCode) and (not rt.node.isReal()) then # NOTE: when t_in/t_out available check them instead of isReal
            rsums := ApplyStrategy(RC(sums), opts.formulaStrategies.rc, UntilDone, opts);
        else
            rsums := sums;
        fi;
        rsums := ApplyStrategy(rsums, opts.formulaStrategies.postProcess, UntilDone, opts);
        rsums := SumsUnification(rsums, opts);
        rsums.ruletree := rt;

        return rsums;
end;

#F See SumsRuleTree
#F
SumsRuleTreeOpts := SumsRuleTree;

RecurseOpts := function(rt, opts)
    local sums, rsums;
    if IsNonTerminal(rt) then rt := RandomRuleTree(rt,opts); fi;
    sums := Recurse(rt, rec(strategy:=opts.formulaStrategies.sigmaSpl, rsteps:=0, opts := opts));
    return sums;
    #return ApplyStrategy(sums, opts.formulaStrategies.postProcess, UntilDone, opts);
end;

######
SumsVerifyRulesForSPL := (S,opts) -> VerifyRules(S,
    (rt, s) -> InfinityNormMat(MatSPL(SumsRuleTree(rt,opts)) - MatSPL(s)) < 1e-11,
    opts);


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

ZeroCoef := function(p)
local i, L;
L:=p.coefficients;
for i in [1..Length(p.coefficients)] do
  if AbsFloat(ReComplex(ComplexAny(p.coefficients[i]))) < FloatString("1e-4") then
    L[i] := p.baseRing.zero;
  fi;
od;
return Polynomial(p.baseRing, L, p.valuation);
end;

# Recursive implementation of the long division for Laurent polynomials h(z) and g(z)
# There are deg(h)-deg(g)+2 different schemes for division
LaurentRemainder := function ( arg )

local a, b, q, q1, r1, q2, r2, Ls, Lt, z,
      deghigha, deglowa, deghighb, deglowb, dega, degb, coefa, coefb,  divs;

if IsList(arg[1]) then arg:=arg[1];fi;
a:=arg[1];
b:=arg[2];

#Print("check");
if not (IsPolynomial(a) and IsPolynomial(a)) then
  Error("<a> and <b> must be polynomials");
elif not a.baseRing.name = b.baseRing.name then
  Error("Both polynomials have to be over the same field");
elif not IsBound(PolynomialRing(a.baseRing).isEuclideanRing) then
  Error("Polymonials must be constructed over a field, not over a ring");
fi;

z:=Indeterminate(a.baseRing);
z.name:="z";
if Length(arg)=3 then q:=arg[3]; else q:=0*z^0; fi;

deglowa := a.valuation;
deglowb := b.valuation;
dega := Length(a.coefficients)-1;
degb := Length(b.coefficients)-1;
deghigha := deglowa + dega;
deghighb := deglowb + degb;

coefa := a.coefficients;
coefb := b.coefficients;

if dega<degb then return [[q,a]];fi;

  q1 := Polynomial(a.baseRing,[coefa[dega+1]/coefb[degb+1]],deghigha-deghighb);
  r1 := a - q1 * b;
  q1 := q1 + q;
#  Print(q1);
#  Print("\n");
#  Print(r1);
#Print("\n");

#  q2 := coefa[1]/coefb[1] * z^(deglowa-deglowb);
  q2 := Polynomial(a.baseRing,[coefa[1]/coefb[1]],deglowa-deglowb);
  r2 := a - q2 * b;
  q2 := q2 + q;
  Ls := [];
  for divs in [[r1,b,q1],[r2,b,q2]] do
     Lt := LaurentRemainder(divs);
     Add(Ls,Lt);
  od;
  Ls := Flat(Ls);
  Ls := List([1..Length(Ls)/2], i->[Ls[2*i-1],Ls[2*i]]);
  Ls := List(Collected(Ls), col->col[1]);
  if not IsList(Ls[1]) then return [Ls];
  else return Ls;
  fi;

end;

#
# Note:
# ! Because we now use double floats, the LiftingScheme generates more LS than
# possible because of finite precision. There will be a bunch of the same lifting
# schemes generated multiple times with coefficients different after 10th
# decimal or so. This requies fixing !
#

#F LiftingScheme( <poly1>, <poly2> )
#F   returns a list of Laurent polynomials that represent
#F   primal and dual lifting steps factorization of the polyphase
#F   matrix of polynomials <poly1> and <poly2>:
#F            P(z)= [[ <poly1_even>, <poly1_odd>],
#F                   [ <poly2_even>, <poly2_odd>]]
#F   In Discrete Wavelet Transform notation, <poly1> represents
#F   a low-pass filter impulse response, whereas <poly2> represents
#F   a high-pass one.
#F
#F   The output is in the following format:
#F   [ b, l1(z), l2(z), .... , ln(z), c1, c2]
#F
#F   b - a binary symbol that determines the type of the
#F       initial lifting step (0 - primal, 1 - dual)
#F   li(z) - polynomial representing i-th lifting step
#F   c1, c2 - scaling constants
#F
#F   Example:
#F
#F   P(z)=  [ c1, 0  ] * [ 1, l2(z)] * [ 1    , 0 ]
#F          [  0, c2 ]   [ 0, 1    ]   [ l1(z), 1 ]
#F          ----------   -----------   ------------
#F           scaling      primal l.s.   dual l.s. (b=1)
#F
#F   Note: Factorization is implemented in Laurent Polynomial
#F   ring over the base field inherited from <poly1>, <poly2>.
#F

LiftingScheme := function ( M )

local he, ho, ge, go, hen, hon, gen, gon, b, QR, z, Le, Lo, L, Lc, Ls, ve, vo, P, factor, scheme;

#Print("Recursion \n");
# Check if h(z) and g(z) in fact represent a perfect reconstruction
# filter bank -> determinant of the polyphase matrix should be
# a monomial
he := M[1][1];ho := M[1][2]; ge := M[2][1]; go := M[2][2];
  if not Length(ZeroCoef(Determinant([[he, ho],[ge, go]])).coefficients) = 1 then
     Error("h(z) and g(z) must satisfy the perfect reconstruction condition P(z)*P'(z)^(-1) = I where P and P' are analysis and synthesis polyphase matrices");
  fi;

  z:=Indeterminate(he.baseRing);
  z.name:="z";
# Factoring Algorithm
  b:=0;

  if LaurentDegree(he) >= LaurentDegree(ho) then b:=1; fi;

  if (ho=0*z^0 or ge=0*z^0) then
    L:=[];
    if ho=0*z^0 then Add(L,ZeroCoef(ge)*
         Polynomial(he.baseRing,[1/ ZeroCoef(go).coefficients[1]], -go.valuation));
    else Add(L,ho/he); fi;
    Add(L, he);
    Add(L, go);
    Add(L, b);
    return [Reversed(L)];
  fi;
  if (he=0*z^0 or go=0*z^0) then
    L:=[];
    if he=0*z^0 then Add(L,ZeroCoef(go)*
         Polynomial(he.baseRing,[1/ ZeroCoef(ge).coefficients[1]], -ge.valuation));
    else Add(L,he/ho); fi;
    Add(L, ho);
    Add(L, ge);
    Add(L, b-2);
    return [Reversed(L)];
  fi;
  L:=[];
    if b=0 then
      QR:=LaurentRemainder(ho,he);
      for factor in QR do
        hon := factor[2];
        hon := ZeroCoef(hon);
        gon := go - ge * factor[1];
        gon := ZeroCoef(gon);

        Lc:=LiftingScheme([[he,hon],[ge,gon]]);
        for scheme in Lc do Add(scheme,factor[1]); od;
        Append(L,Lc);
      od;

    else
      QR:=LaurentRemainder(he,ho);
      for factor in QR do
        hen := factor[2];
        hen := ZeroCoef(hen);
        gen := ge - go * factor[1];
        gen := ZeroCoef(gen);
        Lc:=LiftingScheme([[hen,ho],[gen,go]]);
        for scheme in Lc do Add(scheme,factor[1]); od;
        Append(L,Lc);
      od;

    fi;

return(L);

end;

Lifting := function(M)
local Lt, Ls, L, factor, p, normM, normDiff;

L:=LiftingScheme(M,[]);
Lt:=[];
Ls:=[];

for factor in L do
  Add(Ls, factor);
 if IsInt(factor) then  Add(Lt, Ls); Ls:=[];fi;
od;
return Lt;
end;

VerifyLiftingScheme:=function(P,L2)

local P, M, b, i, l,z,p, normM, normDiff;

  z:=Indeterminate(P[1].baseRing);
  z.name:="z";

#P:=[[P[1],P[2]],[P[3],P[4]]];
M:= [[1, 0],[0,1]];

L2:=Reversed(L2);
l:=Length(L2);
b:=L2[l];
for i in [1 .. l-3] do
 if b=-2 then
   M:= M * [[0*z^0,1*z^0 ], [1*z^0, L2[l-i-2]]]; b:=1;
 elif b=-1 then
   M:= M * [[L2[l-i-2],1*z^0 ], [1*z^0, 0*z^0]]; b:=0;
 elif b=0 then
   M:= M * [[1*z^0, L2[l-i-2]], [0*z^0, 1*z^0]]; b:=1;
 else
   M:= M * [[1*z^0, 0*z^0], [L2[l-i-2], 1*z^0]]; b:=0;
 fi;
od;
M:=[[L2[l-2], 0*z^0],[0*z^0, L2[l-1]]]*M;
#Print([[P[1],P[2]],[P[3],P[4]]]);
if Same(P[1][1].baseRing, Doubles) then
  normM:=Sum(M, l->Sum(l, p->(Sum(p.coefficients,i->i^2))^1/2))^1/2;
  normDiff:=Sum(M-P, l->Sum(l, p->(Sum(p.coefficients,i->i^2))^1/2))^1/2;
  Print("error: ",normDiff/normM,"\n");
  if (normDiff/normM > FloatString("1e-4")) then return false;
  else return true;
  fi;
else return[M,P,M=P];
fi;
end;

ArithmeticCostLS := function(L)
local step, mult, add, i;

mult:=0;
add:=0;
L:=Reversed(L);
for step in [1..Length(L)-3] do
  for i in L[step].coefficients do
     if not i = 0 then mult:=mult+1; add:=add+1;fi;
  od;
od;
mult:=mult+2;
return [add,mult,Length(L)-3];
end;

#F MatPoly := function(h,g)
#F   returns a matrix formed of two filters h(z) and g(z)
#F   also returns largest and the smallest degree of the
#F   coefficients in the matrix.
#F   Format: [<mat>, q, p]

MatPoly := function(h,g)

local hc, gc, hsmall, hlarge,
      gsmall, glarge, q, p, M;

    # Check the degrees of the filters
    hc := h.coefficients;
    gc := g.coefficients;
    hsmall := h.valuation;
    hlarge := Length(hc) + hsmall - 1;
    gsmall := g.valuation;
    glarge := Length(gc) + gsmall - 1;
    q := Minimum([hsmall,gsmall]);
    p := Maximum([hlarge,glarge]);

# Construct the 2 x l matrix of filter coefficients

    M :=
      [
       Concat(List([1..AbsInt(q-hsmall)],i->h.baseRing.zero),
               hc, List([1..p-hlarge],i->h.baseRing.zero)
             ),
       Concat(List([1..AbsInt(q-gsmall)],i->h.baseRing.zero),
               gc, List([1..p-glarge],i->h.baseRing.zero)
             )
      ];

    return ([M, hsmall, hlarge, gsmall, glarge]);
end;

Polyphase := function(p1,p2)

local  he, ho, ge, go, z, h, g;

if not (IsPolynomial(p1) and IsPolynomial(p1)) then
  Error("<p1> and <p2> must be polynomials");
elif not p1.baseRing.name = p2.baseRing.name then
  Error("Both polynomials have to be over the same field");
elif not IsBound(PolynomialRing(p1.baseRing).isEuclideanRing) then
  Error("Polymonials must be constructed over a field, not over a ring");
fi;

  z:=Indeterminate(p1.baseRing);
  z.name:="z";

        h := DownsampleTwo(ListPoly(p1));
        g := DownsampleTwo(ListPoly(p2));
        he := Polynomial(p1.baseRing, h[1][1], h[1][2]);
        ho := Polynomial(p1.baseRing, h[2][1], h[2][2]);
        ge := Polynomial(p1.baseRing, g[1][1], g[1][2]);
        go := Polynomial(p1.baseRing, g[2][1], g[2][2]);

# Find he(z) and ho(z)
#  Lo:=Sublist(p1.coefficients, Filtered([1..Length(p1.coefficients)],IsEvenInt));
#  Le:=Sublist(p1.coefficients, Filtered([1..Length(p1.coefficients)],IsOddInt));
#  if (p1.valuation mod 2 = 1) then L:=Le; Le:=Lo; Lo:=L;
#    ve := Int((p1.valuation + 1)/2); vo := ve;
#  else ve := Int(p1.valuation/2); vo := ve + 1;
#  fi;
#  he:=Polynomial(p1.baseRing, Le, ve);
#  ho:=Polynomial(p1.baseRing, Lo, vo);
#
# Find ge(z) and go(z)

#  Lo:=Sublist(p2.coefficients, Filtered([1..Length(p2.coefficients)],IsEvenInt));
#  Le:=Sublist(p2.coefficients, Filtered([1..Length(p2.coefficients)],IsOddInt));
#  if (p2.valuation mod 2 = 1) then L:=Le; Le:=Lo; Lo:=L;
#    ve := Int((p2.valuation + 1)/2); vo := ve;
#  else ve := Int(p2.valuation/2); vo := ve+1;
#  fi;
#  ge:=Polynomial(p2.baseRing, Le, ve);
#  go:=Polynomial(p2.baseRing, Lo, vo);

return [[he,ho],[ge,go]];
end;

# PolyPolyphaseMat:=function(M)
# Returns a list of 2 polynomials from a polyphase matrix

PolyPolyphaseMat:=function(M)

local he,ho,ge,go,h,g,z;

he:=M[1][1];
ho:=M[1][2];
ge:=M[2][1];
go:=M[2][2];
z:=Indeterminate(he.baseRing);
z.name:="z";

h:=Upsample(he)+z^(1)*Upsample(ho);
g:=Upsample(ge)+z^(1)*Upsample(go);

return([h,g]);
end;


text:
method vfClipTools vfClipTools.mLoad <alias=vfClipTools_mLoad>( )
{
//	this->vForm.mCreateWin()
	ustr empus
	str  emps
	ustr ustmp
	uint comp
	comp as this
		comp.AutoLang=1
		comp.Border=$fbrdSizeToolWin
		comp.Bottom=0
		comp.Caption=empus
		comp.Enabled=1
		comp.FormStyle=$fsPopup
		comp.Height=575
		comp.HelpTopic=empus
		comp.Hint=empus
		comp.HorzAlign=$alhLeft
		comp.IconName=empus
		comp.Left=0
		comp.Name="fClipTools"
		comp.Right=0
		comp.StartPos=$spDesigned
		comp.Style=emps
		comp.TabOrder=0
		comp.Tag=0
		comp.Top=0
		comp.TopMost=1
		comp.VertAlign=$alvTop
		comp.Visible=1
		comp.Width=398
		comp.WindowState=$wsNormal
		comp.OnMouse.Set( this, fClipTools_fm )
		comp.OnCreate.Set( this, fClipTools_Create )
		comp.OnDestroy.Set( this, fClipTools_Destroy )
		comp.OnCloseQuery.Set( this, fClipTools_Close )
		comp as this.Tab0
		comp.Owner = this
			comp.AutoLang=1
			comp.Bottom=0
			comp.Enabled=1
			comp.FixedWidth=0
			comp.Height=530
			comp.HelpTopic=empus
			comp.Hint=empus
			comp.HorzAlign=$alhClient
			comp.ImageList=ustmp.fromutf8("groups")
			comp.Left=0
			comp.Name="Tab0"
			comp.Right=0
			comp.Style=emps
			comp.TabOrder=0
			comp.TabStyle=$tsNone
			comp.Tag=0
			comp.Top=0
			comp.VertAlign=$alvTopBottom
			comp.Visible=1
			comp.Width=380
			comp as this.tiHistory
			comp.Owner = this.Tab0
				comp.AutoLang=1
				comp.Bottom=0
				comp.Caption=ustmp.fromutf8("history")
				comp.Enabled=1
				comp.Height=501
				comp.HelpTopic=empus
				comp.Hint=empus
				comp.HorzAlign=$alhLeft
				comp.ImageId=ustmp.fromutf8("history")
				comp.Left=4
				comp.Name="tiHistory"
				comp.Right=0
				comp.Style=emps
				comp.TabOrder=0
				comp.Tag=0
				comp.Top=25
				comp.VertAlign=$alvTop
				comp.Visible=1
				comp.Width=372
		comp as this.Tab0
		comp.Owner = this
			comp.CurIndex=0
		comp as this.ToolBar0
		comp.Owner = this
			comp.AutoLang=1
			comp.AutoSize=0
			comp.Bottom=5
			comp.Enabled=1
			comp.Height=25
			comp.HelpTopic=empus
			comp.Hint=empus
			comp.HorzAlign=$alhLeft
			comp.ImageList=ustmp.fromutf8("main")
			comp.Left=0
			comp.Name="ToolBar0"
			comp.Right=0
			comp.ShowCaption=$tscNone
			comp.ShowDivider=0
			comp.Style=emps
			comp.TabOrder=1
			comp.Tag=0
			comp.Top=500
			comp.VertAlign=$alvBottom
			comp.Vertical=0
			comp.Visible=0
			comp.Width=75
			comp.Wrapable=0
			comp as this.tbiMainWin
			comp.Owner = this.ToolBar0
				comp.AutoLang=1
				comp.Caption=ustmp.fromutf8("mainwin")
				comp.Checked=0
				comp.Enabled=1
				comp.Hint=empus
				comp.ImageId=ustmp.fromutf8("mclip")
				comp.Index=0
				comp.Name="tbiMainWin"
				comp.ShowCaption=0
				comp.Tag=0
				comp.TBIStyle=$tbsAsCheckBox
				comp.Visible=1
				comp.OnClick.Set( this, fClipTools_MainWin )
			comp as this.ToolBarItem0
			comp.Owner = this.ToolBar0
				comp.AutoLang=1
				comp.Caption=ustmp.fromutf8("ToolBarItem0")
				comp.Checked=0
				comp.Enabled=1
				comp.Hint=empus
				comp.ImageId=empus
				comp.Index=1
				comp.Name="ToolBarItem0"
				comp.ShowCaption=0
				comp.Tag=0
				comp.TBIStyle=$tbsButton
				comp.Visible=1
		comp as this.Tray0
		comp.Owner = this
			comp.AutoLang=1
			comp.Caption=ustmp.fromutf8("progname")
			comp.Image=ustmp.fromutf8("main\\mclip")
			comp.Name="Tray0"
			comp.RBtnPopupMenu=.pmTray
			comp.Tag=0
			comp.Visible=0
			comp.OnMouse.Set( this, fClipTools_Tray )
		comp as this.pmTray
		comp.Owner = this
			comp.AutoLang=1
			comp.Name="pmTray"
			comp.Tag=0
			comp as this.miExit
			comp.Owner = this.pmTray
				comp.AutoCheck=0
				comp.AutoLang=1
				comp.Caption=ustmp.fromutf8("exit")
				comp.Checked=0
				comp.Ellipsis=0
				comp.Enabled=1
				comp.Image=empus
				comp.Name="miExit"
				comp.RadioCheck=0
				comp.Separator=0
				comp.ShortKey=empus
				comp.Tag=0
				comp.Visible=1
				comp.OnClick.Set( this, fClipTools_Exit )
		comp as this.timerRestore
		comp.Owner = this
			comp.AutoLang=1
			comp.Enabled=0
			comp.Interval=1000
			comp.Name="timerRestore"
			comp.Tag=0
			comp.OnTimer.Set( this, fClipTools_Restore )
		comp as this.pmBtns
		comp.Owner = this
			comp.AutoLang=1
			comp.Name="pmBtns"
			comp.Tag=0
			comp.OnBeforeShow.Set( this, fClipTools_PopupBtns )
			comp.OnAfterShow.Set( this, fClipTools_PopupClose )
			comp as this.miCopyToNotices
			comp.Owner = this.pmBtns
				comp.AutoCheck=0
				comp.AutoLang=1
				comp.Caption=ustmp.fromutf8("copytonotices")
				comp.Checked=0
				comp.Ellipsis=0
				comp.Enabled=1
				comp.Image=empus
				comp.Name="miCopyToNotices"
				comp.RadioCheck=0
				comp.Separator=0
				comp.ShortKey=empus
				comp.Tag=0
				comp.Visible=0
			comp as this.miDelete
			comp.Owner = this.pmBtns
				comp.AutoCheck=0
				comp.AutoLang=1
				comp.Caption=ustmp.fromutf8("delete")
				comp.Checked=0
				comp.Ellipsis=0
				comp.Enabled=1
				comp.Image=empus
				comp.Name="miDelete"
				comp.RadioCheck=0
				comp.Separator=0
				comp.ShortKey=empus
				comp.Tag=0
				comp.Visible=1
				comp.OnClick.Set( this, fClipTools_DeleteItem )
			comp as this.MenuItem0
			comp.Owner = this.pmBtns
				comp.AutoCheck=0
				comp.AutoLang=1
				comp.Caption=ustmp.fromutf8("edit")
				comp.Checked=0
				comp.Ellipsis=0
				comp.Enabled=1
				comp.Image=empus
				comp.Name="MenuItem0"
				comp.RadioCheck=0
				comp.Separator=0
				comp.ShortKey=empus
				comp.Tag=0
				comp.Visible=0
				comp.OnClick.Set( this, fClipTools_EditItem )
	comp as this
		comp.ClientHeight=530
		comp.ClientWidth=380

	return this
}

method vfClipTools vfClipTools.init( )
{
   this.pTypeId = vfClipTools
   return this
}
func init_vfClipTools <entry>()
{
   regcomp( vfClipTools, "vfClipTools", vForm, $vForm_last,
      %{ %{$mLoad,     vfClipTools_mLoad}},
      0->collection )

}

max_stars_repo_name: novostrim/macroclip

max_issues_repo_name: novostrim/macroclip

max_forks_repo_name: novostrim/macroclip


text:
#############################################################################
##
#W    matgrp3.gi                                                Karel Dekimpe
#W                                                               Bettina Eick
##
##    This file contains the 3-dimensional almost crystallographic groups
##    as integral matrix groups. There are 17 types of groups.
##

ACDim3Nr01 := function ( k1 )
local a, b, c;
a:=[[1, 0, -k1/2, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]];
b:=[[1, k1/2, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]];
c:=[[1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
return Group( [a , b, c] , IdentityMat(4) );
end;

ACDim3Nr02 := function ( k1, k2, k3, k4)
local a, b, c, alfa;
a:=[[1, 0, -k1/2, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]];
b:=[[1, k1/2, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]];
c:=[[1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
alfa:=[[1, k2, k3, k4/2], [0, -1, 0, 0], [0, 0, -1, 0], [0, 0, 0, 1]];
return Group( [a , b, c, alfa] , IdentityMat(4) );
end;

ACDim3Nr03 := function ( k1, k2 )
local a, b, c, alfa;
a:=[[1, 0, -k1/2, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]];
b:=[[1, k1/2, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]];
c:=[[1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
alfa:=[[-1, -k2, 0, 0], [0, 1, 0, 0], [0, 0, -1, 0], [0, 0, 0, 1]];
return Group( [a , b, c, alfa] , IdentityMat(4) );
end;

ACDim3Nr04 := function ( k1, k2 )
local a, b, c, alfa;
a:=[[1, 0, -k1, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]];
b:=[[1, k1, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]];
c:=[[1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
alfa:=[[-1, 2*k2, k1/2, 0], [0, 1, 0, 1/2], [0, 0, -1, 0], [0, 0, 0, 1]];
return Group( [a , b, c, alfa] , IdentityMat(4) );
end;

ACDim3Nr05 := function ( k1, k2 )
local a, b, c, alfa;
a:=[[1, 0, -k1/2, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]];
b:=[[1, k1/2, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]];
c:=[[1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
alfa:=[[-1, -k2, -k2, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]];
return Group( [a , b, c, alfa ] , IdentityMat(4) );
end;

ACDim3Nr06 := function ( k1, k2, k3, k4)
local a, b, c, alfa, beta;
a:=[[1, 0, -k1/2, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]];
b:=[[1, k1/2, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]];
c:=[[1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
alfa:=[[1, k2, k3, k4/2], [0, -1, 0, 0], [0, 0, -1, 0], [0, 0, 0, 1]];
beta:=[[-1, -k2, 0, 0], [0, 1, 0, 0], [0, 0, -1, 0], [0, 0, 0, 1]];
return Group( [a , b, c, alfa, beta] , IdentityMat(4) );
end;

ACDim3Nr07 := function ( k1, k2, k3, k4)
local a, b, c, alfa, beta;
a:=[[1, 0, -k1, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]];
b:=[[1, k1, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]];
c:=[[1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
alfa:=[[1, k2, -2*(k3 + k4), k3/2], [0, -1, 0, 0], [0, 0, -1, 0], [0, 0, 0, 1]];
beta:=[[-1, k1/2 - k2, 0, 0], [0, 1, 0, 0], [0, 0, -1, 1/2], [0, 0, 0, 1]];
return Group( [a , b, c, alfa, beta] , IdentityMat(4) );
end;

ACDim3Nr08 := function ( k1, k2, k3, k4)
local a, b, c, alfa, beta;
a:=[[1, 0, -k1, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]];
b:=[[1, k1, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]];
c:=[[1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
alfa:=[[1, k1 + 2*k3, -k1 + 2*k2 - 2*k3, k4/2], [0, -1, 0, 0], [0, 0, -1, 0],
 [0, 0, 0, 1]];
beta:=[[-1, -k1/2 - 2*k3, k1/2, 0], [0, 1, 0, 1/2], [0, 0, -1, 1/2], [0, 0, 0, 1]];
return Group( [a , b, c, alfa, beta] , IdentityMat(4) );
end;

ACDim3Nr09 := function ( k1, k2, k3, k4)
local a, b, c, alfa, beta;
a:=[[1, 0, -k1/2, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]];
b:=[[1, k1/2, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]];
c:=[[1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
alfa:=[[1, k2, -k2 + 2*k3, k4/2], [0, -1, 0, 0], [0, 0, -1, 0], [0, 0, 0, 1]];
beta:=[[-1, -k3, -k3, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]];
return Group( [a , b, c, alfa, beta] , IdentityMat(4) );
end;

ACDim3Nr10 := function ( k1, k2, k3, k4)
local a, b, c, alfa;
a:=[[1, 0, -k1/2, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]];
b:=[[1, k1/2, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]];
c:=[[1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
alfa:=[[1, k2, k3, k4/4], [0, 0, -1, 0], [0, 1, 0, 0], [0, 0, 0, 1]];
return Group( [a , b, c, alfa] , IdentityMat(4) );
end;

ACDim3Nr11 := function ( k1, k2, k3, k4)
local a, b, c, alfa, beta;
a:=[[1, 0, -k1/2, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]];
b:=[[1, k1/2, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]];
c:=[[1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
alfa:=[[1, k2, k3, k4/4], [0, 0, -1, 0], [0, 1, 0, 0], [0, 0, 0, 1]];
beta:=[[-1, -k2 - k3, 0,  0], [0, 1, 0, 0], [0, 0, -1, 0], [0, 0, 0, 1]];
return Group( [a , b, c, alfa, beta] , IdentityMat(4) );
end;

ACDim3Nr12 := function ( k1, k2, k3, k4)
local a, b, c, alfa, beta;
a:=[[1, 0, -k1, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]];
b:=[[1, k1, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]];
c:=[[1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
alfa:=[[1, k2, k1 - k2 - 2*k3, k4/4], [0, 0, -1, 0], [0, 1, 0, 0], [0, 0, 0, 1]];
beta:=[[-1, -k1/2 + 2*k3, k1/2, 0], [0, 1, 0, 1/2], [0, 0, -1, 1/2], [0, 0, 0, 1]];
return Group( [a , b, c, alfa, beta] , IdentityMat(4) );
end;

ACDim3Nr13 := function ( k1, k2, k3, k4)
local a, b, c, alfa;
a:=[[1, 0, -k1/2, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]];
b:=[[1, k1/2, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]];
c:=[[1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
alfa:=[[1, k2, -k1/2 + k3, k4/3], [0, 0, -1, 0], [0, 1, -1, 0], [0, 0, 0, 1]];
return Group( [a , b, c, alfa] , IdentityMat(4) );
end;

ACDim3Nr14 := function ( k1, k2, k3, k4)
local a, b, c, alfa, beta;
a:=[[1, 0, -k1/2, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]];
b:=[[1, k1/2, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]];
c:=[[1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
alfa:=[[1, k2, -k1/2 + k3, k4/3], [0, 0, -1, 0], [0, 1, -1, 0], [0, 0, 0, 1]];
beta:=[[-1, -k2 , k2 , 0], [0, 0, -1, 0], [0, -1, 0, 0], [0, 0, 0, 1]];
return Group( [a , b, c, alfa, beta] , IdentityMat(4) );
end;

ACDim3Nr15 := function ( k1, k2, k3, k4)
local a, b, c, alfa, beta;
a:=[[1, 0, -k1/2, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]];
b:=[[1, k1/2, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]];
c:=[[1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
alfa:=[[1, k1 - 2*k3 + 3*k2, -k1/2 + k3, k4/3], [0, 0, -1, 0], [0, 1, -1, 0],
 [0, 0, 0, 1]];
beta:=[[-1, -k2, -k2,  0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]];
return Group( [a , b, c, alfa, beta] , IdentityMat(4) );
end;

ACDim3Nr16 := function ( k1, k2, k3, k4)
local a, b, c, alfa;
a:=[[1, 0, -k1/2, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]];
b:=[[1, k1/2, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]];
c:=[[1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
alfa:=[[1, -k1/2 + k2, k3, k4/6], [0, 1, -1, 0], [0, 1, 0, 0],
 [0, 0, 0, 1]];
return Group( [a , b, c, alfa] , IdentityMat(4) );
end;

ACDim3Nr17 := function ( k1, k2, k3, k4)
local a, b, c, alfa, beta;
a:=[[1, 0, -k1/2, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]];
b:=[[1, k1/2, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 1]];
c:=[[1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]];
alfa:=[[1, -k1/2 + k2, k3, k4/6], [0, 1, -1, 0], [0, 1, 0, 0],
 [0, 0, 0, 1]];
beta:=[[-1, -k3, -k3, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]];
return Group( [a , b, c, alfa, beta] , IdentityMat(4) );
end;

#############################################################################
##
#V ACDim3Funcs
#V ACDim3Param
##
ACDim3Funcs := [ ACDim3Nr01, ACDim3Nr02, ACDim3Nr03, ACDim3Nr04, ACDim3Nr05,
ACDim3Nr06, ACDim3Nr07, ACDim3Nr08, ACDim3Nr09, ACDim3Nr10, ACDim3Nr11,
ACDim3Nr12, ACDim3Nr13, ACDim3Nr14, ACDim3Nr15, ACDim3Nr16, ACDim3Nr17 ];
MakeReadOnlyGlobal( "ACDim3Funcs" );

ACDim3Param := [ 1,4,2,2,2,4,4,4,4,4,4,4,4,4,4,4,4 ];
MakeReadOnlyGlobal( "ACDim3Param" );

max_stars_repo_name: alex-konovalov/aclib

max_issues_repo_name: alex-konovalov/aclib

max_forks_repo_name: alex-konovalov/aclib


text:
originalMatrix := [[1, 1, 1, 1],
                   [2, 4, 8, 16],
                   [3, 9, 27, 81],
                   [4, 16, 64, 256],
                   [5, 25, 125, 625]];
transposedMatrix := TransposedMat(originalMatrix);


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

Declare(IMDCT_Odd, IMDST_Odd);

Class(MDCTBase, DTTBase, rec(
    dims := self >> [ self.params, 2*self.params ],
    SmallRandom := () -> Random([2,4,6,8,12,16,18,24])
));

Class(IMDCTBase, DTTBase, rec(
    dims := self >> [ 2*self.params, self.params ],
    SmallRandom := () -> Random([2,4,6,8,12,16,18,24])
));

#F MDCT_Odd(<n>) - Modified Discrete Cosine Transform non-terminal, (oddly-stacked system)
#F Definition: (n x 2n)-matrix [ cos((2i+1)(2j+1+n)/(4n)) | i = 0..n-1, j = 0..2n-1 ]
#F Note:       MDCT_Odd is the transpose of IMDCT_Odd
#F Example:    MDCT_Odd(8)
#F
Class(MDCT_Odd, MDCTBase, rec(
    terminate := self >> let(n := self.params, Mat(
	List([0 .. n-1], i ->
	    List([0 .. 2*n - 1], j ->
		CosPi((2*j + 1 + n)*(2*i + 1)/(4*n)))))),

    transpose := self >> IMDCT_Odd(self.params),
));

MDCT := MDCT_Odd;

#F MDST_Odd(<n>) - Modified Discrete Sine Transform non-terminal, (oddly-stacked system)
#F Definition: (n x 2n)-matrix [ sin((2i+1)(2j+1+n)/(4n)) | i = 0..n-1, j = 0..2n-1 ]
#F Note:       MDST_Odd is the transpose of IMDST_Odd
#F Example:    MDST_Odd(8)
#F
Class(MDST_Odd, MDCTBase, rec(
    terminate := self >> let(n := self.params, Mat(
	List([0 .. n-1], i ->
	    List([0 .. 2*n - 1], j ->
		SinPi((2*j + 1 + n)*(2*i + 1)/(4*n)))))),

    transpose := self >> IMDST_Odd(self.params),
));

#F MDCT_Even(<n>) - Modified Discrete Cosine Transform non-terminal, (evenly-stacked system)
#F Definition: (n   x 2n)-matrix [ cos((2i)(2j+1+n)/(4n)) | i = 0..n-1, j = 0..2n-1 ], n even
#F             (n+1 x 2n)-matrix [ cos((2i)(2j+1+n)/(4n)) | i = 0..n,   j = 0..2n-1 ], n odd
#F
#F Example:    MDCT_Even(8)
#F
Class(MDCT_Even, MDCTBase, rec(
    dims := self >> [ self.params + When(IsOddInt(self.params), 1, 0), 2*self.params ],

    terminate := self >> let(n := self.params, isodd := When(IsOddInt(n), 1, 0), Mat(
	List([0 .. n-1 + isodd], i ->
	    List([0 .. 2*n - 1], j ->
		CosPi((2*j + 1 + n)*(2*i)/(4*n)))))),

#    transpose := self >> IMDCT_Even(self.params),
));

#F MDST_Even(<n>) - Modified Discrete Sine Transform non-terminal, (evenly-stacked system)
#F Definition: (n   x 2n)-matrix [ sin((2i)(2j+1+n)/(4n)) | i = 1..n,   j = 0..2n-1 ], n even
#F             (n-1 x 2n)-matrix [ sin((2i)(2j+1+n)/(4n)) | i = 1..n-1, j = 0..2n-1 ], n odd
#F
#F Example:    MDST_Even(8)
#F
Class(MDST_Even, MDCTBase, rec(
    dims := self >> [ self.params - When(IsOddInt(self.params), 1, 0), 2*self.params ],

    terminate := self >> let(n := self.params, isodd := When(IsOddInt(n), 1, 0), Mat(
	List([1 .. n-isodd], i ->
	    List([0 .. 2*n - 1], j ->
		SinPi((2*j + 1 + n)*(2*i)/(4*n)))))),

#    transpose := self >> IMDST_Even(self.params),
));

#F IMDCT_Odd(<n>) - Inverse Modified Discrete Cosine Transform non-terminal (oddly-stacked system)
#F Definition: (2n x n)-matrix [ cos((2j+1)(2i+1+n)/(4n)) | i=0..2n-1, j=0..n-1 ]
#F Note:       IMDCT_Odd is the transpose of MDCT_Odd
#F Example:    IMDCT_Odd(8)
#F
Class(IMDCT_Odd, IMDCTBase, rec(
    terminate := self >> let(n:=self.params, Mat(
	List([0 .. 2*n-1], i ->
	    List([0 .. n-1], j ->
		CosPi((2*i + 1 + n)*(2*j + 1)/(4*n)))))),

    transpose := self >> MDCT_Odd(self.params),
));

IMDCT := IMDCT_Odd;

#F IMDST_Odd(<n>) - Inverse Modified Discrete Sine Transform non-terminal (oddly-stacked system)
#F Definition: (2n x n)-matrix [ cos((2j+1)(2i+1+n)/(4n)) | i=0..2n-1, j=0..n-1 ]
#F Note:       IMDST_Odd is the transpose of MDST_Odd
#F Example:    IMDST_Odd(8)
#F
Class(IMDST_Odd, IMDCTBase, rec(
    terminate := self >> let(n:=self.params, Mat(
	List([0 .. 2*n-1], i ->
	    List([0 .. n-1], j ->
		SinPi((2*i + 1 + n)*(2*j + 1)/(4*n)))))),

    transpose := self >> MDST_Odd(self.params),
));

_shiftcut := (transform, n, start, cut, shift, boundary) -> Checked(cut in ["re", "im"],
    let(ofs := When(cut="re", 0, 1),
	N   := Cols(transform),
	i   := Ind(N),

	Gath(H(Rows(transform), n, ofs+start, 2)) *
	transform *
	When(boundary=1, I(N), Diag(Lambda(i, cond(leq(i,shift-1), boundary, 1)))) *
	Z(N, -shift))
);

RulesFor(MDCT_Odd, rec(
    #F MDCT_toDCT4:  MDCT_n = DCT4_n * sums
    #F
    MDCT_Odd_toDCT4 := rec (
	info         := "MDCT_n -> DCT4_n",
	isApplicable := P -> P mod 2 = 0,
	allChildren  := P -> [[ DCT4(P) ]],
	rule := (P, C) ->
	    C[1] *
	    J(P) *
	    DirectSum(Tensor(Mat([[1, -1]]), I(P/2)),
		      Tensor(Mat([[-1, -1]]), I(P/2))) *
	    DirectSum(J(P/2), I(P/2), I(P/2), J(P/2))
    ),

    MDCT_Odd_toPRDFT34 := rec(
	info         := "MDCT_n -> PRDFT4_2n or PRDFT3_2n",
	isApplicable := P -> true,
	allChildren  := P -> When(IsEvenInt(P), [[ PRDFT4(2*P) ]], [[ PRDFT3(2*P) ]]),
	rule := (P, C) -> let(n:=P,
	    _shiftcut(C[1], n, 0, "re", Int((n+1)/2), -1))
    )

));

RulesFor(MDST_Odd, rec(
    #F MDST_toDST4:  MDST_n = DST4_n * sums
    #F
    MDST_Odd_toDST4 := rec (
	info         := "MDST_n -> DST4_n",
	isApplicable := P -> P mod 2 = 0,
        # NOTE: implement rule for n odd, converts to DST3
	allChildren  := P -> [[ DST4(P) ]],
	rule := (P, C) -> let(n:=P,
	    C[1] *
	    J(n) *
	    DirectSum(Tensor(Mat([[1, 1]]), I(n/2)),
		      Tensor(Mat([[1, -1]]), I(n/2))) *
	    DirectSum(J(n/2), I(n/2), I(n/2), J(n/2)))
    ),

    MDST_Odd_toPRDFT34 := rec(
	info         := "MDST_n -> PRDFT4_2n or PRDFT3_2n",
	isApplicable := P -> true,
	allChildren  := P -> When(IsEvenInt(P), [[ PRDFT4(2*P) ]], [[ PRDFT3(2*P) ]]),
	rule := (P, C) -> let(n:=P,
	    _shiftcut(C[1], n, 0, "im", Int((n+1)/2), -1))
    )

));

RulesFor(MDCT_Even, rec(
    MDCT_Even_toPRDFT12 := rec(
	info         := "MDCT_n -> PRDFT1_2n or PRDFT2_2n",
	isApplicable := P -> true,
	allChildren  := P -> When(IsEvenInt(P), [[ PRDFT2(2*P) ]], [[ PRDFT1(2*P) ]]),
	rule := (P, C) -> let(n:=P,
	    _shiftcut(C[1], n+(n mod 2), 0, "re", Int((n+1)/2), 1))
    )
));

RulesFor(MDST_Even, rec(
    MDST_Even_toPRDFT12 := rec(
	info         := "MDST_n -> PRDFT1_2n or PRDFT2_2n",
	isApplicable := P -> true,
	allChildren  := P -> When(IsEvenInt(P), [[ PRDFT2(2*P) ]], [[ PRDFT1(2*P) ]]),
	rule := (P, C) -> let(n:=P,
	    _shiftcut(C[1], n-(n mod 2), 2, "im", Int((n+1)/2), 1))
    )
));


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

Declare(Tensor, I, L);

# ==========================================================================
# Tensor(<spl1>, <spl2>, ...) . . . . . . . . . .  Kronecker product of SPLs
#    Arguments can be lists of spls of any depth. So that the following is
#    legal:
#       Tensor(I(2), [[F(2), [I(2)]]])
#    and is equivalent to
#       Tensor(I(2), F(2), I(2))
# ==========================================================================
Class(Tensor, BaseOperation, rec(
    abbrevs := [ arg -> [Flat(arg)] ],

    new := meth(self, L)
        local scl,mat;
        Constraint(IsList(L) and Length(L) >= 1);
        L:=Filtered(L,x->x<>I(1));
        if (Length(L)=0) then return I(1); fi;
        #DoForAll(L, x -> Constraint(IsSPL(x)));
	if Length(L)=1 then return L[1]; fi;

	#if Length(L)=2 then
	#    scl := Cond(Dimensions(L[1])=[1,1], 1,
	#	        Dimensions(L[2])=[1,1], 2,
	#		false);
	#    if scl <> false then
	#	mat := L [ When(scl=1, 2, 1) ];
	#	scl := MatSPL(L[scl])[1][1];
	#
	#	return When(scl=1, mat, spiral.spl.operators.Scale(scl, mat));
	#    fi;
	#fi;

        return SPL(WithBases( self,
	    rec( _children := L,
		dimensions := [ Product(L, t -> t.dims()[1]),
		                Product(L, t -> t.dims()[2]) ] )));
    end,
    #-----------------------------------------------------------------------
    dims := self >> [ Product(self._children, t -> t.dims()[1]),
	              Product(self._children, t -> t.dims()[2]) ],
    #-----------------------------------------------------------------------
    isPermutation := self >> ForAll(self._children, IsPermutationSPL),
    #-----------------------------------------------------------------------
    isSymmetric := self >> ForAll(self._children, c -> c.isSymmetric()),
    #-----------------------------------------------------------------------
    toAMat := self >> TensorProductAMat(List(self._children, AMatSPL)),
    #-----------------------------------------------------------------------
    transpose := self >>
        CopyFields(self, rec(_children := List(self._children, x->x.transpose()),
		          dimensions := Reversed(self.dimensions))),
    conjTranspose := self >>
        CopyFields(self, rec(_children := List(self._children, x->x.conjTranspose()),
		          dimensions := Reversed(self.dimensions))),
    inverse := self >>
        CopyFields(self, rec(_children := List(self._children, x->x.inverse()),
		          dimensions := Reversed(self.dimensions))),
    #-----------------------------------------------------------------------
    rightBinary  := self >> FoldR1(self._children, (p,x) -> let(base:=self.__bases__[1], base(x, p))),

    leftBinary := self >> FoldL1(self._children, (p,x) -> let(base:=self.__bases__[1], base(p, x))),

    split := meth(self)
        local i, lft, rt, ch, res;
        ch := self.children();
        rt := Product(Drop(ch, 1), Rows);
        lft := 1;
        res := Tensor(ch[1], I(rt));
        for i in [2..self.numChildren()] do
            lft := lft * Cols(ch[i-1]);
            rt := rt / Rows(ch[i]);
            res := res * Tensor(I(lft), ch[i], I(rt));
        od;
        return res;
    end,

    vectorForm := self >> Checked(self.numChildren() = 2,
	let(c1 := self.child(1), c2 := self.child(2),

	    Cond(IsIdentitySPL(c1) and IsIdentitySPL(c2), self,
		 IsIdentitySPL(c1), L(Rows(c2)*Rows(c1), Rows(c1)) *
		                    Tensor(c2, I(Rows(c1))) *
				    L(Cols(c2)*Cols(c1), Cols(c2)),
		 IsIdentitySPL(c2), self,
	         # else
		 Tensor(c1, I(Rows(c2)) * Tensor(I(Cols(c1)), c2).vectorForm()
	    )))),

    parallelForm := self >> Checked(self.numChildren() = 2,
	let(c1 := self.child(1), c2 := self.child(2),

	    Cond(IsIdentitySPL(c1) and IsIdentitySPL(c2), self,
		 IsIdentitySPL(c2), L(Rows(c2)*Rows(c1), Rows(c1)) *
		                    Tensor(I(Rows(c2)), c1) *
				    L(Cols(c2)*Cols(c1), Cols(c2)),
		 IsIdentitySPL(c1), self,
	         # else rc1*rc2 X cc1*cc2
		 Tensor(c1, I(Rows(c2))).parallelForm() * Tensor(I(Cols(c1)), c2)
	    ))),

    #-----------------------------------------------------------------------
    arithmeticCost := meth(self, costMul, costAddMul)
        local t, left, right, cost;
        left  := 1;
	right := self.dimensions[1];
	cost  := costMul(0) - costMul(0); # will work even when costMul(0) <> 0
	for t in self.children() do
	    right := right / t.dimensions[1];
	    cost  := cost + left * right * t.arithmeticCost(costMul, costAddMul);
	    left  := left * t.dimensions[2];
	od;
	return cost;
    end,
    #-----------------------------------------------------------------------
    latexSymbol := "\\tensor"
));


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

#F FreeVars(<o>)
#F   list of free (unbound) variables in an object
#F   Expects a .free() method in <o>.
FreeVars := o -> Cond(
    IsRec(o) and IsBound(o.free), o.free(),
    IsList(o), Set(Union(List(o, FreeVars))),
    Set([])
);


text:
InstallMethod(UnipotentClass,
              "Unipotent conjugacy class in algebraic group",
              [IsString,IsUnipotent,IsList,IsList],
              function(label,u,torus,parabolics)
              local object,
              sys,permute_rep,exp;

              Print(label,"\n");

              object:=Objectify(NewType(NewFamily("UnipotentClassFamily"),
                                        IsAttributeStoringRep and
                                        IsUnipotentClass),
                                rec());

              sys:=chevalleyAdj(u);
              SetalgebraicU(object,sys);
              SetLabel(object,label);
              SetRepresentative(object,u);

#              permute_rep:=function(rep,torus)
#                  local result,signs,pr,perm,perm_matrix;
#
#                  perm:=PositiveTorusPermutation(rootSystem(sys),torus);
#                  perm_matrix:=Inverse(PermutationMatrix(weylGroup(sys),perm));
#                  pr:=positiveRoots(sys);
#
#                  #result:=List(rep,j->perm_matrix*pr[j[1]]);
#                  result:=List(rep,j->pr[j[1]]*perm_matrix);
#                  signs:=List(rep,j->PermutationSign(sys,perm,pr[j[1]]));
#                  #Print(List([1..Length(result)],j->[Position(pr,result[j]),signs[j]*rep[j][2]]));
#                  return List([1..Length(result)],j->[Position(pr,result[j]),signs[j]*rep[j][2]]);
#              end;

              if torus <> [] then
                  SetCocharacter(object,torus);
                  SetBorelPerm(object,PositiveTorusPermutation(rootSystem(sys),torus));
                  #SetBorelRep(object,Unipotent(sys,permute_rep(coefficients(u),torus)));# !!! Verifica asta
                  #SetBorelRep(object,Unipotent(sys,permute_rep(coefficients(u),torus),
                  #                             [1..Length(positiveRoots(sys))]));# !!! Verifica asta
                  SetBorelRep(object,Unipotent(sys,coefficients(ToPositiveBorel(object,u)),
                                               [1..Length(positiveRoots(sys))]));# !!! Verifica asta
                  SetOrder(BorelRep(object),Order(u));
              fi;

              if parabolics <> [[],[]] then
                  SetSupportParabolic(object,ParabolicSystem(rootSystem(sys),parabolics[1]));
                  SetDistinguishedParabolic(object,ParabolicSystem(rootSystem(sys),parabolics[2]));
              fi;

              if HasOrder(u) then SetOrder(object,Order(u)); fi;
#              elif Characteristic(sys)>0 then SetOrder(object,OrderOp(u)); fi; This can take to much time

              SetName(object,Concatenation("<unipotent conjugacy class ",label," for ",type(sys),String(rank(sys)),
                                           " in characteristic ",String(Characteristic(sys)),">"));
              return object;
end);

InstallMethod(ToPositiveBorel,
              "Permutes root to positive Borel subgroup with the cocharacter",
              [IsUnipotentClass,IsPosInt],
              function(orb,root_index)
              local sys,perm,perm_mat,pr,roots;

              sys:=algebraicU(orb);
              perm:=PositiveTorusPermutation(rootSystem(sys),Cocharacter(orb));
              perm_mat:=Inverse(PermutationMatrix(weylGroup(sys),perm));
              pr:=positiveRoots(sys);
              roots:=allRoots(sys);

              #return Position(pr,perm_mat*roots[root_index]);
              return Position(pr,roots[root_index]*perm_mat);
end);

InstallMethod(ToPositiveBorel,
              "Permutes roots (in roots_index) to positive Borel subgroup with the cocharacter",
              [IsUnipotentClass,IsList],
              function(orb,roots_index)
              return List(roots_index,i->ToPositiveBorel(orb,i));
end);

InstallMethod(FromPositiveBorel,
              "Permutes root from positive Borel subgroup with the cocharacter",
              [IsUnipotentClass,IsPosInt],
              function(orb,root_index)
              local sys,perm,perm_mat,pr,roots;

              sys:=algebraicU(orb);
              perm:=PositiveTorusPermutation(rootSystem(sys),Cocharacter(orb));
              perm_mat:=Inverse(PermutationMatrix(weylGroup(sys),Reversed(perm)));
              #perm_mat:=PermutationMatrix(weylGroup(sys),perm);
              pr:=positiveRoots(sys);
              roots:=allRoots(sys);

              #return Position(pr,perm_mat*roots[root_index]);
              return Position(pr,roots[root_index]*perm_mat);
end);

InstallMethod(FromPositiveBorel,
              "Permutes roots (in roots_index) from positive Borel subgroup with the cocharacter",
              [IsUnipotentClass,IsList],
              function(orb,roots_index)
              return List(roots_index,i->FromPositiveBorel(orb,i));
end);

InstallMethod(ToPositiveBorel,
              "Permutes roots (in roots_index) from positive Borel subgroup with the cocharacter",
              [IsUnipotentClass,IsUnipotent],
              function(orb,u)
              local result,coeffs,sys,signs,pr,perm,perm_matrix;

              sys:=chevalleyAdj(u);

              perm:=BorelPerm(orb);
              #perm_matrix:=Inverse(PermutationMatrix(weylGroup(sys),perm));
              perm_matrix:=PermutationMatrix(weylGroup(sys),Reversed(perm));
              pr:=positiveRoots(sys);

              coeffs:=coefficients(u);
              result:=List(coeffs,j->pr[j[1]]*perm_matrix);
              signs:=List(coeffs,j->InversePermutationSign(sys,perm,pr[j[1]]));
              result:=List([1..Length(result)],j->[Position(pr,result[j]),signs[j]*coeffs[j][2]]);
####              result:=List([1..Length(result)],j->[Position(pr,result[j]),coeffs[j][2]]);
#Print(result,"\n");
              if Filtered(result,i->i[1]=fail)<>[] then return fail; fi;

              result:=Unipotent(sys,result);#,Ordering(u));
              if HasOrder(u) then SetOrder(result,Order(u)); fi;

              return result;
end);

InstallMethod(FromPositiveBorel,
              "Permutes roots (in roots_index) from positive Borel subgroup with the cocharacter",
              [IsUnipotentClass,IsUnipotent],
              function(orb,u)
              local result,coeffs,sys,signs,pr,perm,perm_matrix;

              sys:=chevalleyAdj(u);

              perm:=BorelPerm(orb);
              #perm_matrix:=Inverse(PermutationMatrix(weylGroup(sys),perm));
              #perm_matrix:=Inverse(PermutationMatrix(weylGroup(sys),Reversed(perm)));
              perm_matrix:=PermutationMatrix(weylGroup(sys),perm);
              pr:=positiveRoots(sys);

              coeffs:=coefficients(u);
              result:=List(coeffs,j->pr[j[1]]*perm_matrix);
              signs:=List(coeffs,j->PermutationSign(sys,perm,pr[j[1]]));
              result:=List([1..Length(result)],j->[Position(pr,result[j]),signs[j]*coeffs[j][2]]);
####              result:=List([1..Length(result)],j->[Position(pr,result[j]),coeffs[j][2]]);

              if Filtered(result,i->i[1]=fail)<>[] then return fail; fi;

              result:=Unipotent(sys,result);#,Ordering(u));
              if HasOrder(u) then SetOrder(result,Order(u)); fi;

              return result;
end);

#
# ---------- Overclass for classes --------
#

InstallMethod(UnipotentClasses,
              "Over class for unipotent orbits of a simple algebraic group",
              [IsAlgebraicU,IsString],
              function(sys,source)
              local object,
              pr_len,p,
              APR,avarnames,
              file,dir,reps,exps,
              rep,root;

              object:=Objectify(NewType(NewFamily("UnipotentClassesFamily"),
                                        IsAttributeStoringRep and
                                        IsUnipotentClasses),
                                rec());

              SetalgebraicU(object,sys);
              p:=Characteristic(sys);

              dir:=data_dir;
              if type(sys) in ["B","C","D","E","F","G"] and p = 2 or
                 type(sys) in ["E","F","G"] and p = 3 or
                 type(sys) = "E" and rank(sys) = 8 and p = 5 then
#              if type(sys) in ["B","C","D","F"] and p = 2 or
#                 type(sys) = "E" and rank(sys) in [7,8] and p = 2 or
#                 type(sys) in ["E","G"] and rank(sys) in [2,8] p = 3 or
                   file:=Concatenation("data",type(sys),String(rank(sys)),"char",String(p),source,".gi");
              Print(file,"\n");
#                   if source then file:=Concatenation("data",type(sys),String(rank(sys)),"char",String(p),"Mizuno.gi");
#                   else file:=Concatenation("data",type(sys),String(rank(sys)),"char",String(p),".gi"); fi;
              else file:=Concatenation("data",type(sys),String(rank(sys)),".gi"); fi;
              file:=Filename(dir,file);
              file:=ReadAsFunction(file);
              file:=file();

              if file[6]=[] then file[6]:=List([1..Length(file[1])],i->-1); fi;# if Orders not given then compute them
              reps:=List([1..Length(file[1])],i->Unipotent(sys,file[1][i],file[6][i]));
              if file[7]<>[] then
                  for rep in [1..Length(reps)] do
                      SetAduJordanBlocks(reps[rep],file[7][rep]);
                  od;
              fi;

              if p<>0 then
                  exps:=List(reps,i->LogInt(Order(i),p));
                  SetWittStructure(object,Witt(Maximum(exps),sys));
              fi;

              SetAllClasses(object,List([1..Length(file[1])],i->UnipotentClass(#witt_structures[exp[i]],
                  file[3][i],
                  reps[i],
                  file[2][i],
                  [file[4][i],file[5][i]])));

              SetName(object,Concatenation("<unipotent conjugacy classes ",type(sys),String(rank(sys)),
                                           " in characteristic ",String(Characteristic(sys)),">"));

              return object;
end);


text:
# Copyright (c) 2018-2020, Carnegie Mellon University
# See LICENSE for details
#
# Contains the CNOT and SWAP non-terminals and rewrite rules
#

##
#F qSWAPT( <qi>, <qj>, <n>,, <arch> ) Object
##
#F Awaps qi and qj in a system of n qubits, with architecture arch
#F Definition: (2^n x 2^n)-matrix
#F that swaps qubit qi and qubit qj, with the identity operation applied to all others
Class(qSWAPT, TaggedNonTerminal, rec(
  abbrevs   := [ (qi, qj, n, arch) -> Checked(IsPosInt(n), [qi, qj, n, arch]) ],
  dims      := self >> let(size := 2^self.params[3], [size, size]),
  terminate := self >> GenSwapMat(self.params[1], self.params[2], self.params[3]),
  isReal    := self >> true,
  rChildren := self >> self.params,
  from_rChildren := (self, rch) >> self.__bases__[1](rch[1], rch[2], rch[3], rch[4]),
  SmallRandom := () -> Random([2..5]),
  LargeRandom := () -> Random([6..15]),
  normalizedArithCost := self >> Error("ArithCost not implemented"),
  TType := T_Real(64)
));

##
#F CnotTerm( <i>, <j>, <n> )
##
## .terminate() function for the qCNOT Non-terminal. Expands qCNOT non-terminal into a 2^n by 2^n matrix
CnotTerm := function (i, j, n)
    return When(i < j, Tensor(Tensor(I(2^(i)), GenCnotMat(j, 0)), I(2^(n-j))), Tensor(Tensor(I(2^(j)), GenCnotMat(i, 1)), I(2^(n-i)))); # create the intended cnot matrix
end;

##
#F CnotSPL( <i>, <j>, <n> )
##
## Used for the qCNOT_Base breakdown rule. Represents a qCNOT Non-terminal as a tensor of CNOT and I SPL objects
CnotSPL := function (i, j, n)
    return When(i < j, Tensor(Tensor(I(2^(i)), CNOT(j-i, 0)), I(2^(n-1-j))), Tensor(Tensor(I(2^(j)), CNOT(i-j, 1)), I(2^(n-1-i)))); # create the intended cnot matrix
end;

##
#F qCNOT( <qj>, <dir>, <arch> ) Object
##
#F CNOT Non-terminal
#F Definition: (2^n x 2^n)-matrix  that applies a CNOT from qubit 0 to j
Class(qCNOT, TaggedNonTerminal, rec(
  abbrevs   := [ (qj, dir, arch) -> Checked([qj, dir, arch]) ],
  dims      := self >> let(size := 2^(self.params[1]+1), [size, size]),
  terminate := self >> GenCnotMat(self.params[1], self.params[2]),
  isReal    := self >> true,
  groups    := self >> [[2]],
  recursive_def := (self, arch) >> self.__bases__[1](self.params[1], self.params[2], arch),
  rChildren := self >> self.params,
  from_rChildren := (self, rch) >> self.__bases__[1](rch[1], rch[2], rch[3]),
  connected := self >> true,
  SmallRandom := () -> Random([2..5]),
  LargeRandom := () -> Random([6..15]),
  normalizedArithCost := self >> Error("ArithCost not implemented"),
  TType := T_Real(64)
));

##
#F qCNOT Breakdown Rules
##
NewRulesFor(qCNOT, rec(

    #F qCNOT_Base: qCNOT(_) =  CNOT() SPL object
    #F Directly represent as an implementable gate, if a proper edge exists
    qCNOT_Base := rec(
        info             := "qCNOT_(_) -> CNOT(_)",
        forTransposition := false,
        applicable       := (self, nt) >> (HasEdge(0, nt.params[1], nt.params[3]) = 1),
        apply            := (nt, c, cnt) -> CNOT(nt.params[1], nt.params[2]),
    )

));

##
#F CnotSPL( <i>, <j>, <n> )
##
## Used for the qSWAP_Base breakdown rule. Represents a qCNOT Non-terminal as a tensor of CNOT and I SPL objects
CnotSPL := function (i, j, n)
    return When(i < j, Tensor(Tensor(I(2^(i)), CNOT(j-i, 0)), I(2^(n-1-j))), Tensor(Tensor(I(2^(j)), CNOT(i-j, 1)), I(2^(n-1-i)))); # create the intended cnot matrix
end;

##
#F SwapGenChildren( <start>, <ending>, <n>, <arch> )
##
## decompose a qSWAPT transform into smaller qSWAPT transforms
SwapGenChildren := function (start, ending, n, arch)
    local move, children;
    move := Best_Move(start, ending, arch);
    children := [];
    Add(children, [qSWAPT(start, move, n, arch) * qSWAPT(move, ending, n, arch) * qSWAPT(start, move, n, arch)]);
    return children;
end;

##
#F qSWAPT Breakdown Rules
##
NewRulesFor(qSWAPT, rec(

    #F qSWAPT_Base: qSWAPT(_) =  qCNOT() SPL objects
    #F SWAP to CNOT conversion, only if an edge exists
    #F represent a swap as a series of 3 CNOTs with the middle being an alternate direction
    qSWAPT_Base := rec(
        info             := "qSWAPT() -> CNOT * CNOT * CNOT ",
        forTransposition := false,
        applicable       := (self, nt) >> ( (HasEdge(nt.params[1], nt.params[2], nt.params[4]) = 1) and (nt.params[1] <> nt.params[2])),
        children         := nt -> List( [ [ (CnotSPL(nt.params[1], nt.params[2], nt.params[3]))*(CnotSPL(nt.params[2], nt.params[1], nt.params[3]))*(CnotSPL(nt.params[1], nt.params[2], nt.params[3])) ] , [ (CnotSPL(nt.params[2], nt.params[1], nt.params[3]))*(CnotSPL(nt.params[1], nt.params[2], nt.params[3]))*(CnotSPL(nt.params[2], nt.params[1], nt.params[3])) ]  ] ),
        apply            := (nt, c, cnt) -> c[1],
    ),

    #F qSWAPT_ID: qSWAPT(_) =  I() SPL object
    #F SWAP to Identity conversion if the swap is trivial
    qSWAPT_Id := rec(
        info             := "qSWAPT() -> Identity ",
        forTransposition := false,
        applicable       := (self, nt) >> (nt.params[1] = nt.params[2]),
        apply            := (nt, c, cnt) -> I(2^(nt.params[3])),
    ),

    #F qSWAPT_Rec: qSWAPT(_) =   qSWAPT(_) *  qSWAPT(_)
    #F Break down an unimplementable qSWAP into an implementable swap and a swap that is 1 step closer to being implementable
    #F TODO, make this a bit more intelligent, right now SwapGenChildren just does a random walk
    qSWAPT_Rec := rec(
        info             := "qSWAPT() -> qSWAPT() * qSWAPT()",
        forTransposition := false,
        applicable       := (self, nt) >> ((nt.params[1] <> nt.params[2]) and (HasEdge(nt.params[1], nt.params[2], nt.params[4]) = 0)),
        children         := nt -> SwapGenChildren(nt.params[1], nt.params[2], nt.params[3], nt.params[4]),
        apply            := (nt, c, cnt) -> c[1],
    )

));

##
#F ShiftMat( <swap1>, <swap2>, <arch>, <n> )
##
## generates a martix of SWAP operations that swap the qubit groups swap1 and swap2 in parallel
## swap1 is the target position, and swap2 are original positions
ShiftMat := function (swap1, swap2, arch, n)
    local idx, format, backmat, targets, t, positions, e, sw1, sw2, s, i, test, already_swapped;
    idx := 1;
    format := I(2^(n));
    backmat := I(2^(n));
    positions := [];
    targets := [];
    for s in swap1 do
        Add(positions, s);
    od;
    for t in swap2 do
        Add(targets, t);
    od;
    for s in swap1 do
        if ((targets[idx] <> positions[idx])) then # just save the extra trouble, if the swap is trivial
            format := format * qSWAPT(positions[idx], targets[idx], n ,arch);
            backmat := qSWAPT(positions[idx], targets[idx], n , arch) * backmat;
            sw1 := targets[idx];
            sw2 := positions[idx];
            i := 1;
            for e in positions do
                if e = sw1 then
                    positions[i] := sw2;
                fi;
                if e = sw2 then
                    positions[i] := sw1;
                fi;
                i := i + 1;
            od;
        fi;
        idx := idx + 1;
    od;
    return [format, backmat];
end;

max_stars_repo_name: spiral-software/spiral-package-quantum

max_issues_repo_name: spiral-software/spiral-package-quantum

max_forks_repo_name: spiral-software/spiral-package-quantum


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

normCos := normalizeCosine;
dnormCos := x -> When(IsDelay(x) and Type(x)=T_FUNCCALL and x[1] = D(normCos),
	              x, PEV(normCos(x)));

dnormCos := x -> When(IsRat(x), normCos(x), x);

#F SkewDTT(<non-terminal for DCT/DST type 3/4>, <rat>)
#F Parameters:       [ <transform for DCT or dst type 3 or 4>, <rat> ]
#F Definition:       Let t be dct or dst type 3 or 4 as non-terminal spl.
#F                   Transform( "SkewDTT", [t, r] ) represents the (n x n)-matrix
#F                     [ cos(r_k*l*pi/n) | k,l = 0...n-1 ]
#F                     [ sin(r_k*(l+1)*pi/n) | k,l = 0...n-1 ]
#F                     [ cos(r_k*(l+1/2)*pi/n) | k,l = 0...n-1 ]
#F                     [ sin(r_k*(l+1/2)*pi/n) | k,l = 0...n-1 ]
#F                   for dct3, dst3, dct4, dst4, respectively
#F                   where the list L of the r_k arises from the list
#F                     [ (r + 2k)/n | 0 <= k < n ]
#F                   by
#F                     1. normalizing the entries to lie in [0, 1] such
#F 		       that their cosine does not change
#F                     2. ordering the values increasing by size
#F Note:             for r = 1/2, the skew DCT3 reduces to the ordinary DCT3.
#F                   The transpose is provided by the transform TSkewDTT
#F Example:          SkewDTT( DCT3(8), 1/3 )
#F
Class(SkewDTT, NonTerminal, rec(
    _short_print := true,
    abbrevs := [ (nt, r) -> Checked(IsNonTerminal(nt), IsRatSym(r),
	                            ObjId(nt) in [DCT3, DST3, DCT4, DST4],
				    [Copy(nt), dnormCos(r)]) ],

    dims := self >> self.params[1].dimensions,

    terminate := self >> let(
	S := ObjId(self.params[1]),
	N := self.params[1].dimensions[1],
	r := EvalScalar(self.params[2]),
	Cond(S = DCT3, Mat(SkewDCT_IIIunscaled(N, r)),
	     S = DST3, Mat(SkewDST_IIIunscaled(N, r)),
	     S = DCT4, Mat(SkewDCT_IVunscaled(N, r)),
	     S = DST4, Mat(SkewDST_IVunscaled(N, r)),
	     Error("Unrecognized <self.params[1].symbol>"))),

    isReal := True

#    HashId := self >> self.params[1]
));

RulesFor(SkewDTT, rec(
    SkewDTT_Base2 := rec(
        info             := "DTT_2(r) -> F_2",
        isApplicable     := P -> P[1].dimensions[1] = 2,
        allChildren      := P -> [ [ ] ],
        rule := (P,C) ->
	    let(S := ObjId(P[1]),
		r := P[2],
		Cond(
		    S = DCT3,  F(2) * Diag([1, CosPi(r/2)]),

		    S = DST3,  F(2) * Diag([SinPi(r/2), SinPi(r)]),

		    S = DCT4,  Diag([CosPi(r/4), SinPi(r/4)]) * F(2) *
				 Mat([[1, -1], [0, 2*CosPi(r/2)]]),

		    S = DST4,  Diag([SinPi(r/4), CosPi(r/4)]) * F(2) *
				 Mat([[1, 1], [0, 2*CosPi(r/2)]])))
    ),

    SkewDCT3_VarSteidl := rec(
        info             := "SkewDCT3_n -> IterDirectSum(SkewDCT3_n/2)",
	switch           := true,
	forTransposition := false,
        isApplicable     := P -> let(S := ObjId(P[1]), n := Rows(P[1]),
	    S = DCT3 and IsPrimePowerInt(n) and n mod 2 = 0 and n > 2),

	# create one parametrized SkewDTT, use SampledVar() to make timing possible
        allChildren := P -> let(n := Rows(P[1]), r := P[2],
	    [[ SkewDTT(DCT3(n/2), r/2),
	       SkewDTT(DCT3(n/2), 1 - r/2) ]]),

        rule := (P,C) -> let(n := Rows(P[1]), r := P[2],
	    K(n, n/2) *
	    DirectSum(C[1], C[2]) *
	    Tensor(F(2), I(n/2)) *
	    sparse1(n, r/2))
    ),

    # NOTE: We need to be able to time SkewDTT(n, r) where r is free variable
    #        otherwise DP does not work
    SkewDCT3_VarSteidlIterative := rec(
        info             := "SkewDCT3_n -> IterDirectSum(SkewDCT3_n/2)",
	switch           := false,
	forTransposition := false,
        isApplicable     := P -> let(S := ObjId(P[1]), n := Rows(P[1]),
	    S = DCT3 and IsPrimePowerInt(n) and n mod 2 = 0 and n > 2),

	# create one parametrized SkewDTT, use SampledVar() to make timing possible
        allChildren := P -> let(n := Rows(P[1]), r := P[2], i:=Ind(2),
	    [[ SkewDTT(DCT3(n/2), cond(i,1,0) + cond(i,-1,1)*r/2) ]]),

        rule := (P,C, Nonterms) -> let(
	    n := Rows(P[1]),
	    r := P[2],
	    ivar := Nonterms[1].params[2].args[1].args[1],

	    K(n, n/2) *
	    IterDirectSum(ivar, 2, C[1]) *
	    Tensor(F(2), I(n/2)) *
	    sparse1(n, r/2))
    )
));


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

#
# Using these "radix-2" rules will yield optimal 2nlog(n)-4n+6
# arithmetic cost for PRDFT1 (Split-Radix), and an almost equivalent
# cost of 2nlog(n)-4n+8 for PDHT1 (better than the literature).
#

# Each entry is J(2) times conjugated twiddle w* times PDHT3(2)=F(2),
# twiddles are accessed at stride 2.
# (w^*) * F_2 = [[r,i],[-i,r]] * F_2 = [[r+i, r-i], [r-i], [-i-r]]
Class(H3_CasTwid, Sym, rec(
    def := (n, k) ->
	DirectSum(List([0..n/2-1], j -> let(w := E(2*n)^(k*j),
		    J(2) * Mat([[Re(w)+Im(w),  Re(w)-Im(w)],
			        [Re(w)-Im(w), -Re(w)-Im(w)]]))))
));

# These twiddles are combined with PRDFT3(2) (just as for PDHT3 above)
# since PRDFT3(2,1)=I(2), and PRDFT3(2,-1)=Diag(1,-1) we have the
# variable 'm' scaling the second column of rotations.
Class(R3_Twid, Sym, rec(
    def := (n, k) -> let(m := When(k mod 4 = 1, 1, -1),
	DirectSum(List([0..n/2-1], j -> let(w := E(2*n)^(k*j),
		           Mat([[Re(w), -Im(w)*m],
				[Im(w),  Re(w)*m]])))))
));
# since IPRDFT2(2,1)=Diag(1,-1), and IPRDFT2(2,-1)=I(2) we have 'm' again,
# it scales second row
Class(IR3_Twid, Sym, rec(
    def := (n, k) -> let(m := When(k mod 4 = 1, -1, 1),
	DirectSum(List([0..n/2-1], j -> let(w := E(2*n)^(k*j),
		           Mat(2*  [  [Re(w), -Im(w)],
				    m*[Im(w),  Re(w)]])))))
));

RulesFor(PDHT3, rec(
    PDHT3_CT_Radix2 := rec(
	isApplicable := P -> P[1] > 2 and IsEvenInt(P[1]),
	allChildren  := P -> [[ DFT1(P[1]/2, P[2]) ]],
	rule := (P,C) -> let(N := P[1], k := P[2],
	        RC(LIJ(N/2)) *
		DirectSum(Tensor(I(N/4), J(2)), I(N/2)) *
		RC(C[1]) *
		H3_CasTwid(N,k) *
		L(N,N/2)
	))
));
RulesFor(PRDFT3, rec(
    PRDFT3_CT_Radix2 := rec(
	isApplicable := P -> P[1] > 2 and IsEvenInt(P[1]),
	allChildren  := P -> [[ DFT1(P[1]/2, P[2]) ]], # SkewDFT(n, E(4)), and then no twids are necessary
	rule := (P,C) -> let(N := P[1], k := P[2],
	        RC(LIJ(N/2)) *
		Diag(BHD(N/2, 1, -1)) *
		RC(C[1]) *
		R3_Twid(N,k) *
		L(N,N/2)
	))
));
RulesFor(IPRDFT2, rec(
    IPRDFT2_CT_Radix2 := rec(
	isApplicable := P -> P[1] > 2 and IsEvenInt(P[1]),
	allChildren  := P -> [[ DFT1(P[1]/2, P[2]) ]], # SkewDFT(n, E(4)), and then no twids are necessary
	rule := (P,C) -> let(N := P[1], k := P[2],
		L(N,2) *
		IR3_Twid(N,k) *
		RC(C[1]) *
		Diag(BHD(N/2, 1, -1)) *
	        RC(LIJ(N/2).transpose())
	))
));

RulesFor(PDHT1, rec(
    PDHT1_CT_Radix2 := rec(
	isApplicable := P -> P[1] > 2 and IsEvenInt(P[1]),
	allChildren  := P -> [[ PDHT1(P[1]/2, P[2]), PDHT3(P[1]/2, P[2]) ]],
	rule := (P,C) -> let(N := P[1], k := P[2],
	        RC(OddStride(N/2+1, N/4+1)) *
		DirectSum(C[1], C[2]) *
		Tensor(F(2), I(N/2))
	))
));

RulesFor(PRDFT1, rec(
    PRDFT1_CT_Radix2 := rec(
	isApplicable := P -> P[1] > 2 and IsEvenInt(P[1]),
	allChildren  := P -> [[ PRDFT1(P[1]/2, P[2]), PRDFT3(P[1]/2, P[2]) ]],
	rule := (P,C) -> let(N := P[1], k := P[2],
	        RC(L_or_OS(N/2+1, Int(N/4)+1)) *
		DirectSum(C[1], C[2]) *
		Tensor(F(2), I(N/2))
	))
));
RulesFor(IPRDFT1, rec(
    IPRDFT1_CT_Radix2 := rec(
	isApplicable := P -> P[1] > 2 and IsEvenInt(P[1]),
	allChildren  := P -> [[ IPRDFT1(P[1]/2, P[2]), IPRDFT2(P[1]/2, P[2]) ]],
	rule := (P,C) -> let(N := P[1], k := P[2],
		Tensor(F(2), I(N/2)) *
		DirectSum(C[1], C[2]) *
	        RC(OddStride(N/2+1, 2))
	))
));

RulesFor(PRDFT2, rec(
    PRDFT2_CT_Radix2 := rec(
	isApplicable := P -> P[1] > 2 and IsEvenInt(P[1]),
	allChildren  := P -> [[ PRDFT2(P[1]/2, P[2]), PRDFT4(P[1]/2, P[2]) ]],
	rule := (P,C) -> let(N := P[1], k := P[2],
	        RC(OddStride(N/2+1, N/4+1)) *
		DirectSum(C[1], C[2]) *
		Tensor(I(N/2), F(2)) ^ L(N, N/2)
	))
));

# SwitchRules(PRDFT3, [1,3]);
# SwitchRules(PRDFT1, [1,7]);#H3_CasTwid.cost := n -> 4 + (n/2-1)*6;

#DHT1.cost := n -> n + DHT1.cost(n/2) + DHT3.cost(n/2);
#DHT3.cost := n -> H3_CasTwid.cost(n) + DFT.cost(n/2);

#DHT1.cost := n -> n + DHT1.cost(n/2) + H3_CasTwid.cost(n/2) + DFT.cost(n/4);
# DHT1(32) cost = 16*2 (F2) + size 64

# DFT_8 = 56
# DFT_16 = 168
# DFT_32 = 456
# DFT_64 = 1160


text:
#############################################################################
##
#A  util.gi                 toric library          David Joyner
##
##
##  This file contains miscellaneous functions
##
##

TORIC:=rec();

##
##  returns 1 if i=j, 0 otherwise
##
TORIC.krondecker_delta:=function(i, j)
if i = j then return 1; fi;
return 0;
end;

##
##  takes a lists of lists and returns the elements
##
TORIC.flatten:=function(L)
local L0,L1,i;
 L0:=ShallowCopy(L[1]);
 if Size(L) = 1 then return L0; fi;
 for i in [2..Size(L)] do
   Append(L0, L[i]);
 od;
 return Set(L0);
end;

##
##  the ``inner product'' v1^t*A*v2
##
TORIC.inner_product:=function(v1,v2,A) return v1*A*v2; end;

##
##  returns a list of dim-vectors of length dim
##
TORIC.make_basis:=function(dim)
local e,i,j;
 e:=[];
 for i in [1..dim] do
  e[i]:=List(TORIC.kronecker_delta(i,j), j in [1..dim]);
 od;
end;

##
##  Input: L is a list of vectors L:=[v1,...,vn]
##  Output: max abs value of all coords of all vecs vi
##
TORIC.max_vectors:=function(L)
 local n,max0,min0,i,max;
  n:=Size(L);
  max:=0;
  for i in [1..n] do
   max0:=Maximum(L[i]);
   min0:=Minimum(L[i]);
   max:=Maximum([max,AbsoluteValue(max0),AbsoluteValue(min0)]);
  od;
 return max;
end;

##
##  for lists L1, L2, returns Boolean Is L1 subset L2?
##
TORIC.is_sublist:=function(L1, L2)
 local i,j,ans,x,y;
  x:=Set(L1);
  y:=Set(L2);
  if IsSubset(y,x) then return true; fi;
 return false;
end;

##
##      Input: L is a list of n-1 vectors in Q^n
##      Output: a vector v s.t. <v,h>=0, for all h in span(L)
##
TORIC.normal_to_hyperplane:=function(L)
 local normal;
  normal:=NullspaceMat(TransposedMat(L));
 return normal[1];
end;

##
##      Input: L is a list of r vectors
##      Output: true if all vectors have same length
##              false otherwise
##
TORIC.consistent_vectors:=function(L)
 local a,b,i,n;
 n:=Length(L);
 b:=Maximum(List([1..n],i->Length(L[i])));
 a:=Minimum(List([1..n],i->Length(L[i])));
 if a<>b then
  Error("\n The vectors must all have the same length.\n");
 fi;
end;

max_stars_repo_name: jamesjer/toric

max_issues_repo_name: jamesjer/toric

max_forks_repo_name: jamesjer/toric


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

ImportAll(paradigms.smp); # for AParSMP
ImportAll(paradigms.vector); # for VRCLR
Import(approx); # for CeilingRat.

_swrap := (nt) -> let(t:=nt.firstTag(), ScratchWrap(t.size,t.nsgmts,t.linesize));

_nodeFitsTensorDoesnt := (nt) -> let (
    k := nt.getTag(1).size,
    m := nt.params[1].dims()[2],
    n := nt.params[2],

    2*m <= k and m*n >= k
);

_isALStoreTag := (nt) -> nt.isTag(1,ALStore) or nt.isTag(1,ALStoreCx);

NewRulesFor(TTensorI, rec(
#   (In x Am)_{LS(<=k)} -> In x (Am)_{LS(<=k)} for m > k
    IxA_scratch_push := rec(
        forTransposition := false,
        applicable := nt -> nt.hasTags()
            and _isALStoreTag(nt)
            and IsParPar(nt.params)
            and nt.params[1].dims()[2] > nt.getTag(1).size,

        children := nt -> [[ nt.params[1].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> Tensor(I(nt.params[2]), c[1])
    ),

    AxI_scratch_push := rec(
        forTransposition := false,
        applicable := nt -> nt.hasTags()
            and _isALStoreTag(nt)
            and IsVecVec(nt.params)
            and nt.params[1].dims()[2] > nt.getTag(1).size,
        children := nt -> [[nt.params[1].withTags(nt.getTags())]],
        apply := (nt, c, cnt) -> Tensor(c[1], I(nt.params[2]))
    ),

#   ==========================================================
# synchronous rules

#   (In x Am)_{LS(<=k)} -> Imn/k x (Ik/m x Am) for m|k and k|mn
   IxA_scratch := rec(
        forTransposition := false,
        applicable := nt -> nt.hasTags()
            and _isALStoreTag(nt)
	        and IsParPar(nt.params)
            and IsPosInt(nt.getTag(1).size / nt.params[1].dims()[2])
            and IsPosInt(nt.params[1].dims()[2] * nt.params[2] / nt.getTag(1).size),

        children := nt -> [[
            TTensorI(
                nt.params[1],
                nt.getTag(1).size / nt.params[1].dims()[2],
                APar, APar
            ).withTags(Drop(nt.getTags(),1)).setWrap(
                _swrap(nt)
            )
        ]],
        apply := (nt, c, cnt) -> let(k := nt.getTag(1).size,  m := nt.params[1].dims()[2], n:= nt.params[2],
                DMAFence(Tensor(
                    I(m*n/k),
                    LSKernel(c[1], nt.params[1].normalizedArithCost() * k/m)
                )))
    ),

#   (Am x In)_{LS(<=k)} -> (L^m^2n/k_m x Ik/m) (Imn/k x (Am x Ik/m)) (L^m^2n/k_mn/k x Ik/m) for m|k and k|mn
    AxI_scratch := rec(
        forTransposition := false,
        applicable := nt -> nt.hasTags()
                            and _isALStoreTag(nt)
	                        and IsVecVec(nt.params)
                            and IsPosInt(nt.getTag(1).size / nt.params[1].dims()[2])
                            and IsPosInt(nt.params[1].dims()[2]*nt.params[2] / nt.getTag(1).size),

        children := nt -> [[ let(t := nt.firstTag(),
            TTensorI(
                nt.params[1],
                t.size / nt.params[1].dims()[2],
                AVec, AVec
            ).withTags(Drop(nt.getTags(),1)).setWrap(_swrap(nt))
        ) ]],

        apply := (nt, c, cnt) -> let(
            k := nt.getTag(1).size,
            m := nt.params[1].dims()[2],
            n:= nt.params[2],

            DMAFence(Prm(fTensor(L(m^2*n/k, m), fId(k/m))) *
            Tensor(
                I(m*n/k),
                LSKernel(c[1], nt.params[1].normalizedArithCost() * k/m)
            ) *
            Prm(fTensor(L(m^2*n/k, m*n/k), fId(k/m))))
        )
    ),

    IxAL_scratch := rec(
        forTransposition := false,
        applicable := nt -> nt.hasTags()
                            and _isALStoreTag(nt)
                            and IsParVec(nt.params)
                            and IsPosInt(nt.getTag(1).size/ nt.params[1].dims()[2])
                            and IsPosInt(nt.params[1].dims()[2]*nt.params[2]/nt.getTag(1).size),
        children := nt -> let(
            k := nt.firstTag().size,
            m := nt.params[1].dims()[2],

            [[ TTensorI(
                nt.params[1],
                k/m,
                APar,AVec
            ).withTags(Drop(nt.getTags(),1)).setWrap(_swrap(nt))
        ]]),

        apply := (self, nt, c, cnt) >> let(
            k := nt.getTag(1).size,
            m := nt.params[1].dims()[2],
            n := nt.params[2],

            DMAFence(
                Tensor(
                    I(m*n/k),
                    LSKernel(c[1],nt.normalizedArithCost() * k/m)
                ) *
                Prm(fTensor(L(m^2*n/k,m*n/k),fId(k/m)))
            )
        ),
    ),

    L_IxA_scratch :=  rec(
        forTransposition := false,
        applicable := nt -> nt.hasTags()
            and _isALStoreTag(nt)
	        and IsVecPar(nt.params),

        children := nt -> [[
            TTensorI(
                nt.params[1],
                nt.getTag(1).size / nt.params[1].dims()[2],
                AVec, APar
            ).withTags(Drop(nt.getTags(),1)).setWrap(
                _swrap(nt)
            )
        ]],

        apply := (self, nt, c, cnt) >> let(
            k := nt.getTag(1).size,
            m := nt.params[1].dims()[2],
            n:= nt.params[2],

            DMAFence(
                Prm(fTensor(L(m^2*n/k, m*n/k), fId(k/m))
                ) * Tensor(LSKernel(
                    c[1], nt.normalizedArithCost() * k/m
                ), I(m*n/k))
            )
        )
     ),
#====================================================================
# SWP Rules - double buffering is applied
   IxA_scratch_swp := rec(
        forTransposition := false,
        applicable := nt -> nt.hasTags()
            and _isALStoreTag(nt)
	        and IsParPar(nt.params)
            and IsPosInt(nt.getTag(1).size / (2 * nt.params[1].dims()[2]))
            and IsPosInt(2 * nt.params[1].dims()[2] * nt.params[2] / nt.getTag(1).size),

        children := nt -> [[
            TTensorI(
                nt.params[1],
                nt.getTag(1).size / (2 * nt.params[1].dims()[2]),
                APar, APar
            ).withTags(Drop(nt.getTags(),1)).setWrap(
                _swrap(nt)
            )
        ]],
        apply := (nt, c, cnt) -> let(k := nt.getTag(1).size,  m := nt.params[1].dims()[2], n:= nt.params[2],
                DMAFence(Tensor(
                    I(2*m*n/k),
                    LSKernel(c[1], nt.params[1].normalizedArithCost() * k/(2*m))
                )))
    ),

    AxI_scratch_swp := rec(
        forTransposition := false,
        applicable := nt -> nt.hasTags()
                            and _isALStoreTag(nt)
	                        and IsVecVec(nt.params)
                            and IsPosInt(nt.getTag(1).size / (2*nt.params[1].dims()[2]))
                            and IsPosInt(2*nt.params[1].dims()[2]*nt.params[2] / nt.getTag(1).size),

        children := nt -> [[ let(t := nt.firstTag(),
            TTensorI(
                nt.params[1],
                t.size / (2*nt.params[1].dims()[2]),
                AVec, AVec
            ).withTags(Drop(nt.getTags(),1)).setWrap(_swrap(nt))
        ) ]],

        apply := (nt, c, cnt) -> let(
            k := nt.getTag(1).size,
            m := nt.params[1].dims()[2],
            n:= nt.params[2],

            DMAFence(Prm(fTensor(L(2*m^2*n/k, 2*m), fId(k/(2*m)))) *
            Tensor(
                I(2*m*n/k),
                LSKernel(c[1], nt.params[1].normalizedArithCost() * k/(2*m))
            ) *
            Prm(fTensor(L(2*m^2*n/k, 2*m*n/k), fId(k/(2*m)))))
        )
    ),

    IxAL_scratch_swp := rec(
        forTransposition := false,
        applicable := nt -> nt.hasTags()
                            and _isALStoreTag(nt)
                            and IsParVec(nt.params)
                            and IsPosInt(nt.getTag(1).size/ (2*nt.params[1].dims()[2]))
                            and IsPosInt(2*nt.params[1].dims()[2]*nt.params[2]/nt.getTag(1).size),
        children := nt -> let(
            k := nt.firstTag().size,
            m := nt.params[1].dims()[2],

            [[ TTensorI(
                nt.params[1],
                k/(2*m),
                APar,AVec
            ).withTags(Drop(nt.getTags(),1)).setWrap(_swrap(nt))
        ]]),

        apply := (self, nt, c, cnt) >> let(
            k := nt.getTag(1).size,
            m := nt.params[1].dims()[2],
            n := nt.params[2],

            DMAFence(
                Tensor(
                    I(2*m*n/k),
                    LSKernel(c[1],nt.normalizedArithCost() * k/(2*m))
                ) *
                Prm(fTensor(L(2*m^2*n/k,2*m*n/k),fId(k/(2*m))))
            )
        ),
    ),

    L_IxA_scratch_swp :=  rec(
        forTransposition := false,
        applicable := nt -> nt.hasTags()
            and _isALStoreTag(nt)
	        and IsVecPar(nt.params)
            and IsPosInt(nt.getTag(1).size/ (2*nt.params[1].dims()[2]))
            and IsPosInt(2*nt.params[1].dims()[2]*nt.params[2]/nt.getTag(1).size),

        children := nt -> [[
            TTensorI(
                nt.params[1],
                nt.getTag(1).size / (2*nt.params[1].dims()[2]),
                AVec, APar
            ).withTags(Drop(nt.getTags(),1)).setWrap(
                _swrap(nt)
            )
        ]],

        apply := (self, nt, c, cnt) >> let(
            k := nt.getTag(1).size,
            m := nt.params[1].dims()[2],
            n:= nt.params[2],

            DMAFence(
                Prm(fTensor(L(2*m^2*n/k, 2*m*n/k), fId(k/(2*m)))
                ) * Tensor(LSKernel(
                    c[1], nt.normalizedArithCost() * k/(2*m)
                ), I(2*m*n/k))
            )
        )
     )
));