/*
Filename: mylib_list.mac
Contents: private maxima library on the list
Created by Hideo Kodama 2022/12/8
2022/12/14: bugs in selectremove_has are fixed.
2022/12/15: bugs related with the operator "/" are fixed in selectremove
2022/12/16: maxima bug on the confusion of local and global variables is coped with.
2022/12/16: version 1.00
2024/7/17: proc: selectremove_has is modified so that  function names becomes the
  target of selection
2024/9/7: version 1.10
2024/9/28: a serious bug fix of the proc: selectremove
2024/9/29: version 1.20
Last update: 2024/9/29
*/

/******> Top */

/*****>List functions */
/****> proc:sumlist */
/* sumlist( w ) # w::list=root vector in the SRB => level height */
sumlist(w) := block(
    return (apply("+", w))
)$ /* end of proc: sumlist */

/****>proc: slist */
/* slist( L:: list, m, n) => [L[m],..., L[n]] */
slist ( L, m, n) := block ( [i],
  if listp(L) then (
      if ( not integerp(m)) or (not integerp(n) ) then (
          error("slist: the range indeces must be integers")
      ), /* end if*/ 
      if m<0 then ( m: 1),  /*end if*/
      if n> length(L) then ( n: length(L)), /*end if*/
      if m> n then ( 
          return ([])
      ) else (
          return ( makelist(L[i], i, m, n) )
      )  /*end if*/
  ) else (
      error("Usage: slist( <list>, m, n )")
  )
)$

/****> proc: memberp */
/* memberp( <x>, <list>) => a position list [i1,...] */
memberp (x, L) := block([poslist],
  if not member(x, L) then (
      return ([])
  ) else (
      poslist : [],
      for i:1 thru length(L) do (
         if L[i] = x  then (
             poslist : endcons( i, poslist)
         )
      ),
      return (poslist)
   )
)$  /*end proc: memberp */

/****> proc: memberp_has */
/* memberp_has( <x>, <list>) => a position list [i1,...] */
memberp_has (x, L) := block([rs, poslist],
   poslist: [],
   for i:1 thru length(L)  do (
      rs: selectremove_has( x, L[i]),
      if rs[1]#1 then (
         poslist: endcons(i, poslist)
      )  /*end if*/
   ), /*end do*/
   return (poslist)
)$  /*end proc: memberp_has */

/****> proc: selectremove */
/* selectremove( bfunc, eq) =>  [ eq1, eq2 ], eq=eq1*eq2
    eq = "*"(...) or f(*)
    bfunc: bool valued function
    eq1 = part of eq consisting of  elements such that bfunc=true 
    eq2 = part of eq consisting of  elements such that bfunc=false 
*/
selectremove(bfunc, _ex0) := block(
    [ /*list */
        _list0, 
      /* misc */
        _op0, expr, expr1, expr2, xx1, xx2
   ],
   
   expr: expand(_ex0),
   if atom(expr) or atom(-expr)  then (
       if bfunc(expr)  then (
            expr1: expr,  expr2:1
        ) elseif bfunc(-expr) then (
            expr1: -expr, expr2: -1
        ) else (
             expr1:1, expr2: expr
        )  /*end if*/
   ) elseif op(expr) ="+" then (
        error("expr should be a single term equation")
   ) elseif op(expr)="-" then (
        [expr1, xx2] : selectremove(bfunc, -expr),
        expr2: -xx2
   )elseif op(expr)="/" then (
       [expr1,expr2] : selectremove(bfunc, args(expr)[1]),
       [xx1, xx2]: selectremove(bdunc,args(expr)[2]),
       expr1: expr1/xx1, expr2: expr2/xx2
   ) elseif op(expr)="*" then (   
        expr1:1, expr2:1,
        for xx in args(expr) do (
             [xx1,xx2]: selectremove(bfunc, xx),
             expr1: expr1*xx1, expr2: expr2*xx2
        )  /*end do*/
    ) elseif  op(expr) ="^" then (
        if  bfunc(args(expr)[1]) then (
             expr1: expr, expr2:1
        ) else (
             expr1:1, expr2:expr
        ) /*end if*/
    ) elseif bfunc(expr) then ( 
        expr1: expr, expr2:1
    ) else (
        expr1: 1, expr2: expr
    ), /*end if */
    return ([ expr1, expr2])
)$   /*end proc: selectremove*/

/****> proc: selectremove_has */
/* selectremove_has( v, eq) 
    v = name
    eq = "*"(...) or f(*) =>  [ eq1, eq2 ], eq=eq1*eq2
    eq1 = part of eq containing v
    eq2 = part of eq not containing v
*/
selectremove_has(v, _ex0) := block(
    [ /*list */
        _list0, 
      /* misc */
        _op0, _ex, _ex1, _ex2, _xx1,_xx2
   ],
   
   _ex: expand(_ex0),
   if atom(_ex) then (
       if _ex=v then (
            _ex1: v,  _ex2:1
        ) else (
             _ex1:1, _ex2: _ex
        )  /*end if*/
   ) elseif atom(-_ex) then (
       if _ex=-v then (
            _ex1: v,  _ex2: -1
        ) else (
             _ex1:1, _ex2: v
        )  /*end if*/
   ) else (
       _op0 : op(_ex),  _list0: args(_ex),
       if _op0="+" then (
            error("ex should be a single term equation")
       ) elseif _op0="-" then (
           [_xx1,_xx2]: selectremove_has(v, -_ex),
           return ( [_xx1, -_xx2])
       ) elseif string(_op0) = string(v) then (
           return ( [_ex, 1] )
       ) elseif  _op0 ="^" then (
           if  atom(_list0[1]) then (
               if  _list0[1]=v then (
                    _ex1: _ex, _ex2:1
               ) else (
                    _ex1:1, _ex2:_ex
               ) /*end if*/
           ) elseif op(_list0[1])=v then (
               _ex1: _ex, _ex2:1
           ) elseif args(_list0[1])[1]=v or args(_list0[1])[2]=v then (
               _ex1: _ex, _ex2:1
           ) else (
               _ex1:1, _ex2:_ex
           )  /*end if*/
       ) elseif _op0=v then (
           _ex1: _ex, _ex2:1
       )elseif _op0="/" then (
           [_ex1,_ex2] : selectremove_has(v, _list0[1]),
           [_xx1, _xx2]: selectremove_has(v, _list0[2]),
           _ex1: _ex1/_xx1, _ex2: _ex2/_xx2
       ) elseif _op0="*" then (     
           _ex1:1, _ex2:1,
           for _xx in _list0 do (
               [_xx1,_xx2]: selectremove_has(v, _xx),
               _ex1: _ex1*_xx1, _ex2: _ex2*_xx2
           )  /*end do*/
       ) else (
           _ex1: 1, _ex2: _ex,
           for _xx in _list0 do (
               if selectremove_has(v, _xx)[1]#1 then (
                   _ex1: _ex, _ex2:1,
                   break
               )  /*end if*/
           )  /*end do*/
       )
   ), /*end if */
   return ([ _ex1, _ex2])
)$   /*end proc: selectremove_has*/

/****>proc: selectremove_has_sum */
/* selectremove_has_sum( v, eq) 
    v = name
    eq = "+"(...) or f(*) =>  [ eq1, eq2 ], eq=eq1+eq2
    eq1 = the sum  of terms  containing v in eq
    eq2 = the sum  of terms  not  containing v in eq
*/
selectremove_has_sum(v, _ex0) := block( [_ex, _xx1,_xx2, _ex1,_ex2,_term,  _termlist],

    _ex: expand(_ex0),
    if (not atom(_ex)) and op(_ex)="-" then (
        return( - selectremove_has_sum(-_ex) )
    ), /*end if*/
    if  atom(_ex) or ((not atom(_ex)) and op(_ex)#"+" )  then (
        [_xx1,_xx2]: selectremove_has(v, _ex),
        if _xx1#1 then (
            _ex1: _ex, _ex2: 0
        ) else (
            _ex1:0,  _ex2: _ex
        ) /*end if*/
    ) else (         /* "+"(....) */
        _termlist : args(_ex),
        _ex1: 0, _ex2: 0,
        for _term in _termlist do (
            [_xx1,_xx2]: selectremove_has( v, _term),
            if _xx1#1 then (
                _ex1: _ex1 + _term
            ) else (
                _ex2: _ex2 + _term
            )  /*end if*/
        ) /*end do*/
    ),  /*end if */
    return ([ _ex1, _ex2])
)$  /*end proc: selectremove_has_sum */
        


/****>proc:permlist */
/* permlist(indlist) => the list of permutated indices
*/
permlist(indlist) :=block(
    [ /*list*/
       sl, elist, olist, splist,eplist, oplist,
       /*misc*/
       len, x
    ],

    len: length(indlist),
    if len=1 then (
        return ( [[indlist],[]] )
    ) else (
        sl: slist(indlist, 1, (len-1)),
        splist: permlist(sl),
/*  display(splist), */
        eplist: splist[1], 
        oplist: splist[2],
        x:indlist[len],
/*  display(eplist, oplist), */
/*  display(sl, splist, x), */
         elist:[], olist:[],
         for i :  len thru 1 step -1 do (
              if (-1)^(len-i)=1 then (
                  elist: append(elist, map(lambda([y],  append(slist(y, 1, (i-1)),[x],slist(y, i, (len-1)))), eplist)),
                  olist: append(olist, map(lambda([y],  append(slist(y, 1, (i-1)),[x],slist(y, i, (len-1)))), oplist))
              ) else (
                  elist: append(elist, map(lambda([y],  append(slist(y, 1, (i-1)),[x],slist(y, i, (len-1)))), oplist)),
                  olist: append(olist, map(lambda([y],  append(slist(y, 1, (i-1)),[x],slist(y, i, (len-1)))), eplist))
              )  /*end if*/
         ), /*end do*/
         return ( [elist,olist])
    )  /*end if*/
)$  /*end proc: permlist*/

/****>proc: permlist1 */
/* permlist1(indlist) => list of permutated index with signature
*/
permlist1(indlist) := block( [list1,list2],
    [list1,list2]: permlist(indlist),
    list1: map( lambda( [x], [x, +1]), list1),
    list2: map( lambda( [x], [x, -1]), list2),
    return (append(list1,list2))
)$  /*end proc: permlist1 */


/*****>Lattice related functions */

/****> proc:mklattice */
/*  mklattice(rn::integer,[n1,n2]) =>  lattice point list
      rn = lattice dimensions
      [n1,n2] = the range of the integer coordinates of lattice points
*/
mklattice ( rn, rng) :=  block( [sublattice,newlattice,x ],
    if  (not nonnegintegerp(rn)) or ( not listp(rng)  or  length(rng)<2 
         or not integerp(rng[1]) or not integerp( rng[2] ) ) then (
         error("Usage: mklattice(rn::positive integer, [n1::integer,n2::integer])")
    ), /*end if*/
    simp: true,
    
    if rn=0 then (
         newlattice: [[]]
    ) elseif  rn=1 then (
         newlattice :  makelist( [i], i, rng[1], rng[2])
    ) else (
        sublattice : mklattice(rn-1,rng),
        newlattice : [],
        for x in sublattice do (
            newlattice: append( newlattice, makelist( endcons(i, x), i, rng[1], rng[2]))
        )  /*end do*/
    ),  /*end if*/
    return (newlattice)
)$ /* end proc: mklattice */


/****> proc:mklevellattice */
/*  mklevellattice(rn::integer>0, level::integer>=0)=> level lattice point set 
*/
mklevellattice (rn, level) := block( [ x,lattice,  levellattice],
    lattice : mklattice(rn,[0,level]),
    levellattice : [],
    for x in lattice do (
        if sumlist(x)=level then (
            levellattice : endcons(x, levellattice)
        ) /*end if*/
    ), /*end do*/
    return (levellattice)
)$ /* end proc: mklevellattice */

/******>End */
print("loading mylib_list.mac  --- done");