# Limits of functions
# this function is a general limit function and allows (for instance) one to
# take limits of functions, derivatives, integrals, etc.
#FIXME: I need a way to truncate the return value to disregard epsilon < tolerance (so 1/x -> 0 ;-).

# documentation:
# this detect the existence of a limit by a function being almost constant over
# a long enough interval.
# Because of the existence of continuous (even of smooth) function with compact
# support (i.e., that vanish outside of some finite interval [a,b]),
# it is impossible to determine the limit of an arbitrary continuous/smooth
# function as it goes to infinity
# in a finite number of steps, since it is always possible that a function will
# become constant (or close to it), and then suddenly have a step or bulge at
# 10^43 (for instance).
# Thus, this function is necessarily imperfect.
# If you have more information about a function (for instance, it is a
# rational function), then there are other ways of computing its limits.
# 
# 

# Limits:
SetHelp ("NumericalLimitAtInfinity", "calculus", "Attempt to calculate the limit of f(step_fun(i)) as i goes from 1 to N")
function NumericalLimitAtInfinity(_f,step_fun,tolerance,successive_for_success,N) =
# this evaluates _f(step_fun(i)) as i goes from 1 to N.
# If s_f_s of the pairwise successive differences are less than tolerance, then
# it returns the last value; if this doesn't happen, it returns Null.
# INPUT: N had better be at least 1
#FIXME: should have a way of dealing with +infinity, -infinity, and
# bounded oscillation (like sin(x))
(
    current_limit = _f(step_fun(1));
    number_of_consecutive_differences_within_tolerance = 0;
    for i = 2 to N do (
        new_limit = _f(step_fun(i));
        if (|new_limit-current_limit| < tolerance) then (
            number_of_consecutive_differences_within_tolerance = number_of_consecutive_differences_within_tolerance+1
        ) else (
            number_of_consecutive_differences_within_tolerance = 0
        );
        current_limit = new_limit;
        if (number_of_consecutive_differences_within_tolerance >= successive_for_success) then
            return current_limit;
    );
    null
)
protect ("NumericalLimitAtInfinity")


# The following are simple functions to find
# lim_{x -> x_0} f(x)
# ...from the left
# ...from the right
# ...two-sided

SetHelp ("ContinuousTolerance", "parameters", "Tolerance for continuity of functions and for calculating the limit")
parameter ContinuousTolerance=10.0^(-5);
SetHelp ("ContinuousSFS", "parameters", "How many successive steps to be within tolerance for calculation of continuity")
parameter ContinuousSFS=20;
SetHelp ("ContinuousNumberOfTries", "parameters", "How many iterations to try to find the limit for continuity and limits")
parameter ContinuousNumberOfTries=100;

SetHelp ("LeftLimit", "calculus", "Calculate the left limit of a real-valued function at x0")
function LeftLimit(f,x0) =
(
	NumericalLimitAtInfinity (f,
				  `(n)=(x0-2.0^(-n)),
				  ContinuousTolerance,
				  ContinuousSFS,
				  ContinuousNumberOfTries)
)
protect ("LeftLimit")

SetHelp ("RightLimit", "calculus", "Calculate the right limit of a real-valued function at x0")
function RightLimit(f,x0) =
(
	NumericalLimitAtInfinity (f,
				  `(n)=(x0+2.0^(-n)),
				  ContinuousTolerance,
				  ContinuousSFS,
				  ContinuousNumberOfTries)
)
protect ("RightLimit")

SetHelp ("Limit", "calculus", "Calculate the limit of a real-valued function at x0.  Tries to calculate both left and right limits.")
function Limit(f,x0) =
(
	LeftLim = LeftLimit(f,x0);
	RightLim = RightLimit(f,x0);
	if ( not IsNull(LeftLim) and
	     not IsNull(RightLim) and
	     |LeftLim-RightLim| < 2*ContinuousTolerance ) then
		(LeftLim+RightLim)/2
)
protect ("Limit")

SetHelp ("IsContinuous", "calculus", "Try and see if a real-valued function is continuous at x0 by calculating the limit there")
function IsContinuous(f,x0) =
(
    l = Limit(f,x0);
    not IsNull(l) and |l-f(x0)| < ContinuousTolerance
)
protect ("IsContinuous")