#morscripts 08/11/1999 19:24 
with(linalg):
with(grobner):
readlib(mtaylor):

fixorder:=proc(base,vars)local i,s,base_:
base_:=base:
s:=sum(base_['i'],'i'=1..nops(base_)):
[]:
for i to nops(base_)do
lt(s,vars,plex):
[op(""),"]:
s:=s-"":""
od:
base_:=":
end:

mdc:=
proc()
local i,list_or_set:
  if nargs=1 then list_or_set:=args[1]:
  else list_or_set:={seq(args[i],i=1..nargs)}
  fi:
  list_or_set[1]:
	for i from 2 to nops(list_or_set) do gcd(list_or_set[i],")
	od:
	RETURN("):
end:    #mdc

omit:=
proc()
local i,j,lista,x:
    if nargs<2 or not type(args[1],list)then 
    ERROR(`USAGE omit(list1,list2) or omit(list1,x1,x2,...)`)else
    lista := args[1]:
    if type(args[2],list) or type(args[2],set) then x := args[2]
    else x := [seq(args[i],i = 2 .. nargs)]
    fi:
    j := []:
    for i to nops(lista) do  if not member(lista[i],x) then j := [op(j),lista[i]] fi od:
    RETURN(j)
    fi
end:    

intersec:=proc(f,g,vars) local s_,i,g_,gg_:
    if nargs=2 then
    [s_,op(indets(f)union indets(g))]:
    elif nargs=3 then[s_,op(args[3])]
    fi:
    g_:=gbasis({seq(f[i]*s_,i=1..nops(f)),seq(g[i]*(1-s_),i=1..nops(g))}
	,",plex):
    gg_:={}:
    for i to nops(g_) do
	collect(op(i,g_),s_):
	degree(",s_):
	if "=0 then gg_:=gg_ union{simplify("")} fi:
    od:
RETURN(gg_)
end:

######
sature:=proc(g,t,fn) local i,s_,g_,gg_,g0_,j:
    if type(g,list) or type(g,set) then
	gg_:=convert(g,set):
	g0_:={}:
	if nargs=3 then j:=0:fi:
	while gg_<>g0_ do
	    if nargs=3 then j:=j+1:lprint(`starting round #`,j)fi:
	    g0_:=gg_:
	     {seq(g0_[i]*s_,i=1..nops(")),t*(1-s_)}:
	     [s_,op(indets(g0_)union indets(t))]:
	    g_:=gbasis("",",plex):
	    gg_:={}:
	    for i to nops(g_) do
	       collect(op(i,g_),s_):
	       degree(",s_):
	       if "=0 then gg_:=gg_ union{simplify(""/t)} fi:
	    od:
	    if nargs=3 then
		lprint(`writing output of round #`,j,`to `,fn):
		interface(quiet=true):
		appendto(fn):
		lprint(`#output of round #`,j):
		lprint(gg_,`:`):
		writeto(terminal):
		interface(quiet=false):
	    fi:
	od:
    subs(t=0,gg_):
    lprint(`lim`,t,` ->0 =`):
    lprint("):
    RETURN([gg_,"]):
    else ERROR(`WRONG TYPE`)
    fi:
end:#sature

##

si:=
proc()
local i,j,ro:
    if type(args[1],matrix) then
        ro := eval(args[1]):
        for i to rowdim(ro) do  for j to coldim(ro) do  ro[i,j] := simplify(ro[i,j]) od od
    elif type(args[1],vector) then
        ro := eval(args[1]):
        #print(vector,op(2,op(2,eval(ro)))):
        for i to op(2,op(2,eval(ro))) do  ro[i] := simplify(ro[i]) od
    else ro := simplify(args[1])
    fi:
    RETURN(eval(ro))
end:
##############

choosevar:=proc(a,vars)
##############
indets(a):
if"<>{}then grobner[leadmon](a,vars)[2]fi:
RETURN("):
end:

##############  #17:14  22/05/94

cleanupdown:=proc(m,r,c,rows)
##############
local temp,i,j,test:
   j := copy(m):
   if type(j,matrix) and type(r,integer) and type(c,integer)then
	if c <=coldim(j) and c >=1 then
	    j[r,c]:=si(j[r,c]):
	    if j[r,c] <> 0 then test:=0:
		for i in rows do
		    if i <> r then
			temp := si(j[i,c]/j[r,c]):
			if temp <> 0 then
			  if type(denom(temp),integer)or 
				type(j[r,c],numeric) then
			    j := addrow(j,r,i,-temp):
			  else
				if test=0 then
		      lprint(r,c,`WARNING:ROW `,i,UNCHANGED,denom(temp)):
				test:=1
				else
		      lprint(r,c,`WARNING:ROW `,i,UNCHANGED):
				fi
			  fi
			fi
		    fi
		od
	    fi
	fi
    fi:
    RETURN(eval(j))
end:

########

coefmon:=proc() #f,vars,monomio
########
local n,f0,i,v,monomio:
    if nargs = 2 then monomio := args[2]: v := indets(monomio)
    elif nargs = 3 then monomio := args[3]: v := args[2]
    fi:
    n := nops(v):
    degree(monomio,v[1]):
    f0 := coeff(collect(args[1],v[1]),v[1],"):
    for i from 2 to n while f0 <> 0 do
	degree(monomio,v[i]): f0 := coeff(collect(f0,v[i]),v[i],")
    od:
    RETURN(f0)
end:

###

eqs:=
###
proc(m0,rows_no_pivo_m0,vars)
local w,eqs_0,eqs_0_1,rows_with_eqn_0,position_eqn_m0,i,j,rows_no_pivo_m:
    lprint(`DOING EQS IN`,'m0'):
    eqs_0 := []:
    eqs_0_1 := []:
    rows_with_eqn_0 := {}:
    position_eqn_m0 := {}:
    rows_no_pivo_m := [op(rows_no_pivo_m0)]:
    for i to nops(rows_no_pivo_m) do
    	if i=1 then w:=`1st` elif i=2 then w:=`2nd` elif i=3 then 
		w:=`3rd` elif i>3 then w:=cat(i,`th`)
	fi:
        lprint(`searching row`,rows_no_pivo_m[i],w,` out of `,
            nops(rows_no_pivo_m)):
        for j to coldim(m0) do
            mt(m0[rows_no_pivo_m[i],j],vars):
            if " <> 0 then
                choosevar(",vars):
                if " <> {} then
                    position_eqn_m0 :=
                        [op(position_eqn_m0),[rows_no_pivo_m[i],j]]:
                    "":
                    rows_with_eqn_0 :=
                        {rows_no_pivo_m[i]} union rows_with_eqn_0:
                    "":
                    eqs_0 := [op(eqs_0),m0[rows_no_pivo_m[i],j]]:
                    eqs_0_1 := [op(eqs_0_1),""]:
                    lprint(`found `.(nops(eqs_0_1)).` equations sofar`)
                fi
            fi
        od
    od:
    RETURN([rows_with_eqn_0,position_eqn_m0,expand(eqs_0),eqs_0_1])
end
:#of eqs

#####

find:=proc()    #expr,matrix,row
#####
    local j:
    for j to coldim(args[2])while expand(args[2][args[3],j]-args[1])<>0do
    od:
    if j>coldim(args[2])then ERROR(`NO SUCH ENTRY`)else RETURN(j)fi
end:


##########
local_triv:=proc(m1,rows_with_eqn_0,eqs01)local i,j:
##########
lprint(`DOING LOCAL_TRIV IN`,'m1'):
    for i in rows_with_eqn_0 do
      for j to coldim(m1)do lprint(`SIMPLIFICANDO i,j=`,i,j):
      si(m1[i,j]/eqs01):
#DIVISIBILITY ExPECTED BECAUSE  IDEAL SPANNED BY THAT ROW BECOMES <eqs_01>
	if type(denom("),integer)then m1[i,j]:=" else
	ERROR(`MISTAKE AT `,`m1`,`?`,i,j,eqs01)
	fi
      od:
    od:
0:0:0:
end:

lt:=proc()
if nargs=2 then leadmon(args[1],args[2])
 else leadmon(args[1],args[2],args[3])
fi: "[1]*"[2] 
end:


#####
mon:=
#####
proc()
local k,d,vars,degs:
  if nargs<2 then ERROR(`mon(d,vars), mon(d,vars,degs)`)else
    d:=args[1]:vars:=args[2]:
    if nargs=2 then degs:=seq(1,k=1..nops(vars))else degs:=args[3]fi:
    product('1/(1-vars[k]*t^degs[k])','k' = 1 .. nops(vars)):
    expand(convert(taylor(",t,d+1),polynom)):
    collect(",t):
    coeff(",t,d): 
    if type(",`+`)then convert(",list):else ["] fi:
    fixorder(",vars):
    RETURN(")
  fi
end:

####
mt:=proc()local vars,i:
####
    if nargs = 1 then RETURN(mtaylor(args[1],indets(args[1]),2))
    elif nargs = 2 then 
      if type(args[2],set)or type(args[2],list)then 
	vars:=args[2]:i:=2 
      else vars:=indets(args[1]):i:=args[2]:
      fi:
      RETURN(mtaylor(args[1],vars,i))
    fi
end:

######
origin:=proc()local i,vars,expr:
######
#lprint(`computing at 0`):
    if nargs=2 then vars:=args[1]:expr:=eval(args[2])
    else expr:=args[1]:vars:=indets(convert(expr,set))
    fi:#print(vars):
RETURN(subs({seq(vars[i]=0,i=1..nops(vars))},copy(eval(expr))))
end:

##############
pivo:=              #pivo(matrix,row,vars)   #14/4/95
##############
proc()
local m,j,vars,jj,i:
    m := eval(args[1]):
    i := args[2]:
    if type(m,matrix) and type(i,integer) then
	if nargs = 3 then vars := args[3]
	else vars := indets(convert(m,set))
	fi:
	for j to coldim(m) while si(mtaylor(m[i,j],vars,1)) = 0 do
	od:
	#print(j):
	if coldim(m) < j then j := -1
	elif type(si(m[i,j]),numeric) = false then
	    print(`NONCONST. PIVO, TRY AGAIN`):
	    #vars := j:
	    #j := j+1:
	    for jj from j to coldim(m) while type(si(m[i,jj]),numeric) = false or
		si(m[i,jj]) = 0 do

	    od:
	    if jj <= coldim(m) and si(m[i,jj]) <> 0 then
		j:=jj:
		print(`2nd trial successful`)
	    else
		lprint(`#WARNING :NONCONSTANT PIVO`,[i,j],"):
		lprint(`#CHECK HONESTY OF NONCONSTANT PIVO`,[i,j],"):
	    fi
	fi:
	RETURN(j)
    fi
end:    #of pivo

######

prep:=proc()    #m0,position_eqn_m0,eqs_0,eqs_0_1,basis,den #,ver
######
local i,is,j,base,eqs_,m,i1,j1,q,d,M:
    if type(args[1],string)then
    M:=convert(args[1],name) else M:=`matrix`fi:
    lprint(`PREP EQS IN`,M):
    if nargs<5 then 
ERROR(`USAGE prep(m0,position_eqn_m0,eqs_0,eqs_0_1,basis)`,
`OR prep(m0,position_eqn_m0,eqs_0,eqs_0_1,basis,den)`, 
`OR prep(m0,position_eqn_m0,eqs_0,eqs_0_1,basis,den,ver)`
)
    else
    m:=eval(args[1]):
    eqs_ := args[3]:
    base := args[5]:
    j1:=0:
    for i to nops(args[3]) do   
	is:=collect(eqs_[i],args[4][i]):
	d:=coeff(is,args[4][i],1):
	if degree(is,args[4][i])=1 and 
		(type(d,numeric)
		or
(nargs>=6 and
{sqrfree(numer(d))[2][1][1],sqrfree(denom(d))[2][1][1]}minus
{1,args[6]}={}
)               )
		then
		#is:=coeff(eqs_[i],args[4][i],1):
lprint(d,`CLEAR `,args[4][i],` OFF OTHER EQS`):
	for j to nops(args[3]) do
	    if j <> i then
		rem(eqs_[j],eqs_[i],args[4][i],'q'):
		if q<>0 then lprint(`PREP'S CHANGING `,
			M,`AND BASIS`):
		j1:=1:  
		m := addcol(m,args[2][i][2],args[2][j][2],-q):
		[seq(base[i1],i1=1..(args[2][i][2]-1)),
			base[args[2][i][2]]+
				q*base[args[2][j][2]],
		  seq(base[i1],i1=args[2][i][2]+1..nops(base))]:        
		base :=":q:='q':
		eqs_ := [seq(eqs_[i1],i1=1..j-1),
				si(m[op(args[2][j])]),
			 seq(eqs_[i1],i1=j+1..nops(eqs_))
			]:
		fi:#  if "<>0 
	    fi:    #if j <> i
	od:
	elif type(d,numeric)=false and nargs<6 then     
	print(`EQN `,i ,` IS NOT MONIC WRT `,args[4][i],d,nargs)
	#elif
	fi:        #if type(is,numeric)...
    od:         #for i to nops(args[3])
if nargs=8 then
if j1=1 then
lprint(`verificando estrago na base`):
for i1 to rowdim(m)do
sum(m[i1,'j1']*base['j1'],'j1'=1..nops(base)):
sum(args[1][i1,'j1']*args[5]['j1'],'j1'=1..nops(base)):
if si("-"")<>0 then ERROR(`DIVERGES AT i1 = `,i1,`antes=`,",
				`depois=`,"")
fi
od:    
fi:     #if j1=1
fi:     #if nargs=8
    RETURN([expand(eqs_),base,eval(m)])
fi:     #if nargs<5 
end:    #of prep
#prep(mzx,position_eqn_mzx,eqs_zx,eqs_zx_1,z3x2y1,den):"[1]:

#######
relexcn:=
#######
proc(eqs_0,eqs_0_1,newvar,i_)
local i,j,eqs_,eqs1,rel_,others:
    eqs_ := [eqs_0[i_],op(omit(eqs_0,eqs_0[i_]))]:
    eqs1 := [eqs_0_1[i_],op(omit(eqs_0_1,eqs_0_1[i_]))]:
    rel_ := solve(
    {seq(eqs_[i]-newvar[i]*eqs_[1],i = 2 .. nops(eqs_0))},
    {seq(eqs1[i],i = 2 .. nops(eqs_0))}
    ):
    RETURN([rel_,eqs_,eqs1])
end:

   #relexcn




####
ran := proc()   #vars,m
####
local i,vars,m:
    if nargs=2 then vars:=args[1]:m:=eval(args[2])
    else m:=args[1]:vars:=indets(convert(args[1],set))
    fi:if vars<>{}then
    randvector(nops(vars)):
    [seq(vars[i] = "[i],i = 1 .. nops(vars))]:fi:
    subs(",eval(m)):
    RETURN(rank("))
end:

#############
rows_no_pivo:=
#############
proc()
local m0,i,j,pivo_m0,rows_no_pivo_m0,test,vars,rows,M:
    if nargs=0 then 
    ERROR(`usage rows_no_pivo(matrix) or rows_no_pivo(matrix,vars)\
 or rows_no_pivo(matrix,vars,eco)`)
    else
    m0 := args[1]:if type(args[1],string)then
    M:=convert(args[1],name) else M:=`matrix`fi:
    #lprint(`FINDING ROWS NO PIVO IN`,M):
    if nargs = 2 then vars := args[2] else 
    	vars := indets(convert(m0,set))
    fi:
    rows_no_pivo_m0 := {}:
    pivo_m0 := {}:
    if nargs = 3 then 
     lprint(antes,rows_no_pivo_m0,pivo_m0):
    fi:
    rows:={$1..rowdim(m0)}:
    for i to rowdim(m0) do
	rows:=rows minus{i}:    
	if nargs = 3 then lprint(`rows=`,"):fi:
	j := pivo(m0,i,vars):
	if nargs = 3 then print(i,j):fi:
	if j <> -1 then
	    m0 := cleanupdown(m0,i,",rows):
	    pivo_m0 := {[i,j]} union pivo_m0
	else    rows_no_pivo_m0 := {i} union rows_no_pivo_m0:
		rows:=rows union":
		    if nargs = 3 then lprint(`rows=`,"):fi:
	fi
    od:
    if nargs = 3 then lprint(agora,rows_no_pivo_m0,pivo_m0):fi:
    test := {}:
    for i in pivo_m0 do 
    if nargs = 3 then lprint(`testing pivo `,i):fi:
	m0[op(i)] := si(m0[op(i)]):
	if " <> 0 then test := test union {i}
	else rows_no_pivo_m0 := rows_no_pivo_m0 union {i[1]}
	fi
    od:
    pivo_m0 := test:
    #if 2 < 1 then
	for i in test do
	    if type(m0[op(i)],numeric) = false then
		lprint(`CHECK HONESTY OF PIVO`,i,` IN`,M):
		convert(row(m0,i[1]),set):
		if not solvable(",vars) then
		    lprint(PIVO,i,`IS HONEST!`)
		else
		    lprint(`PIVO `,i,`IS not HONEST!`):
		    lprint(EXCLUDE,i):
		    pivo_m0 := pivo_m0 minus {i}:
		    rows_no_pivo_m0 := {i[1]} union rows_no_pivo_m0
		fi
	    fi
	od
    #fi:        #if 2 < 1 then
    :RETURN([rows_no_pivo_m0,pivo_m0,eval(m0)])
fi:    
end:    #rows_no_pivo



###
sol:=proc()
###
if nargs=1 then solve(convert(args[1],set))elif
nargs=2 then solve(convert(args[1],set),convert(args[2],set))
fi:end:

#######
verifiq:=       #proc(m0,rows_with_eqn_0,eqs_0,eqs_0_1) <===OLD
#######

		 proc(m0,rows_no_pivo_,eqs_0,eqs_0_1)
local i,j:
    print(VERIFYING):
    j := sol(eqs_0,eqs_0_1):
    for i in rows_no_pivo_ do
	ec(subs(j,row(m0,i)),set):
	if " = {0} then
	    lprint(`FOUND 0 HENCE GIVEN EQNS `,
		` DO SPAN EQNS FROM ROW `,i,` OF rows_no_pivo_`.m0)
	else lprint(`#MISTAKE AT `,m0,`?`,i)
	fi
    od
end:    #of verif


###

zof:=proc(a)
###
local i,j,b,r:
if type(a,array)=false then ERROR(`WRONG TYPE:`)else
  r:=[op(2,eval(a))]:
  b:=array(sparse,op(r)):
  if nops(r)=2 then
    for i in $r[1] do for j in$r[2] do a[i,j]:
      if"<>0 then b[i,j]:=si(")fi od od:
  elif nops(r)=1 then
    for i in $r[1] do a[i]:
    if"<>0 then b[i]:=si(")fi od:
  fi
fi:
RETURN(eval(b))end:

supertex :=
proc()
local j,w0,l,i,w,teste_,x,s:
    if nargs = 1 then x := args[1]: s := 65
    elif nargs = 2 then x := args[1]: s := args[2]
    fi:
    latex(x,teste_):
    j := 1:
    w0 := ``:
    while j <> 0 do
        j := readline(teste_):
        l := length(j):
        i := 0:
        w := ` `:
        while i < l do
            i := i+1:
            if substring(j,i .. i+4) <> `\\left` then
                if substring(j,i .. i+5) <> `\\right` then
                    if substring(j,i .. i+1) <> `\\,` then
                        if not member(substring(j,i .. i),[`{`,`}`])
                             then
                            if not member(
                                substring(j,i .. i+1),[`\\{`,`\\}`])
                                 then
                                w := cat(w,substring(j,i .. i))
                            else w := cat(w,substring(j,i .. i+1))
                            fi
                        fi
                    else i := i+1
                    fi
                else i := i+5
                fi
            else i := i+4
            fi
        od:
        w0 := cat(w0,w)
    od:
    l := length(w0):
    if l <= s then lprint(w0)
    else
        x := irem(l,s,'j'):
        w := 0:
        for i to j do
            if substring(w0,i*s-1 .. i*s) = `\\i` then
                lprint(substring(w0,w+1+(i-1)*s .. i*s+1)): w := 1
            elif substring(w0,i*s .. i*s) = `\\` then
                lprint(substring(w0,w+1+(i-1)*s .. i*s+2)): w := 2
            else lprint(substring(w0,w+1+(i-1)*s .. i*s)): w := 0
            fi
        od:
        lprint(substring(w0,1+j*s+w .. j*s+x))
    fi:
    RETURN(NULL)
end: