/* Comparison operators. They call the internal comparison routines when
 * both arguments are numbers. The value Infinity is also understood.
*/

// Undefined is a very special case as we return False for everything
1 # Undefined <  _x  <--  False;
1 # Undefined <= _x  <--  False;
1 # Undefined >  _x  <--  False;
1 # Undefined >= _x  <--  False;
1 # _x <  Undefined  <--  False;
1 # _x <= Undefined  <--  False;
1 # _x >  Undefined  <--  False;
1 # _x >= Undefined  <--  False;


// If n and m are numbers, use the standard LessThan function immediately
5 # (n_IsNumber < m_IsNumber) <-- LessThan(n-m,0);


// If n and m are symbolic after a single evaluation, see if they can be coerced in to a real-valued number. 
LocalSymbols(nNum,mNum)
[
  10 # (_n < _m)_[nNum:=N(Eval(n)); mNum:=N(Eval(m));IsNumber(nNum) And IsNumber(mNum);] <-- LessThan(nNum-mNum,0);
];

// Deal with Infinity
20 #  (Infinity < _n)_(Not(IsInfinity(n)))  <-- False;
20 #  (-Infinity < _n)_(Not(IsInfinity(n))) <-- True;
20 #  (_n < Infinity)_(Not(IsInfinity(n)))  <-- True;
20 #  (_n < -Infinity)_(Not(IsInfinity(n))) <-- False;

// Lots of known identities go here
30 # (_n1/_n2) < 0  <--  (n1 < 0) != (n2 < 0);
30 # (_n1*_n2) < 0  <--  (n1 < 0) != (n2 < 0);

// This doesn't sadly cover the case where a and b have opposite signs
30 # ((_n1+_n2) < 0)_((n1 < 0) And (n2 < 0))  <--  True;
30 # ((_n1+_n2) < 0)_((n1 > 0) And (n2 > 0))  <--  False;
30 #  _x^a_IsOdd  < 0  <--  x < 0;
30 #  _x^a_IsEven < 0  <--  False; // This is wrong for complex x

// Add other functions here!  Everything we can compare to 0 should be here.
40 # (Sqrt(_x))_(x > 0) < 0          <--  False;

40 # (Sin(_x) < 0)_(Not(IsEven(N(x/Pi))) And IsEven(N(Floor(x/Pi)))) <-- False;
40 # (Sin(_x) < 0)_(Not(IsOdd (N(x/Pi))) And IsOdd (N(Floor(x/Pi)))) <-- True;

40 # Cos(_x) < 0 <-- Sin(Pi/2-x) < 0;

40 # (Tan(_x) < 0)_(Not(IsEven(N(2*x/Pi))) And IsEven(N(Floor(2*x/Pi)))) <-- False;
40 # (Tan(_x) < 0)_(Not(IsOdd (N(2*x/Pi))) And IsOdd (N(Floor(2*x/Pi)))) <-- True;

// Functions that need special treatment with more than one of the comparison
// operators as they always return true or false.  For these we must define
// both the `<' and `>=' operators.
40 # (Complex(_a,_b) <  0)_(b!=0) <--  False;
40 # (Complex(_a,_b) >= 0)_(b!=0) <--  False;
40 # (Sqrt(_x))_(x < 0) <  0      <--  False;
40 # (Sqrt(_x))_(x < 0) >= 0      <--  False;

// Deal with negated terms
50 # -(_x) < 0 <-- Not((x<0) Or (x=0));

// Define each of {>,<=,>=} in terms of <
50 # _n >  _m <-- m < n;
50 # _n <= _m <-- m >= n;
50 # _n >= _m <-- Not(n<m);


Function("!=",{aLeft,aRight}) Not(aLeft=aRight);

/* Shifting operators */

n_IsInteger << m_IsInteger <-- ShiftLeft(n,m);
n_IsInteger >> m_IsInteger <-- ShiftRight(n,m);

0 # Sqrt(0) <-- 0;
0 # Sqrt(Infinity) <--  Infinity;
0 # Sqrt(-Infinity) <-- Complex(0,Infinity);
0 # Sqrt(Undefined) <--  Undefined;
1 # Sqrt(x_IsPositiveInteger)_(IsInteger(SqrtN(x))) <-- SqrtN(x);
2 # Sqrt(x_IsPositiveNumber)_InNumericMode() <-- SqrtN(x);
2 # Sqrt(x_IsNegativeNumber) <-- Complex(0,Sqrt(-x));
/* 3 # Sqrt(x_IsNumber/y_IsNumber) <-- Sqrt(x)/Sqrt(y); */
3 # Sqrt(x_IsComplex)_InNumericMode() <-- x^(1/2);
/* Threading  */
Sqrt(xlist_IsList) <-- MapSingle("Sqrt",xlist);

90 # (Sqrt(x_IsConstant))_(IsNegativeNumber(N(x))) <-- Complex(0,Sqrt(-x));

400 # x_IsRationalOrNumber * Sqrt(y_IsRationalOrNumber)  <-- Sign(x)*Sqrt(x^2*y);
400 # Sqrt(y_IsRationalOrNumber) * x_IsRationalOrNumber  <-- Sign(x)*Sqrt(x^2*y);
400 # x_IsRationalOrNumber / Sqrt(y_IsRationalOrNumber)  <-- Sign(x)*Sqrt(x^2/y);
400 # Sqrt(y_IsRationalOrNumber) / x_IsRationalOrNumber  <-- Sign(x)*Sqrt(y/(x^2));
400 # Sqrt(y_IsRationalOrNumber) / Sqrt(x_IsRationalOrNumber)  <-- Sqrt(y/x);
400 # Sqrt(y_IsRationalOrNumber) * Sqrt(x_IsRationalOrNumber)  <-- Sqrt(y*x);
400 # Sqrt(x_IsInteger)_IsInteger(SqrtN(x)) <-- SqrtN(x);
400 # Sqrt(x_IsInteger/y_IsInteger)_(IsInteger(SqrtN(x)) And IsInteger(SqrtN(y))) <-- SqrtN(x)/SqrtN(y);

/* Integer divisions */
0 # Div(n_IsInteger,m_IsInteger) <-- DivN(n,m);
1 # Div(0  ,_m) <-- 0;
2 # Div(n_IsRationalOrNumber,m_IsRationalOrNumber) <--
[
  Local(n1,n2,m1,m2,sgn1,sgn2);
  n1:=Numer(n);
  n2:=Denom(n);
  m1:=Numer(m);
  m2:=Denom(m);
  sgn1 := Sign(n1*m2);
  sgn2 := Sign(m1*n2);
  sgn1*sgn2*Floor(DivideN(sgn1*n1*m2,sgn2*m1*n2));
];
30 # Div(n_CanBeUni,m_CanBeUni)_(Length(VarList(n*m))=1) <--
[

  Local(vars,nl,ml);
  vars:=VarList(n*m);
  nl := MakeUni(n,vars);
  ml := MakeUni(m,vars);
  NormalForm(Div(nl,ml));
];

0 # Mod(_n,m_IsRationalOrNumber)_(m<0) <-- `Hold(Mod(@n,@m));

1 # Mod(n_IsNegativeInteger,m_IsPositiveInteger) <--
[
  Local(result);
  result := ModN(n,m);
  If (result < 0,result := result + m);
  result;
];
1 # Mod(n_IsPositiveInteger,m_IsPositiveInteger) <-- ModN(n,m);
2 # Mod(0,_m) <-- 0;
2 # Mod(n_IsPositiveInteger,Infinity) <-- n;
3 # Mod(n_IsInteger,m_IsInteger) <-- ModN(n,m);
4 # Mod(n_IsNumber,m_IsNumber) <-- NonN(Mod(Rationalize(n),Rationalize(m)));

5 # Mod(n_IsRationalOrNumber,m_IsRationalOrNumber)/*_(n>0 And m>0)*/ <--
[
  Local(n1,n2,m1,m2);
  n1:=Numer(n);
  n2:=Denom(n);
  m1:=Numer(m);
  m2:=Denom(m);
  Mod(n1*m2,m1*n2)/(n2*m2);
];

6 # Mod(n_IsList,m_IsList) <-- Map("Mod",{n,m});
7 # Mod(n_IsList,_m) <-- Map("Mod",{n,FillList(m,Length(n))});


30 # Mod(n_CanBeUni,m_CanBeUni) <--
[
  Local(vars);
  vars:=VarList(n+m);
  NormalForm(Mod(MakeUni(n,vars),MakeUni(m,vars)));
];

0 # Rem(n_IsNumber,m_IsNumber) <-- n-m*Div(n,m);
30 # Rem(n_CanBeUni,m_CanBeUni) <-- Mod(n,m);


0 # Gcd(0,0) <-- 1;
1 # Gcd(0,_m) <-- Abs(m);
1 # Gcd(_n,0) <-- Abs(n);
1 # Gcd(_m,_m) <-- Abs(m);
2 # Gcd(_n,1) <-- 1;
2 # Gcd(1,_m) <-- 1;
2 # Gcd(n_IsInteger,m_IsInteger) <-- GcdN(n,m);
3 # Gcd(_n,_m)_(IsGaussianInteger(m) And IsGaussianInteger(n) )<-- GaussianGcd(n,m);

4 # Gcd(-(_n), (_m)) <-- Gcd(n,m);
4 # Gcd( (_n),-(_m)) <-- Gcd(n,m);
4 # Gcd(Sqrt(n_IsInteger),Sqrt(m_IsInteger)) <-- Sqrt(Gcd(n,m));
4 # Gcd(Sqrt(n_IsInteger),m_IsInteger) <-- Sqrt(Gcd(n,m^2));
4 # Gcd(n_IsInteger,Sqrt(m_IsInteger)) <-- Sqrt(Gcd(n^2,m));

5 # Gcd(n_IsRational,m_IsRational) <--
[
  Gcd(Numer(n),Numer(m))/Lcm(Denom(n),Denom(m));
];


10 # Gcd(list_IsList)_(Length(list)>2) <--
    [
      Local(first);
      first:=Gcd(list[1],list[2]);
      Gcd(first:Tail(Tail(list)));
    ];
14 # Gcd({0}) <-- 1;
15 # Gcd({_head}) <-- head;

20 # Gcd(list_IsList)_(Length(list)=2) <-- Gcd(list[1],list[2]);

30 # Gcd(n_CanBeUni,m_CanBeUni)_(Length(VarList(n*m))=1) <--
[
  Local(vars);
  vars:=VarList(n+m);
  NormalForm(Gcd(MakeUni(n,vars),MakeUni(m,vars)));
];

100 # Gcd(n_IsConstant,m_IsConstant) <-- 1;
110 # Gcd(_m,_n) <--
[
  Echo("Not simplified");
];


/* Least common multiple */
5  # Lcm(a_IsInteger,b_IsInteger) <-- Div(a*b,Gcd(a,b));

10 # Lcm(list_IsList)_(Length(list)>2) <--
[
	Local(first);
	first:=Lcm(list[1],list[2]);
	Lcm(first:Tail(Tail(list)));
];

10 # Lcm(list_IsList)_(Length(list)=2) <-- Lcm(list[1],list[2]);


/* Expand expands polynomials.
 */
10 # Expand(expr_CanBeUni) <-- NormalForm(MakeUni(expr));
20 # Expand(_expr) <-- expr;

10 # Expand(expr_CanBeUni(var),_var) <-- NormalForm(MakeUni(expr,var));
20 # Expand(_expr,_var) <-- expr;



RuleBase("Object",{pred,x});
Rule("Object",2,0,Apply(pred,{x})=True) x;

10 # Abs(n_IsNumber) <-- AbsN(n);
10 # Abs(n_IsPositiveNumber/m_IsPositiveNumber) <-- n/m;
10 # Abs(n_IsNegativeNumber/m_IsPositiveNumber) <-- (-n)/m;
10 # Abs(n_IsPositiveNumber/m_IsNegativeNumber) <-- n/(-m);
10 # Abs( Sqrt(_x)) <-- Sqrt(x);
10 # Abs(-Sqrt(_x)) <-- Sqrt(x);
10 # Abs(Complex(_r,_i)) <-- Sqrt(r^2 + i^2);
10 # Abs(n_IsInfinity) <-- Infinity;
10 # Abs(Undefined) <-- Undefined;
20 # Abs(n_IsList) <-- MapSingle("Abs",n);

100 # Abs(_a^_n) <-- Abs(a)^n;
100 # Abs(_a)^n_IsEven <-- a^n;
100 # Abs(_a)^n_IsOdd <-- Sign(a)*a^n;

10 # Sign(n_IsPositiveNumber) <-- 1;
10 # Sign(n_IsZero) <-- 0;
20 # Sign(n_IsNumber) <-- -1;
15 # Sign(n_IsInfinity)_(n < 0) <-- -1;
15 # Sign(n_IsInfinity)_(n > 0) <-- 1;
15 # Sign(n_IsNumber/m_IsNumber) <-- Sign(n)*Sign(m);
20 # Sign(n_IsList) <-- MapSingle("Sign",n);

100 # Sign(_a)^n_IsEven <-- 1;
100 # Sign(_a)^n_IsOdd <-- Sign(a);

5 # Floor(Infinity) <-- Infinity;
5 # Floor(-Infinity) <-- -Infinity;
5 # Floor(Undefined) <-- Undefined;
5 # Ceil(Infinity) <-- Infinity;
5 # Ceil(-Infinity) <-- -Infinity;
5 # Ceil(Undefined) <-- Undefined;
5 # Round(Infinity) <-- Infinity;
5 # Round(-Infinity) <-- -Infinity;
5 # Round(Undefined) <-- Undefined;

/* Changed by Nobbi before redefinition of Rational
10 # Floor(x_IsNumber) <-- FloorN(x);
10 # Ceil (x_IsNumber) <-- CeilN (x);
10 # Round(x_IsNumber) <-- FloorN(x+0.5);

20 # Floor(x_IsRational) _ (IsNumber(Numer(x)) And IsNumber(Denom(x))) <-- FloorN(N(x));
20 # Ceil (x_IsRational) _ (IsNumber(Numer(x)) And IsNumber(Denom(x))) <-- CeilN (N(x));
20 # Round(x_IsRational) _ (IsNumber(Numer(x)) And IsNumber(Denom(x))) <-- FloorN(N(x+0.5));
*/

10 # Floor(x_IsRationalOrNumber) 
   <-- 
   [
     x:=N(Eval(x));
//Echo("x = ",x);
     Local(prec,result,n);
     Set(prec,BuiltinPrecisionGet());
     If(IsZero(x),
       Set(n,2),
       If(x>0,
         Set(n,2+FloorN(N(FastLog(x)/FastLog(10)))),
         Set(n,2+FloorN(N(FastLog(-x)/FastLog(10))))
       ));
     If(n>prec,BuiltinPrecisionSet(n));
//Echo("Before");
     Set(result,FloorN(x));
//Echo("After");
     BuiltinPrecisionSet(prec);
     result;
   ];

//     FloorN(N(x));

10 # Ceil (x_IsRationalOrNumber) 
   <-- 
   [
     x:=N(x);
     Local(prec,result,n);
     Set(prec,BuiltinPrecisionGet());
     If(IsZero(x),Set(n,2),
     If(x>0,
       Set(n,2+FloorN(N(FastLog(x)/FastLog(10)))),
       Set(n,2+FloorN(N(FastLog(-x)/FastLog(10))))
       ));
     If(n>prec,BuiltinPrecisionSet(n));
     Set(result,CeilN(x));
     BuiltinPrecisionSet(prec);
     result;
   ];
//   CeilN (N(x));
10 # Round(x_IsRationalOrNumber) <-- FloorN(N(x+0.5));
10 # Round(x_IsList) <-- MapSingle("Round",x);

20 # Round(x_IsComplex)  _ (IsRationalOrNumber(Re(x)) And IsRationalOrNumber(Im(x)) )
		<-- FloorN(N(Re(x)+0.5)) + FloorN(N(Im(x)+0.5))*I;


// Canonicalise an expression so its terms are grouped to the right
// ie a+(b+(c+d))
// This doesn't preserve order of terms, when doing this would cause more
// subtractions and nested parentheses than necessary.
1 # CanonicalAdd((_a+_b)+_c) <-- CanonicalAdd(CanonicalAdd(a)+
                                              CanonicalAdd(CanonicalAdd(b)+
                                                           CanonicalAdd(c)));
1 # CanonicalAdd((_a-_b)+_c) <-- CanonicalAdd(CanonicalAdd(a)+
                                              CanonicalAdd(CanonicalAdd(c)-
                                                           CanonicalAdd(b)));
1 # CanonicalAdd((_a+_b)-_c) <-- CanonicalAdd(CanonicalAdd(a)+
                                              CanonicalAdd(CanonicalAdd(b)-
                                                           CanonicalAdd(c)));
1 # CanonicalAdd((_a-_b)-_c) <-- CanonicalAdd(CanonicalAdd(a)-
                                              CanonicalAdd(CanonicalAdd(b)+
                                                           CanonicalAdd(c)));
2 # CanonicalAdd(_a)         <-- a;

////////////////////// Log rules stuff //////////////////////

// LnExpand
1 # LnExpand(Ln(x_IsInteger))
                            <-- Add(Map({{n,m},m*Ln(n)},Transpose(Factors(x))));
1 # LnExpand(Ln(_a*_b))     <-- LnExpand(Ln(a))+LnExpand(Ln(b));
1 # LnExpand(Ln(_a/_b))     <-- LnExpand(Ln(a))-LnExpand(Ln(b));
1 # LnExpand(Ln(_a^_n))     <-- LnExpand(Ln(a))*n;
2 # LnExpand(_a)            <-- a;

// LnCombine is nice and simple now
LnCombine(_a) <-- DoLnCombine(CanonicalAdd(a));

// Combine single terms.  This can always be done without a recursive call.
1 # DoLnCombine(Ln(_a))              <-- Ln(a);
1 # DoLnCombine(Ln(_a)*_b)           <-- Ln(a^b);
1 # DoLnCombine(_b*Ln(_a))           <-- Ln(a^b);

// Deal with the first two terms so they are both simple logs if at all
// possible.  This involves converting a*Ln(b) to Ln(b^a) and moving log terms
// to the start of expressions.  One of either of these operations always takes
// us to a strictly simpler form than we started in, so we can get away with
// calling DoLnCombine again with the partly simplified argument.

// TODO: Make this deal with division everywhere it deals with multiplication

// first term is a log multiplied by something
2 # DoLnCombine(Ln(_a)*_b+_c)        <-- DoLnCombine(Ln(a^b)+c);
2 # DoLnCombine(Ln(_a)*_b-_c)        <-- DoLnCombine(Ln(a^b)-c);
2 # DoLnCombine(_b*Ln(_a)+_c)        <-- DoLnCombine(Ln(a^b)+c);
2 # DoLnCombine(_b*Ln(_a)-_c)        <-- DoLnCombine(Ln(a^b)-c);

// second term of a two-term expression is a log multiplied by something
2 # DoLnCombine(_a+(_c*Ln(_b)))      <-- DoLnCombine(a+Ln(b^c));
2 # DoLnCombine(_a-(_c*Ln(_b)))      <-- DoLnCombine(a-Ln(b^c));
2 # DoLnCombine(_a+(Ln(_b)*_c))      <-- DoLnCombine(a+Ln(b^c));
2 # DoLnCombine(_a-(Ln(_b)*_c))      <-- DoLnCombine(a-Ln(b^c));

// second term of a three-term expression is a log multiplied by something
2 # DoLnCombine(_a+((Ln(_b)*_c)+_d)) <-- DoLnCombine(a+(Ln(b^c)+d));
2 # DoLnCombine(_a+((Ln(_b)*_c)-_d)) <-- DoLnCombine(a+(Ln(b^c)-d));
2 # DoLnCombine(_a-((Ln(_b)*_c)+_d)) <-- DoLnCombine(a-(Ln(b^c)+d));
2 # DoLnCombine(_a-((Ln(_b)*_c)-_d)) <-- DoLnCombine(a-(Ln(b^c)-d));

2 # DoLnCombine(_a+((_c*Ln(_b))+_d)) <-- DoLnCombine(a+(Ln(b^c)+d));
2 # DoLnCombine(_a+((_c*Ln(_b))-_d)) <-- DoLnCombine(a+(Ln(b^c)-d));
2 # DoLnCombine(_a-((_c*Ln(_b))+_d)) <-- DoLnCombine(a-(Ln(b^c)+d));
2 # DoLnCombine(_a-((_c*Ln(_b))-_d)) <-- DoLnCombine(a-(Ln(b^c)-d));

// Combine the first two terms if they are logs, otherwise move one or both to
// the front, then recurse on the remaining possibly-log-containing portion.
// (the code makes more sense than this comment)
3 # DoLnCombine(Ln(_a)+Ln(_b))       <-- Ln(a*b);
3 # DoLnCombine(Ln(_a)-Ln(_b))       <-- Ln(a/b);
3 # DoLnCombine(Ln(_a)+(Ln(_b)+_c))  <-- DoLnCombine(Ln(a*b)+c);
3 # DoLnCombine(Ln(_a)+(Ln(_b)-_c))  <-- DoLnCombine(Ln(a*b)-c);
3 # DoLnCombine(Ln(_a)-(Ln(_b)+_c))  <-- DoLnCombine(Ln(a/b)-c);
3 # DoLnCombine(Ln(_a)-(Ln(_b)-_c))  <-- DoLnCombine(Ln(a/b)+c);

// We know that at least one of the first two terms isn't a log
4 # DoLnCombine(Ln(_a)+(_b+_c))      <-- b+DoLnCombine(Ln(a)+c);
4 # DoLnCombine(Ln(_a)+(_b-_c))      <-- b+DoLnCombine(Ln(a)-c);
4 # DoLnCombine(Ln(_a)-(_b+_c))      <-- DoLnCombine(Ln(a)-c)-b;
4 # DoLnCombine(Ln(_a)-(_b-_c))      <-- DoLnCombine(Ln(a)+c)-b;

4 # DoLnCombine(_a+(Ln(_b)+_c))      <-- a+DoLnCombine(Ln(b)+c);
4 # DoLnCombine(_a+(Ln(_b)-_c))      <-- a+DoLnCombine(Ln(b)-c);
4 # DoLnCombine(_a-(Ln(_b)+_c))      <-- a-DoLnCombine(Ln(b)+c);
4 # DoLnCombine(_a-(Ln(_b)-_c))      <-- a-DoLnCombine(Ln(b)-c);

// If we get here we know that neither of the first two terms is a log
5 # DoLnCombine(_a+(_b+_c))          <-- a+(b+DoLnCombine(c));

// Finished
6 # DoLnCombine(_a)                  <-- a;