(* elementary operations *)

isRelQ[x_] := MatchQ[Head[x], relHeads]

notRel[x__] := And @@ (Not[isRelQ[#]]& /@ {x})

dividedBy[x_][y_] := y/x /; notRel[x, y]

dividedBy[x_][y_Equal] := #/x& /@ y /; notRel[x]

dividedBy[x_Equal, y_Equal] := x[[1]]/y[[1]] == x[[2]]/y[[2]] 

divideBy = dividedBy;

times[x_][y_] := y x /; notRel[x, y]

times[x_][y_Equal] := #*x& /@ y /; notRel[x]

times[x_Equal, y_Equal] := x[[1]]*y[[1]] == x[[2]]*y[[2]] 

multipliedBy = times;

add := plus;

plus[x_][y_] := y + x /; notRel[x, y]

plus[x_][y_?isRelQ] := #+x& /@ y /; notRel[x]

plus[x_?isRelQ, y_?isRelQ] := 
 Module[{hx, hy, hh},
  hx = Head[x]; hy = Head[y];
  Switch[{hx, hy},
   {Equal, _}, hy[x[[1]] + y[[1]], x[[2]] + y[[2]]],
   {_, Equal}, hx[x[[1]] + y[[1]], x[[2]] + y[[2]]],
   {_hh, _hh}, If[hh =!= Unequal, hx[x[[1]] + y[[1]], x[[2]] + y[[2]]],
                   unresolved[x[[1]] + y[[1]], x[[2]] + y[[2]]]],
   {GreaterEqual, Greater} | {Greater, GreaterEqual},
     Greater[x[[1]] + y[[1]], x[[2]] + y[[2]]],
   {LessEqual, Less} | {Less, LessEqual},
     Less[x[[1]] + y[[1]], x[[2]] + y[[2]]],
   {_, _}, unresolved[x[[1]] + y[[1]], x[[2]] + y[[2]]]]]

increasedBy = plus;

minus[x_][y_] := y - x /; notRel[x, y]

minus[x_][y_?isRelQ] := #-x& /@ y /; notRel[x]

minus[x_?isRelQ, y_?isRelQ] :=
 Module[{hx, hy, hh},
  hx = Head[x]; hy = Head[y];
  Switch[{hx, hy},
   {Equal, _}, hy[x[[1]] - y[[1]], x[[2]] - y[[2]]],
   {_, Equal}, hx[x[[1]] - y[[1]], x[[2]] - y[[2]]],
   {_hh, _hh}, If[hh =!= Unequal, hx[x[[1]] - y[[1]], x[[2]] - y[[2]]], 
                   unresolved[x[[1]] - y[[1]], x[[2]] - y[[2]]]],
   {GreaterEqual, Greater} | {Greater, GreaterEqual},
     Greater[x[[1]] - y[[1]], x[[2]] - y[[2]]],
   {LessEqual, Less} | {Less, LessEqual},
     Less[x[[1]] - y[[1]], x[[2]] - y[[2]]],
   {_, _}, unresolved[x[[1]] - y[[1]], x[[2]] - y[[2]]]]]

decreasedBy = minus;

squared[x_] := x^2 /; notRel[x, y]

squared[x_?isRelQ] :=
 Module[{hx},
  hx = Head[x]; 
  Return[If[MatchQ[hx, Greater | GreaterEqual | Equal], hx, unresolved] @@ 
                                                      ({x[[1]], x[[2]]} ^2)]]

inverted[s_] := 1/s


(* relationships *)

reverse[(op_)?rhm[lhs_, rhs_]] :=
 If[!MatchQ[Head[rhs], Which | If],
     reverseOp[op][rhs, lhs],
     Print[" "]; Print["warning -- reversing conditional relationship"];
     reverseOp[op][rhs, lhs]]

reverse[If[c_, rel_]] := If @@ {c, reverse[rel]}
 
reverseOp[op_] := Switch[op, Greater, LessEqual, GreaterEqual, Less,
 Equal, Equal, Less, GreaterEqual, LessEqual, Greater, Unequal, Unequal]

solveLinear[v_][r_] :=
 Module[{h, step1, step2, step3, step4, step5, a, b, $},
  h = Head[r];
  If[Not[MatchQ[h, relHeads]] || Length[r] != 2, Return[invalid]];
  step1 = r[[1]] - r[[2]];
  step2 = step1 /. v :> $;
  step3 = step2 ~ Collect ~ $;
  If[Not[MatchQ[step3, _.$ + _.]] ||
       MatchQ[step3, _. $^_ + _.], Return[invalid]];
  step4 = workOutwards[Plus][collect[$]][step2];
  step5 = step4 /. a_. $ + b_ :> - b/a;
  h = If[TrueQ[Negative[step4 /. a_. $ + b_ :> a]],
       Switch[h, Greater, Less, GreaterEqual, LessEqual,
                 LessEqual, GreaterEqual, Less, Greater, _, h], h];
  Return[h[v, step5]]]

clearFractions[s_Equal] :=
 s // times[Times @@ Flatten[subfactor[Denominator /@ List @@ s]]]

zeroRight[s_] :=
 Module[{h},
  h = Head[s];
  Return[If[Not[MatchQ[h, relHeads]] || Length[s] != 2, invalid,
          h[s[[1]] - s[[2]], 0]]]]


(*                collection                                   *)

instancesOf[pattern_][expression_] :=
 part[expression, #]& /@ Position[expression, pattern] // Union

collect[][s_] := s

collect[p1_, pr___][s_] :=
 Module[{prList, localInstances, interim, currentInstance, i, hold, result},
print["entering collect, s = ", s];
print["p1 = ", p1];
 prList = {pr};
print["prList = ", prList];
 If[Head[s] =!= Plus || {p1, pr} === {}, Return[s]];
 localInstances = instancesOf[p1][s];
print["localInstances = ", localInstances];
 interim = s;
 If[Length[localInstances] > 0,
 Do[
print["i = ", i];
  currentInstance = localInstances[[i]];
print["currentInstance = ", currentInstance];
print["interim = ", interim];
    test = interim // toEachTermThat[matches[_. currentInstance]][hold];
    If[Length[Cases[test, hold[_]]] > 1, 
    interim =
     collectivelyToTermsThat[match[_. currentInstance]][
      (#/currentInstance &) /@ # &,
       If[prList === {}, Identity, collect @@ prList],
       # currentInstance &][interim]], {i, Length[localInstances]}]];
print["interim = ", interim];
print["prList = ", prList];
 result = If[prList === {}, Identity, collect @@ prList][interim];
print["result = ", result];
 Return[result]]

collectProductsOf := collect

collectPowersOf[v___] := collect @@ (# ^(_.) &) /@ {v}

collectByArguments[f_][s_] := 
 (collect @@ 
   Select[
    Flatten[
     If[Head[#] === Times, List @@ #, #]& /@ List @@ s, 1],
   MemberQ[FixedPointList[Head, #], f[___]^_.]&])[s]

collectFunctionsOf[v_][s_] := 
 Module[{asList, step1, step2, step3, step4, step5, step6},
  If[Head[s] =!= Plus, Return[s]];
  asList = List @@ s;
  step1 = Position[asList, v];
  step2 = DeleteCases[step1, {___, 0}];
  step3 = #[[1]]& /@ step2 // Union;
  step4 = asList[[#]]& /@ step3;
  step5 = 
   toEachFactorThat[doesNotContain[v]][ReplacePart[#, 1, {{}}]&] /@ step4;
  step6 = (collect @@ Union[step5])[s];
  Return[step6]]

collectFunctionsOf[u_, v__][s_] := 
 collectFunctionsOf[v][collectFunctionsOf[u][s]]

(* distribution *)

takeSubseq[s_][{index___}] := Head[s] @@ (part[s, #]& /@ {index})

indexList[s_] := Range[Length[s]]

distribution[qspec___][rspec___][jspec___][s_] :=
 Module[{localLength, null, jjspec, jPath, jj, qqList, qPath, qList, rrList,
         t, rList, rPath, multiplier, hold1, hold2, step1, step2, step3, step4},
print["{qspec} = ", {qspec}];
print["{rspec} = ", {rspec}];
print["{jspec} = ", {jspec}]; print["s = ", s];
 If[Head[s] =!= Times, Return[s]];
 localLength = Length[s];
print["localLength = ", localLength];
(*                                                                *)
 jjspec = If[Length[{jspec}] == 0, null, jspec];
print["jjspec = ", jjspec];
 Switch[
  {jjspec},
  {null}, jPath = Position[s, Plus, 2];
      If[jPath === {}, Return[s], jj = jPath[[1, 1]]],
  {_Integer},
      If[1 <= Abs[jspec] <= localLength,
         jj = If[# > 0, #, localLength + # +1]&[jspec], Return[s]],
  {_Plus},
      jPath = Position[s, jspec, 1];
      If[jPath === {}, Return[s], jj = jPath[[1, 1]]],
  {_}, jPath = Position[jspec /@ List @@ s, True];
       jj = If[jPath === {}, Return[s], jj = jPath[[1, 1]]]];
print["jj = ", jj];
(*                                                                *)
  qqList = DeleteCases[Range[localLength], jj];
print["qqList = ", qqList];
  Switch[
   {qspec},
   {}, qList = qqList,
   {_Symbol}, qList = Position[s, qspec, 1][[1]],
   {__Integer},
       qList = If[# > 0, #, localLength + # + 1]& /@ {qspec};
       If[Complement[qList, qqList] =!= {}, Return[s]],
   {_Times},
       qPath = Position[s, #, 1]& /@ {qspec};
       If[Cases[qPath, {}] =!= {}, Return[s], qList = #[[1]]& /@ qPath],
   {_}, qPath = Position[qspec /@ List @@ s, True];
        If[qPath =!= {}, qList = Flatten[qPath], Return[s]]];
print["qList = ", qList];
(*                                                                *)
   t = s[[jj]];
   rrList = indexList[t];
print["rrList = ", rrList];
   Switch[
    {rspec},
    {}, rList = rrList,
    {_Symbol}, rList = Position[t, rspec, 1][[1]],
    {_Integer},
        rPath = Position[t, rspec, 1];
        If[rPath =!= {}, rList = rPath[[1, 1]],
            If[MemberQ[rrList, Abs[rspec]], rList = {rspec}, Return[s]]],
    {__Integer},
        rList = {rspec};
        If[Complement[Abs /@ rList, rrList] =!= {}, Return[s]],
    {_Plus},
        rPath = Position[t, #, 1]& /@ {rspec};
        If[Cases[rPath, {}] =!= {}, Return[s], rList = #[[1]]& /@ rPath],
    {_}, rPath = Position[rspec /@ List @@ t, True];
         If[rPath =!= {}, rList = Flatten[rPath], Return[s]]];
print["rList = ", rList];
(*                                                                *)
  multiplier = takeSubseq[s][qList];
print["multiplier = ", multiplier];
  step1 = s // toFactor[jj][hold1];
print["step1 = ", step1];
  step2 = step1 // toThe[hold1][toTerms[rList][hold2]];
print["step2 = ", step2];
  step3 = step2 // 
   toThe[hold1][collectivelyToTerms[notHeadedBy[hold2]][times[multiplier]]];
print["step3 = ", step3];
  step4 = step3 // toThe[hold1][toEach[hold2][times[multiplier]]];
print["step4 = ", step4];
  Return[step4/multiplier //. {hold1 -> Identity, hold2 -> Identity}]]

leftDistribution[s_] := 
  Module[{step1, hold, pathList}, 
   If[Head[s] =!= Times, Return[s]]; 
    step1 = s // toEachFactor[hold];
    pathList = 
     (#[[1]] & /@ Position[step1, Plus, {3}]) ~ DeleteCases ~ 1;
    If[pathList === {}, Return[s]];
    Return[distribution[1][][pathList[[1]] ][s]]]

rightDistribution[s_] :=
  Module[{step1, hold, pathList}, 
   If[Head[s] =!= Times, Return[s]]; 
    step1 = s // toEachFactor[hold];
    pathList = 
     (#[[1]] & /@ Position[step1, Plus, {3}]) ~ DeleteCases ~ Length[s];
    If[pathList === {}, Return[s]];
    Return[distribution[-1][][pathList[[-1]] ][s]]]

topDistribution[s_] := 
 If[Head[Numerator[s]] === Plus && Denominator[s] =!= 1, 
   #/Denominator[s]& /@ Numerator[s], s]

expandLeft[s_] :=
 Module[{pathList},
  If[Head[s] =!= Times, Return[s]];
  pathList = Position[s, Plus, {2}];
  If[pathList === {}, Return[s]];
  Return[distribution[][][pathList[[1,1]]][s]]]

expandRight[s_] :=
 Module[{pathList},
  If[Head[s] =!= Times, Return[s]];
  pathList = Position[s, Plus, {2}];
  If[pathList === {}, Return[s]];
  Return[distribution[][][pathList[[-1,1]]][s]]]


(* factoring *)

factorOut[x_][s_] := If[Head[s] =!= Plus, s, x * (#/x& /@ s)]

factorOut[-1][s_] :=
 Module[{dummy, a, b},
  Return[Switch[Head[s],
          Times, s /. a_ b_Plus  :> -a Hold[-b] // ReleaseHold,
          Plus, -HoldForm @@ {-s}, _, s]]]

factor[s_] :=
 Module[{asList, nPair, dPair},
 If[Head[s] =!= Plus, Return[s]];
 asList = List @@ s;
 nPair = subfactor[Numerator /@ asList];
 dPair = subfactor[Denominator /@ asList];
 Return[nPair[[1]] / dPair[[1]] *
          Plus @@ Thread[#1/#2&[nPair[[2]], dPair[[2]]]]]]

subfactor[asList_List] :=
 Module[{nl2, nl3, nl4, nl5, nl6, nl7, nl8, nl9, nD, result},
 nl2 = If[Head[#] === Times, List @@ #, {#}]& /@ asList;
 nl3 = If[MatchQ[#[[1]], _?NumberQ], #, # ~ Prepend ~ 1]& /@ nl2;
 nl4 = Switch[#, _?NumberQ, #, _Power, List @@ #, _, {#, 1}]& /@ # & /@ nl3;
 nl5 = #[[1]]& /@ nl4;
 nl6 = Drop[#, 1]& /@ nl4;
 nl7 = Intersection @@ (#[[1]]& /@ # & /@ nl6);
 nl8 = Cases[Flatten[nl6, 1], {#, _}]& /@ nl7;
 nl9 = {#[[1, 1]]& /@ nl8, Min /@ (#[[2]]& /@ #& /@ nl8)} // Transpose;
 nD = (GCD @@ nl5) Times @@ (Power @@ #& /@ nl9);
 result = {nD, asList/nD};
 Return[result]]

workOutwards[h_][f_][s_] := 
 Module[{levelList, $, i},
  levelList = Union[Length /@ Position[s, h] //. {$___, 0} :> {$}] - 1;
  interim = s;
  Do[interim = to[h][level[levelList[[i]]]][f][interim],
       {i, Length[levelList], 1, -1}];
  Return[interim]]  

factorIfPossible[s_] :=
  Which[
   Head[s] === Plus && Head[Factor[s]] === Times, Factor[s],
   Head[s] === Times, Apply[Times, Map[factorIfPossible, Apply[List, s]]],
   Head[s] === Power, (factorIfPossible[s[[1]]])^s[[2]],
   True, s]

(* aggregates *)

nullRangeRule =
 op_[m_, mlo_, mhi_][v_] /; mlo > mhi :>
  Switch[op, sum, 0, prod, 1,
        orAgg, False, andAgg, True, seq, Null] 

binop[op_] :=
 Switch[op, sum, Plus, prod, Times, orAgg, Or, andAgg, And, 
         seq, Sequence, agg[_], op[[1]]]

inverseBinop[op_] :=
 Switch[op, sum, (#2 - #1)&, prod, (#2/#1)&, orAgg, Or, andAgg, And]

leftExpand[op_[n_, nlo_, nhi_][s_]] :=
  (s /. n -> nlo) ~ binop[op] ~ (op[n, nlo+1, nhi][s] /. nullRangeRule)

rightExpand[op_[n_, nlo_, nhi_][s_]] :=
  (op[n, nlo, nhi-1][s] /. nullRangeRule) ~ binop[op] ~ (s /. n -> nhi)

fullExpand[op_[n_, nlo_, nhi_][s_]] :=
 If[IntegerQ[nhi - nlo],
     binop[op] @@ Table[s, {n, nlo, nhi}], op[n, nlo, nhi][s]] 

leftExtend[op_[n_, nlo_, nhi_][s_]] :=
 (s /. n -> nlo-1) ~ 
   Switch[op, sum, (#2 - #1)&, prod, (#2/#1)&, orAgg, Or, andAgg, And] ~ 
     op[n, nlo-1, nhi][s]

rightExtend[op_[n_, nlo_, nhi_][s_]] :=
 op[n, nlo, nhi+1][s] ~ 
  Switch[op, sum, (#1 - #2)&, prod, (#1/#2)&, orAgg, Or, andAgg, And] ~ 
    (s /. n -> nhi+1)

splitTheRange[ns_][op_[n_, nlo_, nhi_][s_]] :=
 If[TrueQ[ns - nlo < 0] || TrueQ[nhi - ns < 0], 
  op[n, nlo, nhi][s],
  op[n, nlo, Expand[ns]][s] ~ binop[op] ~ op[n, Expand[ns+1], nhi][s] ]

joinTheRanges[
  op_[n_, nlo1_, nhi1_][s_] ~ inop_ ~ op_[n_, nlo2_, nhi2_][s_] /; 
    inop === binop[op] && nlo2-nhi1-1 == 0] := op[n, nlo1, nhi2][s]

grule[joinTheRanges] =
 (op_[n_, nlo1_, nhi1_][s_] ~ inop_ ~ op_[n_, nlo2_, nhi2_][s_] /;
    inop === binop[op] && nlo2-nhi1-1 == 0) :> op[n, nlo1, nhi2][s]

joinTheRanges[
 f_[seq[n_, nlo1_, nhi1_], seq[n_, nlo2_, nhi2_]][s_] /; 
    nlo2-nhi1-1 == 0] := f[seq[n, nlo1, nhi2]][s] 

splitTheOperands[op_[n_, nlo_, nhi_][s_]] :=
 If[Head[s] === binop[op], 
    op[n, nlo, nhi] /@ s, sum[n, nlo, nhi][s]]

grule[joinTheOperands] =
 Composition[Apply[RuleDelayed, #]&,
  {#[n_, nlo_, nhi_][s1_] ~ binop[#] ~ #[n_, nlo_, nhi_][s2_], 
   #[n, nlo, nhi][s1 ~ binop[#] ~ s2]}&] /@ 
    {sum, prod, orAgg, andAgg};

moveCoefficientRight[a_. sum[n_, nlo_, nhi_][s_]] := sum[n, nlo, nhi][a s]

moveCoefficientRight[a_. prod[n_, nlo_, nhi_][s_]] := 
 prod[n, nlo, nhi][a ^(1/(nhi-nlo+1)) s]

moveCoefficientLeft[sum[n_, nlo_, nhi_][s_]] := 
 Module[{coefficient, $},
 If[Head[s] =!= Times, Return[sum[n, nlo, nhi][s]]];
 coefficient = s // toEachFactorThat[contains[n]][$_ :> 1];
 Return[coefficient * sum[n, nlo, nhi][s/coefficient]]] 

moveCoefficientLeft[prod[n_, nlo_, nhi_][s_]] := 
 Module[{coefficient, $},
 If[Head[s] =!= Times, Return[prod[n, nlo, nhi][s]]];
 coefficient = s // toEachFactorThat[contains[n]][$_ :> 1];
 Return[coefficient^(nhi-nlo+1) * prod[n, nlo, nhi][s/coefficient]]]

reindex[inxnew_, itsval_][op_[n_, nlo_, nhi_][s_]] :=
 Module[{inx, inverse, floc, newlo, newhi, newhold, result, expanded},
 inverse = 
  (solveLinear[n][inx == itsval] /. inx :> inxnew // Expand)[[2]]; 
  floc = Apply[Function, {n, itsval}];
  newlo = Expand[floc[nlo]] /. grule[infinityAdd];
  newhi = Expand[floc[nhi]] /. grule[infinityAdd];
  If[newlo > newhi || newhi - newlo + nhi - nlo == 0,
      newhold = newhi; newhi = newlo; newlo = newhold];
  result = op[inxnew, newlo, newhi] [s /. n :> inverse];
  expanded = result // toEach[-(inverse)][Expand];
 Return[expanded]]

switchOrderOfSummation[
 sum[n1_, n1lo_, n1hi_][sum[n2_, n2lo_, n2hi_][s_]]] :=
  sum[n2, n2lo, n2hi][sum[n1, n1lo, n1hi][s]]

grule[infinityAdd] =
 Module[{a}, {infinity + (a_) -> infinity, -infinity + (a_) -> -infinity}];


(*   solution  *)

solve[{eqn__}, {v__}] :=
 Module[{fullRuleList, ruleList, eqnList},
  fullRuleList = Solve[{eqn}, {v}];
  If[Length[fullRuleList] == 0, Return["unsolved"]];
  ruleList = fullRuleList[[1]];
  eqnList = Map[(#[[1]] == #[[2]])&, ruleList];
  Return[eqnList]]