climatrep.mac

/*******************************
Clifford algebra
a lightweight package for performing Geometric Algebra calculations

Clifford Algebra Matrix representation code + FVS inverse algorithm

version 	
			2.5.7. 3 Jun 2023
			- fadlevicg2cp
			- fadlevicg1cp
			2.5.6 10 May 2022
			- print formatting
			2.1.5 5 Dec 2021
			- multtrace
			2.1 10 Jul 2021
			2.0 29 Jan 2021
			1.0 11 May 2016

@depends 'clifford

**********************************
 * @license This library is free software	 you can redistribute it and/or
 *      modify it under the terms of the GNU Lesser General Public
 *      License as published by the Free Software Foundation	 either
 *      version 2.1 of the License, or (at your option) any later version.
 *
 *      This library is distributed in the hope that it will be useful,
 *      but WITHOUT ANY WARRANTY	 without even the implied warranty of
 *      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *       Lesser General Public License for more details.
 *
 *      You should have received a copy of the GNU Lesser General Public
 *      License along with this library	 if not, write to the Free Software
 *      Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */


if get('clifford,'version)=false then (	
	err:errcatch(
		load("clifford.mac")
	)
);

/* move to clifford ?*/
if get('climatrep,'version)=false then (
	remove("@", operator),
	"inner product",
	infix ("@", 130, 129),
	"@"(a, b):=  scprod(a,b),
	print("warning: redefining @")

);

/**
 multiplication table, trace
*/
multtrace():=block([s:[1], lst],
	lst:elements(all),
	for i:1 thru length(lst) do 
		s:endcons(dotsimpc(lst[i].lst[i]), s),
	s
);

/*
diagnoal transposition of a matrix
*/
dtranspose(M):= map(reverse, M);

diag(MM):=block([n],
	if not matrixp(MM) then error("matrix required"),
	n:length(MM),
	genmatrix( lambda([i,j],  if i=j then MM[i,j] else 0 ), n)
);

adiag(MM):=block([n],
	if not matrixp(MM) then error("matrix required"),
	n:length(MM),
	genmatrix( lambda([i,j],  if i=j then MM[n-i +1, j] else 0 ), n)
);

diaglist(MM):=block([n, LL],
	if not matrixp(MM) then error("matrix required"),
	n:length(MM),
	local(LL),
	LL:genmatrix( lambda([i,j],  if i=j then MM[i,j] else 0 ), n),
	list_matrix_entries(matrix(makelist(1, n)).LL)
);

mdiag(MM):=block([n],
	if not listp(MM) then error("list required"),
	n:length(MM),
	genmatrix( lambda([i,j],  if i=j then MM[i] else 0 ), n)
);
 
 
diaglist(M):=block([n],
	if not matrixp(M) then error("matrix required"),
	n:length(M),
	makelist(M[i,i], i, n)
);


elements([cpl]):=block([c:1, lst:[], u, v],
	for i:1 thru ndim do c:c.(1+ asymbol[i] ),
	if emptyp(cpl) then cpl:0
	else cpl:first(cpl),
	c:expand(c),
	lst:inargs(c),
	if cpl=0 then pop(lst)
	else if cpl=1 then
		lst:sublist(lst, lambda([u], maxgrade(u)<2^(ndim/2-1)))
	else if cpl=2 then
		lst:sublist(lst, lambda([u], maxgrade(u)>=2^(ndim/2-1))),
	sort(lst, lambda([u,v], 
			ordergreatp(u,v) and if not (numberp(u) or numberp(v)) then length(u)< length(v)
			)
	)
);



climatrep3(vv):=block([q, lst: elements(all), lsta, lstb, nn, le],
	nn: ceiling(ndim/2),
 
	[lsta, lstb]:oppart(lst, lambda([q], maxgrade(q)<nn)),
	if oddp(length(lstb)) then  (
		le:pop(lstb),
		if vv=le then vv: -%i*clihdual(le),
		if _debug1=true then
			display(lsta, lstb, vv),
		climatrep21(vv, lstb)
	) else
		climatrep21(vv, lstb)
);


/*
computes matrix representation of the algebra;
experimental code
*/
climatrep(vv):=block([n, q, AA, lst: elements(all), nn, dd, G, EE],
	local(AA, G, EE, q),
	n:length(lst),
	/* multiplication table of the algebra */
	AA:genmatrix( lambda([i,j], dotsimpc(lst[i] . lst[j] ) ), n),
	/* signature of the algebra */
	G:diag(AA),
	ratvars(q),
	nn:expand(ratnum(vv)), dd:ratdenom(vv),
	display(nn, dd),
	if dd #1 then
		EE:matrixmap(lambda([q],  scalarpart(ratcoeff (q, nn))/dd), AA)
	else 
		EE:matrixmap(lambda([q],  scalarpart(ratcoeff (q, nn))), AA),
	G.EE
);

climatrep2(vv, lst):=block([n, q, AA, nn, dd, G, EE],
	local(AA, G, EE, q),
	n:length(lst),
	/* multiplication table of the algebra */
	AA:genmatrix( lambda([i,j], dotsimpc(lst[i] . lst[j] ) ), n),
	/* signature of the algebra */
	G:diag(AA),
	ratvars(q),
	nn:expand(ratnum(vv)), dd:ratdenom(vv),
	display(nn, dd, AA),
	EE:matrixmap(lambda([q],  scalarpart(ratcoeff (q, nn))/dd), AA),
	display(EE, G),
	G.EE
);


climatrep21(vv, lst, [ex]):=block([n, q, AA, nn, dd, G, EE, vi, dcase:1],
	local(AA, G, EE, q),
	n:length(lst),
	/* multiplication table of the algebra */
	AA:genmatrix( lambda([i,j], dualrep(dotsimpc(lst[i] . lst[j] )) ), n),
	if not emptyp(ex) then ex:first(ex),
	
	if maxgrade(vv)>ndim/2 and vv#ex then
		vv: dualrep(vv),
	vi:ratcoeff(vv, %i),
	
	if vi#0 then ( 
		vv:vi, 
		dcase:-%i, 
		if _debug1=true then
		display(vv, vi)
	),
	/* signature of the algebra */
	G:diag(AA),
	ratvars(q),
	nn:expand(ratnum(vv)), dd:ratdenom(vv),
	if _debug1=true then
		display( nn, dd, AA),
	EE:matrixmap(lambda([q],   scalarpart(ratcoeff (q, nn))/dd), AA),
	if _debug1=true then
		display(EE, G),
	ratsimp(G.EE*dcase)
);

climatrep21(vv, lst ):=block([n, q, AA, nn, dd, G, EE, vi, l, r, dcase:1],
	local(AA, G, EE, q),
	n:length(lst),
	/* multiplication table of the algebra */
	AA:genmatrix( lambda([i,j], dualrep(dotsimpc(lst[i] . lst[j] )) ), n),
  	if maxgrade(vv)>ndim/2   then
		vv: dualrep(vv),
	vi:ratcoeff(vv, %i),
	
	if vi#0 then ( 
		vv:vi, 
		dcase:-%i, 
		if _debug1=true then
		display(vv, vi)
	),
	[l, r]: subst(nil=1, oppart(vv, lambda([u], freeof(asymbol, u)))),
	/* signature of the algebra */
	G:diag(AA),
	ratvars(q),
	nn:expand(ratnum(r)), dd:ratdenom(r),
	if _debug1=true then
		display( nn, dd, AA),
	EE:matrixmap(lambda([q],   scalarpart(ratcoeff ( q, nn))/dd), AA),
	if _debug1=true then
		display(EE, G),
	ratsimp(l*G.EE*dcase)
); 

dualrep(ee):=block([la, lb, u ],
	if maxgrade(ee) <ndim/2 then ee 
	else(
		if listp(ee) then
			map(dualrep, ee)
		elseif inop(ee)="." then
			%i*clihdual(ee)
		elseif inop(ee)="*" then (
			[la, lb]:subst(nil=1, oppart(ee, lambda([u], freeof(asymbol, u)  ))),
			la*%i*clihdual(lb)
		)else (
			[la, lb]:subst(nil=0, oppart(ee, lambda([u], maxgrade(u)<ndim/2 ))),
			la+%i*clihdual(lb)
		)
	)
);


nelements([cpl]):=block([c:1, lst:[] ],
	if emptyp(cpl) then ( 
		for i:1 thru ndim do c:c.(1+ asymbol[i] ),
		c:expand(c),
		c:dotsimpc(c),
		lst:inargs(c)
	)else (
		cpl:cpl[1],
		for i:1 thru length(cpl) do (
			c:c  & cpl[i] ,
			push(c, lst)
		)
	),
	lst:sort(lst, lambda([u,v], 
				ordergreatp(u,v) and if not (numberp(u) or numberp(v)) then length(u)< length(v)
				)
			),
	lst
);

climatrep1(vv):=block([n, AA, lst, EE, l,r, q, opsubst :false, gs:gensym(string(asymbol)), MM ],
	local(AA, EE),
	if emptyp(%elements) then lst:elements(all)
	elseif %elements[1] #1 then lst: cons(1, %elements) else lst: %elements,	
	n:length(lst),
	/* special case */
	/* multiplication table of the algebra */
	AA:genmatrix( lambda([i,j], dotsimpc(lst[i] . lst[j] ) ), n),
	AA:subst( asymbol[gs]=asymbol[1], subst(1=gs, AA)), 
	vv:subst( asymbol[gs]=asymbol[1], subst(1=gs, vv)), 
	/*display(vv, gs),*/
	[l, r]: oppart(vv, lambda([u], freeof(asymbol, u))),
	/*display(AA),*/
	if l='nil then l:1,
	if r='nil then r:gs,
	display(l,r),
	EE:matrixmap(lambda([q], ratcoeff (q, r)), AA),
	display(EE),
	/*  twiddle to get the signs right*/
	MM:genmatrix( lambda([i,j],   dotsimpc( lst[i] . lst[i]) * EE[i,j]), n),
	subst( gs=1,l*MM)
);


/* matrix FVS algorithm */
fadlevim(A):=block([n, m, c, M, K, I, i, cq, Ai],
	if not matrixp(A) then error ("not a matrix ", A),
	[n,m]:matrix_size(A),
	if n#m then error ("not a square matrix"),
	c: makelist(1,n+1),
	M:ident(n),
	I:ident(n),
	A:rat(A),
	for i:1 thru n do (
		K: rat(A.M),
		c[i+1]: -1/i* mtrace(K),
		M:(K + c[i+1]*I),
		if i=n-1 then  Ai:M
	), 
	c:factor(c),
	cq:last(c),
	if cq=0 then cq:1,
	display(M),
	M:rat(-Ai/cq),
	[ M, c]
);

/* 
matrix version
characteristic polynomial Faddeev -Le Verrier- Sourieau 
*/
fadlevcpm(A, v):=block([n, m, c, M, K, I, ss:0], 
	local(c),
	if not matrixp(A) then error ("not a matrix ", A),
	[n,m]:matrix_size(A),
	if n#m then error ("not a square matrix"),
	array(c,n+1), 
	for r:0 thru n+1 do c[r]:1,
	ratvars(v),
	A:rat(A),
	M:ident(n),
	I:ident(n),
	ss:c[1]*v^n,
	for i:1 thru n do (
		K: A.M,
		c[i+1]: -1/i* mtrace(K),
		M:(K + c[i+1]*I),
		ss:ss+c[i+1]*v^(n-i)
	),
	/*[ ss , M]*/
	factor(ss)
);

/*
Generalized inverse
naive implementation full  dimensions
*/
fadlevic(A):=block([c, M:1, K, i, n:2^ndim, cq],
	c: makelist(1,n+1),
	for i:1 thru n do (
		K: dotsimpc(expand(A.M)),
		c[i+1]: -n/i* scalarpart(K),
		if i<n then  
			M:ratsimp(K + c[i+1])
	), 
	c:factor(c),
	cq:last(c),
	if cq=0 then cq:1,
	M:ratsimp(-(M)/cq),
   [M, c]
);	


/*
Generalized inverse + characteristic polynomial
 full grade  dimensions
*/
fadlevicg(A, v):=block([c, M:1, K, i, n, cq,  n: 2^ndim, ss:0],
	if _debug1=true then display(n),
	array(c,n), 
	c[0]:1,
	A:rat(A),
	ratvar(v),
	for i:1 thru n do (
		K:  (dotsimpc(expand (A.M))),
		c[i]: -n/i* scalarpart(K),
		if _debug1=true then
			print(i, " ", c[i]),
		/*print(i, " ", c[i], " ", K),*/
		if i<n then  
			M:rat(K + c[i])
	), 
	c:factor(c),
	cq: c[n],
	if cq=0 then (
		if _debug1=true then print("determinant ", cq),
		M:0
	) else
		M:factor(-(M)/cq),
	ss:sum(c[n-i]*v^(i),i,0, n),
	ss:factor(ss),
	[M, ss]
);

fadlevicp(A, v):=block([ M:1, K, i:1, n, k, cq, c, ss],
	n: 2^ceiling(ndim/2),
	if _debug1=true then print("n=",n),
	array(c,n+1), 
	for r:0 thru n+1 do c[r]:1,
	A:rat(A),
	ss:c[1]*v^^n,
	while i<n and K#0 do (
		K: dotsimpc(expand (A.M)),
		cq: -n/i* scalarpart(K),
		if _debug1=true then print(i, " ", cq),
		if _debug1=all then print("t_{",i,"}=", cq, " m_{",i,"}=", (K+cq), "\\\\"),
		if K#0 then
			M: rat(K + cq), 
		c[i+1]: cq,
		ss:ss+c[i+1]*v^^(n-i),
		i:i+1
	), 
	 
	K: dotsimpc(expand(A.M)),
	cq: -n/i* scalarpart(K),
	if _debug1=all then print("t_{",i,"}=", cq, " m_{",i,"}=", (K+cq), "\\\\"),
	ss:ss+cq,
	if cq=0 then cq:1,
	M:factor(-(M)/cq),
	[M, ss]
);

/*
Generalized inverse 
 reduced dimensions
*/
fadlevicg2(A):=block([ M:1, K, i:1, n, sc:1, cq],
	n: 2^ceiling(ndim/2),
	if _debug1=all then print("n=",n),
	A:rat(A),
	while i<n and K#0 do (
		K: dotsimpc(expand (A.M)),
		cq: n/i* scalarpart(K),
		if _debug1=true then
			print(i, " ", cq),
		if _debug1=all then
			print("t[",i,"]=", cq, " m[",i,"]=", K),
		if K#0 then
			M: rat(K - cq), 
		i:i+1
	), 
	 
	K: dotsimpc(expand(A.M)),
	cq: -n/i* scalarpart(K),
	if _debug1=all then print("t_{",i,"}=", cq, " m_{",i,"}=", (K-cq), "\\\\"),
	if cq=0 then cq:1,
	M:factor(-(M)/cq)
);


/* 
pseudoinverse : for n idempotent it returns the complement zero divisor 
full grade dimensions
*/
fadlevicg1(A):=block([ M:1, K, i:1, n,  k: length(clv(A)),  cq],
	n: 2^k,
	if _debug1=true then print("n=",n),
	A:rat(A),
	while i<n and K#0 do (
		K: dotsimpc(expand (A.M)),
		cq: -n/i* scalarpart(K),
		if _debug1=true then print(i, " ", cq),
		if _debug1=all then print("t_{",i,"}=", cq, " m_{",i,"}=", K, "\\\\"),
		if K#0 then
			M: rat(K + cq), 
		i:i+1
	), 
	K: dotsimpc(expand(A.M)),
	cq: -n/i* scalarpart(K),
	if _debug1=all then print("t_{",i,"}=", cq, " m_{",i,"}=", (K+cq), "\\\\"),
	cq:factor(cq),
	if cq=0 then cq:1,
	M:factor(-(M)/cq)
	
);


/* inverse +  full grade, characteristic polynomial  FVS */
fadlevicg1cp(A, v):=block([ M:1, K, i:1, n, k: length(clv(A)),  cq, c, ss],
	n: 2^k,
	if _debug1=true then print("n=",n),
	array(c,n+1), 
	for r:0 thru n+1 do c[r]:1,
	A:rat(A),
	ss:c[1]*v^^n,
	while i<n and K#0 do (
		K: dotsimpc(expand (A.M)),
		cq: -n/i* scalarpart(K),
		if _debug1=true then print(i, " ", cq),
		if _debug1=all then print("t_{",i,"}=", cq, " m_{",i,"}=", (K+cq), "\\\\"),
		if K#0 then
			M: rat(K + cq), 
		c[i+1]: cq,
		ss:ss+c[i+1]*v^^(n-i),
		i:i+1
	), 
	 
	K: dotsimpc(expand(A.M)),
	cq: -n/i* scalarpart(K),
	if _debug1=all then print("t_{",i,"}=", cq, " m_{",i,"}=", (K+cq), "\\\\"),
	ss:ss+cq,
	if cq=0 then cq:1,
	M:factor(-(M)/cq),
	[M, ss]
);

/*
minimal polynomial of the multivector
*/
minpoly(P, x):=block([dP:rat(P), md:0, %a, Q, Q1, n, k, cc:[], sol,  m, qq],
	local(%a),
	n:hipow(P,x),
	if _debug1=all then  display ( n),
	while gcd(dP, rat(P))#1 do (
		dP:diff(dP,x),
		md:md+1
	),
	if md=1 then (
	 print (P, "is minimal"),
	 return(P)
	),
	if _debug1=all then  display (md),
	m: n/md,
	if not integerp (m) then return (false),
	Q:sum(%a[i]*v^i,i,0, m),

	Q1:expand(P-Q^(md)),
	 
	for k:n thru 1 step -1 do (
		qq:ratcoeff(Q1, x, k),
		qq: subst(cc, qq),
		sol:solve(qq, %a[m]),
		if _debug1=all then display(qq, k, m, sol),
		if listp(sol) then
				cc: push(last(sol), cc)
		else return( [ subst(cc,Q), cc] ),
		m:m-1,
		if m <0 then return( [ subst(cc,Q), cc] )
		),
	[ subst(cc, Q), cc] 
);


/*
Generalized inverse + characteristic polynomial
 reduced grade dimensions
*/
fadlevicg2cp(A, v):=block([ M:1, K, i:1, n, k: length(clv(A)),  cq, c, ss],
	n: 2^ ceiling( k/2),
	if _debug1=all then print("n=",n),
	array(c,n+1), 
	for r:0 thru n+1 do c[r]:1,
	A:rat(A),
	ss:c[1]*v^^n,
	while i<n and K#0 do (
		K: dotsimpc(expand (A.M)),
		cq: -n/i* scalarpart(K),
		if _debug1=true then print(i, " ", cq),
		if _debug1=all then print("t_{",i,"}=", cq, " m_{",i,"}=", K, "\\\\"),
		if K#0 then
			M: rat(K + cq), 
		c[i+1]: cq,
		ss:ss+c[i+1]*v^^(n-i),
		i:i+1
	), 
	 
	K: dotsimpc(expand(A.M)),
	cq: -n/i* scalarpart(K),
	if _debug1=all then print("t_{",i,"}=", cq, " m_{",i,"}=", (K+cq), "\\\\"),
	ss:ss+cq,
	if cq=0 then cq:1,
	M:factor(-(M)/cq),
	[M, ss]
);

fadlevicg3cp(A, v, n):=block([ M:1, K, i:1,   cq, c, ss],
	if not integerp(n) then error("n is not integer ", n),
	array(c,n+1), 
	for r:0 thru n+1 do c[r]:1,
	A:rat(A),
	ss:c[1]*v^^n,
	while i<n and K#0 do (
		K: dotsimpc(expand (A.M)),
		cq: -n/i* scalarpart(K),
		if _debug1=true then print(i, " ", cq),
		if _debug1=all then print("t_{",i,"}=", cq, " m_{",i,"}=", K, "\\\\"),
		if K#0 then
			M: rat(K + cq), 
		c[i+1]: cq,
		ss:ss+c[i+1]*v^^(n-i),
		i:i+1
	), 
	 
	K: dotsimpc(expand(A.M)),
	cq: -n/i* scalarpart(K),
	if _debug1=all then print("t_{",i,"}=", cq, " m_{",i,"}=", (K+cq), "\\\\"),
	ss:ss+cq,
	if cq=0 then cq:1,
	M:factor(-(M)/cq),
	[M, ss]
);

/*
guessgrade(A):=block([ A2:A, k:1, ipdecomp:true],
	while A2#0 do (
		A2: A2 | A2,
		display(A2),
		k:k*2
	),
	k/2
);
*/

/* maximal rank*/
mgr(A):= block([k:length(clv(A))],
	2^ceiling(k/2)
);

mpowers(A):=block([powers:[], m: mgr(A), i:0, Ai:1],
	while i<= m  do (
		push(Ai,powers),
		Ai:dotsimpc(expand( Ai . A)), 
		i:i+1
	),
	reverse(powers)
);

guessgrade2(A):=block([M, Ls:[], u],
	Ls: mpowers(A),
	/*Ls: map(lambda([u], list_matrix_entries(elem2matr(u))) ,Ls),  */
	Ls: map(lambda([u], elem2lst(u)), Ls),  
    M: apply(matrix, Ls),
	display(M),
	rank(M)
);


trimzerocols(M):=block([MT: transpose(M), L:[], m,n, i:1, Z],
	[m,n]:matrix_size(M),
	Z:zeromatrix(1,m),
	while i<=n do (
		if row(MT,i)#Z then
			L:cons( list_matrix_entries(row(MT,i)), L),
		i:i+1
	),
	L:reverse(L),
	transpose(apply(matrix, L))
);

guessgrade21(A):=block([M, Ls:[], u, L],
	Ls: mpowers(A),
	Ls: map(lambda([u], list_matrix_entries(elem2matr(u))) ,Ls),  
    M: apply(matrix, Ls),
	M:trimzerocols(M),
	display(M),
	rank(M)
);

guessgrade22(A):=block([M, Ls:[], u, L],
	Ls: mpowers(A),
	Ls: map(lambda([u], elem2lst(u)) ,Ls),  
    M: apply(matrix, Ls),
	M:trimzerocols(M),
	display(M),
	rank(M)
);

/*
Generalized inverse  float
 reduced grade dimensions
*/
fadlevicg2f(A):=block([ M:1, K, i:1, n, k: length(clv(A)), sc:1, cq],
	n: 2^ ceiling( k/2),
	if _debug1=all then print("n=",n),
	A:float(A),
	while i<n and K#0 do (
		K: dotsimpc(expand (A.M)),
		cq: n/i* scalarpart(K),
		if _debug1=true then
			print(i, " ", cq),
		if _debug1=all then
			print("t[",i,"]=", cq, " m[",i,"]=", K),
		if K#0 then
			M: float(K - cq), 
		i:i+1
	), 
	 
	K: dotsimpc(expand(A.M)),
	cq: -n/i* scalarpart(K),
	if _debug1=all then print("t_{",i,"}=", cq, " m_{",i,"}=", (K+cq), "\\\\"),
	if cq=0 then cq:1,
	M:factor(-(M)/cq)
	
);




elem2mat3(expr, [ulst] ):=block(
	if emptyp(ulst) then ulst:false 
	else ulst:true,
	if mapatom(expr) then return(climatrep3(expr)),
	if ulst=true then
		maplist(climatrep3, expr)
	else
		map(climatrep3, expr)
);

coeffmap(vv):=block([n, AA, lst, EE],
	local(AA,  EE),
	lst: elements(all),
	n:length(lst),
	/* multiplication table of the algebra */
	AA:genmatrix( lambda([i,j], dotsimpc(lst[i] . lst[j] ) ), n),
	if vv=1 then (AA:subst(1=asymbol[0], AA), vv:asymbol[0]),
	EE:matrixmap(lambda([q], ratcoeff (vv, q)), AA)
);

/*
identifies the abstract element from a matrix representation
*/
matelem(MM):=block([LL],
	local(LL),
	if not matrixp(MM) then error("matrix required"),
	LL:matrix(elements(all)),
	LL:MM.LL,
	LL[1,1]
);

matelem1(MM):=block([LL],
	local(LL),
	if not matrixp(MM) then error("matrix required"),
	LL: (matrix(elements(all))),
	LL:row(MM,1).LL
);


matclimin(MM):=block([n, q:0, LL,  z:0, s],
	if not matrixp(MM) then error("matrix required"),
	local(LL),
	n:length(MM),
	while n>1 do (
		n:n/2,
		q:q+1
	),
	LL:diag(MM),
	matrixmap(lambda([q], if q>0 then z:z+1), LL),
	n:length(MM),
	s:n-z,
	while n>q do (
		if s>1 then s:s-1, 
		if z>1 then z:z-1,
		n:s+z
	),
	[z,s]
);

lclimin(lst):=block([n, q:0,  z:0, s],
	if not listp(lst) then error("list required"),
	n:length(lst),
	while n>1 do (
		n:n/2,
		q:q+1
	),
	map(lambda([q], if q>0 then z:z+1), lst),
	n:length(lst),
	s:n-z,
	while n>q do (
		if s>1 then s:s-1, 
		if z>1 then z:z-1,
		n:s+z
	),
	[z,s]
);

/*  checks if M is an orthogonal matrix*/
morthogonalp(MM):=block([n, m, LL  ],
	if not matrixp(MM) then error("matrix required"),
	[m,n]:matrix_size(MM),
	if n#m then error("matrix non square"),
	local(LL),
	LL:transpose(MM),
	LL:MM.LL - ident(n),
	is(zeromatrix (n,n) = LL)
);

idemps():=block([n, ee, ilst, rilst, iilst, ret],
	ee: elements(all),
	n:length(ee),
	ilst:makelist(dotsimpc(ee[i] . ee[i] ), i, n ),
	rilst:(sublist_indices(ilst, lambda([u], u>0))),
	iilst:(sublist_indices(ilst, lambda([u], u<0))),
	ret:makelist((ee[i] +1)/2, i, rilst),
	append(ret, makelist((ee[i]*%i +1)/2, i, iilst))
);

/*
elem2mat(expr):=block([L:[], ei, ret, E:[], qq],
	if emptyp(%elements) then E:elements(all)
	elseif %elements[1] #1 then E: cons(1, %elements) else E: %elements,
	
	ret:expand(expr),
	L:map(lambda([u], 
		qq: u.expr, dotsimpc(expand(qq)) 
		), E),
	genmatrix( lambda([i,j],  scprod( L[i] ,  E[j]) * scprod( E[j] ,  E[j])), length(L))
);
*/

elementsr(A):=block([c:1, lst:clv(A), u, v, n ],
	n: length(lst),
	for i:1 thru n do c:c.(1+ lst[i] ),
	c:expand(c),
	lst:inargs(c),
	sort(lst, lambda([u,v], 
			ordergreatp(u,v) and if not (numberp(u) or numberp(v)) then length(u)< length(v)
			)
	)
);

elem2mat(expr):=block([L:[], ei, ret, E:[], u, qq],
	if emptyp(%elements) then E:elements(all)
	elseif %elements[1] #1 then E: cons(1, %elements) else E: %elements,	
	ret:expand(expr),
	ei:maplist(lambda([u],  dotsimpc(u.u).u  ), E),
	L:map(lambda([u], 
		qq: u.expr, dotsimpc(expand(qq)) 
		), E),
	genmatrix( lambda([i,j],  L[i] @ ei[j] ), length(L))
);


/* computes only the 1st row*/
elem2matr(expr):=block([L:[], ei, ret, E:[], u, qq],
	E:elementsr(A), 
	ret:expand(expr),
	ei:maplist(lambda([u],  dotsimpc(u.u).u  ), E),
	L:map(lambda([u], 
		qq: u.expr, dotsimpc(expand(qq)) 
		), E),
	genmatrix( lambda([i,j],   L[i] @ ei[j] ), 1, length(L))
);

elem2lst(expr):=block([L:[], ei, ret, E:[], u, qq],
	E:elementsr(A), 
	ret:expand(expr),
	ei:maplist(lambda([u],  dotsimpc(u.u).u  ), E),
	 makelist(  expr @ ei[j], j, 1, length(E))
);


permatrixsym(lsta, lstb):= block([n,m],
	if not (listp(lsta) or listp(lstb)) then error("list required"),
	n:length(lsta), m:length(lstb),
	genmatrix( lambda([i,j], coeff(lsta[i], lstb[j] )), n, m)
);

permatrixnum(lsta, lstb):= block([n,m, MM, ss:0],
	if not (listp(lsta) or listp(lstb)) then error("list required"),
	n:length(lsta), m:length(lstb),
	local(MM),
	MM:genmatrix( lambda([i,j], if lsta[i] =lstb[j] then 1 else 0), n, m),
	for i:1 thru n do (
		ss:sum(MM[i, q], q, 1, m),
		if ss#0 then MM[i]:MM[i]/ss
	),
	MM
);

mtrace(MM):=block([m, n],
	if not matrixp(MM) then error("matrix required"),
	n:length(MM), m:length(MM[1]),
	if n#m then error("matrix not square"),
	sum(MM[i,i], i, 1, n)
);

permcheck(lsta, lstb):=block([PP, H, E],
	local(PP, H, E),
	PP:permatrixnum(lsta, lstb),
	/* permutation matrix should be pseudo-orthogonal*/
	H:transpose(PP).PP,
	/* its square is idempotent */
	E:H.H-H,
	[is(mtrace(E)=0), E, H, PP]
);

/*  random stuff*/

lrandomrat(n, m):=block([rs, rn, rd],
 makelist( (rs:random(2), rn:random(m), rd:subst(0=1,random(m)),  
 (if oddp(rs) then -1 else 1))*rn/rd, n) 
);


lrandomrats(n, m):=block([rs, rn, rd],
 makelist( (rs:random(3)-1, rn:random(m), rd:subst(0=1,random(m)),  
 (if oddp(rn) then 0 else 1)*rs*rn/rd), n) 
);

/*
sparse random multivector
*/
rancelem(m):=block([n:2^ndim, lst],
 lst:lrandomrats(n, m),
 celem(lst)
);

celemr( rng):=block([ee, ss:0, n, cc, qq],
	if not numberp(rng) then return (false),
	ee:elements(all),
	n:length(ee),
	rng:float(rng), 
	ss:sum (random(rng)*ee[i], i, 1, n),
	ss
);

put('climatrep, 'v25,'version);
put('climatrep, "Dimiter Prodanov", 'author);
put('climatrep, "(C) - Dimiter Prodanov, 2016 - 2023", 'copyright);


print("package name: climatrep.mac");
print("author: ", get('climatrep,'author));
print("version:", get('climatrep,'version));
print("Recommended location: share/contrib");
print("last update: 01 Aug 2023");