Yes, I wrote all this "for fun" on a 320x240 screen with a non qwerty keyboard. Terry would be so proud.

apm.xcas

cis:= (theta_)→
  cos(theta_)+i*sin(theta_);

// Find complex solutions to equations
// by isolating the imaginary and real
// parts
gcsolve:= (eqs_, vars_)→
  solve(concat({}, 
               im(eqs_), 
               re(eqs_)), 
        vars_);

comp2cart:= (expr_)→
  [re(expr_), im(expr_)];

cart2comp:= (x_, y_)→
  x_+i*y_;

withcis:= (comp_)→
  abs(comp_)*cis(arg(comp_));

compexp(z_, exp_):= begin
  if !realn(exp_) then
    return eval(ret_me);
  end;
  exp_:= exact(exp_);
  z_:= (z_)^(numer(exp_));
  return gapply(
    gprogram(['k_'],[0],
      ((abs(z_))^(1/denom(exp_)))*
        cis((arg(z_)+2*π*'k_')/
            denom(exp_))
    ),
    {(k_)$(k_=(0..(denom(exp_)-1)))} 
  );
end;

compexp2:= (u_,w_,k_)->(e^(
    re(w_)*ln(abs(u_))-
    im(w_)*(arg(u_)+2*π*k_))
  *cis(im(w_)*ln(abs(u_))+
       re(w_)*(arg(u_)+2*π*k_)));

cln:= (z_, k_)→ (
  ln(abs(z_))+*(arg(z_)+2*π*k_)
);

cgte:= (f_, z_, z0_, n_)→ (
  subst(
    ((2/((n_-1)!))*
      diff(f_, z_, n_-1)*π*),
  z_=z0_)
);

is_e_harmonic:= (expr_, v1_, v2_)→(
  eq(diff(expr_, v1_, 2) +
     diff(expr_, v2_, 2), 0)
);

e_harmonic_conj:= (expr_, v1_, v2_)→(
potential({-(diff(expr_ ,v2_)),
             diff(expr_ ,v1_)},
          {v1_, v2_})
);

lineint_acheck:= (vectorf_,
           vars_,
           curve_,
           param_,
           limits_)→ (
  ( type(vectorf_) == type([]) and
    type(vars_) == type([]) and 
    type(curve_) == type([]) and
    type(limits_) == type([])) and 
  type(param_) == type('x') and
  ( size(vectorf_) == size(vars_) ==
    size(curve_)) and
  size(limits_) == 2 and
  gand(apply((x_)→type(x_)==type('x'), 
             vars_))
);

lineint:= (vectorf_,
           vars_,
           curve_,
           param_,
           limits_)→ (
  lineint_acheck(op(tail(args))) ? (
    gint(subst(vectorf_, vars_ = curve_)*
          diff(curve_, param_),
          param_, op(limits_))
  ) : (eval(ret_me))
);

fluxint:= (vectorf_,
           vars_,
           curve_,
           param_,
           limits_)→ (
  lineint_acheck(op(tail(args))) ? (
    gint(subst(vectorf_, vars_ = curve_)*
          (((v_)→([v_[2], -v_[1]]))(
            diff(curve_, param_))),
          param_, op(limits_))
  ) : (eval(ret_me))
);

compint(vectorf_,
        vars_,
        curve_,
        param_,
        limits_):= begin
  if !lineint_acheck(
       op(tail(args))) or
     size(vectorf_) != 2 then
    return eval(ret_me);
  end;
  vectorf_:= [vectorf_[1], -vectorf_[2]];
  return lineint(op(tail(args)))+
         *fluxint(op(tail(args))) 
end;

de01(coeffs_, vars_, var_):= begin
  (e.^(exact(proot(coeffs_)).*var_))*vars_;
end;

mdomain:= (expr_, vars_) →(
  (type(vars_) == DOM_LIST and
  gand(op(apply(
         (v_)→(type(v_) == DOM_IDENT), 
         vars_)))
  ) ? (
    (
      op(apply(
           (v_)→domain(expr_, v_),
           vars_))
    )
  ) : eval(ret_me)
);

sinsum:= (x1_, x2_)→ (
  sin(x1_)*cos(x2_)+
  cos(x1_)*sin(x2_)
);

cossum:= (x1_, x2_)→ (
  cos(x1_)*cos(x2_)-
  sin(x1_)*sin(x2_)
);

wronskian:= (exprs_, var_)→ (
  apply((n_)→ (
          apply(
            (f_)→diff(f_, var_, n_),
            exprs_)
        ),
        {seq(t, t=(0..(size(exprs_)-1)))})
);

vpar_yp(yhs_, var_, fun_):= begin
  local wrons_:= wronskian(yhs_, var_);
  local dw_, swap_, i_, (ans_:= []);
  dw_:= det(wrons_);
  for i_ from 1 to size(yhs_) do
    swap_:= transpose(wrons_);
    swap_[i_]:=
      [seq(0,t=1..size(yhs_)-1), fun_]; 
    ans_[i_]:= 
      gint(det(transpose(swap_))/dw_,
           var_);
  end;
  return ans_;
end;

hs2:= (x1_, x2_, var_)→ (
  Heaviside(var_-x1_)-
    Heaviside(var_-x2_) 
);

hs:='Heaviside';

laplace_diff:=
(fun_, var_, lvar_, ord_)→(
  laplace(fun_, var_, lvar_)*
    lvar_^ord_ -
  sum(apply(
    (k_)→diff(fun_, var_, k_)(var_=0)*
           lvar_^(ord_-k_-1),
    range(ord_)))
);

applib.xcas

ag_active:= (NULL)→ (
  type(Instruction)==DOM_STRING
);

ag_vars():= begin
  local (vars_:= {}), (i_:= 1);
  if !(ag_active()) then
    return eval(ret_me); end;
  if Instruction() == "" then
    return vars_; end;
  while (Instruction(i_-1, 2) != 
        Instruction(i_, 2)) or 
        (i_ == 1) do
    vars_:= append(vars_, 
                   Instruction(i_, 2));
    i_:=i_+1;
  end;
  return vars_;
end;

ag_instruction:= (n_, k_, set_='nop', 
                  force_=false)→ (
  (ag_active() and
   type(n_) == DOM_INT and
   type(k_) == DOM_INT) ? (
    (set_ != 'nop' or force_) ? (
      expr("Instruction("+
              n_+","+k_+"):= set_")
    ) : (
      expr("Instruction("+
              n_+","+k_+")")
    )
  ) : eval(ret_me)
);

cachelib.xcas

is_cache:= (cache_)→ (
  type(cache_) == DOM_MAP 
);

if !is_cache(cache_store) then
  cache_store:= table();
end;

// preferible for "pure" functions
cache(getter_, targs_, force_=false,
      extra_key_=NULL):= begin
  local key_, value_, result_;
  key_:= table('getter'=getter_, 
               'targs'=targs_,
               'extra_key'=[extra_key_]);
  if cache_store[key_] == 0 or force_ then
    result_:= gtry2((getter_)@top, targs_);
    if result_[1] then
      throw(
        "cache: Exception raised when " +
        "calling to cache value with " +
        "message: " + result_[3]); return;
    end;
    value_:= op(result_[2]);
    cache_store[key_]:= table(
      'result'=value_);
  end;
  return cache_store[key_]['result'];
end;

// an order of magnitude faster
cache2(getter_, arg_, force_=false,
       extra_key_=NULL):= begin
  local (key_:= folder[]),
        value_,
        (result_:=folder[true,0,0]),
        (error_msg_:= 0),
        (tkey_:=table());
  key_[1]:= getter_;
  key_[2]:= arg_;
  key_[3]:= extra_key_;
  tkey_:= table(1=key_);
  value_:= cache_store[tkey_];
  if !(value_ in [0]) and
     !force_ then
    return value_[2];
  end;
  gtry_catch('result_[2]:= getter_(arg_)',
             'error_msg_', 'nop',
             'result_[1]:= false');
  result_[3]:= error_msg_;
  if result_[1] then
    throw(
      "cache: Exception raised when " +
      "calling to cache value with " +
      "message: " + result_[3]); return;
  end;
  value_:= folder[1];
  value_[2]:= result_[2];
  cache_store[tkey_]:= value_;
  return value_[2];
end;

cache2_orig:= cache2;
cache2:= prog_tmpify(
  cache2_orig,
  gminus(lname(cache2_orig),
    ['cache_store','gtry_catch','throw'])
)();

clear_cache(NULL):= begin
  local old_cache_:= cache_store;
  cache_store:= table();
  return old_cache_;
end;

calc.xcas

idiff(function_, vars_):= begin
  local result_, index_;
  result_ := function_;
  for index_ from 1 to dim(vars_) do
    result_ := diff(result_, vars_[index_]);
  end;
  return result_;
end;

fcurl:=(unapply@top)@(
  (tab_)→(
    tab_[1]:= (curl@top)(tab_)
))@(sym_env);

fdiff:= (fun_)→(
  fun_[4]:= [diff(fun_[4][3], 
              vectorize(fun_[2]))]
);

riint(function_, vars_, limits_):= begin
  return iint(function_,
              vars_,
              gapply(rtolist, limits_)); 
end;

giint:= (expr_, vars_)→ (
  reduce((accum_, var_)→
    gint(accum_, left(var_),
           left(right(var_)),
           right(right(var_))),
    vars_, expr_)
);

iint(function_, vars_, limits_):= begin
  local gint_, index_, size_, result_;
  gint_ := 'int';
  result_ := function_;
  size_ := min(size(vars_), size(limits_));
  for index_ from 1 to size_ do
    result_ := gint_(result_,
                     vars_[index_], 
                     limits_[index_](1), 
                     limits_[index_](2));
  end;
  return result_;
end;

iintf:= (expr_, vars_, limits_) →
  reduce((accum_, el_)→
      gint(accum_, el_[1],
           el_[2][1], el_[2][2]),
    gzip(pair, vars_, limits_),
    expr_
  )
;

dmean(function_, vars_, limits_):= begin
  local index_, size_, result_, total_;
  result_ := {};
  total_ := iint(function_, vars_, limits_);
  size_ := min(size(vars_), size(limits_));
  for index_ from 1 to size_ do
    result_[index_] := 
      iint(function_ * vars_[index_],
           vars_, limits_)/total_;
  end;
  return result_;
end;

dmeanf:= (expr_, vars_, limits_) →
  gapply((v_)→iintf(expr_ * v_,
      vars_, limits_), vars_) ./ 
    iintf(expr_, vars_, limits_)
;

frenet_serret(vector_, par_, var_):= begin
  local param_s_, param_t_, gint_ := 'int';
  local result_ := {};
  param_s_ := 
    solve(
      gint_(vabs(diff(vector_, par_)), par_, 0, par_) = var_,
      par_)[1];
  result_[1] := diff(subst(vector_, par_ = param_s_), var_);
  result_[2] := diff(result_[1], var_)/vabs(diff(result_[1], var_));
  param_t_ := solve(param_s_=par_, var_);
  result_[3] := cross(result_[1], result_[2]);
  return result_, param_t_;
end;

// WIP
frenet_serretf:= (vec_, par_, var_) →
  gapply((s_)→
    reduce(
      (state_, elem_)→(elem_(state_)),
      [
        (state_)→(state_['t']:= solve(state_['s']=par_, var_)),
        (state_)→(state_['result']:=append(gclear(vec_), 
          diff(subst(vec_, par_ = state_['s']), var_))),
        (state_)→(state_['result']:=append(state_['result'],
          diff(state_['result'][1], var_)/vabs(diff(state_['result'][1], var_)))),
        (state_)→(state_['result']:=append(state_['result'],
          cross(state_['result'][1], state_['result'][2])))
      ],
      table(
        's' = s_
      )
    ),
    solve(gint(vabs(diff(vec_, par_)), par_, 0, par_) = var_, par_)
  )
;

fourier_series(expr_, var_, T_, n_, a_, save_ = true):= begin
  local args_, i_, result_, an_, bn_, saved_b_, saved_a_, ni_;
  args_ := [expr_, var_, T_, nvar_, a_];
  if save_ then
    saved_a_ := fourier_an(op(args_));
    saved_b_ := fourier_bn(op(args_));
    an_ := 'limit(saved_a_, nvar_, ni_)';
    bn_ := 'limit(saved_b_, nvar_, ni_)';
  else
    an_ := 'fourier_an(op(subst(args_, nvar_=ni_)))'
    bn_ := 'fourier_bn(op(subst(args_, nvar_=ni_)))'
  end;
  ni_ := 0;
  result_ := (1/2) * eval(an_);
  for ni_ from 1 to n_ do
    result_ += eval(an_) * cos(2*ni_*π*var_/T_) + eval(bn_) * sin(2*ni_*π*var_/T_);
  end;
  return result_;
end;

upolar:= (r_, t_) →
  [r_*cos(t_), r_*sin(t_)];

cylindric:= (r_, t_, z_) →
  [op(upolar(r_, t_)), z_];

spheric:= (r_, t_, a_) →
  [r_*sin(t_)*cos(a_),
   r_*sin(t_)*sin(a_),
   r_*cos(t_)]; 

jacobian(functions_, vars_):= begin
  return diff(functions_, vars_);
end;

vabs := (vector_) →
  sqrt(dot(vector_, vector_));

vconvert := (vector_, unit_, depth_ = 1) → 
  dapply(((x_)→convert(x_, unit_)), vector_, depth_);

vusimplify := (vector_, depth_ = 1) →
  dapply(((x_)→usimplify(x_)),vector_, depth_);

// WIP
iltransform(vars_, limits_, new_vars_, trns_):= begin
  local lower_limits_ := {};
  local upper_limits_ := {};
  local limit_, index_;
  local scoped_vars_ := {};
  local head_ := (list_, up_to_) -> 
    suppress(list_, MAKELIST(x_, x_, up_to_ + 1, size(list_) + 1));
  for index_ from 1 to size(vars_) do
    lower_limits_[index_] := limits_[index_, 1];
    upper_limits_[index_] := limits_[index_, 2];
  end;
  for index_ from 1 to size(vars_) do
    upper_limits_[index_] := tolist(solve(prepend(
      head_(trns_ = makelist(0, 1, size(trns_)), index_ - 1),
        subst(vars_[index_] = limits_[index_, 1], vars_=trns_)),
      head_(new_vars_, index_)));
  end;
  return upper_limits_;
end;

ciran.xcas

si:= (args_)→(inv(
  sum(apply(inv, 
            vectorize(args_)))
  )
);

star:= (r1_, r2_, r3_)→(
  r1_+r2_,r2_+r3_,r3_+r1_
);

tri:= (ra_, rb_, rc_)→(
  si(rc_,ra_+rb_),
  si(ra_,rb_+rc_),
  si(rb_,rc_+ra_) 
);

tusimp:= (expr_)→(
  subst(expr_, {sin=a})
);

tf2zp:= (numer_, denom_, var_)→ (
  id(gfroot(numer_, var_),
     gfroot(denom_, var_),
     lcoeff(numer_, var_)/
     lcoeff(denom_, var_)) 
);

// returns a mapping of type
// order=location for each root
// and pole of a rational polynomial
// expr_ of variable var_
gfroot:= gfcomp((zeros_poles_)→ (
  reduce((curr_, pair_)→ (
      (curr_)(pair_[1]=pair_[2])
    ), 
    zeros_poles_,
    fcurry(table, size(zeros_poles_))
  )()
),(expr_, var_)→ (
  chunked(froot(expr_, var_), 2)
));

//gfroot:= (expr_, var_)→ (
//  ((zeros_poles_)→ (
//    reduce(pprint((curr_, pair_)→(
//        (curr_)(pair_[1]=pair_[2])
//      ),
//      zeros_poles_,
//      fcurry(table, size(zeros_poles_)) 
//    ))()
//  ))(chunked(froot(expr_, var_), 2))
//);

gfcoeff:=(roots_, var_)→(
  fcoeff(reduce(
    (a_,e_)→(
      append(a_,left(e_),right(e_))),
    table2list(roots_),[]), var_)
);

bodize:= (H_w_, w_, x_)→ (
  apply(
    (y_)→20*log10(subst(y_, 
      (w_=10^(x_)))),
    [abs(H_w_), arg(H_w_)])
);

mcoeff(inner_, outer_, ms_):= begin
  local loi_;
  local coeffs_:=seq([1, size(inner_)]);
  for lo_ in 1..size(outer_) do
    loi_:= lo_ + size(inner_) - 2 
    for li_ in 1..size(inner_) do
      loi_ += 1;
      sign(inner_[li_])*outer_[lo_]*ms_[loi_] 
    end;
  end;
end;

linterp:= (p1_, p2_, r_)→ (
  p1_+(p2_ - p1_)*r_
)

cprod_strat.xcas

cprods_v1(vars_, limits_):= begin
  local gand_, gprogram_, conds_, size_, index_;
  gand_ := 'and'
  gprogram_ := 'program';
  conds_ := {};
  size_ := min(size(vars_), size(limits_));
  for index_ from 1 to size_ step 1 do
    conds_[index_] := 'limits_'[index_,1] ≤ 'vars_'[index_] ≤ 'limits_'[index_,2];
  end;
  return gprogram_(vars_, nop, gand_(conds_));
end;

//cprods_v2(vars_, limits_):= begin

cprod(variable_, limits_):= begin
  local gprogram_;
  local lower_limit_, upper_limit_;
  gprogram_ := 'program';
  lower_limit_:=limits_[1];
  upper_limit_:=limits_[2];
  return gprogram_(variable_, 
    nop, 
    'lower_limit_ ≤ variable_ ≤ upper_limit_');
end;

cprod2d(variables_, limits_):= begin
  local gprogram_, index_;
  local lower_limit_1_, upper_limit_1_, variable_1_;
  local lower_limit_2_, upper_limit_2_, variable_2_;
  gprogram_ := 'program';
  variable_1_:=variables_[1];
  lower_limit_1_:=limits_[1,1];
  upper_limit_1_:=limits_[1,2];
  variable_2_:=variables_[2];
  lower_limit_2_:=limits_[2,1];
  upper_limit_2_:=limits_[2,2];
  return gprogram_({variable_1_, variable_2_}, 
    nop, 
    'lower_limit_1_ ≤ variable_1_ ≤ upper_limit_1_ AND lower_limit_2_ ≤ variable_2_ ≤ upper_limit_2_');
end;

make_cprod(nvars_ = 1):= begin
  local gprogram_, index_, cond_;
  local var_name_, ll_name_, ul_name, vars_arg_, lims_arg_;
  local body_, var_list_;
  vars_arg_:="variables_";
  lims_arg_:="limits_";
  var_name_:=(index_)->"variable_" + index_ + "_";
  ll_name_:=(index_)->"lower_limit_" + index_ + "_";
  ul_name_:=(index_)->"upper_limit_" + index_ + "_";
  cond_ := "";
  var_list_ := "{ ";
  body_ := 
"begin
  local gprogram_ := 'program';";
  for index_ from 1 to nvars_ do
    body_ := body_ + "
  local " + var_name_(index_) + " := " + vars_arg_ + "[" + index_ + "];
  local " + ll_name_(index_) + " := " + lims_arg_ + "[" + index_ + ", 1];
  local " + ul_name_(index_) + " := " + lims_arg_ + "[" + index_ + ", 2];";
    if index_ > 1 then
      cond_ := cond_ + " AND ";
      var_list_ := var_list_ + ", ";
    end;
    var_list_ := var_list_ + var_name_(index_);
    cond_ := cond_ + ll_name_(index_) + " ≤ " + var_name_(index_) + " ≤ " + ul_name_(index_);
  end;
  var_list_ := var_list_ + " }";
  body_ := body_ + "
  return gprogram_(" + var_list_ + ", 'nop', '" + cond_ + "');
end;";
  return ("(" + vars_arg_ + ", " + lims_arg_ + ")->" + body_);
end;


cprods(variables_, limits_):= begin
  local size_, index_, var_name_, ll_name_, ul_name_;
  local var_list_, cond_, result_, gprogram_, tmp_mark_;
  size_ := min(size(variables_), size(limits_));
  // Since we're defining temp globals we should
  // randomize the names to avoid name collisions
  tmp_mark_ := rand(100000000);
  gprogram_ := 'program';
  var_name_:=(index_)->"variable_" + index_ + "_" + tmp_mark_ + "_";
  ll_name_:=(index_)->"lower_limit_" + index_ + "_" + tmp_mark_ + "_";
  ul_name_:=(index_)->"upper_limit_" + index_ + "_" + tmp_mark_ + "_";
  cond_ := "";
  var_list_ := "{ "
  for index_ from 1 to size_ step 1 do
    (#(var_name_(index_))) := variables_[index_];
    (#(ll_name_(index_))) := limits_[index_,1];
    (#(ul_name_(index_))) := limits_[index_,2];
    if index_ > 1 then
      cond_ := cond_ + " AND "
      var_list_ := var_list_ + ", ";
    end;
    var_list_ := var_list_ + var_name_(index_);
    cond_ := cond_ + ll_name_(index_) + " ≤ " + var_name_(index_) + " ≤ " + ul_name_(index_);
  end;
  var_list_ := var_list_ + " }";
  result_ := expr("gprogram_("+ var_list_ +", nop, '" + cond_ + "')");
  for index_ from 1 to size_ step 1 do
    cas("purge(" + var_name_(index_) + ");");
    cas("purge(" + ll_name_(index_) + ");");
    cas("purge(" + ul_name_(index_) + ");");
  end;
  return result_;
end;

deprecated.xcas

// Programs that served a purpose until
// I found less hacky alternatives

// Replacement: ``func_(op(args_))``


// Replacement: ``op(vec_)``
// Convers an iterable to a sequence
toplist(vec_):= begin
  local func_sym_;
  purge(func_sym_);
  return vacall(func_sym_, vec_)[3];
end;

// Replacement: ``{op(vector_)}``
tolist(vector_, size_ = -1):= begin
  local index_;
  local result_ := {};
  if size_ == -1 then size_ := size(vector_); end;
  for index_ from 1 to size_ do
    result_[index_] := vector_[index_];
  end;
  result result_;
end;

drawlib.xcas

// H/V
screen_ratio:= 159/109;

sframe:= (min_,
          offset_=0,
          ratio_=screen_ratio)→ (
  nop 
);

// plotlist
segments:= (items_)→(
  apply(segment@op,
        adjacent_pairs(
          vectorize(items_)))
);

clamp:= (value_, range_)→(
  is_range(range_) ? ( 
    ((value_ MOD
      (right(range_) - left(range_)))
    ) + left(range_)) : (
    eval(ret_me))
);

// circumcircle
circle3:= (p1_, p2_, p3_)→
  ((r_)→(circle(r_, p2_-r_)))(
    single_inter(
      perpen_bisector(p1_, p2_), 
      perpen_bisector(p2_, p3_))
);

plotcont:= (expr_, v1_, v2_, vals_)→ (
  (type(vals_)==DOM_LIST and
  type(expr_)==DOM_SYMBOLIC) ? (
    apply(
      (k_)→plotimplicit(expr_=k_,[v1_,v2_]),
      vals_)
    ) : ( eval(ret_me)
  )
);

xsum:= (objs_)→ (
  plotimplicit(sum(apply(
    (obj_)→one(
      solve(gequation(obj_),'x')),
  vectorize(objs_)))='x')
);

ysum:= (objs_)→ (
  plotimplicit(sum(apply(
    (obj_)→one(
      solve(gequation(obj_),'y')),
  vectorize(objs_)))='y')  
);

emag.xcas

mu0:=1.2566370614ᴇ−6_(H/m);

exprlib.xcas

// walks expressions
opply(f_, expr_):= begin
  local i_, tmp_expr_, opd_, nested_;
  //print(f_, expr_);
  if type(expr_) == DOM_LIST then
    for i_ from 1 to size(expr_) do
      expr_[i_]:= f_(
        seq_nest_symb('expr_[i_]')());
    end;
    return expr_;
  end;
  if type(expr_) == DOM_SYMBOLIC then
    if gsommet(expr_) == at then
      opd_:= op(expr_);
      expr_[1]:= f_(at);
      expr_[2]:= f_(
        seq_nest_symb('opd_[1]')());
      expr_[3]:= f_(
        seq_nest_symb('opd_[2]')());
      return expr_;
    end;
    for i_ from 1 to size({op(expr_)})+1 do
      print(i_);
      print(expr_[i_]);
      expr_[i_]:= f_(
        seq_nest_symb('expr_[i_]')());
    end;
    return expr_;
  end;
  if type(expr_) == DOM_FUNC then
    if gsommet(expr_) != 'program' then
      return f_(expr_);
    end;
    for i_ from 1 to size({op(expr_)})+1 do
      expr_[i_]:= f_(
        seq_nest_symb('expr_[i_]')());
    end;
    return expr_;
  end;
  return f_(expr_);
end;

is_atomic:= (expr_) →(
  type(expr_) != DOM_LIST and
  type(expr_) != DOM_SYMBOLIC and
  (type(expr_) != DOM_FUNC or 
   gsommet(expr_) != 'program') 
);

quote_symb:= (symb_, n_)→ (
  expr("(seq[])→ "+
         "quote("*n_+"symb_"+")"*n_
  )
);

subst_eq:= (expr_, symbol_) →(
  subst(ffget(
    folder_build(expr_)(geq)(symbol_)
  ))
);

// attempt at making programs with
// local variable names randomized
prog_tmpify(fun_, temp_vars_):= begin
  local tmp_fn_, locals_, interm_,
        locals_expr_, temp_symbols_;
  if !is_program(fun_) or
     size(op(fun_[4])) != 2 or
     fun_[4][1] != glocal then
    return eval(ret_me); end;
  tmp_fn_:= expr(temp_var());
  locals_:= [op(set[
    op(lname(fun_)),
    op(temp_vars_:=lname(temp_vars_))
  ])];
  locals_expr_:= folder[];
  temp_symbols_:= apply(
    (var_)→str2ident(
      temp_var(ident2str(var_))),
    temp_vars_
  );
  locals_[size(locals_)+1]:= tmp_fn_;
  locals_expr_[1]:= op(locals_);
  locals_expr_[2]:= seq[];
  //print(temp_symbols_[3][1]);
  interm_(NULL):= begin
    local locals_;
    purge(locals_);
    temp_vars_:= temp_symbols_;
    tmp_fn_:= (NULL)→ NULL;
    // subst gets messed up for default
    // values set to seq[]
    tmp_fn_[2]:= [];
    tmp_fn_[4][1]:= glocal;
  end;
  interm_[4][2]:= [op(locals_expr_)];
  interm_[4][3][1][2]:= op(locals_);
  interm_[4][3][2][3]:= op(temp_vars_);
  interm_[4][3][3][2]:= fun_;
  interm_[4][3][3][2][2]:=seq[];
  // hack to avoid symb_ = seq[] being
  // substituted by seq[symb_] = seq[],
  // might break with seqs of lengths >1 
  interm_[4][3][4][2]:= subst(ffget(
    folder_build(fun_[2])(geq)('geq')
  ));
  interm_[4][3][3][2][4][1]:= quote;
  interm_
end;

deep_subst(symb_, val_, expr_):= begin
  local nested_:= seq_nest_symb('expr_')();
  local self_args_, opply_args_;
  //print(symb_, val_, expr_);
  if is_atomic(nested_) then
    print("atomic");
    return subst(expr_, symb_, val_);
  end;
  //print("not atomic");
  self_args_:= folder[];
  self_args_[1]:= symb_;
  self_args_[2]:= val_; 
  opply_args_:= folder[];
  opply_args_[1]:= (oexpr__)→ (
    gquote('deep_subst'(op(
      ((args__=self_args_)→ (
        args__[3]:= oexpr__
      ))()
    )))
  );
  opply_args_[2]:= expr_;
  opply(op(opply_args_))
end;

gsubst(eqs_):= begin
  local eqs_subst_:= folder[]; 
  local args_reduce_:= folder[0,0,0];
  local e_, (i_:= 1), eq_subst_;
  local reduce; purge(reduce);
  for e_ in eqs_ do
    eq_subst_:= folder[seq[]];
    eq_subst_[2]:= left(e_);
    eq_subst_[3]:= right(e_); 
    eqs_subst_[i_]:= eq_subst_;
    i_++; 
  end;
  args_reduce_[2]:= eqs_subst_;
  args_reduce_[1]:= (accum_, eq_)→ (
    subst(op(eq_[1]:=gquote(accum_)))
  );
  (expr_)→ (
    reduce(op(
      ((tmp_arg_=args_reduce_)→
         (tmp_arg_[3]:=expr_))()
    ))
  )
end;

seq_nest_symb(
    symb_,
    temp_var_=expr(temp_var())):= begin
  local (temp_expr_:= ((seq[])→0)),
        doer_;
  temp_expr_[3]:= seq[0];
  doer_:= (temp_var_=NULL)→ (
    (temp_var_[3][1]:= symb_)[3]);
  doer_[2][1][3]:= temp_expr_;
  doer_
end;

sto_symb(left_, right_):= begin
  (NULL)→ (left_:= right_)
end;

sto_quoted(right_, left_):= begin
  local func_:= (NULL)→ 'left_:=right_';
  func_[4][2][3]:= left_;
  func_[4][2][2]:= right_;
  func_()
end;

csubst(expr_, eqs_):= begin
  local my_sto_, subster1_, subster2_;
  local subst_args_:= folder[];
  subster1_:= gsubst(eqs_);
  subster2_:= gsubst(eqs_);
  subst_args_[1]:= gquote(expr_);
  subst_args_[2]:= sto;
  subst_args_[3]:= sto_quoted@subster1_@quote;
  subster2_(gquote(subst(op(subst_args_))))
end;

// binary ops
bin_symb(op_):= begin
  local tmp_:= expr(temp_var());
  local result_:= 
  (left_, right_)→ begin
    local expr_;
    expr_:= (NULL)→ 'right_:= left_';
    expr_[4][2][2]:= left_;
    expr_[4][2][3]:= right_;
  end;
  result_[4][3][1][2][4][2][1]:= op_;
  result_
end;

eq_symb:= bin_symb('=');

explode(expr_):= begin
  local root_, leaves_;
  if is_leaf(expr_) then
    return expr_;
  end;
  [root_, leaves_]:= parts(expr_);
  return [root_, apply(explode, leaves_)];
end;

str2ident(str_):= begin
  local result_;
  if type(str_) != type("") then
    return eval(ret_me); end;
  if contains(str_, "`") then
    warn("str2ident: " +
         "identifiers cannot contain `.");
    return eval(ret_me);
  end; 
  result_:= expr("'`" + str_ + "`'");
  if type(result_) != DOM_IDENT then
    warn("str2ident: " +
         "resulting expression is not " +
         "an identifier");
  end;
  return result_;
end;

gsommet(expr_):= begin
  if type(expr_) == DOM_SYMBOLIC then
    if part(expr_, 0) == "at" then
      return at;
    end;
    return expr_[1];
  end;
  if type(expr_) == DOM_FUNC and
     size({op(expr_)}) == 3 then
    return expr_[1];
  end;
  return id;
end;

ident2str(ident_):= begin
  local result_;
  if type(ident_) != DOM_IDENT then
    return eval(ret_me);
  end;
  result_:= str(ident_);
  if contains(result_, "`") then
    warn("ident2str: " +
         "identifier contains `.");
  end;
  if type(result_) != type("") then
    warn("ident2str: " +
         "result is not a string.");
  end
  return result_; 
end;

is_leaf:= (expr_)→ (
  part(expr_, 0) ==
    str(part(expr_, 1))
);

//is_atomic(expr_):= begin
//  return 0;
//end;

is_call:= (expr_)→ (
  (type(expr_) == DOM_SYMBOLIC) and
  (part(expr_) == 2) and
  (part(expr_, 0) == "of") and
  (expr_[1] == 'of')
);

evalp:= (call_)→ (
  is_call(call_) ? 
    gpart(call_, 1, eval(part(call_, 1)))
);

gpart(expr_, n_=NULL, set_=NULL):= begin
  local show_size_:=size(args) < 3;
  local set_value_:=size(args) ≥ 4;
  //print(expr_, type(expr_));
  //print(n_, type(n_));
  //print(set_, type(set_));
  if type(expr_) == DOM_LIST then
    if show_size_ then
      return size(expr_); end;
    if n_ == 0 and !set_value_ then
      return str(gclear(expr_)); end;
    if 1 ≤ n_ ≤ size(expr_) then
      return expr_[n_]; end;
    return part(0, 2);  // error out
  end;
  if type(expr_) == DOM_RAT then
    if show_size_ then return 2; end;
    if n_ == 0 and !set_value_ then
      return "/"; end;
    if n_ == 1 then
      return (set_value_) ? (
                set_/denom(expr_)
             ) : ( 
                numer(expr_)
             ); end;
    if n_ == 2 then
      return (set_value_) ? (
                numer(expr_)/set_
             ) : ( 
                denom(expr_)
             ); end;
    return part(0, 2);
  end;
  if set_value_ then
    (part(expr_) ≥ n_ > 0) ? (
      return expr_[n_ + 1] := (set_)) : (
      part(0, 2)); 
  end;
  return show_size_ ?
           part(expr_) :
           part(expr_, n_);
end;

parts(expr_):= begin
  local parts_:= {};
  local result_;
  local i_:= 1;
  expr_:= head({expr_});
  while true do
    result_:= gtry2(
      part@op,
      [expr_, i_]);
    pprint(result_[1]) ? break;
    parts_[i_]:=result_[2];
    i_++;
  end;
  return part(expr_, 0), parts_;
end;

folderlib.xcas

// Folders act like any other arrays
// except for a few notable cases_
//   - contents do not get substituted
//     or evaluated, even with curry
//     substitutions
//   - can be crated using member
//     notation

ddp:= 'a::b'[1];

is_folder:= (obj_) →(
  type(obj_) == DOM_LIST and
  str(gclear(obj_)) == "folder[]" 
);

mkfolder(rows_=NULL):= begin
  local (i_:=1), row_,
        (result_:=folder[]);
  for row_ in rows_ do
    result_[i_]:= row_; 
    i_++;
  end;
  result_
end;

folder_build(elem_=seq[]):= begin
  local (state_:= folder[]),
        (me_:= folder_build);
  me_[4][2][1][1][2]:= ( 
    state_[size(state_)+1]:= elem_
  )
end;

fbuild:= folder_build;

build_get:= (builder_) →
  builder_[4][2][1][1][2];

ffget:= (builder_) →
  op(builder_[4][2][1][1][2]);

fun.xcas

// Functional programming
indices:= (sizes_, offset_=1)→ (
  (size(sizes_)<1 and 
   type(offset_) == DOM_INT) ? 
    set[] : (
    gprod(
      apply((size_)→ 
              set[op(range(size_).+
                       offset_)],
             mkfolder(sizes_)))
  )
);

// function composition
gfcomp(funs_):= begin
  local args_, result_, i_;
  if type(funs_) != DOM_LIST
     or size(funs_) < 1 then
    return eval(ret_me); end;
  if size(funs_) == 1 then
    return funs_[1]; end;
  args_:= expr(temp_var("args_"));
  result_:= gunapply(
    (sin@@(size(funs_)))(args_),
    args_);
  // substituting in reverse order
  // so we don't pull the rug from
  // under our feet
  for i_ from (size(funs_)) downto 1 do
    expr("result_[4]"+"[2]"*(i_-1)+
         "[1]:= funs_[i_]");
  end;
  return result_;
end;

reduce(fun_, vec_,
       id_=NULL):= begin
  local (accum_:= id_), elem_, args_;
  if (size((args))≤2) then
    accum_:= head(vec_);
    vec_:= tail(vec_);
  end;
  for elem_ in vec_ do
    args_:= seq[0, 0];
    args_[1]:= accum_;
    args_[2]:= elem_;
    accum_ := fun_(args_);
  end;
  return accum_;
end;

reduce_orig:= reduce;
reduce:= cache2(
  ((fun_, vars_)→(
    prog_tmpify(fun_, vars_)()))@op,
  [reduce_orig, lname(reduce_orig)]
);

// deprecated
gmzip(args_):= begin
  local (func_:= head([args_])),
        (vecs_:= tail([args_])),
        lowest_, result_, index_;
  purge(index_);
  lowest_:= min(apply(size, vecs_));
  result_:= []; 
  for index_ in 1..lowest_ do
    result_:= append(result_, 
      func_(op(apply(((v_)→v_[index_]), 
                     vecs_)))
    );
  end;
  return result_;
end;

gmzip:= vfscramble_simple(gmzip,
          ['func_', 'vecs_', 'index_',
           'lowest_', 'result_']);

// deprecated
ggmzip(args_):= begin
  local vecs_:= tail([args_]);
  local func_:= head([args]);
  if size(vecs_) < 1 then
    return eval(ret_me); end;
  return(
        reduce(append, tail(vecs_),
          apply(mkvec, head(vecs_)))); 
end;

// deprecated
enumerate:= (vec_)→ (
  zip(mkvec,
      vectorize(vec_),
      range(1,size(vectorize(vec_))+1))
);
enumerate:= 
  vfscramble_simple(enumerate);

// deprecated
reduce_init(fun_, vec_, id_=undef):= begin
  local accum_ := id_;
  local first_ := true;
  for elem_ in vec_ do
    accum_ := fun_(accum_, elem_, first_);
    first_ := false;
  end;
  return accum_;
end;

group_by(vec_, key_='id'):= begin
  local elem_, (result_:= table()),
        class_vec_, class_;
  for elem_ in vec_ do
    class_:= key_(elem_);
    class_vec_:= result_[class_]; 
    if class_vec_ in [0] then
      class_vec_:= folder[] 
    end;
    result_[class_]:= (
      class_vec_[size(class_vec_)+1]:=
        elem_);
  end;
  result_;
end;

group_by_orig:= group_by;
group_by:= cache2(
  ((fun_, vars_)→(
    prog_tmpify(fun_, vars_)()))@op,
  [group_by_orig,
   lname(group_by_orig)]
);

inclusive_scan(fun_, vec_, id_):= begin
  local result_:= gclear(vec_);
  // last, current
  local pair_:= folder[];
  local i_;
  pair_[1]:= id_;
  for i_ from 1 to size(vec_) do
    pair_[2]:= vec_[i_]; 
    pair_[1]:= fun_(op(pair_)); 
    result_[i_]:= pair_[1];
  end;
  result_
end;

inclusive_scan_orig:= inclusive_scan;
inclusive_scan:= cache2(
  ((fun_, vars_)→(
    prog_tmpify(fun_, vars_)()))@op,
  [inclusive_scan_orig,
   gminus(lname(inclusive_scan_orig),
     ['gclear'])]
);

//inclusive_scan:= 
//  vfscramble_simple(inclusive_scan, 
//                    ['result_', 'elem_']);

// unsafe
adjacent_pairs:= (vec_)→(
  size(vec_) < 1 ? gclear(vec_) :
  inclusive_scan(
    (last_,elem_)->(id({last_[2],elem_})),
    tail(vec_),
    {0,head(vec_)})
);

chunked(vec_, n_):= begin
  local (result_:= []), (i_:=0);
  if type(vec_) != DOM_LIST or
     type(n_) != DOM_INT or
     n_ < 1 then
    return eval(ret_me); end;
  for i_ from n_ to size(vec_) step n_ do
    result_[i_/n_]:= vec_[i_-n_+1..i_]; 
  end;
  if (i_-n_) < size(vec_) then
    result_[(i_/n_)]:=
      vec_[i_-n_+1..size(vec_)];
  end;
  return result_;
end;

// it'd be neat to have different
// functions for different stages
// like some sort of state machine
// args:
// - func, nargs
// - nth_arg
// - ...
// - NULL
fcurry(args_=NULL):= begin
  local (i_:=1), (max_:=nop),
        (table_:=table()), (func_:=id),
        (self_:=fcurry);
  local (teq_:='1=1');
  if type(max_) != type(1) then
    if !is_seq(args_) or
       size(args_) != 2 then
      return eval(ret_me);
    end; 
    self_[4][2][1][2][2]:= args_[2];
    self_[4][2][1][4][2]:= args_[1];
    return self_;
  end;
  if i_ ≤ max_ then
    teq_[2]:= i_;
    self_[4][2][1][3][2]:=
      table_ + table(teq_[3]:=args_);
    self_[4][2][1][1][2]:= i_+1;
    self_[4][2][1][2][2]:= max_;
    self_[4][2][1][4][2]:= func_;
    return self_;
  end;
  if !is_null(args_) then
    throw(
      "fcurry: expected no arguments " +
      "to confirm end of the table."); 
  end;
  return (func_@top)(table_); 
end;

fcurry2(fun_, nargs_=NULL):= begin
  local fcopy_;
  if type(nargs_) == DOM_INT then
    if nargs_<1 then eval(ret_me); end;
    fcopy_:= fcurry2_sized_currier;
    fcopy_[4][2][1][1][2]:= fun_;
    fcopy_[4][2][1][5][2]:= nargs_;
    return fcopy_;
  end;
  fcopy_:= fcurry2_free_currier;
  fcopy_[4][2][1][1][2]:= fun_;
  return fcopy_;
end;

make_curre(NULL):= begin
  local curried_:= 'curried_';
  local i_:= 'curried_[4][2][1][4][2]';
  local f_:= 'curried_[4][2][1][1][2]';
  local args_:= 'curried_[4][2][1][2][2]';
  local sns_:= 'seq_nest_symb';
  (curried_)→ (
    (i_<2) ? (
      f_(sns_('op(args_)'))
    ) : (f_(op(args_)))
  ) 
end;

fcurre:= make_curre();

fcurry2_sized_currier(arg_=NULL):= begin
  local func_:= id;
  local args_:= folder[];
  local base_:= fcurry2_sized_currier;
  local i_:= 1;
  local max_:= 0;
  // edit below this line
  if i_>max_ then
    if is_null(arg_) then
      if i_≤2 then
        args_:= op(args_);
        return func_(
          seq_nest_symb('args_')());
      end; 
      return func_(op(args_));
    end;
    throw(
      "fcurry2: expected no arguments " +
      "to confirm end of the call."); 
  end;
  base_[4][2][1][1][2]:=func_;
  base_[4][2][1][5][2]:=max_;
  args_[i_]:= arg_;
  base_[4][2][1][2][2]:= args_;
  base_[4][2][1][4][2]:= i_+1;
  return base_;
end;

// WIP
fcurry2_auto_currier(arg_=NULL):= begin
  local func_:= id;
  local args_:= folder[];
  local base_:= fcurry2_auto_currier;
  local i_:= 1;
  local max_:= 1;
  // edit below this line
  args_[i_]:= arg_;
  if i_++>max_ then
    return func_(op(args_)); 
  end;
  base_[4][2][1][1][2]:= func_;
  base_[4][2][1][2][2]:= args_;
  base_[4][2][1][4][2]:= i_;
  base_[4][2][1][5][2]:= max_;
  return base_; 
end;

fcurry2_free_currier(arg_=NULL):= begin
  local func_:= id;
  local args_:= folder[];
  local base_:= fcurry2_free_currier;
  local i_:= 1;
  base_[4][2][1][1][2]:= func_;
  // edit below this tine
  args_[i_]:= arg_;
  base_[4][2][1][2][2]:= args_;
  base_[4][2][1][4][2]:= i_+1;
  return base_;
end;

funcurry(fun_):= begin
  local uncurrier_:= funcurry_args; 
  uncurrier_[4][2][1][1][2]:= fun_
end;

funcurry_args(args_):= begin
  local accum_:= id;
  for arg_ in args_ do
    accum_:= accum_(arg_);
  end;
  accum_
end;

teed_inner(arg_):= begin
  nop(arg_);
  return arg_;
end;

// basically an identity function
// with a side-effecty function
// called on the input
teed(fun_):= begin
  local copy_:= teed_inner; 
  return copy_[4][1][1]:= fun_;
end;

// ugly af
lincomb(operator_, mat1_, mat2_):= begin
  local result_, index1_, index2_, temp_;
  // transposing twice will turn any
  // horizontal vectors into matrices
  mat1_ := (transpose@@2)(mat1_);
  mat2_ := (transpose@@2)(mat2_);
  result_ := [];
  for index1_ from 1 to dim(mat1_)[1] do
    temp_ := [];
    for index2_ from 1 to dim(mat2_)[2] do
      temp_[index2_] := 
        operator_(mat1_[index1_],
                  transpose(mat2_)[index2_]);
    end;
    result_[index1_] := temp_;
  end;
  return result_;
end;

functools.xcas

indices:= (vec_)→ (
  op(apply((id_)→(id_-1),
           vectorize(vec_)))
);

// WIP
vfscramble_simple(fun_, extras_=[]):= begin
  local (substs_:= {}), varn_,
        (gbsubst:=('a|b'[1])), i_;
  if !(is_program(fun_)) or
     type(extras_)!=DOM_LIST then
    return eval(ret_me); end;
  for i_ from 1 to size(fun_[2]) do
    if type(fun_[2][i_]) != DOM_IDENT then
      return eval(ret_me)
    end;
    varn_:= expr(rnd_name(
          str(fun_[2][i_])+"_")); 
    substs_:= append(substs_,
      (fun_[2][i_] = varn_));
    fun_[2][i_]:= varn_;
  end;
  for i_ from 1 to size(extras_) do
    if type(extras_[i_]) != DOM_IDENT
      then return eval(ret_me) end;
    varn_:= expr(rnd_name(
          str(extras_[i_])+"_")); 
    substs_:= append(substs_,
      (extras_[i_] = varn_));
  end; 
  return gbsubst(fun_, substs_);
end;

// sometimes evaluates the contents of
// program(...), we'll need an 
// expression walker for this
vfscramble_subst:= (f_)→
  gbsubst(f_,
    program = fcomp(vfscramble_simple,
                     program, op));

glocal:=((args)-> begin local a; return a; end)[4][1];

// Calls a function passing a veriable number
// of arguments, contained in the iterable ``args_``
vacall(func_, args_):= begin
  local args_sum_, i_;
  local args_sym_;
  purge(args_sym_);
  args_sum_ := 0;
  for i_ from 1 to size(args_) do
    args_sum_ := args_sum_ + args_sym_[i_];
  end;
  args_sum_[1] := func_;
  args_sym_ := args_;
  return eval(args_sum_);
end;

// Evaluates a function where
// every argument is assigned
// its symbolic value
// (no default values allowed) 
sym_env(func_):= begin
  local params_:= func_[2];
  if !is_program(func_) then
    warn("sym_env: func_ is not a " +
         "program.")
    return eval(ret_me);
  end;
  return tcomp(2)(
   func_(params_))(
   [params_])();
end;

sym_env2(func_):= begin
  local params_, f_;
  params_:= mkfolder([func_[2]]);
  f_:= (NULL)→ begin
    purge(params_);
    func_(
      op(apply(
        (var_)→expr(str(var_)),
        params_
    )))
  end;
  f_[2]:= func_[2]
end;

tcall(func_, args_):= begin
  local args_sum_, (i_:= 1);
  local args_sym_, keys_;
  if !is_table(args_) then
    return eval(ret_me); end;
  purge(args_sym_);
  args_sum_ := 0;
  keys_:= tkeys(args_);
  while i_ in keys_ do
    args_sum_ += args_sym_[i_];
    i_++;
  end;
  if i_ ≤ 1 then
    return func_(); end;
  if i_ ≤ 2 then
    args_sum_ := inv(args_sym_[1]); end;
  args_sum_[1] := func_;
  args_sym_ := args_;
  return eval(args_sum_);
end;

twrap:=(func_)→(
  (args_)→(
    (('tnop')@('tcall'))(func_,args_)
  )
);

rnd_name:= (prefix_="temp_")-> (
  prefix_ + str(randint(exact(1e9)))
);

// WIP
// Variables 'wrapper__' and 'in__' are declared
// by this function, use different names for
// your variables in func_, args_ and vars_
// (func_, vars_, args_='auto')
vfcapture(in__):= begin
  local wrapper__;
  in__:=table("args" = [in__]);
  // Input sanitization
  if size(in__["args"]) < 2 or
     size(in__["args"]) > 3 or
     type(in__["args"][1]) != DOM_FUNC or
     gor(op(apply((x_)→(type(x_) != DOM_IDENT),
              in__["args"][2]))) then
    return eval(ret_me);
  end;
  if size(in__["args"]) < 3 then
    if is_program(in__["args"][1]) then
      in__["args"]:=append(in__["args"], [in__["args"][1][2]]);
    else return eval(ret_me); end;
  else
    if gor(op(apply((x_)→(type(x_) != DOM_IDENT),
                in__["args"][3]))) then
    return eval(ret_me); end;
  end;
  // Build the thing
  in__["callee"]:= rnd_name();
  in__["assignments"]:= size(in__["args"][2]) > 0 ? ("local " + (
      reduce_init((accum_, e_, first_)->(first_ ? e_ : (accum_ + ", " + e_)),
        apply((e_)->(
          "(" + str(e_) + ":=0)"
        ), in__["args"][2])
      )
    ) + ";") : "";
  in__["wrapper_args"]:= size(in__["args"][3]) > 0 ? 
    str(op(in__["args"][3])) : "NULL";
  in__["wrapper_params"]:= str(op(in__["args"][3]));
  in__["wrapper_str"]:= ("\
  \
  (" + in__["wrapper_args"] + ")-> begin
    local " + in__["callee"] + ":=0;
    " + in__["assignments"] + "
    return (" + in__["callee"] + ")(" + in__["wrapper_params"] + ");
  end;");
  // Do the thing
  wrapper__:= expr(in__["wrapper_str"])
  wrapper__[4,2,1,1,2]:=in__["args"][1];
  in__["i"]:=1;
  while (in__["i"] <= size(in__["args"][2])) do
    wrapper__[4,2,1,in__["i"]+1,2]:=
      eval(eval(in__["args"][2][in__["i"]]));
    wrapper__[4,2,1,in__["i"]+1,3]:=
      in__["args"][2][in__["i"]];
    (in__["i"])++;
  end;
  return wrapper__;
end;

partial(func_, vars_):= begin 
  local wrapper_, var1_, var2_;
  purge(var1_, var2_);
  wrapper_(args_=NULL):= begin
    local (var1_:=vars_), (var2_:=func_);
    return var2_(var1_, args_);
  end;
  return wrapper_;
end;

xyfsubst(expr_, x_, y_):= begin
  return pfsubst(
    expr_, 
    [x_, y_], ['X', 'Y']);
  //return ((X, Y)→(
  //  gprogram(
  //    [X, Y],
  //    [0, 0],
  //    subst(expr_, {x_=X, y_=Y}))
  //))('X', 'Y');
end;

// Creates a function using expr_ with
// the variables in orig_ substituted
// by the ones in substs_ as the body, 
// in an eval safe environment for
// said variables, all this because
// hp decided to use the same variables
// for graphing and cursor position
pfsubst(expr_, orig_, substs_):= begin
  if size(substs_) != size(orig_) or
    type(orig_) != DOM_LIST or
    type(substs_) != DOM_LIST or 
    size(orig_) < 1 or
    sum(
      apply(
        (x_)→(type(x_) != DOM_IDENT),
        concat({}, orig_, substs_)
      )) > 0 then return eval(ret_me); end;
  
  orig_, substs_:= op(transpose(
    zip((x_,y_)→{x_, y_}, 
        orig_, substs_)))
  return expr("("+str(substs_)+"→(
    purge("+str(substs_)[2..-1]+"),
    gprogram(
      " + str(substs_) + ",
      " + str([(0)$(k=1..size(substs_))]) + ",
      subst(expr_, orig_=substs_))
  )[0])(" + 
      reduce((str_, e_)→
          (str_ + ", quote(" + str(e_)) + ")",
        tail(substs_), 
        "quote(" + str(head(substs_)) + ")") 
       + ")");
end;

genlib.xcas

// func_= (n_, prev_, op(args))
// 
generate(func_, l_, args_={},
         prev_='nop', n_=1):= begin 
  local cont_, value_;
  local next_me:='((args)[1])(
               func_, l_, args_, value_, n_+1)';
  if !(is_program(func_) and 
     ( (size(func_[2]) == 
          size(args_) + 2) or 
       (size(func_[2]) ≤ 1)) and
     type(args_) == DOM_LIST and
     type(l_) == DOM_INT and
     type(n_) == DOM_INT) then
    return eval(ret_me) end;

  (cont_, value_):= 
    gen_cont((x_)→func_(n_, prev_, op(args_)));
  if !(cont_) then (return NULL) end;
  if n_ == l_ then
    return value_;
  else
    if n_ < l_ then
      return value_,
             eval(next_me);
    else
      return value_,
             eval(next_me); 
    end;
  end; 
end;

gen_end:= (NULL)→(
  throw('gen_end__error')
);

is_gen:= (gen_)→(
  part(gen_, 0) == "of" and
  part(gen_, 1) == 'generate' 
);

gen_pack:= (gen_)→ (
  is_gen({gen_}[0]) ? (
    {gen_}[0], 
    reverse(tail(reverse({gen_})))
  ) : ({},{gen_})
);

gen_collect(gen_):= begin
  local next_, prevs_;
  (next_, prevs_):= gen_pack(gen_);
  while is_gen(next_) do
    (next_, prevs_):= 
      gen_pack(op(prevs_) ,eval(next_));
  end;
  return prevs_;
end;

gen_cont(next_):= begin
  local res_;
  res_:= gtry(next_, 0);
  if part(res_, 0) == "of" then
    return eval(ret_me)
  end;
  if res_[1] then
    if (gtry((x_)->gen_end(),0)[3] == 
       res_[3]) then
      return false, nop;
    else
      throw(res_[3]);
    end; 
  else
    return true, res_[2];
  end; 
end;

glib.xcas

// This library contains functions and
// symbols declared "sane" or "good"
// alternatives to the built-in ones
// by me after using CAS for a long
// time

// most ef these do not work correctly
// when "Recursive Replacement" is set
// over 1 on CAS settings


// Unaccesible operators made usable
// as plain functions, all arguments
// will be evaluated

newline:="
";
throw:='error';
gfor:=('for a in [1] do print(a); END')[1];
// fails,error_var,do_if_true,do_if_false
gtry_catch:=('IFERR nop THEN nop ELSE nop END')[1];
glocal:=((args)-> begin local a; return a; end)[4][1];
gprogram:='program';
gquote:=quote;
gpiecewise:='piecewise';
gunapply:='unapply';
gor:='or';
gand:='and';
gof:='of';
gmod:='0 MOD 0'[1];
gbsubst:='0|0'[1];
gne:='0≠0'[1];
geq:='=';
gprod:='*';
gassume:='assume';
gminus:=('a minus b')[1];
glog:='log';
fcomp:='@';
fncomp:='@@';
gcopy:= (('testa():= begin
  local x:=[];
  return x;
end')[2][4][2][1][1][2][1]);
// I have no idea of how I came
// up with this
at:=(((x_)->(return x_))(
  [][undef]))[1];
usum:=('([a,b]) .+ b')[1];
gpurge:='purge';

gindex:=(n_=NULL)→ (
  (type(n_)==DOM_INT and
   (n_==1 or n_==0)) ?
    head([gindex(), (index:= n_)]) : (
    -op(expr("'tmp_[0]'"))[2])
);

//geq(args_= seq[]):= begin
//  local of_args_:= folder[];
//  of_args_[1]:= 'eq';
//  of_args_[2]:= args_;
//  return (gof@op)(of_args_);
//end;

gequation:= (obj_,x_='x',y_='y')→ (
  equation(obj_)({'x','y'}={x_,y_})
);

eqinvert:=(eq_,x_)→ (
  one(solve(eq_,x_)) = x_
);

invsum:=(eqs_,x_,y_)→ (
  eqinvert(sum(
    apply(
      (eq_)→one(solve(eq_, x_)),
      eqs_))=x,
  y_)
);

targs(args_):= begin
  local eq_:='1=1';
  table(eq_[3]:=(args_))
end;

warn:= (msg_)→
  head([true, 
        print("warning: " + msg_)]);

mkvec:=(
  (args_=seq[])→[args_]
);

pprint:= teed(print);

// WIP
gtry2(test_, arg_, def_=0):= begin
  // {"error?", {"ret", ..}, "msg"}
  local err_msg_:= 0;
  local result_:= 
    gtry_catch('def_:= test_(arg_)',
               'err_msg_',
               true,false);
  return tcomp(3)(result_)(
    {def_})(err_msg_)();
end;

gtry(test_, arg_, def_=0):= begin
  local (result_:= folder[true,0,0]),
        (error_msg_:= 0);
  result_[2]:= def_;
  gtry_catch('result_[2]:= test_(arg_)',
             'error_msg_', 'nop',
             'result_[1]:= false');
  result_[3]:= error_msg_;
  result_
end;

gtry_orig:= gtry;
gtry:= cache2(
  ((fun_, vars_)→(
    prog_tmpify(fun_, vars_)()))@op,
  [gtry_orig,
   gminus(lname(gtry_orig),
          ['gtry_catch'])]
);

// Misc

gint:='int';

// Does not work correctly for
// fuction calls with no arguments.
// With functions that take lambdas
// as arguments, they will repr as
// mangled, but in practice they 
// should still be usable
ret_me:='((args)[1] != undef) ? 
((args)[1])(
  op(tail(args))
) : throw("ret_me: cannot return " +
          "self on temp function")';
me:=(NULL)→(eval(ret_me));

is_program:= (prog_)->(
  is_func(prog_) and
  size({op(prog_)}) == 3 and
  part(prog_, 0) == "program"
);

is_eq_expr:= (expr_)→ (
  type(expr_) == type(1 = 1) and
  size([op(expr_)]) == 2 and
  expr_[1] == '='
);

is_func:= (func_)→(
  type(func_) == DOM_FUNC
);

temp_var:= (prefix_="temp_")-> (
  prefix_ + str(randint(exact(1e9)))
);

expr_size:= (expr_)→(
  head(expr_)
);

// Numeric
realn:= (expr_)→(
  gor(apply((type_)→(type(expr_)==type_),
  [type(1), type(1/2), type(1.2)]))
);

engn:= (num_)→(
  (realn(num_)) ? (
    {mant(num_), xpon(num_)}
  ) : ((has_unit(num_) and 
        realn(guval(num_))) ? (
    {mant(guval(num_))*gupart(num_), 
      xpon(guval(num_))} 
  ) : (eval(ret_me)))
);

late_mod:= (v1_, v2_)→ (
  realn(v1_) and realn(v2_) ? (
    v1_ MOD v2_) : 
    eval(ret_me)
);

// Vector helper functions

is_seq:= (args_=NULL)→(
  (type(args_) == type([])) and
  (size([args_])!=1)
);

is_null:= (expr_=NULL)→(
  is_seq(expr_) and
  size([expr_]) == 0
);

seq_nest(expr_, n_):= begin
  local (temp_expr_:= ((seq[])→ 0)),
         i_, (result_:='1=1');
  if type(n_) != type(1) or 
     n_ < 0 then
    return eval(ret_me); end;
  for i_ from 0 to n_-1 do
    expr("temp_expr_[3]"+"[1]"*i_+
         ":=seq[1]");
  end;
  expr("temp_expr_[3]"+"[1]"*n_+
       ":=expr_");
  return table(result_[3]:=
                (op(temp_expr_)[2]));
end;

vectorize:= (vec_)→(
  (is_seq(vec_) or
  type(vec_) != DOM_LIST) ? 
    ([vec_]) :
    (vec_)
);

one:= (vec_)→ (
  size(vec_) != 1 ? 
    throw("one: expected only one " +
          "element, got " + size(vec_)):
    head(vec_)
);

gclear:= (vec_=NULL)→ (
  // I don't have the documentation of ``clear``
  // on hand, but I guess it does what we need
  size(vec_) < 1 ? 
    vec_ : 
    clear(vec_)
);

atolist:= (args_) → {args_};
atovec := (args_) → [args_];
pair:= (e1_, e2_) → {e1_, e2_};

// Same as apply but without the list->vector
// conversion
gapply(fun_, vec_):= begin
  return concat(gclear(vec_),
                apply(fun_, vec_));
end;

gzip(fun_, v1_, v2_):= begin
  return concat(gclear(v1_),
                zip(fun_, v1_, v2_));
end;

// Very clean hack imo
gsorted:= (vec_, key_='id')→ (
  (type(vec_) == DOM_LIST and
   type(key_) == DOM_FUNC) ? (
    reduce((accum_, e_)→(
        append(accum_, e_[3])),
      sort(
        {seq({key_(vec_[i_]), i_, vec_[i_]},
             i_=1..size(vec_))}
      ),
      gclear(vec_)
    )
  ) : eval(ret_me)
);

dapply(fun_, vec_, depth_=1):= begin
  local index_;
  local myself_ := args[1];
  if depth_ <= 0 then
    // At depth 0 (dimention 0) 
    // we expect single elements
    return function_(vec_);
  end;
  return gapply(
    (v_)->myself_(fun_, v_, depth_ - 1), 
    vec_);
end;

last:= (vec_, id_='NULL')→(
  type(vec_) == DOM_LIST ? (
    size(vec_) < 1 ? id_ : vec_[0]
  ) : eval(ret_me)
);

// Sets

mkset:= (elems_)→ set[op(elems_)];

// Ranges

range_op:= (0..0)[1];

is_range:=(expr_)→(
  type(expr_) == DOM_SYMBOLIC and
  size({op(expr_)}) == 2 and 
  expr_[1] == range_op
);

//rdist:= (range_)→(
//  
//)

// linear interpolation
rinterp:= (range_, prop_)→(
  is_range(range_) ? (
    left(range_)+
    (right(range_) - left(range_))*prop_
  ) : eval(ret_me) 
); 

rcentre:= (range_)→ rinterp(range_, 1/2);

in_range:= (num_, range_, inc_=true)→(
  is_range(range_) ? (
    inc_ ? ( 
      left(range_) ≤ num_ ≤ right(range_)
      ) : (
      left(range_) < num_ < right(range_) 
    ) 
  ) : eval(ret_me)
);

rtolist(range_):= begin
  return {left(range_), right(range_)};
end;

// Iterables

cenumerate(iterable_):= begin
  local (elements_:= folder[]),
        element_, (i_:= 1), pair_;
  for element_ in iterable_ do
    pair_:= folder[];
    elements_[i_]:= (pair_[2]:= element_);
    i_++;
  end;
  elements_
end;

// Tables

is_table:= (table_)→ (
  type(table_) == DOM_MAP
);

tkeys(table_):= begin
  local key_, el_, keys_;
  keys_ := folder[];
  if !is_table(table_) then
    return eval(ret_me);
  end;
  for el_ in cenumerate(table_) do
    keys_[el_[1]]:= el_[2] .+ 1;
  end;
  return keys_;
end;

tvalues(table_):= begin
  local el_, values_;
  values_ := folder[];
  if !is_table(table_) then
    return eval(ret_me);
  end;
  for el_ in cenumerate(tkeys(table_)) do
    values_[el_[1]] := table_[el_[2]];
  end;
  return values_; 
end;

tsize(table_):= begin
  local (i_:=0), e_;
  if !is_table(table_) then
    return eval(ret_me); end;
  for e_ in table_ do
    i_++;
  end;
  return i_; 
end;

table2list(table_):= begin
  return tkeys(table_) = tvalues(table_);
end;

ttrans(table_):= begin
  return table(
    tvalues(table_) = tkeys(table_));
end;

// also sets for values like seq[0, 0]
//tset(table_, key_, value_=NULL):= begin
//  local expr_:='1=1';
//  expr_[2]:= key_;
//  return table(op(
//           table2list(table_)),
//           (expr_[3]:=value_)); 
//end;

// Updates a table dest_ with indices
// in orig_. Like tset, this does not
// remove "null" elements
tup(dest_, orig_):= begin
  return table(op(table2list(dest_)),
               op(table2list(orig_)))
end;

// accumulates a fixed amount of
// elements into a table,
// ofthen called "currying"
tcomp(args_=NULL):= begin
  local (i_:=1), (max_:=0),
        (table_:=nop);
  local (self_:= tcomp), (teq_:= '1=1');
  if !is_table(table_) then
    if type(args_) == type(1) 
       and args_ > 0 then
      self_[4][2][1][2][2]:= args_;
      self_[4][2][1][3][2]:= table();
      return self_;
    end;
    warn("tcomp: expected number of " +
         "arguments, got: " + str(args_))
    return eval(ret_me);
  end;
  if i_ ≤ max_ then
    teq_[2]:= i_;
    self_[4][2][1][3][2]:= 
      table_ + table(teq_[3]:=args_);
    self_[4][2][1][1][2]:= i_+1;
    self_[4][2][1][2][2]:= max_;
    return self_;
  end;
  if !is_null(args_) then
    throw(
      "tcomp: expected no arguments " +
      "to confirm end of the table.");
    return;
  end;
  return table_;
end;

list2table(list_):= begin
  local result_, expr_, i_;
  result_:= table(); 
  if type(list_) != DOM_LIST then
    return eval(ret_me); end;
  for i_ in 1..size(list_) do
    expr_:= '1=1';
    expr_[2]:= i_;
    result_+= table(
      expr_[3]:= list_[i_]);
  end;
  return result_;
end;

test_list2table(func_=((x,y)→x+y)):= begin
  return (list2table@op)(func_);
end;

tnop:= gfcomp('list2table', 'mkvec');

top(table_):= begin
  local (i_:= 1), keys_,
        (result_:= []);
  if type(table_) != DOM_MAP then
    return eval(ret_me); end;
  keys_:= tkeys(table_);
  while i_ in keys_ do
    result_[i_]:= table_[i_];
    i_++; end;
  if i_ ≤ 2 and is_seq(table_[1]) then
    warn("top: single element sequences reduce to one element");
  end;
  return op(result_);
end;

probz.xcas

vec_acc:= (vec_) →
  inclusive_scan('+', vec_, 0);

count_unique:= (vec_) →
  reduce((tab_,y_)→tab_[y_]++,vec_,table());

group_adjacent_if:= (vec_, cond_, fun_='id')->
  concat(gclear(vec_), op(gapply(
    (x_)->((x_[1]) ? {fun_(x_[2])} : x_[2]),
    group_by(vec_,cond_)
  )));

hypergeometric:= (x_, N_, k_, n_) →
  comb(k_, x_)*comb(N_-k_, n_-x_)/
    comb(N_, n_);

exp_pf:= (l_, x_) →
  l_*e^(-l*x_);

cuantile_interp(list_, coef_):= begin
  //local pos_ := coef_ * size(list_) + 1;
  local pos_, lower_, upper_;
  pos_ := coef_ * (size(list_) - 1) + 1;
  list_ := append(list_, list_[0]);
  lower_ := list_[floor(pos_)];
  upper_ := list_[ceiling(pos_)];
  return (upper_ - lower_) * (pos_ - floor(pos_)) + lower_;
end;

pvar1:=(X_,S_,N_,err_)->
  ((X_-t_*S_/(√N_)) .. (X_+t_*S_/(√N_))
    )(t_ = (student_icdf(N_,1-err_/2)));

pvar2:=(X_,s_,N_,err_)->
  ((X_-z_*s_/(√N_)) .. (X_+z_*s_/(√N_))
    )(z_ = (normald_icdf(1-err_/2)));

pvar3:=(p_,n_,err_)->
  ((p_-z_*sqrt(p_*(1-p_)/n_)) .. (p_+z_*sqrt(p_*(1-p_)/n_))
    )(z_ = (normald_icdf(1-err_/2)));

pvar4:=(s1_,s2_,n1_,n2_,err_)->
  ((fisher_icdf(n1_-1,n2_-1,1-err_/2)*s2_^2/s1_^2)^((-1)/2)) .. 
  ((fisher_icdf(n1_-1,n2_-1,err_/2)*s2_^2/s1_^2)^((-1)/2));

interval2diff:=(in_)->
  (right(in_) - left(in_))/2;

pvar01:= [
  ['nop',    "no rechazar", "rechazar"],
  ["cierta", "1-alpha",     "alpha"   ],
  ["falsa",  "beta (II)",   "1-beta"  ] 
];

// observed, expected
pvarU:=(O_, e_)→(((O_-e_)^2)/e_);

// cells, est_pars, signif
pvarXb:=(k_, r_, a_)→
  (CHISQUARE_ICDF(k_-r_-1,1-a_));

puvlib.xcas

ep0:=(8.854187817ᴇ−12_(F/m));
epr_si:=(11.68);
qe:=(1.6021766208ᴇ−19_C);
ni_si:=(1.45ᴇ10_(cm^(-3)));
kb:=(1.38064852ᴇ−23_(J/K));
nu_si:=2;
nu_ge:=1;

puvlib01:= (q1_, r_) →
  (1/(4*π*ep0))*(q1_)/((r_)^2); 

puvlib02:= (q1_, q2_, r_) →
  puvlib01(q1_, r_)*q2_;

p2c:= (length_,angle_) →
  [COS(angle_), SIN(angle_)] * length_;

vnorm:=(v_)->v_/(√(v_*v_));

lrc_c:= (e0_,l_,r_,w_,c_,t_) → 
  p2c(e0_/(√(r_^2+(w_*l_-1/(w_*c_))^2)),
      w_*t_-atan(usimplify(
        limit((w_*l_-1/(w_*c_))/r__,r__,r_,1))
      ));

pw_cont:= (exprs_, lims_, var_)→ (
  (type(exprs_) == DOM_LIST and
  type(lims_) == DOM_LIST and
  size(exprs_) - 1 == (size(lims_)) and
  type(var_) == DOM_IDENT) ? (
    prepend(zip(
      (f_, l_)→ (
        f_[2]-
        subst(f_[2], var_=l_)+
        subst(f_[1], var_=l_)
      ),
      adjacent_pairs(exprs_), lims_)
    , head(exprs_))
  ) : eval(ret_me)
);

phys_Lh:= (Dh_, th_)→ ((Dh_*th_)^(1/2));

phys_J:= (muh_, mue_, n_, p_, E_)→ (
  (n_*mue_+p_*muh_)*qe*E_
);

phys_Je0:= (De_, npo_, Le_, V_, T_)→ (
  ((qe*De_*npo_)/Le_)*(e^((qe*V_)/(kb*T_))-1)
);

phys_J0:= (ni_, De_, Le_, Na_, Dh, Lh_, Nd_)→ (
  qe*(ni_)^2*(De_/(Le_*Na_)+Dh_/(Lh_*Nd_))
);

phys_dJ:= (J0_, V_, T_)→ (
  J0_*(e^((qe*V_)/(kb*T_))-1) 
);

phys_dI:= (I0_, V_, T_, nu_)→ (
  I0_*(exp((q*V_)/(nu_*kb*T_))-1) 
);

phys_vn:= (Nd_, ep_, x_, xn_)→ (
  ((qe*Nd_)/(ep_))*(x_-xn_)
);

phys_En:= (Nd_, ep_, x_, xn_)→ (
  ((qe*Nd_)/(ep_))*(x_-xn_)
);

phys_V0:= (T_, Na_, Nd_, ni_)→ (
  ((T_*kb)/qe)*ln((Na_*Nd_)/(ni_^2))
);

phys_Vn:= (Nd_, ep_, xn_, x_)→ (
  ((qe*Nd_)/(2*ep_))*((x_)^2-2*xn_*x_)
);

phys_xn:= (V0_, ep_, Na_, Nd_, qe_=qe)→ (
  √((2*ep_*V0_/(qe_*Nd_))*(Na_/(Nd_+Na_))) 
);

phys_xp:= (V0_, ep_, Na_, Nd_, qe_=qe)→ (
  √((2*ep_*V0_/(qe_*Na_))*(Nd_/(Nd_+Na_))) 
);

phys_Rh:= (Ey_, Bz_, Jx_)→ (
  Ey_/(Bz_*Jz_)
);

phys_Rhnp:= (np0_, t_='n')→ (
  1/(qe*np0_)*(t_=='n' ? -1 : 1)
);

phys_muhn:= (cond_, Rn_)→ (
  -cond_*Rn_ 
);

signsys.xcas

eodec:= (f_,t_)->( 
  [(f_(t_)+f_(-t_))/2,
   (f_(t_)-f_(-t_))/2]
);

defv:= (vec_, def_)→(
  (idx_)→(
    1≤idx_≤(size(vec_)) ? (
      vec_[idx_]) : (def_)
  )
);

convsum:= (s1_, s2_)→(
  zip((i_, s_)→(
        mid(s1_, i_, s_)*
        revlist(mid(s2_, i_, s_))
      ),((size_)->(
          {seq(1,k_ = (1 .. (size_-1))),
           seq(1 .. size_)},
          {seq(1 .. (size_-1)),
           seq(size_,k_ = (1 .. size_))})
        )(dim([s1_, s2_])[2])
  )
);

convint:= (f1_,
           f2_, 
           var_,
           lim_=(-∞)..∞,
           tau_=τ)→(
  (((f1_)(var_=tau_))*
       ((f2_)(var_=(var_-tau_))),
      tau_, left(lim_), right(lim_))
);

convlt:= (f1_,
          f2_,
          var_,
          ll_,
          s_=s)→(
  gilaplace( 
    reduce('*',
      apply((f_)→(
              glaplace(
                ((f_)(var_=var_+ll_))*
                  hs(var_),
                var_, s_)
            ),
            vectorize(f1_, f2_)),1),
  s_, var_)(var_=var_-ll_*2)*hs(var_)
);

laplaced:= (f_, t_, s_, ord_)→ (
  (s_^ord_)*glaplace(f_(t_), t_, s_)-
  sum(
    zip((n_, d_)→(
          s_^(n_-1)*
          (fncomp('function_diff', d_))(
            gunapply(f_(t_), t_))(0)
        ),
        range(ord_,0,-1),
        range(ord_)
    )
  )
);

laplacedc:= (f_, t_, s_, conds_)→ (
  (s_^size(conds_))*
  glaplace(f_(t_), t_, s_)-
  sum(zip((n_, fdn_)→(s_^(n_-1)*fdn_),
          range(size(conds_),0,-1),
          conds_))
);

// Discards assumptions leaked by
// 'laplace'
glaplace(expr_, t_=t, s_=s):= begin
  local conds_:= apply(about, [t_, s_]);
  local result_:= laplace(expr_, t_, s_);
  apply(purge, [t_, s_]);
  zip('sto', conds_, [t_, s_]);
  return result_;
end;

gilaplace(expr_, s_=s, t_=t):= begin
  local conds_:= apply(about, [t_, s_]);
  local result_:= ilaplace(expr_, s_, t_);
  apply(purge, [t_, s_]);
  zip('sto', conds_, [t_, s_]);
  return result_;
end;

laplacep:= (expr_, t_, s_, T_)→ (
  glaplace(expr_*hs(t_-T_), t_, s_)/
    (1-e^(-s_*T_))
);

sinc:= (arg_)→ (
  piecewise(eq(arg_, 0), 1, sin(arg_)/arg_)
);

fourier_cns:= (expr_, t_, T_, n_, a_)→ (
  (fourier_cn(expr_, t_, T_, n_, a_))*
  e^(*n_*(2*π/T_)*t_)
);

fourier_c(expr_, t_, T_, r_, a_):= begin
  local k_, kth_;
  purge(k_);
  kth_:= fourier_cns(
    expr_, t_, T_, k_, a_);
  return sum(kth_+(kth_)(k_=-k_), k_, 1, r_)+
         fourier_cns(expr_, t_, T_, 0, a_);
end;

hs2p:= (ll_, ul_, var_)→(
  hsp(ll_, var_)-hsp(ul_, var_)
);

hsp:= (lim_, var_)→(
  piecewise(var_<lim_,0,1)
);

plotdiscr:= (list_, origin_=1, step_=1)→(
  zip((i_, y_)→(
        point([(i_-origin_)*step_, y_])
      ), 
      [seq(1..size(list_))], list_) 
);

switch.xcas

sw_test(int_):= begin
  switch (int_) {
    case 1: {1;};
  };
end;

tablelib.xcas

// the index order the code was
// evaluated with
midx_:= gindex();

tisolate(table_, key_):= begin
  local eq_:='1=1';
  if !is_table(table_) then
    return eval(ret_me); end;
  eq_[midx_+1]:= key_;
  return table(
    eq_[midx_+2]:=table_[key_]); 
end;

tparts(table_):= begin
  local key_, eq_, index_, i_,
        (result_:=[]);
  if !is_table(table_) then
    return eval(ret_me); end;
  index_:= gindex(0);
  i_:= midx_;
  for key_ in table_ do
    eq_:= '1=1';
    eq_[midx_+1]:= key_;
    result_[i_]:=
      table(eq_[midx_+2]:=table_[key_]);
    i_++;
    //pprint(result_);
  end;
  gindex(index_);
  return result_;
end;

tkeys2(table_,
       key_offset_=gindex()):= begin
  local key_, eq_, index_, i_,
        (result_:=[]);
  if !is_table(table_) or
     type(key_offset_) != DOM_INT then
    return eval(ret_me); end;
  index_:= gindex(0);
  i_:= midx_;
  eq_:= (0, key_offset_);
  for key_ in table_ do
    eq_[midx_]:= key_;
    result_[i_]:= usum(eq_);
    i_++;
  end;
  gindex(index_);
  return result_;
end;

tselect(table_,
        keys_):= begin
  local (result_:=[]), i_, key_;
  if type(table_) != DOM_MAP or
     type(keys_) != DOM_LIST then
    return eval(ret_me); end;
  i_:= midx_;
  for key_ in keys_ do
    result_[i_]:= table_[key_];
    i_++;
  end;
  return result_;
end;

thas_key:= (table_, key_)→
  !(table_[key_] in [0]) or
  key_ in tkeys2(table_)
;

tvolatile_value(value_):= begin
  local (temp_:= table()), key_;
  temp_[1]:= value_;
  for key_ in temp_ do
    return false;
  end;
  return true;
end;

tget(table_, key_, default_=0):= begin
  local result_:= table_[key_]; 
  if result_ in [0] and 
     !(key_ in tkeys2(table_)) then
    return default_;
  end;
  return result_;
end;

tsuppress(table_, key_):= begin
  local pair_:= 0=0;
  pair_[2]:= key_;
  pair_[3]:= undef;
  table_+= table(pair_);
  gpurge('table_'[key_]);
  return table_;
end;

tsingle(key_=1):= begin
  local next_:= tsingle_work;
  if type(key_) == DOM_INT then
    next_:= tsingle_work_int
    return (
      next_[4][2][1][1][2][2]:= key_)
  end;
  next_[4][3][1][2][2]:= key_
end;

tsingle_work_any(value_):= begin
  local result_:= '0=0';
  result_[2]:= '0'; 
  table(result_[3]:= value_)
end;

// can break for keys like 'undef',
// slightly faster 
tsingle_work_int(value_):= begin
  local result_:= 0=0;
  table(result_[3]:= value_)
end;

tsingle_one:= tsingle(1);

tsingle_two:= tsingle(2);

tsingle_three:= tsingle(3);

// about 4 times slower than doing
// simple assignments, but preserves
// volatile values
tset(table_, key_, value_):= begin
  local pair_:= 0=0;
  pair_[2]:= key_;
  pair_[3]:= undef;
  table_+= table(pair_);
  gpurge('table_'[key_]);
  pair_[3]:= value_;
  return table_ + table(pair_);
end;

te.xcas

purge(te);

te::rcolor:=[
  ["black",  0,     1ᴇ0,    'nop'  ],
  ["brown",  1,     1ᴇ1,    .1/100 ],
  ["red",    2,     1ᴇ2,    .2/100 ],
  ["orange", 3,     1ᴇ3,    'nop'  ],
  ["yellow", 4,     1ᴇ4,    'nop'  ],
  ["green",  5,     1ᴇ5,    .5/100 ],
  ["blue",   6,     1ᴇ6,    .25/100],
  ["violet", 7,     'nop',  .01/100],
  ["grey",   8,     'nop',  'nop'  ],
  ["white",  9,     'nop',  'nop'  ],
  ["gold",   'nop', 1ᴇ-1,   5./100 ],
  ["silver", 'nop', 1ᴇ-2,   10./100]
];


// WIP
sres_near:= (val_, series_, opts_)→ (
  (has_unit(val_) and ( 
    sto(
      gtry((x_)→convert(x_, 1_Ohm), 
           val_, 'failed')[2],
        val_) != 'failed') and (
    series_ in transpose(te_sres)[1]
  ) and type(opts_) == DOM_INT) ? (
    val_
  ) : (
    eval(ret_me) 
  ) 
);

sres_near(val_, series_, opts_):= begin
  val_:= gtry((x_)→convert(x_, 1_Ohm), 
              val_, 'failed')[2]; 
  if !(val_ != 'failed') or
     type(opts_) != DOM_INT or
     !(series_ in transpose(te_sres)[1])
       then return eval(ret_me) end;
  val_:= ((x_)→({mant(x_), xpon(x_)}))(
    guval(val_));
 
end;

CopyVar(me, cparallel);
cparallel_eval:= (zn_)→ (
  sum(apply((zn_)→zn_^(-1),
        vectorize(zn_)))^(-1)
);


CopyVar(me, cparallel);
cparallel_eval:= (zn_)→ (
  sum(apply((zn_)→zn_^(-1),
        vectorize(zn_)))^(-1)
);

CopyVar(me, cseries);
cseries_eval:= (zn_)→ (
  sum(vectorize(zn_)) 
);

imp_eval:= (expr_)→ (
  eval(
    subst(expr_, 
      {'cseries'= 'cseries_eval',
       'cparallel'= 'cparallel_eval'}))
);

te::se1:=
[ 1.0 ];

te::se3:=
[ 1.0, 2.2, 4.7 ]
;

te::se6:=
[ 1.0, 1.5, 2.2, 3.3, 4.7, 6.8 ]
;

te::se12:=
[ 1.0, 1.2, 1.5, 1.8, 2.2, 2.7, 3.3, 3.9, 4.7,
  5.6, 6.8, 8.2 ]
;

te_sres:= [
  ["E1",   'te::se1',   'nop'],
  ["E3",   'te::se3',   .20],
  ["E6",   'te::se6',   .20],
  ["E12",  'te::se12',  .10]
];

note1:= 
"
izmin:= .05 * iz
izmed:= 1/4 or 1/3 * izmax
izmax:= Pz/Vz";

unitlib.xcas

unit_op:= ((0)_(0))[1]; 
unit_operator:=unit_op;

unorm(exprs_):= begin
  local ref_:= upart(exprs_[1]);
  local i_;
  for i_ from 2 to size(exprs_) do
    exprs_[i_]:= convert(exprs_[i_],
                         ref_);
  end;
  return exprs_;
end;

guconvert:=(expr_, unit_)→ (
  uval(expr_)*
    convert(upart(expr_), unit_) 
);

try_usimplify:= (expr_)→ (
  (
    (res_)→ (res_[1] ? 
                expr_ : 
                op(res_[2])
      )
  )(
  gtry2((x_)→usimplify(op(x_)),
        {expr_})
  )
);

has_unit:= (symbol_)→ (
  (type(symbol_) == type((0)_(0))) and
  (size({op(symbol_)}) == 2) and
  (symbol_[1] == unit_operator)
);

guval:= (symbol_)→ (
  has_unit(symbol_) ? (
    symbol_[2] 
  ) : symbol_
);

gupart:= (symbol_)→ (
  has_unit(symbol_) ? (
    symbol_[2]:=1
  ) : 1
);

gusymb:= (symbol_)→ (
  has_unit(symbol_) ? (
    symbol_[3]
  ) : 1
); 

// MKSA definitions
u_m_:= 'u_m_';
u_kg_:= 'u_kg_';
u_s_:= 'u_s_';
u_A_:= 'u_A_';

u_conversions:= table(
  ((1_(m))[3])=u_m_,
  ((1_(kg))[3])=u_kg_,
  ((1_(s))[3])=u_s_,
  ((1_(A))[3])=u_A_
);

uredef_op:= (coeff_, unit_)→ (
  coeff_*guval(
      unit_:=exact(mksa(unit_)))*
    gusymb(
      subst(unit_,
            table2list(u_conversions)))
);

uextract:= (expr_)→ (
  fcurry(subst, 2)(expr_)(
    unit_op=(uredef_op@op)
  )()
);

urestore:= (expr_)→ (
  fcurry(subst, 2)(expr_)(
    table2list(ttrans(u_conversions))
  )()
);

gusimp:= (expr_)→ (
  urestore(simplify(uextract(expr_))) 
);