create table div;   # initialise package table

create table div.aux

div.initstring :=
   'div package v0.3.0 as of May 12, 2013\n';

# Metatable for adding and subtracting fractions
div.div_mt := [
   '__add' ~ proc(a, b) is
                return div.add(a, b)
             end,
   '__sub' ~ proc(a, b) is
                return div.subtract(a, b)
             end,
   '__mul' ~ proc(a, b) is
                return div.multiply(a, b)
             end,
   '__div' ~ proc(a, b) is
                return div.divide(a, b)
             end,
   '__pow' ~ proc(a, b) is
                return div.pow(a, b)
             end,
   '__unm' ~ proc(a) is
                return div.multiply(-1, a)
             end,
   '__le' ~  proc(a :: div, b :: div) is
                return div.todec(a) <= div.todec(b)
             end,
   '__lt' ~  proc(a :: div, b :: div) is
                return div.todec(a) < div.todec(b)
             end,
   '__eq' ~  proc(a :: div, b :: div) is
                return div.todec(a) = div.todec(b)
             end,
   '__abs' ~ proc(a :: div) is
                return div.div(abs(a[1]), a[2])
             end,
   '__sqrt' ~ proc(a :: div) is
                return div.sqroot(a)
             end,
   '__ln' ~  proc(a :: div) is
                return div.loge(a)
             end,
   '__exp' ~ proc(a :: div) is
                return div.expe(a)
             end,
   '__sin' ~ proc(a :: div) is
                return div.sine(a)
             end,
   '__cos' ~ proc(a :: div) is
                return div.cosine(a)
             end,
   '__tan' ~ proc(a :: div) is
                return div.tangent(a)
             end,
   '__arctan' ~ proc(a :: div) is
                return div.atangent(a)
             end,
   '__sign' ~ proc(a :: div) is
                return sign(a[1])
             end,
   '__tostring' ~ proc(a) is
                if a[1] = 0 then return 0
                elif a[2] = 1 then return a[1]
                else return 'div(' & a[1] & ', ' & a[2] & ')'
                fi
             end
];


div.aux.partition := proc(x, y, procname) is  # split two fractions into its numerators and denominators
   if x :: number then
      x := div.div(x, 1)
   elif x :- 'div' then
      argerror(x, procname, 'expected a number or `div` value (fraction) as first argument')
   fi;
   if y :: number then
      y := div.div(y, 1)
   elif y :- 'div' then
      argerror(y, procname, 'expected a number or `div` value (fraction) as second argument')
   fi;
   return bottom(x), top(x), bottom(y), top(y)
end;


# define a fraction
div.div := proc(?) is
   local r, item, d;
   create seq r(3);
   if nargs < 1 or nargs > 3 then
      error('Error in `div.div`: wrong number of arguments.')
   fi;
   item := varargs[1];
   if item :: pair then  # handle div format [hh:]mm:ss (as pair or nested pair)
      if right(item) :: pair then
         insert left(item), left(right(item)), right(right(item)) into r
      else
         insert left(item), right(item) into r
      fi;
      if not checktype(r, sequence, number) then
         error('Error in `div.div`: number expected in pair.')
      fi
   else
      local v;
      for i to nargs do
         v := varargs[i];
         if v :- number then
            argerror(v, 'div.div', 'number expected')
         else
            insert v into r
         fi
      od
   fi;
   if size r = 3 then  # handle improper fractions
      local i := bottom(r);
      pop bottom from r;
      return i + div.div(r[1], r[2])
   fi;
   if size(r) = 1 then
      r[2] := 1
   elif r[2] = 0 then
      #error('Error in `div.div`: denominator is zero.')
      return undefined
   fi;
   if r[1] = 0 then r[2] := 1 fi;
   if r[2] < 0 then r[1] := -r[1]; r[2] := -r[2] fi;
   d := math.gcd(r[1], r[2]);
   r[1] := r[1]/d;
   r[2] := r[2]/d;
   settype(r, 'div');
   setmetatable(r, div.div_mt);
   return r
end;


# adds two values of type div
div.add := proc(x, y) is
   local a, b, c, d := div.aux.partition(x, y, 'div.add');
   return div.div(a*d + c*b, b*d)
end;


# subtracts two values of type div
div.subtract := proc(x, y) is
   local a, b, c, d := div.aux.partition(x, y, 'div.subtract');
   return div.div(a*d - c*b, b*d)
end;


# multiply two values of type div
div.multiply := proc(x, y) is
   local a, b, c, d := div.aux.partition(x, y, 'div.multiply');
   return div.div(a*c, b*d)
end;


# divide two values of type div
div.divide := proc(x, y) is
   local a, b, c, d := div.aux.partition(x, y, 'div.divide');
   return div.div(a*d, b*c)
end;


div.pow := proc(x, y) is
   local a, b, c, d;
   a, b, c, d := div.aux.partition(x, y, 'div.pow');
   if d = 1 then
      return div.div(a^c, b^c)
   else
      local f, p, q, eps, fl;
      f := c/d;
      fl := a^f / b^f;
      p, q, eps := math.fraction(fl);
      if eps = 0 then
         return div.div(p, q)
      else
         return fl
      fi
   fi;
end;


div.sqroot := proc(x) is
   local f := sqrt(x[1]) / sqrt(x[2]);
   if nan(f) then  # changed to `nan` operator, 2.1.2
      return undefined
   else
      local p, q, eps := math.fraction(f);
      if eps = 0 then
         return div.div(p, q)
      else
         return f
      fi
   fi
end;


div.loge := proc(x::div) is
   local f := ln(x[1]) - ln(x[2]);
   if nan(f) then  # changed to `nan` operator, 2.1.2
      return undefined
   else
      local p, q, eps := math.fraction(f);
      if eps = 0 then
         return div.div(p, q)
      else
         return f
      fi
   fi
end;


div.expe := proc(x::div) is
   local f, p, q, eps;
   f := exp(x[1]/x[2]);
   p, q, eps := math.fraction(f);
   if eps = 0 then
      return div.div(p, q)
   else
      return f
   fi
end;


div.sine := proc(x::div) is
   local f, p, q, eps;
   f := sin(x[1]/x[2]);
   p, q, eps := math.fraction(f);
   if eps = 0 then
      return div.div(p, q)
   else
      return f
   fi
end;


div.cosine := proc(x::div) is
   local f, p, q, eps;
   f := cos(x[1]/x[2]);
   p, q, eps := math.fraction(f);
   if eps = 0 then
      return div.div(p, q)
   else
      return f
   fi
end;


div.tangent := proc(x::div) is
   local f, p, q, eps;
   f := tan(x[1]/x[2]);
   p, q, eps := math.fraction(f);
   if eps = 0 then
      return div.div(p, q)
   else
      return f
   fi
end;


div.atangent := proc(x::div) is
   local f;
   f := arctan(x[1]/x[2]);
   if nan(f) then  # changed to `nan` operator, 2.1.2
      return undefined
   else
      local p, q, eps := math.fraction(f);
      if eps = 0 then
         return div.div(p, q)
      else
         return f
      fi
   fi
end;


div.todec := proc(a::div) is
   return a[1]/a[2]
end;


div.todiv := proc(a::number) is
   local n, d, err := math.fraction(a);
   return div.div(n, d), err
end;


div.numer := proc(a::div) is  # 1.11.3
   return a[1]
end;


div.denom := proc(a::div) is  # 1.11.3
   return a[2]
end;

