Commit a32af779 authored by Jiri (George) Lebl's avatar Jiri (George) Lebl

Thu Jul 23 12:26:47 2009 Jiri (George) Lebl <jirka@5z.com>

	* src/funclib.c: Add UndefineAll, ProtectAll, make Undefine an alias
	  for undefine.

	* src/dict.[ch], src/eval.c: allow setting of protected parameters,
	  they just can't be deleted or changed

	* lib/*.gel: Use the ProtectAll function

	* src/testscope.gel: add extra tests
parent ef35f07e
Thu Jul 23 12:26:47 2009 Jiri (George) Lebl <jirka@5z.com>
* src/funclib.c: Add UndefineAll, ProtectAll, make Undefine an alias
for undefine.
* src/dict.[ch], src/eval.c: allow setting of protected parameters,
they just can't be deleted or changed
* lib/*.gel: Use the ProtectAll function
* src/testscope.gel: add extra tests
Wed Jul 22 14:44:10 2009 Jiri (George) Lebl <jirka@5z.com>
* configure.in: raise version, require slightly newer glib (2.10)
......
......@@ -44,7 +44,6 @@ function OneSidedThreePointFormula(f,x0,h) =
(-3*f(x0)+4*(f call (x0+h))-(f call (x0+2*h)))/(2*h)
)
SetHelp("OneSidedThreePointFormula","calculus","Compute one-sided derivative using three-point formula");
protect("OneSidedThreePointFormula");
# Two-sided three-point formula, Section 4.1, Formula 4.5, p. 161
function TwoSidedThreePointFormula(f,x0,h) =
......@@ -63,7 +62,6 @@ function TwoSidedThreePointFormula(f,x0,h) =
((f call (x0+h))-(f call (x0-h)))/(2*h)
)
SetHelp("TwoSidedThreePointFormula","calculus","Compute two-sided derivative using three-point formula");
protect("TwoSidedThreePointFormula");
# One-sided five-point formula, Section 4.1, Formula 4.7, p. 161
function OneSidedFivePointFormula(f,x0,h) =
......@@ -82,7 +80,6 @@ function OneSidedFivePointFormula(f,x0,h) =
(-25*(f call (x0))+48*(f call (x0+h))-36*(f call (x0+2*h))+16*(f call (x0+3*h))-3*(f call (x0+4*h)))/(12*h)
)
SetHelp("OneSidedFivePointFormula","calculus","Compute one-sided derivative using five point formula");
protect("OneSidedFivePointFormula");
# Two-sided five-point formula, Section 4.1, Formula 4.6, p. 161
function TwoSidedFivePointFormula(f,x0,h) =
......@@ -101,7 +98,6 @@ function TwoSidedFivePointFormula(f,x0,h) =
((f call (x0-2*h))-8*(f call (x0-h))+8*(f call (x0+h))-(f call (x0+2*h)))/(12*h)
)
SetHelp("TwoSidedFivePointFormula","calculus","Compute two-sided derivative using five-point formula");
protect("TwoSidedFivePointFormula");
SetHelp ("DerivativeTolerance", "parameters", "Tolerance for calculating the derivatives of functions")
parameter DerivativeTolerance=10.0^(-5);
......@@ -119,7 +115,6 @@ function IsDifferentiable(f,x0) =
| NumericalLeftDerivative (f,x0) - NumericalRightDerivative (f,x0) | < 2*DerivativeTolerance
)
SetHelp("IsDifferentiable","calculus","Test for differentiability by approximating the left and right limits and comparing");
protect("IsDifferentiable");
# Simple derivative function
#FIXME!!!! BROKEN!!!! ???
......@@ -134,10 +129,8 @@ function NumericalDerivative(f,x0) =
DerivativeNumberOfTries)
)
SetHelp("NumericalDerivative","calculus","Attempt to calculate numerical derivative");
protect("NumericalDerivative");
SetHelpAlias ("NumericalDerivative", "NDerivative");
NDerivative = NumericalDerivative
protect("NDerivative");
function NumericalLeftDerivative(f,x0) =
(
......@@ -148,7 +141,6 @@ function NumericalLeftDerivative(f,x0) =
DerivativeNumberOfTries)
)
SetHelp("NumericalLeftDerivative","calculus","Attempt to calculate numerical left derivative");
protect("NumericalLeftDerivative");
function NumericalRightDerivative(f,x0) =
(
......@@ -159,7 +151,6 @@ function NumericalRightDerivative(f,x0) =
DerivativeNumberOfTries)
)
SetHelp("NumericalRightDerivative","calculus","Attempt to calculate numerical right derivative");
protect("NumericalRightDerivative");
function Derivative(f,x0) =
(
......@@ -170,4 +161,3 @@ function Derivative(f,x0) =
df(x0)
)
SetHelp("Derivative","calculus","Attempt to calculate derivative by trying first symbolically and then numerically");
protect("Derivative");
......@@ -15,7 +15,6 @@ function NumericalFourierSeriesFunction(f,L,N) =
FourierSeriesFunction(c@(1),c@(2),L)
)
SetHelp("NumericalFourierSeriesFunction","calculus","FIXME");
protect("NumericalFourierSeriesFunction");
function FourierSeriesFunction(a,b,L) =
(
......@@ -42,7 +41,6 @@ function FourierSeriesFunction(a,b,L) =
)
)
SetHelp("FourierSeries","calculus","FIXME");
protect("FourierSeries");
function NumericalFourierSeriesCoefficients(f,L,N) =
(
......@@ -67,7 +65,6 @@ function NumericalFourierSeriesCoefficients(f,L,N) =
`[a,b]
)
SetHelp("NumericalFourierSeriesCoefficients","calculus","FIXME");
protect("NumericalFourierSeriesCoefficients");
function NumericalFourierSineSeriesCoefficients(f,L,N) =
(
......@@ -89,7 +86,6 @@ function NumericalFourierSineSeriesCoefficients(f,L,N) =
b
)
SetHelp("NumericalFourierSineSeriesCoefficients","calculus","FIXME");
protect("NumericalFourierSineSeriesCoefficients");
function NumericalFourierCosineSeriesCoefficients(f,L,N) =
(
......@@ -113,7 +109,6 @@ function NumericalFourierCosineSeriesCoefficients(f,L,N) =
a
)
SetHelp("NumericalFourierCosineSeriesCoefficients","calculus","FIXME");
protect("NumericalFourierCosineSeriesCoefficients");
function PeriodicExtension(f,a,b) =
(
......@@ -131,7 +126,6 @@ function PeriodicExtension(f,a,b) =
)
)
SetHelp("PeriodicExtension","calculus","FIXME");
protect("PeriodicExtension");
function EvenPeriodicExtension(f,L) =
(
......@@ -150,7 +144,6 @@ function EvenPeriodicExtension(f,L) =
)
)
SetHelp("EvenPeriodicExtension","calculus","FIXME");
protect("EvenPeriodicExtension");
function OddPeriodicExtension(f,L) =
(
......@@ -169,4 +162,3 @@ function OddPeriodicExtension(f,L) =
)
)
SetHelp("OddPeriodicExtension","calculus","FIXME");
protect("OddPeriodicExtension");
......@@ -36,7 +36,6 @@ function CompositeSimpsonsRuleTolerance(f,a,b,FourthDerivativeBound,Tolerance) =
CompositeSimpsonsRule (f, a, b, n)
)
protect ("CompositeSimpsonsRuleTolerance")
# FIXME: reference, though this is really stupid
SetHelp ("MidpointRule", "calculus", "Integration by midpoint rule")
......@@ -56,7 +55,6 @@ function MidpointRule(f,a,b,n) =
s = sum i=1 to n do float(f(a+(len*(i-0.5))/n));
(s*len)/n
)
protect ("MidpointRule")
SetHelp ("NumericalIntegralSteps", "parameters", "Steps to perform in NumericalIntegral")
parameter NumericalIntegralSteps = 1000
......@@ -66,4 +64,3 @@ parameter NumericalIntegralFunction = `CompositeSimpsonsRule
SetHelp ("NumericalIntegral", "calculus", "Integration by rule set in NumericalIntegralFunction of f from a to b using NumericalIntegralSteps steps")
function NumericalIntegral(f,a,b) = (NumericalIntegralFunction call (f,a,b,NumericalIntegralSteps))
protect ("NumericalIntegral")
......@@ -44,7 +44,6 @@ function NumericalLimitAtInfinity(_f,step_fun,tolerance,successive_for_success,N
);
null
)
protect ("NumericalLimitAtInfinity")
# The following are simple functions to find
......@@ -69,7 +68,6 @@ function LeftLimit(f,x0) =
ContinuousSFS,
ContinuousNumberOfTries)
)
protect ("LeftLimit")
SetHelp ("RightLimit", "calculus", "Calculate the right limit of a real-valued function at x0")
function RightLimit(f,x0) =
......@@ -80,7 +78,6 @@ function RightLimit(f,x0) =
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) =
......@@ -92,7 +89,6 @@ function Limit(f,x0) =
|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) =
......@@ -100,4 +96,3 @@ function IsContinuous(f,x0) =
l = Limit(f,x0);
not IsNull(l) and |l-f(x0)| < ContinuousTolerance
)
protect ("IsContinuous")
......@@ -32,7 +32,6 @@ function InfiniteSum (func, start, inc) = (
);
null
);
protect("InfiniteSum")
#calculate an infinite sum until the new terms stop making a difference
SetHelp("InfiniteSum2","calculus","Try to calculate an infinite sum for a double parameter function with func(arg,n)");
......@@ -58,7 +57,6 @@ function InfiniteSum2(func,arg,start,inc) = (
);
null
);
protect("InfiniteSum2")
#calculate an infinite product until the new terms stop making a difference
SetHelp("InfiniteProduct","calculus","Try to calculate an infinite product for a single parameter function");
......@@ -84,7 +82,6 @@ function InfiniteProduct (func, start, inc) = (
);
null
);
protect("InfiniteProduct")
#calculate an infinite product until the new terms stop making a difference
SetHelp("InfiniteProduct2","calculus","Try to calculate an infinite product for a double parameter function with func(arg,n)");
......@@ -110,5 +107,4 @@ function InfiniteProduct2(func,arg,start,inc) = (
);
null
);
protect("InfiniteProduct2")
......@@ -10,7 +10,6 @@ function Subfactorial(n) = (
(error("Subfactorial: argument not an integer >= 0");bailout);
(n!) * sum k=1 to n do ((-1)^k)/(k!)
)
protect("Subfactorial")
SetHelp("Catalan","combinatorics","Get n'th catalan number");
function Catalan(n) = (
......@@ -20,7 +19,6 @@ function Catalan(n) = (
(error("Catalan: argument not an integer >= 0");bailout);
nCr(2*n,n)/(n+1)
);
protect("Catalan")
## Factorial
## Defined by n! = n(n-1)(n-2)...
......@@ -30,7 +28,6 @@ function Factorial(n) = (
(error("Factorial: argument not an integer >= 0");bailout);
n!
)
protect("Factorial")
## Double Factorial
## Defined by n!! = n(n-2)(n-4)...
......@@ -40,7 +37,6 @@ function DoubleFactorial(n) = (
(error("DoubleFactorial: argument not an integer >= 0");bailout);
n!!
)
protect("DoubleFactorial")
SetHelp ("RisingFactorial", "combinatorics", "(Pochhammer) Rising factorial: (n)_k = n(n+1)...(n+(k-1))")
## (Pochhammer) Rising factorial: (n)_k = n(n+1)...(n+(k-1))
......@@ -53,10 +49,8 @@ function RisingFactorial(n,k) =
prod i=0 to (k-1) do (n+i)
);
)
protect("RisingFactorial")
SetHelpAlias ("RisingFactorial", "Pochhammer")
Pochhammer = RisingFactorial
protect("Pochhammer")
# Falling factorial: (n)_k = n(n-1)...(n-(k-1))
SetHelp ("FallingFactorial", "combinatorics", "Falling factorial: (n)_k = n(n-1)...(n-(k-1))")
......@@ -69,7 +63,6 @@ function FallingFactorial(n,k) =
prod i=0 to (k-1) do (n-i)
)
)
protect("FallingFactorial");
## Binomial Coefficients
......@@ -87,7 +80,6 @@ function nPr(n,r) = (
prod i=0 to (r-1) do (n-i)
)
);
protect("nPr")
# nCr is built in now!
......@@ -106,7 +98,6 @@ function Multinomial(v,arg...) = (
MatrixSum(m)!/MatrixProduct(ApplyOverMatrix(m,`(x)=x!))
)
protect("Multinomial");
SetHelp("Pascal","combinatorics","Get the pascal's triangle as a matrix");
function Pascal(i) = (
......@@ -120,7 +111,6 @@ function Pascal(i) = (
);
m
);
protect("Pascal")
SetHelp("Triangular", "combinatorics", "Calculate the nth triangular number");
function Triangular(nth) = (
......@@ -130,4 +120,3 @@ function Triangular(nth) = (
(error("Trianglular: argument not an integer larger than or equal to 0");bailout);
(nth*(nth+1))/2
);
protect("Triangular")
......@@ -23,7 +23,6 @@ function Hofstadter(n) = (
_h (n)
)
protect("Hofstadter");
# return the Fibonacci number, calculated using an iterative method
SetHelp("Fibonacci", "combinatorics", "Calculate n'th Fibonacci number");
......@@ -36,10 +35,8 @@ function Fibonacci(x) = (
(error("Fibonacci: argument less than zero");bailout)
else (([1,1;1,0]^x)@(1,2))
);
protect("Fibonacci")
SetHelpAlias ("Fibonacci", "fib");
fib=Fibonacci
protect("fib")
## Harmonic Numbers
## H_n^(r) (the nth harmonic number of order r)
......@@ -50,10 +47,8 @@ function HarmonicNumber(n,r) = (
return ApplyOverMatrix2 (m, n, HarmonicNumber);
sum x=1 to n do x^(-r)
)
protect("HarmonicNumber");
SetHelpAlias ("HarmonicNumber", "HarmonicH");
HarmonicH = HarmonicNumber;
protect("HarmonicH");
# return the Frobenius number
......@@ -94,7 +89,6 @@ function FrobeniusNumber(v,arg...) = (
);
k
);
protect("FrobeniusNumber");
SetHelp("GreedyAlgorithm", "combinatorics", "Use greedy algorithm to find c, for c . v = n. (v must be sorted)");
function GreedyAlgorithm(n,v) = (
......@@ -125,7 +119,6 @@ function GreedyAlgorithm(n,v) = (
);
null
);
protect("GreedyAlgorithm");
SetHelp("StirlingNumberFirst", "combinatorics", "Stirling number of the first kind");
function StirlingNumberFirst(n,m) = (
......@@ -143,10 +136,8 @@ function StirlingNumberFirst(n,m) = (
StirlingNumberFirst(n-1,m-1) - n * StirlingNumberFirst(n-1,m)
)
)
protect("StirlingNumberFirst");
SetHelpAlias ("StirlingNumberFirst", "StirlingS1");
StirlingS1 = StirlingNumberFirst;
protect("StirlingS1");
SetHelp("StirlingNumberSecond", "combinatorics", "Stirling number of the second kind");
function StirlingNumberSecond(n,m) = (
......@@ -157,10 +148,8 @@ function StirlingNumberSecond(n,m) = (
1/(m!) * (sum j=0 to m do ((-1)^j * Binomial(m,j)*(m-j)^n))
)
protect("StirlingNumberSecond");
SetHelpAlias ("StirlingNumberSecond", "StirlingS2");
StirlingS2 = StirlingNumberSecond;
protect("StirlingS2");
## More: BernoulliPolynomial,
## EulerNumber, EulerPolynomial
......
......@@ -9,7 +9,6 @@ SetHelp ("GaloisMatrix", "combinatorics", "Galois matrix given a linear combinin
function GaloisMatrix(combining_rule) = (
[[0;I(columns(combining_rule)-1)],combining_rule.']
)
protect ("GaloisMatrix")
# Linear recursive sequence
SetHelp ("LinearRecursiveSequence", "combinatorics", "Compute linear recursive sequence using galois stepping")
......@@ -28,4 +27,3 @@ function LinearRecursiveSequence(seed_values,combining_rule,n) =
else null # (if sequence is reversible)
)
)
protect ("LinearRecursiveSequence")
......@@ -21,7 +21,6 @@ function EulersMethod(f,x0,y0,x1,n) = (
);
y
)
protect ("EulersMethod")
# See Handbook of Mathematics and Computational Science,
# John W.Harris, Horst Stocker
......@@ -55,4 +54,3 @@ function RungeKutta(f,x0,y0,x1,n) = (
);
y
)
protect ("RungeKutta")
......@@ -49,7 +49,6 @@ function FindRootBisection(f,a,b,TOL,N) =
);
[0,p,iteration]
)
protect ("FindRootBisection")
#function root_find_newton
#FIXME: I can't take derivatives, so this doesn't work
......@@ -89,7 +88,6 @@ function FindRootSecant(f,a,b,TOL,N) =
);
[0,p,iteration]
)
protect ("FindRootSecant")
# The Method of False Position, Section 2.3, Algorithm 2.5, p. 64
SetHelp ("FindRootFalsePosition", "equation_solving",
......@@ -125,7 +123,6 @@ function FindRootFalsePosition(f,a,b,TOL,N) =
);
[0,p,iteration]
)
protect ("FindRootFalsePosition")
# FIXME: This doesn't work right
# The M"uller's Method, Section 2.6, Algorithm 2.8, p. 88
......@@ -171,4 +168,3 @@ function FindRootMullersMethod(f,x1,x2,x3,TOL,N) =
);
[0,p,iteration]
)
protect ("FindRootMullersMethod")
......@@ -44,7 +44,6 @@ function CubicFormula(pol) = (
p/(3*u*omega1) - u*omega1 - a/3
p/(3*u*omega2) - u*omega2 - a/3]
)
protect ("CubicFormula")
## see planetmath
SetHelp ("QuarticFormula", "equation_solving",
......@@ -86,7 +85,6 @@ function QuarticFormula(pol) = (
(r3pr4 + A2)/2
(r3pr4 - A2)/2]
)
protect ("QuarticFormula")
SetHelp ("PolynomialRoots", "equation_solving",
"Find roots of a polynomial (given as vector of coefficients)")
......@@ -106,4 +104,3 @@ function PolynomialRoots(p) = (
else # if elements(p) == 5 then
QuarticFormula (p)
)
protect ("PolynomialRoots")
......@@ -9,9 +9,6 @@ arg=Argument;
SetHelp("Argument","functions","argument (angle) of complex number");
SetHelpAlias ("Argument", "arg");
SetHelpAlias ("Argument", "Arg");
protect("Argument");
protect("Arg");
protect("arg");
#
......@@ -23,35 +20,29 @@ SetHelp("MoebiusDiskMapping","functions","Moebius mapping of the disk to itself
function MoebiusDiskMapping(a,z) = (
(z-a)/(1-conj(a)*z)
)
protect("MoebiusDiskMapping");
SetHelp("MoebiusMapping","functions","Moebius mapping using the cross ratio taking z2,z3,z4 to 1,0, and infinity respectively");
function MoebiusMapping(z,z2,z3,z4) = (
((z-z3)/(z-z4))/((z2-z3)/(z2-z4))
)
protect("MoebiusMapping");
SetHelp("MoebiusMappingInftyToOne","functions","Moebius mapping using the cross ratio taking infinity to 1 and z3,z4 to 0 and infinity respectively");
function MoebiusMappingInftyToOne(z,z3,z4) = (
(z-z3)/(z-z4)
)
protect("MoebiusMappingInftyToOne");
SetHelp("MoebiusMappingInftyToZero","functions","Moebius mapping using the cross ratio taking infinity to 0 and z2,z4 to 1 and infinity respectively");
function MoebiusMappingInftyToZero(z,z2,z4) = (
(z2-z4)/(z-z4)
)
protect("MoebiusMappingInftyToZero");
SetHelp("MoebiusMappingInftyToInfty","functions","Moebius mapping using the cross ratio taking infinity to infinity and z2,z3 to 1 and 0 respectively");
function MoebiusMappingInftyToInfty(z,z2,z3) = (
(z-z3)/(z2-z3)
)
protect("MoebiusMappingInftyToInfty");
SetHelp("cis","functions","The cis function, that is cos(x)+i*sin(x)");
function cis(x) = cos(x)+(sin(x))i
protect("cis");
# FIXME: this is too stupid, must have better implementation
# to be at all useful
......@@ -79,7 +70,6 @@ protect("cis");
#(error("ZetaFunction: Not yet defined for Re z < 0");bailout);
#)
#)
#protect("ZetaFunction");
#SetHelp("GammaFunctionTolerance", "parameters", "Tolerance for the GammaFunction");
#parameter GammaFunctionTolerance = 10.0^(-10);
......@@ -88,4 +78,3 @@ protect("cis");
#function GammaFunction(z) = (
# FIXME: implement
#)
#protect("GammaFunction");
......@@ -10,7 +10,6 @@
## FIXME: this works fine for vectors, but I need to cast matrices
# as vectors to get it to work for matrices.
function DiscreteDelta(v) = KroneckerDelta([0,v])
protect("DiscreteDelta");
SetHelp("DiscreteDelta","functions","Returns 1 iff all elements are zero");
## Kronecker Delta
......@@ -24,12 +23,10 @@ function KroneckerDelta(v) =
);
1
)
protect("KroneckerDelta");
SetHelp("KroneckerDelta","functions","Returns 1 iff all elements are equal");
## Unit Step Function
## The integral of the Dirac Delta function
## FIXME: should have option to make UnitStep(0) be undefined
function UnitStep(x) = if (x<0) then 0 else 1
protect("UnitStep");
SetHelp("UnitStep","functions","The unit step function = 0 for x<0, 1 otherwise. This is the integral of the Dirac Delta function.");
......@@ -13,7 +13,6 @@ function rad2deg(x) = (
(error("rad2deg: argument not a value");bailout);
(x*180)/pi
);
protect("rad2deg")
SetHelp("deg2rad", "functions", "Convert degrees to radians");
function deg2rad(x) = (
......@@ -23,7 +22,6 @@ function deg2rad(x) = (
(error("deg2rad: argument not a value");bailout);
(x*pi)/180
);
protect("deg2rad")
#FIXME: these may not deal well with zero values.
......@@ -36,7 +34,6 @@ function asin(x) = (
if (x==1) then pi/2 else if (x==-1) then -pi/2
else atan(x/sqrt(1-x^2))
);
protect("asin")
arcsin = asin
SetHelpAlias ("asin", "arcsin");
......@@ -48,7 +45,6 @@ function asinh(x) = (
(error("asinh: argument not a value");bailout);
ln(x+sqrt((x^2)+1))
);
protect("asinh")
arcsinh = asinh
SetHelpAlias ("asinh", "arcsinh");
......@@ -61,7 +57,6 @@ function acos(x) = (
if (x==0) then pi/2
else atan(sqrt(1-x^2)/x)+(if x>0 then 0 else pi)
);
protect("acos")
arccos = acos
SetHelpAlias ("acos", "arccos");
......@@ -73,7 +68,6 @@ function acosh(x) = (
(error("acosh: argument not a value");bailout);
ln(x+sqrt((x^2)-1))
);
protect("acosh")
arccosh = acosh
SetHelpAlias ("acosh", "arccosh");
......@@ -85,7 +79,6 @@ function cot(x) = (
(error("cot: argument not a value");bailout);
1/tan(x)
);
protect("cot")
SetHelp("coth","trigonometry","The hyperbolic cotangent function");
function coth(x) = (
......@@ -95,7 +88,6 @@ function coth(x) = (
(error("coth: argument not a value");bailout);
1/tanh(x)
);
protect("coth")
SetHelp("acot","trigonometry","The arccot (inverse cot) function");
function acot(x) = (
......@@ -105,7 +97,6 @@ function acot(x) = (
(error("acot: argument not a value");bailout);
atan(1/x)
);
protect("acot")
arccot = acot
SetHelpAlias ("acot", "arccot");
......@@ -117,7 +108,6 @@ function acoth(x) = (
(error("acoth: argument not a value");bailout);
atanh(1/x)
);
protect("acoth")
arccoth = acoth
SetHelpAlias ("acoth", "arccoth");
......@@ -129,7 +119,6 @@ function tanh(x) = (
(error("tanh: argument not a value");bailout);
sinh(x)/cosh(x)
);
protect("tanh")
SetHelp("atanh","trigonometry","The arctanh (inverse tanh) function");
function atanh(x) = (
......@@ -139,7 +128,6 @@ function atanh(x) = (
(error("atanh: argument not a value");bailout);
ln((1+x)/(1-x))/2
);
protect("atanh")
arctanh = atanh
SetHelpAlias ("atanh", "arctanh");
......@@ -151,7 +139,6 @@ function csc(x) = (
(error("csc: argument not a value");bailout);
1/sin(x)
);
protect("csc")
SetHelp("csch","trigonometry","The hyperbolic cosecant function");
function csch(x) = (
......@@ -161,7 +148,6 @@ function csch(x) = (
(error("csch: argument not a value");bailout);
1/sinh(x)
);
protect("csch")
SetHelp("acsc","trigonometry","The inverse cosecant function");
function acsc(x) = (
......@@ -171,7 +157,6 @@ function acsc(x) = (
(error("acsch: argument not a value");bailout);
asin(1/x)
);
protect("acsc")
arccsc = acsc
SetHelpAlias ("acsc", "arccsc");
......@@ -183,7 +168,6 @@ function acsch(x) = (
(error("acsc: argument not a value");bailout);
asinh(1/x)
);
protect("acsch")
arccsch = acsch
SetHelpAlias ("acsch", "arccsch");
......@@ -195,7 +179,6 @@ function sec(x) = (
(error("sec: argument not a value");bailout);
1/cos(x)
);
protect("sec")
SetHelp("sech","trigonometry","The hyperbolic secant function");
function sech(x) = (
......@@ -205,7 +188,6 @@ function sech(x) = (
(error("sech: argument not a value");bailout);
1/cosh(x)
);
protect("sech")
SetHelp("asec","trigonometry","The inverse secant function");
function asec(x) = (
......@@ -215,7 +197,6 @@ function asec(x) = (
(error("asec: argument not a value");bailout);
acos(1/x)
);
protect("asec")
arcsec = asec
SetHelpAlias ("asec", "arcsec");
......@@ -227,7 +208,6 @@ function asech(x) = (
(error("asech: argument not a value");bailout);
acosh(1/x)
);
protect("asech")
arcsech = asech
SetHelpAlias ("asech", "arcsech");
......@@ -252,7 +232,6 @@ function log(x,b...) = (
(error("log: arguments not values");bailout);
ln(x)/ln(base)
);
protect("log")
parameter ErrorFunctionTolerance = 10.0^(-10);
SetHelp ("ErrorFunctionTolerance", "parameters", "Tolerance of the ErrorFunction")
......@@ -280,7 +259,6 @@ function ErrorFunction (x) = (
) while (|t| > ErrorFunctionTolerance);
s
);
protect("ErrorFunction")
erf = ErrorFunction
SetHelpAlias ("ErrorFunction", "erf");
......@@ -303,4 +281,3 @@ function NewtonsMethodPoly(poly,guess,epsilon,maxn) = (
);
null
)
protect("NewtonsMethodPoly")
......@@ -8,13 +8,11 @@ SetHelp("PoissonKernel","functions","Poisson kernel on D(0,1) (not normalized to
function PoissonKernel(r,sigma) = (
(1 - r^2) / (1 - 2*r*cos(sigma) + r^2)
)
protect("PoissonKernel")
SetHelp("PoissonKernelRadius","functions","Poisson kernel on D(0,R) (not normalized to 1)");
function PoissonKernelRadius(r,sigma) = (
(R^2 - r^2) / (R^2 - 2*R*r*cos(sigma) + r^2)
)
protect("PoissonKernelRadius")
# See Planetmath http://planetmath.org/encyclopedia/DirchletKernel.html
SetHelp("DirichletKernel","functions","Dirichlet kernel of order n");
......@@ -25,7 +23,6 @@ function DirichletKernel(n,t) = (
( sin ((n+0.5)*t) ) / ( sin(t/2))
)
)
protect("DirichletKernel")
# See Planetmath http://planetmath.org/encyclopedia/FejerKernel.html
SetHelp("FejerKernel","functions","Fejer kernel of order n");
......@@ -37,4 +34,3 @@ function FejerKernel(n,t) = (
( ( sin ((n*t)/2) ) / ( sin(t/2)) )^2
)
)
protect("FejerKernel")
......@@ -5,10 +5,8 @@
#just for somewhat of "source compatibility"
SetHelp ("AbsoluteValue", "numeric", "Absolute value");
function AbsoluteValue (x) = |x|
protect ("AbsoluteValue")
SetHelpAlias ("AbsoluteValue", "abs")
abs = AbsoluteValue
protect ("abs")
SetHelp("Sign","numeric","Return the sign (-1,0,1)");
function Sign (x) = (
......@@ -24,10 +22,8 @@ function Sign (x) = (
-1
) else 0
);
protect("Sign")
SetHelpAlias ("Sign", "sign")
sign = Sign
protect("sign")
#-----
......@@ -39,7 +35,6 @@ function FractionalPart(x) =
else
x - IntegerPart (x)
)
protect ("FractionalPart")
# Chop replaces very small number with zero
parameter ChopTolerance = 10.0^(-10);
......@@ -55,7 +50,6 @@ function Chop(x) =
else
x
)
protect ("Chop")
#-----
# Mod (built-in)
......@@ -64,7 +58,6 @@ protect ("Chop")
# Division w/o remainder
SetHelp ("IntegerQuotient", "numeric", "Division w/o remainder")
function IntegerQuotient(m,n) = floor(m/n)
protect ("IntegerQuotient")
# IntergerQuotient w/offset (such that d <= m-r*n < d+n
# IntegerDigits = (convert interger to its list of digits, base b, of length len)
# IntegerExponent = (number of trailing zeros in base b expansion = heighest power
......
......@@ -5,13 +5,11 @@
#FIXME: check parameters
SetHelp ("BilinearForm", "linear_algebra", "Evaluate (v,w) with respect to the bilinear form given by the matrix A")
function BilinearForm(v,A,w) = (v.'*A*w)@(1)
protect ("BilinearForm")
# Evaluate (v,w) with respect to the sesquilinear form given by the matrix A.
#FIXME: check parameters
SetHelp ("SesquilinearForm", "linear_algebra", "Evaluate (v,w) with respect to the sesquilinear form given by the matrix A")
function SesquilinearForm(v,A,w) = (v'*A*w)@(1)
protect ("SesquilinearForm")
# The angle of two vectors, given an inner product
#FIXME: default inner product!
......@@ -29,16 +27,13 @@ function VectorAngle(v,w,B...) =
acos(P(v,w)/(P(v,v)*P(w,w)))
)
protect ("VectorAngle")
#FIXME: check parameters
SetHelp ("BilinearFormFunction", "linear_algebra", "Return a function that evaluates two vectors with respect to the bilinear form given by A")
function BilinearFormFunction(A)=(`(v,w)=(v.'*A*w)@(1))
protect ("BilinearFormFunction")
SetHelp ("SesquilinearFormFunction", "linear_algebra", "Return a function that evaluates two vectors with respect to the sesquilinear form given by A")
function SesquilinearFormFunction(A)=(`(v,w)=(v'*A*w)@(1))
protect ("SesquilinearFormFunction")
# Projection onto a vector space
# Projection of vector v onto subspace W given a sesquilinear form B
......@@ -57,7 +52,6 @@ function Projection(v,W,B...) =
(P(v,c)/P(c,c))*c
)
)
protect ("Projection")
# Gram-Schmidt Orthogonalization
# Described, for instance, in Hoffman & Kunze, ``Linear Algebra'' p.~280
......@@ -85,4 +79,3 @@ function GramSchmidt(v,B...) =
);