# The Calculus package, initiated June 30, 2008

package.checkclib('calc');  # has C dynamic link library been found ?

calc.initstring := 'calc v1.1.3 as of July 16, 2014\n';

# Return all roots of a function in one variable on the interval [a, b]; the function must be
# differentiable on [a, b].
# 0.9.1  - January 21, 2008 / July 02, 2008; uses a modified Regula Falsi algorithm
# 0.22.3 - June 12, 2009: corrected and extended check whether calculated roots are really zeros of f.

# Arguments:
#    f - the function
#    a - left border
#    b - right border
#    stepsize - step size for subintervals at which a change of sign has been detected
#    accuracy - accuracy for regula falsi method used and final determination of possible zeros
calc.zero := proc(f :: procedure, a :: number, b :: number, stepsize, accuracy) is
   local step, eps, poszero, fdiff, fxpdiff, zeros;
   create seq zeros;
   step := if unassigned(stepsize) then 0.1 else stepsize fi;
   eps := if unassigned(accuracy) then Eps else accuracy fi;
   if step :- number or eps :- number then
      error('Error in `calc.zero`: fourth and fifth argument must be numbers.')
   fi;
   # the following calls to regulafalsi, sections, and assignment of f1 take 97 % of total
   # computation time; the loops/inserts are as fast as calls to select
   # get roots in function f and check wehether they are really zeros (in case regula falsi gets out of bounds)
   for i in calc.sections(f, a, b, step) do
      poszero := calc.regulafalsi(f, left(i), right(i), eps);
      if assigned(poszero) and abs(f(poszero)) < eps then
         insert poszero into zeros
      fi
   od;
   # check derivative of function f
   fdiff := << x -> calc.diff(f, x, eps) >>;
   fxpdiff := << x -> calc.xpdiff(f, x, eps) >>;
   # in case the zero is also an extremum, which regula falsi will miss: get possible roots in derivative and
   # check whether they are zeros of f
   for j in calc.sections(fdiff, a, b, step) do
      poszero := calc.regulafalsi(fdiff, left(j), right(j), eps);
      if assigned(poszero) and abs(f(poszero)) < eps then
         insert poszero into zeros
         skip
      fi
      poszero := calc.regulafalsi(fxpdiff, left(j), right(j), eps);
      if assigned(poszero) and abs(f(poszero)) < eps then
         insert poszero into zeros
      fi
   od;
   sort(zeros);
   # substitute values very close to zero with zero
   return map(math.zeroin, unique(zeros), eps);
end;


# calc.minimum: calculate the local minimum locations of a function that is differentiable over
# [a, b]. Returns a sequence of _possible_ minimum locations.
# Try minimum(<< x -> x^5 >>, -10, 10) and minimum(<< x -> 1 >>, -10, 10)).
# 0.26.4, 20.08.2009; patched 28.12.2009; improved 01.01.2010
calc.minimum := proc(f :: procedure, a :: number, b :: number, stepsize, accuracy) is
   local step, eps, p, fdiff, fxpdiff, s, mins, mins2;
   create seq mins, seq mins2;  # mins2 will include all results that have been determined with less precision.
   step := if unassigned(stepsize) then 0.1 else stepsize fi;
   eps := if unassigned(accuracy) then Eps else accuracy fi;
   if step :- number or eps :- number then
      error('Error in `calc.minimum`: fourth and fifth argument must be numbers.')
   fi;
   fdiff := << x -> calc.diff(f, x, eps) >>;           # Conte and de Boor's algorithm
   fxpdiff := << x -> calc.xpdiff(f, x, eps, eps) >>;  # try Richardson's extrapolation algorithm
   # generate intervals
   s := nseq(a, b, step);
   s := << x -> seq(x, x+step) >> @ s;
   pop top from s;  # remove excess interval
   for i in s do
      p := calc.fminbr(f, i[1], i[2], eps);  # possible minimum location (x value)
      if approx(fdiff(p), 0, eps) or approx(fxpdiff(p), 0, eps) then
         # f'(p) = 0 ? Note that there are functions that have a minimum at p, but where f''(p) = 0
         # (e.g. x^4), so do not check for `not(approx(f2(p), 0))`.
         insert p into mins;
      elif eps <= Eps and (approx(fdiff(p), 0, 1e-5) or approx(fxpdiff(p), 0, 1e-5)) then
         # if p did not meet the test, check with a lower precision and if successful, also insert it into
         # the second sequence with all minimums that did not meet the first test
         insert p into mins;
         insert p into mins2
      fi
   od;
   # substitute values very close to zero with zero
   return unique(map(math.zeroin, mins, eps)), unique(map(math.zeroin, mins2, eps))
end;


# calc.maximum: calculate the local maximum locations of a function that is differentiable over
# [a, b]. Returns a sequence of _possible_ maximum locations
# 0.29.0, 23.11.2009
calc.maximum := proc(f :: procedure, a :: number, b :: number, stepsize, accuracy) is
   local g := << x -> -f(x) >>;
   return calc.minimum(g, a, b, stepsize, accuracy)
end;


# wrapper around calc.intde, calc.intdei, and calc.intdeo, 0.32.2, 23.05.2010
calc.integral := proc(f :: procedure, a :: number, b :: number, eps, omega) is
   if eps = null then
      eps := 1e-15
   elif eps :- number then
      error('Error in `calc.integral`: fourth argument is not a number.')
   fi;
   if omega = null then
      omega := 1
   elif omega :- number or omega <= 0 then
      error('Error in `calc.integral`: fifth argument is not a positive number.')
   fi;
   if b = infinity then
      local r, err;
      # try non-oscillatory algorithm
      r, err := calc.intdei(f, a, eps);
      if r = fail then  # try oscillatory algorithm
         r, err := calc.intdeo(f, a, omega, eps);
      fi;
      return r, err
   else
      return calc.intde(f, a, b, eps)
   fi
end;


##############################################################################################
# Backward compatibility to Agena 1.8.6 and earlier                                          #
##############################################################################################

calc.newton := calc.interp;

