[genius] Thu Jul 23 12:26:47 2009 Jiri (George) Lebl <jirka 5z com>
- From: George Lebl <jirka src gnome org>
- To: svn-commits-list gnome org
- Subject: [genius] Thu Jul 23 12:26:47 2009 Jiri (George) Lebl <jirka 5z com>
- Date: Thu, 23 Jul 2009 17:32:14 +0000 (UTC)
commit a32af7797b56284617d2061b5994c577dc1318a1
Author: Jiri (George) Lebl <jirka 5z com>
Date: Thu Jul 23 12:28:31 2009 -0500
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
ChangeLog | 12 +++++
lib/calculus/differentiation.gel | 10 ----
lib/calculus/fourier.gel | 8 ---
lib/calculus/integration.gel | 3 -
lib/calculus/limits.gel | 5 --
lib/calculus/sums_products.gel | 4 --
lib/combinatorics/factorial.gel | 11 ----
lib/combinatorics/misc.gel | 11 ----
lib/combinatorics/recursive_sequences.gel | 2 -
lib/equation_solving/diffeqs.gel | 2 -
lib/equation_solving/find_root.gel | 4 --
lib/equation_solving/formulas.gel | 3 -
lib/functions/complex_numbers.gel | 11 ----
lib/functions/delta.gel | 3 -
lib/functions/elementary.gel | 23 ---------
lib/functions/kernels.gel | 4 --
lib/functions/numerical.gel | 7 ---
lib/linear_algebra/bilinear_forms.gel | 7 ---
lib/linear_algebra/linear_algebra.gel | 42 ----------------
lib/linear_algebra/misc.gel | 27 ----------
lib/linear_algebra/special_matrices.gel | 14 -----
lib/linear_algebra/subspaces.gel | 13 -----
lib/make_loader_gel.sh | 4 +-
lib/misc/misc.gel | 3 -
lib/number_theory/factoring.gel | 5 --
lib/number_theory/misc.gel | 12 -----
lib/number_theory/modulus.gel | 10 ----
lib/number_theory/primes.gel | 6 --
lib/sets/basic.gel | 2 -
lib/statistics/basic.gel | 19 -------
lib/symbolic/differentiation.gel | 3 -
src/dict.c | 37 +++++++++++++-
src/dict.h | 1 +
src/eval.c | 11 ++--
src/funclib.c | 47 +++++++++++++++++-
src/testscope.gel | 75 +++++++++++++++++++++++++++++
36 files changed, 174 insertions(+), 287 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index a31e1d9..8bdd64c 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,15 @@
+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)
diff --git a/lib/calculus/differentiation.gel b/lib/calculus/differentiation.gel
index 751eabd..37730f3 100644
--- a/lib/calculus/differentiation.gel
+++ b/lib/calculus/differentiation.gel
@@ -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");
diff --git a/lib/calculus/fourier.gel b/lib/calculus/fourier.gel
index a8d3fb4..f65a15f 100644
--- a/lib/calculus/fourier.gel
+++ b/lib/calculus/fourier.gel
@@ -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");
diff --git a/lib/calculus/integration.gel b/lib/calculus/integration.gel
index 7e5c9ed..abcc489 100644
--- a/lib/calculus/integration.gel
+++ b/lib/calculus/integration.gel
@@ -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")
diff --git a/lib/calculus/limits.gel b/lib/calculus/limits.gel
index aeb40c2..899d024 100644
--- a/lib/calculus/limits.gel
+++ b/lib/calculus/limits.gel
@@ -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")
diff --git a/lib/calculus/sums_products.gel b/lib/calculus/sums_products.gel
index 2f54c35..0fc1844 100644
--- a/lib/calculus/sums_products.gel
+++ b/lib/calculus/sums_products.gel
@@ -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")
diff --git a/lib/combinatorics/factorial.gel b/lib/combinatorics/factorial.gel
index 97e3810..934976e 100644
--- a/lib/combinatorics/factorial.gel
+++ b/lib/combinatorics/factorial.gel
@@ -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")
diff --git a/lib/combinatorics/misc.gel b/lib/combinatorics/misc.gel
index 492671e..15439f0 100644
--- a/lib/combinatorics/misc.gel
+++ b/lib/combinatorics/misc.gel
@@ -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
diff --git a/lib/combinatorics/recursive_sequences.gel b/lib/combinatorics/recursive_sequences.gel
index 503e54f..67992b4 100644
--- a/lib/combinatorics/recursive_sequences.gel
+++ b/lib/combinatorics/recursive_sequences.gel
@@ -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")
diff --git a/lib/equation_solving/diffeqs.gel b/lib/equation_solving/diffeqs.gel
index 66690fb..2e65089 100644
--- a/lib/equation_solving/diffeqs.gel
+++ b/lib/equation_solving/diffeqs.gel
@@ -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")
diff --git a/lib/equation_solving/find_root.gel b/lib/equation_solving/find_root.gel
index 5266add..8e59bf9 100644
--- a/lib/equation_solving/find_root.gel
+++ b/lib/equation_solving/find_root.gel
@@ -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")
diff --git a/lib/equation_solving/formulas.gel b/lib/equation_solving/formulas.gel
index fefed5c..eb9c576 100644
--- a/lib/equation_solving/formulas.gel
+++ b/lib/equation_solving/formulas.gel
@@ -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")
diff --git a/lib/functions/complex_numbers.gel b/lib/functions/complex_numbers.gel
index cf1a5f5..897cd2f 100644
--- a/lib/functions/complex_numbers.gel
+++ b/lib/functions/complex_numbers.gel
@@ -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");
diff --git a/lib/functions/delta.gel b/lib/functions/delta.gel
index b2cafef..1ccef22 100644
--- a/lib/functions/delta.gel
+++ b/lib/functions/delta.gel
@@ -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.");
diff --git a/lib/functions/elementary.gel b/lib/functions/elementary.gel
index 1f4fc27..7282cc0 100644
--- a/lib/functions/elementary.gel
+++ b/lib/functions/elementary.gel
@@ -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")
diff --git a/lib/functions/kernels.gel b/lib/functions/kernels.gel
index 602f14c..02e60cb 100644
--- a/lib/functions/kernels.gel
+++ b/lib/functions/kernels.gel
@@ -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")
diff --git a/lib/functions/numerical.gel b/lib/functions/numerical.gel
index 94a420e..6174d92 100644
--- a/lib/functions/numerical.gel
+++ b/lib/functions/numerical.gel
@@ -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
diff --git a/lib/linear_algebra/bilinear_forms.gel b/lib/linear_algebra/bilinear_forms.gel
index d28258d..686cb50 100644
--- a/lib/linear_algebra/bilinear_forms.gel
+++ b/lib/linear_algebra/bilinear_forms.gel
@@ -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...) =
);
w
)
-protect ("GramSchmidt")
diff --git a/lib/linear_algebra/linear_algebra.gel b/lib/linear_algebra/linear_algebra.gel
index 7695528..b76bb5e 100644
--- a/lib/linear_algebra/linear_algebra.gel
+++ b/lib/linear_algebra/linear_algebra.gel
@@ -3,12 +3,10 @@
# Transpose of a matrix
SetHelp ("Transpose", "linear_algebra", "Transpose of a matrix")
function Transpose(M)=M.'
-protect ("Transpose")
# Conjugate Transpose of a matrix
SetHelp ("ConjugateTranspose", "linear_algebra", "Conjugate transpose of a matrix (adjoint)")
function ConjugateTranspose(M)=M'
-protect ("Conjugate Transpose")
#FIXME: We'll need more sophisticated tools for Tensors
@@ -25,7 +23,6 @@ function NonPivotColumns(M) =
else
non_pivot
)
-protect ("NonPivotColumns")
# Delete a row from a matrix
# Shamelessly copied from DeleteRowColumn
@@ -37,7 +34,6 @@ function DeleteRow(M,row) = (
else if row==row_dim_M then M@(1:row_dim_M-1,)
else [M@(1:row-1,);M@(row+1:row_dim_M,)]
)
-protect ("DeleteRow")
# Delete a column from a matrix
# Also shamelessly copied from DeleteRowColumn
@@ -49,18 +45,14 @@ function DeleteColumn(M,col) = (
else if col==col_dim_M then M@(,1:col_dim_M-1)
else [M@(,1:col-1),M@(,col+1:col_dim_M)]
)
-protect ("DeleteColumn")
# Column Reduced Echelon Form
SetHelp("cref", "linear_algebra", "Compute the Column Reduced Echelon Form")
function cref(M)=(rref(M.').')
-protect ("cref")
SetHelpAlias ("cref", "CREF")
CREF = cref
-protect ("CREF")
SetHelpAlias ("cref", "ColumnReducedEchelonForm")
ColumnReducedEchelonForm = cref
-protect ("ColumnReducedEchelonForm")
SetHelp ("StripZeroRows", "matrix", "Removes any all-zero rows of M")
function StripZeroRows (M) =
@@ -70,7 +62,6 @@ function StripZeroRows (M) =
(error("StripZeroColumns: argument not a matrix");bailout);
StripZeroColumns (M.').'
)
-protect ("StripZeroRows")
# Cross product of two vectors (in R^3)
##<refentry id="fun-linear-algebra-cross-product">
@@ -118,7 +109,6 @@ function CrossProduct(v,w) =
M=[v;w;1,1,1];
[Minor(M,3,1),Minor(M,3,2),Minor(M,3,3)]
)
-protect ("CrossProduct")
# Direct sum of a vector of matrices
SetHelp ("DirectSumMatrixVector", "linear_algebra", "Direct sum of a vector of matrices")
@@ -132,7 +122,6 @@ function DirectSumMatrixVector(v) = (
else
[v@(1),0;0,DirectSumMatrixVector(v@(1,2:columns(v)))]
)
-protect ("DirectSumMatrixVector");
# Direct sum of matrices
SetHelp ("DirectSum", "linear_algebra", "Direct sum of matrices")
@@ -144,7 +133,6 @@ function DirectSum(M,N...) = (
else
[M,0;0,DirectSumMatrixVector(N)]
)
-protect ("DirectSum");
# (Kronecker) Tensor product of two matrices
#FIXME!
@@ -166,7 +154,6 @@ function IsInvertible(n) =
# determinant is a unit
error ("IsInvertible: argument not a value or matrix");bailout
)
-protect ("IsInvertible")
SetHelp ("IsInvertibleField", "linear_algebra", "Is a matrix (or number) invertible over a field")
function IsInvertibleField(n) =
@@ -178,7 +165,6 @@ function IsInvertibleField(n) =
# determinant is a unit
error ("IsInvertibleField: argument not a value or matrix");bailout
)
-protect ("IsInvertibleField")
# SmithNormal Form for fields (will end up with 1's on the diagonal)
SetHelp("SmithNormalFormField", "linear_algebra", "Smith Normal Form for fields (will end up with 1's on the diagonal)")
@@ -186,7 +172,6 @@ function SmithNormalFormField(A) =
(
ref (ref (A).').'
)
-protect ("SmithNormalFormField")
SetHelp("SmithNormalFormInteger", "linear_algebra", "Smith Normal Form for square integer matrices (not its characteristic)")
function SmithNormalFormInteger(M) =
@@ -200,7 +185,6 @@ function SmithNormalFormInteger(M) =
S@(i,i) = d@(i)/d@(i-1);
S
)
-protect ("SmithNormalFormInteger")
# AuxilliaryUnitMatrix
SetHelp("AuxilliaryUnitMatrix", "linear_algebra", "Get the auxilliary unit matrix of size n")
@@ -211,7 +195,6 @@ function AuxilliaryUnitMatrix(n) =
else
[0,I(n-1);0,0]
)
-protect ("AuxilliaryUnitMatrix")
# AuxilliaryUnitMatrix
SetHelp("JordanBlock", "linear_algebra", "Get the jordan block corresponding to lambda and n")
@@ -226,10 +209,8 @@ function JordanBlock(n,lambda) =
else
lambda * I(n) + AuxilliaryUnitMatrix(n)
)
-protect ("JordanBlock")
SetHelpAlias("JordanBlock", "J")
J = JordanBlock
-protect ("J")
# Dot product
SetHelp("DotProduct", "matrix", "Get the dot product of two vectors (no conjugates)")
@@ -246,7 +227,6 @@ function DotProduct(u,v) =
else # if columns (u) == 1 and rows (v) == 1 then
( ( u .' ) * (v.') )@(1)
)
-protect ("DotProduct")
# Outer product
SetHelp("OuterProduct", "matrix", "Get the outer product of two vectors")
@@ -263,7 +243,6 @@ function OuterProduct(u,v) =
else # if columns (u) == 1 and rows (v) == 1 then
u * conj(v)
)
-protect ("OuterProduct")
SetHelp("InfNorm", "linear_algebra", "Get the Inf Norm of a vector");
function InfNorm(v) = (
@@ -272,7 +251,6 @@ function InfNorm(v) = (
else
max (|v|)
)
-protect("InfNorm")
SetHelp("Norm", "linear_algebra", "Get the p Norm (or 2 Norm if no p is supplied) of a vector");
function Norm(v,p...) = (
@@ -285,10 +263,8 @@ function Norm(v,p...) = (
else
(sum k in v do (|k|^p@(1)))^(1/p@(1))
)
-protect("Norm")
SetHelpAlias ("Norm", "norm")
norm = Norm
-protect("norm")
SetHelp ("CharacteristicPolynomial", "linear_algebra", "Get the characteristic polynomial as a vector")
function CharacteristicPolynomial(M) = (
@@ -303,14 +279,11 @@ function CharacteristicPolynomial(M) = (
v@(n+1) = (-1) ^ n;
v
)
-protect("CharacteristicPolynomial")
SetHelpAlias ("CharacteristicPolynomial", "CharPoly")
CharPoly = CharacteristicPolynomial
-protect("CharPoly")
SetHelp ("CharacteristicPolynomialFunction", "linear_algebra", "Get the characteristic polynomial as a function")
function CharacteristicPolynomialFunction(M) = PolyToFunction (CharacteristicPolynomial (M))
-protect("CharacteristicPolynomialFunction")
SetHelp ("DeterminantalDivisorsInteger", "linear_algebra", "Get the determinantal divisors of an integer matrix (not its characteristic)")
function DeterminantalDivisorsInteger(M) = (
@@ -328,7 +301,6 @@ function DeterminantalDivisorsInteger(M) = (
v@(n) = |det (M)|;
v
)
-protect("DeterminantalDivisorsInteger")
SetHelp("InvariantFactorsInteger", "linear_algebra", "Get the invariant factors of a square integer matrix (not its characteristic)")
function InvariantFactorsInteger(M) =
@@ -341,7 +313,6 @@ function InvariantFactorsInteger(M) =
H@(i) = d@(i)/d@(i-1);
H
)
-protect ("InvariantFactorsInteger")
SetHelp("IsNormal", "linear_algebra", "Is a matrix normal")
function IsNormal(M) = (
@@ -349,7 +320,6 @@ function IsNormal(M) = (
(error("IsNormal: argument not a square matrix");bailout);
M' * M == M * M'
)
-protect ("IsNormal")
SetHelp("IsUnitary", "linear_algebra", "Is a matrix unitary")
function IsUnitary(M) = (
@@ -357,7 +327,6 @@ function IsUnitary(M) = (
(error("IsUnitary: argument not a square matrix");bailout);
M' * M == I(columns(M))
)
-protect ("IsUnitary")
SetHelp("IsHermitian", "linear_algebra", "Is a matrix hermitian")
function IsHermitian(M) = (
@@ -365,7 +334,6 @@ function IsHermitian(M) = (
(error("IsHermitian: argument not a square matrix");bailout);
M' == M
)
-protect ("IsHermitian")
SetHelp("IsSkewHermitian", "linear_algebra", "Is a matrix skew-hermitian")
function IsSkewHermitian(M) = (
@@ -373,7 +341,6 @@ function IsSkewHermitian(M) = (
(error("IsSkewHermitian: argument not a square matrix");bailout);
M' == -M
)
-protect ("IsSkewHermitian")
SetHelp("IsPositiveDefinite", "linear_algebra", "Is a matrix positive definite")
function IsPositiveDefinite(M) = (
@@ -391,7 +358,6 @@ function IsPositiveDefinite(M) = (
) else
false
)
-protect ("IsPositiveDefinite")
SetHelp("IsPositiveSemidefinite", "linear_algebra", "Is a matrix positive semidefinite")
function IsPositiveSemidefinite(M) = (
@@ -409,7 +375,6 @@ function IsPositiveSemidefinite(M) = (
) else
false
)
-protect ("IsPositiveSemidefinite")
SetHelp("RayleighQuotient", "linear_algebra", "Return the Rayleigh quotient of a matrix and a vector")
function RayleighQuotient(A,x) = (
@@ -417,7 +382,6 @@ function RayleighQuotient(A,x) = (
if columns (x) > 1 then x = x.';
( ( x' * A * x ) / ( x' * x ) ) @(1)
)
-protect ("RayleighQuotient")
SetHelp("Eigenvalues", "linear_algebra", "Get the eigenvalues of a matrix (Currently only for up to 4x4 or triangular matrices)")
function Eigenvalues(M) = (
@@ -435,10 +399,8 @@ function Eigenvalues(M) = (
) else
(error("Eigenvalues: Currently only implemented for up to 4x4 or triangular matrices");bailout)
)
-protect ("Eigenvalues")
eig = Eigenvalues;
SetHelpAlias ("Eigenvalues", "eig")
-protect ("eig")
SetHelp("Eigenvectors", "linear_algebra", "Get the eigenvalues and eigenvectors of a matrix (Currently only for up to 2x2 matrices)")
function Eigenvectors(M, evsmults...) = (
@@ -501,7 +463,6 @@ function Eigenvectors(M, evsmults...) = (
) else
(error("Eigenvectors: Currently only implemented for up to 2x2");bailout)
)
-protect ("Eigenvectors")
# Written by David W. Hutchison, dahutchi indiana edu for Genius 0.7.1
# Copyright (C) 2004 David W. Hutchison (GPL)
@@ -549,7 +510,6 @@ function LUDecomposition(A,L,U) = (
*U = upper;
true
);
-protect("LUDecomposition");
# FIXME: Should use more stable algorithm
SetHelp("QRDecomposition", "linear_algebra", "Get the QR decomposition of A, returns R and Q can be a reference")
@@ -565,7 +525,6 @@ function QRDecomposition(A,Q) = (
UpperTriangular (QQ' * A)
)
-protect("QRDecomposition")
# See for example: http://www.cs.utk.edu/~dongarra/etemplates/node97.html
SetHelp("RayleighQuotientIteration", "linear_algebra", "Compute an eigenvalue using the Rayleigh Quotient Iteration method until we are epsilon from eigenvalue or for maxiter iterations")
@@ -605,4 +564,3 @@ function RayleighQuotientIteration(A,x,epsilon,maxiter,vecref) = (
);
null
)
-protect ("RayleighQuotientIteration")
diff --git a/lib/linear_algebra/misc.gel b/lib/linear_algebra/misc.gel
index 6770ba8..8386143 100644
--- a/lib/linear_algebra/misc.gel
+++ b/lib/linear_algebra/misc.gel
@@ -12,7 +12,6 @@ function ApplyOverMatrix(a,func) = (
);
r
);
-protect("ApplyOverMatrix")
SetHelp("ApplyOverMatrix2", "matrix", "Apply a function over all entries of 2 matrices (or 1 value and 1 matrix) and return a matrix of the results")
function ApplyOverMatrix2(a,b,func) = (
@@ -49,7 +48,6 @@ function ApplyOverMatrix2(a,b,func) = (
);
r
);
-protect("ApplyOverMatrix2")
#calculate a trace function
SetHelp("Trace", "linear_algebra","Calculate the trace of a matrix");
@@ -60,10 +58,8 @@ function Trace(m) = (
(error("Trace: matrix not a square");bailout);
sum i = 1 to rows(m) do m@(i,i)
);
-protect("Trace")
SetHelpAlias("Trace", "trace")
trace = Trace
-protect("trace")
#calculate convolution of two horizontal vectors
SetHelp("Convolution", "linear_algebra","Calculate convolution of two horizontal vectors");
@@ -77,10 +73,8 @@ function Convolution(a,b) = (
ca = columns(a);
sum i = 1 to ca do a@(1,i)*b@(1,ca-i+1)
);
-protect("Convolution")
SetHelpAlias("Convolution", "convol")
convol = Convolution
-protect("convol")
#calculate convolution of two horizontal vectors and return the result
#not added together but in a vector
@@ -99,7 +93,6 @@ function ConvolutionVector(a,b) = (
);
r
);
-protect("ConvolutionVector")
#calculate the sum of all elements in a matrix
SetHelp("MatrixSum", "matrix","Calculate the sum of all elements in a matrix");
@@ -108,7 +101,6 @@ function MatrixSum(a) = (
(error("MatrixSum: argument not a value only matrix");bailout);
sum n in a do n
);
-protect("MatrixSum")
SetHelp("MatrixSumSquares", "matrix","Calculate the sum of squares of all elements in a matrix");
function MatrixSumSquares(a) = (
@@ -116,7 +108,6 @@ function MatrixSumSquares(a) = (
(error("MatrixSumSquares: argument not a value only matrix");bailout);
sum n in a do n^2
);
-protect("MatrixSumSquares")
#calculate the product of all elements in a matrix
SetHelp("MatrixProduct","matrix", "Calculate the product of all elements in a matrix")
@@ -125,21 +116,17 @@ function MatrixProduct(a) = (
(error("matprod: argument not a value only matrix");bailout);
prod n in a do n
);
-protect("MatrixProduct")
SetHelp("Submatrix", "matrix", "Return column(s) and row(s) from a matrix")
function Submatrix(m,r,c) = [m@(r,c)]
-protect ("Submatrix")
SetHelp("ComplementSubmatrix", "matrix", "Remove column(s) and row(s) from a matrix");
function ComplementSubmatrix(m,r,c) = [m@(IndexComplement(r, rows(m)), IndexComplement (c, columns (m)))]
-protect ("ComplementSubmatrix")
# Minor of a matrix (determinant of a submatrix given by deleting
# one row and one column)
SetHelp("Minor","linear_algebra", "Get the i-j minor of a matrix")
function Minor(M,i,j) = det (ComplementSubmatrix (M, i, j))
-protect("Minor");
#classical adjoint (adjugate) of a matrix
SetHelp("adj","linear_algebra", "Get the classical adjoint (adjugate) of a matrix");
@@ -159,10 +146,8 @@ function adj(m) = (
);
a
);
-protect("adj")
SetHelpAlias ("adj", "Adjugate")
Adjugate = adj
-protect("Adjugate")
SetHelp("MinimizeFunction","functions","Find the first value where f(x)=0");
function MinimizeFunction(func,x,incr) = (
@@ -173,7 +158,6 @@ function MinimizeFunction(func,x,incr) = (
while(func(x)>0) do x=x+incr;
x
);
-protect("MinimizeFunction")
SetHelp("MakeDiagonal","matrix","Make diagonal matrix from a vector");
function MakeDiagonal(v,arg...) = (
@@ -191,10 +175,8 @@ function MakeDiagonal(v,arg...) = (
);
r
);
-protect("MakeDiagonal")
SetHelpAlias("MakeDiagonal","diag")
diag = MakeDiagonal
-protect("diag")
SetHelp("SwapRows","matrix","Swap two rows in a matrix");
function SwapRows(m,row1,row2) = (
@@ -210,7 +192,6 @@ function SwapRows(m,row1,row2) = (
);
m
);
-protect("SwapRows")
SetHelp("RowSum","matrix","Calculate sum of each row in a matrix");
function RowSum(m) = (
@@ -223,7 +204,6 @@ function RowSum(m) = (
);
r
);
-protect("RowSum")
SetHelp("RowSumSquares","matrix","Calculate sum of squares of each row in a matrix");
function RowSumSquares(m) = (
@@ -236,7 +216,6 @@ function RowSumSquares(m) = (
);
r
);
-protect("RowSumSquares")
#sort a horizontal vector
SetHelp("SortVector","matrix","Sort vector elements");
@@ -258,7 +237,6 @@ function SortVector(v) = (
) while not sorted;
v
);
-protect("SortVector")
SetHelp("ReverseVector","matrix","Reverse elements in a vector");
function ReverseVector(v) = (
@@ -270,7 +248,6 @@ function ReverseVector(v) = (
r@(i) = v@(ev-i+1);
r
);
-protect("ReverseVector")
SetHelp("UpperTriangular", "matrix", "Zero out entries below the diagonal")
function UpperTriangular(M) = (
@@ -283,7 +260,6 @@ function UpperTriangular(M) = (
);
M
)
-protect ("UpperTriangular")
SetHelp("LowerTriangular", "matrix", "Zero out entries above the diagonal")
function LowerTriangular(M) = (
@@ -291,7 +267,6 @@ function LowerTriangular(M) = (
(error("LowerTriangular: argument not a square matrix");bailout);
UpperTriangular (M.').'
)
-protect ("LowerTriangular")
SetHelp("CompoundMatrix", "matrix", "Calculate the kth compund matrix of A")
function CompoundMatrix(k,A) = (
@@ -305,7 +280,6 @@ function CompoundMatrix(k,A) = (
C@(i,j) = det (A@(gamma@(i),omega@(j)));
C
)
-protect ("CompoundMatrix")
SetHelp("MakeVector", "matrix", "Make column vector out of matrix by putting columns above each other")
function MakeVector(A) = (
@@ -318,4 +292,3 @@ function MakeVector(A) = (
);
r
)
-protect ("MakeVector")
diff --git a/lib/linear_algebra/special_matrices.gel b/lib/linear_algebra/special_matrices.gel
index 315f178..a9fb83b 100644
--- a/lib/linear_algebra/special_matrices.gel
+++ b/lib/linear_algebra/special_matrices.gel
@@ -15,7 +15,6 @@ function CompanionMatrix(p)= (
p = p.';
[ [0;I(elements(p)-2)] , -p@(1 : elements(p)-1,1) ]
)
-protect ("CompanionMatrix")
SetHelp ("HankelMatrix", "linear_algebra", "Hankel matrix")
function HankelMatrix(c,r)= (
@@ -25,7 +24,6 @@ function HankelMatrix(c,r)= (
H@(i,j) = if i+j-1 <= elements(c) then c@(i+j-1) else r@(i+j-elements(c));
H
)
-protect ("HankelMatrix")
SetHelp ("HilbertMatrix", "linear_algebra", "Hilbert matrix of order n")
function HilbertMatrix(n)= (
@@ -35,11 +33,9 @@ function HilbertMatrix(n)= (
H@(i,j) = 1/(i+j-1);
H
)
-protect ("HilbertMatrix")
SetHelp ("InverseHilbertMatrix", "linear_algebra", "Inverse Hilbert matrix of order n")
function InverseHilbertMatrix(n) = HilbertMatrix(n)^(-1)
-protect ("InverseHilbertMatrix")
# magic square (exists for n != 2 -- I know how to make them for n odd)
# pascal matrix
# toeplitz matrix
@@ -63,10 +59,8 @@ function VandermondeMatrix(v) =
);
vandermonde
)
-protect ("VandermondeMatrix")
SetHelpAlias ("VandermondeMatrix", "vander")
vander = VandermondeMatrix
-protect ("vander")
SetHelp ("RosserMatrix", "linear_algebra", "Rosser matrix, a classic symmetric eigenvalue test problem")
RosserMatrix = [611,196,-192,407,-8,-52,-49,29
@@ -77,7 +71,6 @@ RosserMatrix = [611,196,-192,407,-8,-52,-49,29
-52,-43,49,44,-599,411,208,208
-49,-8,8,59,208,208,99,-911
29,-44,52,-23,208,208,-911,99];
-protect("RosserMatrix")
#Rotation matrices
# See http://mathworld.wolfram.com/RotationMatrix.html
@@ -87,10 +80,8 @@ function Rotation2D(angle) = (
[cos(angle),sin(angle)
-sin(angle),cos(angle)]
)
-protect("Rotation2D")
SetHelpAlias ("Rotation2D", "RotationMatrix")
RotationMatrix = Rotation2D
-protect ("RotationMatrix")
SetHelp ("Rotation3DX", "linear_algebra", "Rotation around origin in R^3 about the x-axis")
function Rotation3DX(angle) = (
@@ -98,7 +89,6 @@ function Rotation3DX(angle) = (
0,cos(angle),sin(angle)
0,-sin(angle),cos(angle)]
)
-protect("Rotation3DX")
SetHelp ("Rotation3DY", "linear_algebra", "Rotation around origin in R^3 about the y-axis")
function Rotation3DY(angle) = (
@@ -106,7 +96,6 @@ function Rotation3DY(angle) = (
0,1,0
sin(angle),0,cos(angle)]
)
-protect("Rotation3DY")
SetHelp ("Rotation3DZ", "linear_algebra", "Rotation around origin in R^3 about the z-axis")
function Rotation3DZ(angle) = (
@@ -114,7 +103,6 @@ function Rotation3DZ(angle) = (
-sin(angle),cos(angle),0
0,0,1]
)
-protect("Rotation3DZ")
# Ported from octave
SetHelp ("CommutationMatrix", "linear_algebra", "Return the commutation matrix K(m,n) which is the unique m*n by m*n matrix such that K(m,n) * MakeVector(A) = MakeVector(A.') for all m by n matrices A.")
@@ -131,7 +119,6 @@ function CommutationMatrix(m,n) = (
);
k
)
-protect("CommutationMatrix")
#ported from octave
SetHelp ("ToeplitzMatrix", "linear_algebra", "Return the Toeplitz matrix constructed given the first column c and (optionally) the first row r.")
@@ -178,4 +165,3 @@ function ToeplitzMatrix (c, r...) = (
retval
)
-protect("ToeplitzMatrix")
diff --git a/lib/linear_algebra/subspaces.gel b/lib/linear_algebra/subspaces.gel
index 1ea58c7..ee3f1d4 100644
--- a/lib/linear_algebra/subspaces.gel
+++ b/lib/linear_algebra/subspaces.gel
@@ -16,7 +16,6 @@ function ColumnSpace(M) =
(error("ColumnSpace: argument not a matrix");bailout);
StripZeroColumns (cref (M))
)
-protect ("ColumnSpace")
# Row Space of a matrix
SetHelp ("RowSpace", "linear_algebra", "Get a basis matrix for the rowspace of a matrix")
@@ -26,7 +25,6 @@ function RowSpace(M) = (
(error("RowSpace: argument not a matrix");bailout);
ColumnSpace (M.')
)
-protect ("RowSpace")
# Rank of a matrix
SetHelp ("Rank", "linear_algebra", "Get the rank of a matrix")
@@ -37,10 +35,8 @@ function Rank(M) =
(error("Rank: argument not a matrix");bailout);
columns (M) - CountZeroColumns(cref(M))
)
-protect ("Rank")
SetHelpAlias ("Rank", "rank");
rank = Rank
-protect ("rank")
# Nullity of a matrix
SetHelp ("Nullity", "linear_algebra", "Get the nullity of a matrix")
@@ -51,35 +47,29 @@ function Nullity (M) =
(error("Nullity: argument not a matrix");bailout);
CountZeroColumns(cref(M))
)
-protect ("Nullity")
SetHelpAlias ("Nullity", "nullity");
nullity = Nullity
-protect ("nullity")
# Image of a linear transform (=column space of the corresponding matrix)
SetHelp ("Image", "linear_algebra", "Get the image (columnspace) of a linear transform")
Image = ColumnSpace
-protect ("Image")
# Null space/kernel of a linear transform
# NullSpace is now built-in for speed
SetHelp ("Kernel", "linear_algebra", "Get the kernel (nullspace) of a linear transform")
function Kernel(M) = NullSpace(M)
-protect ("Kernel")
## Vector Subspace operations
# Given W1, W2 subspaces of V, returns W1 + W2
# = {w | w=w1+w2, w1 in W1, w2 in W2}
SetHelp ("VectorSubspaceSum", "linear_algebra", "The sum of the vector spaces M and N, that is {w | w=m+n, m in M, n in N}")
function VectorSubspaceSum(M,N)=ColumnSpace([M,N])
-protect ("VectorSubspaceSum")
# Given vector spaces V1, V2, return V1 (+) V2 (the direct sum)
# (given by the column space of the direct sum of matrices)
SetHelp ("VectorSpaceDirectSum", "linear_algebra", "The direct sum of the vector spaces M and N")
function VectorSpaceDirectSum(M,N)=ColumnSpace(DirectSum(M,N))
-protect ("VectorSpaceDirectSum")
# Given vector spaces V1, V2, return V1 (x) V2 (the tensor product)
# (given by the column space of the tensor product of matrices)
@@ -92,7 +82,6 @@ protect ("VectorSpaceDirectSum")
# W1 cap W2 = (W1^perp+W2^perp)^perp
SetHelp ("VectorSubspaceIntersection", "linear_algebra", "Intersection of the subspaces given by M and N")
function VectorSubspaceIntersection(M,N)= OrthogonalComplement(VectorSubspaceSum(OrthogonalComplement(M),OrthogonalComplement(N)))
-protect ("VectorSubspaceIntersection")
# Given a vector subspace of an Inner Product Space (or maybe just
# a space with a (symmetric?) bilinear form (check lang),
@@ -101,7 +90,6 @@ protect ("VectorSubspaceIntersection")
# inner product on
SetHelp ("OrthogonalComplement", "linear_algebra", "Get the orthogonal complement of the columnspace")
function OrthogonalComplement(M)=NullSpace(M')
-protect ("OrthogonalComplement")
## (Vector space + a vector) operations
@@ -115,5 +103,4 @@ protect ("OrthogonalComplement")
#....)
SetHelp ("IsInSubspace", "linear_algebra", "Test if a vector is in a subspace")
function IsInSubspace(v,W) = (v==Projection(v,W,.))
-protect ("IsInSubspace")
diff --git a/lib/make_loader_gel.sh b/lib/make_loader_gel.sh
index 83a2264..d5ab11c 100755
--- a/lib/make_loader_gel.sh
+++ b/lib/make_loader_gel.sh
@@ -2,5 +2,7 @@
echo "# Automatically generated loader, don't touch"
echo
for n in "$@"; do
- echo load $n
+ echo "load $n"
done
+
+echo "ProtectAll()"
diff --git a/lib/misc/misc.gel b/lib/misc/misc.gel
index 906801a..4064d36 100644
--- a/lib/misc/misc.gel
+++ b/lib/misc/misc.gel
@@ -3,11 +3,9 @@
#make something a string
SetHelp("string","basic","Make a string");
function string(s) = s + ""
-protect("string")
SetHelp("Compose","basic","Compose two functions")
function Compose(f,g) = `(x)=f(g(x))
-protect("Compose")
SetHelp("ComposePower","basic","Compose a function with itself n times, passing x as argument, and returning x if n == 0")
function ComposePower(f,n,x) = (
@@ -15,4 +13,3 @@ function ComposePower(f,n,x) = (
x = f call (x);
x
)
-protect("ComposePower")
diff --git a/lib/number_theory/factoring.gel b/lib/number_theory/factoring.gel
index 9fbcde8..0eeffc5 100644
--- a/lib/number_theory/factoring.gel
+++ b/lib/number_theory/factoring.gel
@@ -14,7 +14,6 @@ function PrimeFactors (n) = (
fac@(1,2:columns(fac))
)
-protect("PrimeFactors")
SetHelp("MaximalPrimePowerFactors", "number_theory", "Return all maximal prime power factors of a number")
function MaximalPrimePowerFactors (n) = (
@@ -30,7 +29,6 @@ function MaximalPrimePowerFactors (n) = (
fac@(1,2:columns(fac))
)
-protect("MaximalPrimePowerFactors")
# given two factorizations, returns factorization of the product
# (useful for factoring rational numbers, for instance)
@@ -67,7 +65,6 @@ function CombineFactorizations(a,b) = (
);
combined
)
-protect ("CombineFactorizations")
# Returns all factors of n, i.e., all numbers between 1 and n that divide
# n.
@@ -103,7 +100,6 @@ function Factors(n) =
r = [front_list,back_list];
r@(2:elements(r)-1)
)
-protect("Factors")
SetHelp("FermatFactorization", "number_theory", "Attempt fermat factorization of n into (t-s)*(t+s), returns t and s as a vector if possible, null otherwise")
function FermatFactorization(n,tries) = (
@@ -118,4 +114,3 @@ function FermatFactorization(n,tries) = (
);
null
)
-protect("FermatFactorization")
diff --git a/lib/number_theory/misc.gel b/lib/number_theory/misc.gel
index 24f0b49..facee8d 100644
--- a/lib/number_theory/misc.gel
+++ b/lib/number_theory/misc.gel
@@ -21,7 +21,6 @@ function IsNthPower(m,n) =
(IsNthPower(Numerator(m),n)) and IsNthPower(Denominator(m),n)
)
SetHelp("IsNthPower","number_theory","Tests if a rational number is a perfect power");
-protect("IsNthPower");
# Returns the highest power of p that divides n, where n is rational.
# FIXME: This is only a valuation if p is prime, so this is an abuse
@@ -47,7 +46,6 @@ function PadicValuation(n,p) =
valuation
)
SetHelp("PadicValuation","number_theory","Returns the padic valuation (number of trailing zeros in base p).");
-protect("PadicValuation");
function RemoveFactor(n,m) =
(
@@ -65,11 +63,9 @@ function RemoveFactor(n,m) =
n
)
SetHelp("RemoveFactor","number_theory","Removes all instances of the factor m from the number n");
-protect("RemoveFactor");
function PowerMod(a,b,m) = a^b mod m
SetHelp("PowerMod","number_theory","Compute a^b mod m");
-protect("PowerMod");
function AreRelativelyPrime(a,b) = (
if IsMatrix (a) or IsMatrix (b) then
@@ -77,7 +73,6 @@ function AreRelativelyPrime(a,b) = (
gcd(a,b) == 1
)
SetHelp("AreRelativelyPrime","number_theory","Are a and b relatively prime?");
-protect("AreRelativelyPrime");
function LeastAbsoluteResidue(a,n) = (
if not IsInteger(a) or not IsPositiveInteger(n) then
@@ -92,7 +87,6 @@ function LeastAbsoluteResidue(a,n) = (
b
)
SetHelp("LeastAbsoluteResidue", "number_theory", "Return the residue of a mod n with the least absolute value (in the interval -n/2 to n/2)");
-protect("LeastAbsoluteResidue");
SetHelp("ChineseRemainder", "number_theory", "Find the x that solves the system given by the vector a and modulo the elements of m, using the Chinese Remainder Theorem");
function ChineseRemainder(a,m) = (
@@ -103,10 +97,8 @@ function ChineseRemainder(a,m) = (
)
) % M
)
-protect("ChineseRemainder");
SetHelpAlias ("ChineseRemainder", "CRT")
CRT = ChineseRemainder
-protect("CRT")
SetHelp("ConvertToBase", "number_theory", "Convert a number to a vector of powers for elements in base b");
function ConvertToBase(n,b) = (
@@ -123,7 +115,6 @@ function ConvertToBase(n,b) = (
);
ret
)
-protect("ConvertToBase");
SetHelp("ConvertFromBase", "number_theory", "Convert a vector of values indicating powers of b to a number");
function ConvertFromBase(v,b) = (
@@ -137,7 +128,6 @@ function ConvertFromBase(v,b) = (
n = n + v@(i)*b^(i-1);
n
)
-protect("ConvertFromBase");
# Eric W. Weisstein. "Bernoulli Number." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/BernoulliNumber.html
SetHelp("BernoulliNumber", "number_theory", "Return the nth Bernoulli number");
@@ -158,7 +148,6 @@ function BernoulliNumber(n) = (
)
)
)
-protect("BernoulliNumber");
SetHelp("MoebiusMu", "number_theory", "Return the Moebius mu function evaluated in n");
@@ -181,5 +170,4 @@ function MoebiusMu(n) = (
);
moeb
)
-protect("MoebiusMu");
diff --git a/lib/number_theory/modulus.gel b/lib/number_theory/modulus.gel
index 710cb9e..268be23 100644
--- a/lib/number_theory/modulus.gel
+++ b/lib/number_theory/modulus.gel
@@ -11,7 +11,6 @@ function EulerPhi(n) =
(error("EulerPhi: argument not an integer larger than 0");bailout);
n * (prod p in PrimeFactors(n) do (1-(1/p)))
)
-protect("EulerPhi");
# Silver-Pohlig-Hellman as described in
# Neil Koblitz, A Course in Number Theory and Cryptography, Springer, 1987
@@ -61,7 +60,6 @@ function SilverPohligHellmanWithFactorization(n,b,q,f) = (
# Now return the chinese remainder of the system
ChineseRemainder(x,m)
)
-protect ("SilverPohligHellmanWithFactorization")
# compute the discrete log using the SilverPohligHellman alg
SetHelp("DiscreteLog", "number_theory", "Find discrete log of n base b in F_q where q is a prime using the Silver-Pohlig-Hellman algoritm");
@@ -74,7 +72,6 @@ function DiscreteLog(n,b,q) = (
f = Factorize(q-1);
SilverPohligHellmanWithFactorization(n,b,q,f)
)
-protect ("DiscreteLog")
# precalculates the table of logarithms for the
# factor base mentioned in S for the generator b mod q
@@ -139,7 +136,6 @@ function IndexCalculusPrecalculation (b, q, S) = (
# column the computed logarithms
[S@(,1) , mat^-1 * v mod q-1]
);
-protect ("IndexCalculusPrecalculation")
# outputs log_b n mod q (using factor base S),
# first column of S is the factor base and second column
@@ -180,7 +176,6 @@ function IndexCalculus (n, b, q, S) = (
) while f > 1;
l mod q-1
);
-protect ("IndexCalculus")
# g is primitive mod p iff g^((g-1)/p) == 1 mod q for all prime
# factors p of q-1
@@ -197,7 +192,6 @@ function IsPrimitiveModWithPrimeFactors(g,q,f) = (
);
true
)
-protect ("IsPrimitiveModWithPrimeFactors")
# g is primitive mod p iff g^((g-1)/p) == 1 mod q for all prime
# factors p of q-1
@@ -213,7 +207,6 @@ function IsPrimitiveMod(g,q) = (
);
true
)
-protect ("IsPrimitiveMod")
# Find the first primitive element
SetHelp("FindPrimitiveElementMod", "number_theory", "Find the first primitive element in F_q (q must be a prime)");
@@ -227,7 +220,6 @@ function FindPrimitiveElementMod(q) = (
if g == q then null else g
);
-protect ("FindPrimitiveElementMod")
# Find a random primitive element
SetHelp("FindRandomPrimitiveElementMod", "number_theory", "Find a random primitive element in F_q (q must be a prime)");
@@ -237,7 +229,6 @@ function FindRandomPrimitiveElementMod(q) = (
g = randint(q)
while not IsPrimitiveModWithPrimeFactors(g,q,f)
);
-protect ("FindRandomPrimitiveElementMod")
#ref:
# Alfred J. Menezes, Paul C. van Oorschot, Scott A. Vanstone,
@@ -289,4 +280,3 @@ function SqrtModPrime(n,p) = (
[r,-r] mod p
)
)
-protect ("SqrtModPrime")
diff --git a/lib/number_theory/primes.gel b/lib/number_theory/primes.gel
index 008e4e6..97d41cd 100644
--- a/lib/number_theory/primes.gel
+++ b/lib/number_theory/primes.gel
@@ -17,7 +17,6 @@ function PseudoprimeTest(n,b) =
(
b^(|n|-1) == 1 mod |n|
)
-protect("PseudoprimeTest");
# An integer is said to be a PseudoPrime with base b iff it is:
# odd, composite, and satisfies b^(n-1) == 1 (mod n).
@@ -28,14 +27,12 @@ function IsPseudoprime(n,b) =
IsOdd(n) and PseudoprimeTest(n,b) and not(IsPrime(n))
)
SetHelp("IsPseudoprime","number_theory","If n is a pseudoprime base b but not a prime, that is if b^(n-1) == 1 mod n");
-protect("IsPseudoprime");
function IsStrongPseudoprime(n,b) =
(
IsOdd(n) and gcd(n,b)==1 and StrongPseudoprimeTest(n,b) and not (IsPrime(n))
)
SetHelp("IsStrongPseudoprime","number_theory","Test if n is a strong pseudoprime to base b but not a prime");
-protect("IsStrongPseudoprime");
#######################################
############## Misc ###################
@@ -58,7 +55,6 @@ function LucasLehmer(p) = (
);
S == 0
);
-protect("LucasLehmer")
# See http://primes.utm.edu/mersenne/index.html
# Or http://www.mersenne.org
@@ -110,7 +106,6 @@ MersennePrimeExponents = [2
37156667
42643801
43112609];
-protect("MersennePrimeExponents");
# See http://www.mersenne.org
SetHelp("IsMersennePrimeExponent","number_theory","Test if Mp is a Mersenne prime using a table");
@@ -132,4 +127,3 @@ function IsMersennePrimeExponent(p) = (
error("IsMersennePrimeExponent: Number too large (known values up to: " + 18297089 + ")");
bailout
);
-protect("IsMersennePrimeExponent")
diff --git a/lib/sets/basic.gel b/lib/sets/basic.gel
index 87cb9b4..0dac732 100644
--- a/lib/sets/basic.gel
+++ b/lib/sets/basic.gel
@@ -9,7 +9,6 @@ function Union(X,Y) =
if not IsIn(x,Y) then Y=[Y,x];
Y
)
-protect("Union");
SetHelp ("MakeSet", "sets", "Returns a set where every element of X appears only once");
function MakeSet(X) =
@@ -19,4 +18,3 @@ function MakeSet(X) =
if not IsIn(x,S) then S=[S,x];
S
)
-protect("MakeSet");
diff --git a/lib/statistics/basic.gel b/lib/statistics/basic.gel
index 7820fad..b25f592 100644
--- a/lib/statistics/basic.gel
+++ b/lib/statistics/basic.gel
@@ -19,7 +19,6 @@ function RowMedian(m) = (
);
r
);
-protect("RowMedian")
SetHelp("Median","statistics","Calculate median of an entire matrix");
function Median(m) = (
@@ -37,10 +36,8 @@ function Median(m) = (
(s@(1,columns(s)/2) +
s@(1,(columns(s)/2)+1))/2
);
-protect("Median")
SetHelpAlias ("Median", "median")
median = Median
-protect("median")
SetHelp("RowAverage","statistics","Calculate average of each row in a matrix");
function RowAverage(m) = (
@@ -54,10 +51,8 @@ function RowAverage(m) = (
);
r
);
-protect("RowAverage")
SetHelpAlias("RowAverage", "RowMean")
RowMean = RowAverage
-protect("RowMean")
SetHelp("Average","statistics","Calculate average of an entire matrix");
function Average(m) = (
@@ -65,16 +60,12 @@ function Average(m) = (
(error("Average: argument not value-only matrix");bailout);
MatrixSum (m) / elements (m)
);
-protect("Average")
SetHelpAlias ("Average", "average")
average = Average
-protect("average")
SetHelpAlias ("Average", "Mean")
Mean = Average
-protect("Mean")
SetHelpAlias ("Average", "mean")
mean = Average
-protect("mean")
SetHelp("RowStandardDeviation", "statistics", "Calculate the standard deviations of rows of a matrix and return a vertical vector")
function RowStandardDeviation(m) = (
@@ -92,10 +83,8 @@ function RowStandardDeviation(m) = (
);
r
);
-protect("RowStandardDeviation")
SetHelpAlias ("RowStandardDeviation", "rowstdev")
rowstdev = RowStandardDeviation
-protect("rowstdev")
SetHelp("RowPopulationStandardDeviation", "statistics", "Calculate the population standard deviations of rows of a matrix and return a vertical vector")
function RowPopulationStandardDeviation(m) = (
@@ -111,10 +100,8 @@ function RowPopulationStandardDeviation(m) = (
);
r
);
-protect("RowPopulationStandardDeviation")
SetHelpAlias ("RowPopulationStandardDeviation", "rowstdevp")
rowstdevp = RowPopulationStandardDeviation
-protect("rowstdevp")
SetHelp("StandardDeviation", "statistics", "Calculate the standard deviation of a whole matrix")
function StandardDeviation(m) = (
@@ -128,10 +115,8 @@ function StandardDeviation(m) = (
rr = rr + (k-r)^2;
sqrt(rr/(elements(m)-1))
);
-protect("StandardDeviation")
SetHelpAlias ("StandardDeviation", "stdev")
stdev = StandardDeviation
-protect("stdev")
SetHelp("PopulationStandardDeviation", "statistics", "Calculate the population standard deviation of a whole matrix")
function PopulationStandardDeviation(m) = (
@@ -143,17 +128,14 @@ function PopulationStandardDeviation(m) = (
rr = rr + (k-r)^2;
sqrt(rr/elements(m))
);
-protect("PopulationStandardDeviation")
SetHelpAlias ("PopulationStandardDeviation", "stdevp")
stdevp = PopulationStandardDeviation
-protect("stdevp")
SetHelp("GaussFunction", "statistics", "The normalized Gauss distribution function (the normal curve)")
function GaussFunction(x,sigma) = (
(1/(sigma*sqrt(2*pi))) *
exp(-x^2 / (2*sigma^2))
)
-protect("GaussFunction")
SetHelp ("GaussDistributionTolerance", "parameters", "Tolerance of the GaussDistribution function")
parameter GaussDistributionTolerance = 10.0^(-10)
@@ -166,4 +148,3 @@ function GaussDistribution(x,sigma) = (
3/(sqrt(2*pi)*sigma^5),
GaussDistributionTolerance)
)
-protect("GaussDistribution")
diff --git a/lib/symbolic/differentiation.gel b/lib/symbolic/differentiation.gel
index 642bf01..9c5f0df 100644
--- a/lib/symbolic/differentiation.gel
+++ b/lib/symbolic/differentiation.gel
@@ -15,7 +15,6 @@ function SymbolicNthDerivative(f,n) =
df
)
SetHelp("SymbolicNthDerivative","symbolic","Attempt to symbolically differentiate a function n times");
-protect("SymbolicNthDerivative");
function SymbolicNthDerivativeTry(f,n) =
(
@@ -34,7 +33,6 @@ function SymbolicNthDerivativeTry(f,n) =
df
)
SetHelp("SymbolicNthDerivativeTry","symbolic","Attempt to symbolically differentiate a function n times quietly and return null on failure");
-protect("SymbolicNthDerivativeTry");
function SymbolicTaylorApproximationFunction(f,x0,n) =
(
@@ -60,4 +58,3 @@ function SymbolicTaylorApproximationFunction(f,x0,n) =
`(x) = (tp call (x - x0))
)
SetHelp("SymbolicTaylorApproximationFunction","symbolic","Attempt to construct the taylor approximation function around x0 to the nth degree.");
-protect("SymbolicTaylorApproximationFunction");
diff --git a/src/dict.c b/src/dict.c
index b1d7e5a..d8ff824 100644
--- a/src/dict.c
+++ b/src/dict.c
@@ -508,7 +508,8 @@ d_intern (const char *id)
/* this may be inefficient as it also goes through global,
* but we don't assume we do this kind of thing often. Only
- * done in d_delete which is only done on Undefine. */
+ * done in d_delete and d_delete_global which are only done
+ * on Undefine and UndefineAll. */
static void
whack_from_all_contexts (GelEFunc *func)
{
@@ -523,6 +524,8 @@ d_delete(GelToken *id)
{
GSList *li, *list;
+ g_return_val_if_fail (id != NULL && id->token != NULL, FALSE);
+
id->protected_ = 0;
id->parameter = 0;
id->built_in_parameter = 0;
@@ -540,6 +543,35 @@ d_delete(GelToken *id)
return TRUE;
}
+gboolean
+d_delete_global(GelToken *id)
+{
+ GSList *list;
+ GelEFunc *f;
+
+ g_return_val_if_fail (id != NULL && id->token != NULL, FALSE);
+
+ printf ("FOO3 %s\n", id->token);
+
+ id->protected_ = 0;
+ id->parameter = 0;
+ id->built_in_parameter = 0;
+
+ list = g_slist_last (id->refs);
+ f = list->data;
+
+ if (f->context == 0) {
+ /* inefficient, but it need not be */
+ id->refs = g_slist_remove_link (id->refs, list);
+
+ f->id = NULL;
+ whack_from_all_contexts (f);
+ d_freefunc (f);
+ }
+
+ return TRUE;
+}
+
/*clear all context dictionaries and pop out all the contexts except
the global one
also init the context stack if it hasn't been done*/
@@ -1012,8 +1044,7 @@ d_protect_all(void)
GelEFunc *func = li->data;
if(!func->id || strcmp(func->id->token,"Ans")==0)
continue;
- if ( ! func->id->parameter)
- func->id->protected_ = 1;
+ func->id->protected_ = 1;
}
}
diff --git a/src/dict.h b/src/dict.h
index 758a61a..7c2a7ae 100644
--- a/src/dict.h
+++ b/src/dict.h
@@ -93,6 +93,7 @@ GelEFunc * d_lookup_global (GelToken *id);
GelToken * d_intern (const char *id);
gboolean d_delete(GelToken *id);
+gboolean d_delete_global(GelToken *id);
/*clear all context dictionaries and pop out all the contexts except
the global one
diff --git a/src/eval.c b/src/eval.c
index 9c36fee..528c3a0 100644
--- a/src/eval.c
+++ b/src/eval.c
@@ -5444,17 +5444,16 @@ iter_equalsop(GelETree *n)
}
if(l->type == GEL_IDENTIFIER_NODE) {
- if G_UNLIKELY (d_curcontext() == 0 &&
- l->id.id->protected_) {
- gel_errorout (_("Trying to set a protected id '%s'"),
- l->id.id->token);
- return;
- }
if (l->id.id->parameter) {
GelETree *ret = set_parameter (l->id.id, r);
if (ret != NULL)
replacenode (n, ret);
return;
+ } else if G_UNLIKELY (d_curcontext() == 0 &&
+ l->id.id->protected_) {
+ gel_errorout (_("Trying to set a protected id '%s'"),
+ l->id.id->token);
+ return;
} else if(r->type == GEL_FUNCTION_NODE) {
d_addfunc (d_makerealfunc (r->func.func,
l->id.id,
diff --git a/src/funclib.c b/src/funclib.c
index 4b5a5fd..d1f1867 100644
--- a/src/funclib.c
+++ b/src/funclib.c
@@ -274,6 +274,44 @@ undefine_op (GelCtx *ctx, GelETree * * a, gboolean *exception)
}
static GelETree *
+UndefineAll_op (GelCtx *ctx, GelETree * * a, gboolean *exception)
+{
+ GSList *li;
+ GSList *list;
+
+ list = g_slist_copy (d_getcontext_global ());
+
+ for (li = list;
+ li != NULL;
+ li = li->next) {
+ GelEFunc *f = li->data;
+ GelToken *tok = f->id;
+ if (tok->parameter) {
+ printf ("FOO %s %d %d %p\n",
+ tok->token,
+ tok->parameter,
+ tok->protected_,
+ f);
+ }
+ if ( ! tok->protected_ &&
+ strcmp (tok->token, "Ans") != 0) {
+ printf ("FOO2 %s\n", tok->token);
+ d_delete_global (tok);
+ }
+ }
+
+ return gel_makenum_null ();
+}
+
+static GelETree *
+ProtectAll_op (GelCtx *ctx, GelETree * * a, gboolean *exception)
+{
+ d_protect_all ();
+
+ return gel_makenum_null ();
+}
+
+static GelETree *
true_op (GelCtx *ctx, GelETree * * a, gboolean *exception)
{
return gel_makenum_bool (1);
@@ -6339,12 +6377,15 @@ gel_funclib_addall(void)
FUNC (StringToAlphabet, 2, "str,alphabet", "misc", N_("Convert a string to a vector of 0-based alphabet values (positions in the alphabet string), -1's for unknown letters"));
FUNC (AlphabetToString, 2, "vec,alphabet", "misc", N_("Convert a vector of 0-based alphabet values (positions in the alphabet string) to a string"));
- FUNC (protect, 1, "id", "basic", N_("Protect a variable from being modified"));
- FUNC (unprotect, 1, "id", "basic", N_("Unprotect a variable from being modified"));
+ FUNC (protect, 1, "id", "basic", N_("Protect a variable from being modified. It will be treated as a system defined variable from now on. Protected parameters can still be modified."));
+ FUNC (unprotect, 1, "id", "basic", N_("Unprotect a variable from being modified. It will be treated as a user defined variable from now on."));
VFUNC (SetFunctionFlags, 2, "id,flags", "basic", N_("Set flags for a function, currently \"PropagateMod\" and \"NoModuloArguments\""));
FUNC (GetCurrentModulo, 0, "", "basic", N_("Get current modulo from the context outside the function"));
FUNC (IsDefined, 1, "id", "basic", N_("Check if a variable or function is defined"));
- FUNC (undefine, 1, "id", "basic", N_("Undefine a variable (including locals and globals)"));
+ FUNC (undefine, 1, "id", "basic", N_("Undefine a variable (including all locals and globals of the same name)"));
+ ALIAS (Undefine, 1, undefine);
+ FUNC (UndefineAll, 0, "", "basic", N_("Undefine all unprotected (user defined) global variables and parameters. Does not reset or change protected (system) parameters."));
+ FUNC (ProtectAll, 0, "", "basic", N_("Mark all currently defined variables as protected. They will be treated as system defined variables from now on."));
FUNC (Parse, 1, "str", "basic", N_("Parse a string (but do not execute)"));
FUNC (Evaluate, 1, "str", "basic", N_("Parse and evaluate a string"));
diff --git a/src/testscope.gel b/src/testscope.gel
index b11e2ba..948de86 100644
--- a/src/testscope.gel
+++ b/src/testscope.gel
@@ -46,5 +46,80 @@ function f(x) = (
);
if f(0) != 1+5 then (error("global lookup over all local failed");exit());
+A = 1;
+x = 9;
+function f(x) = (
+ local *;
+ A = 8;
+ function g(y) = A+5+x;
+ g(x)
+);
+if f(0) != 1+5+9 then (error("global lookup over all local (2) failed");exit());
+
+A = 1;
+function f(x) = (
+ A = 8;
+ function g(y) = A+5;
+ g
+);
+h = f(0);
+function f2(g) = (
+ A=-9;
+ g(0)
+);
+if f2(h) != 8+5 then (error("subst pass up 1 failed");exit());
+
+A = 1;
+function f(x) = (
+ A = 8;
+ function g(y) = A+5;
+ g
+);
+h = f(0);
+function f2(g) = (
+ A=-9;
+ g(5);
+ g
+);
+h2=f2(h);
+if h2(0) != 8+5 then (error("subst pass up 2 failed");exit());
+
+function f(x) = (
+ local *;
+ function g(y) = 5;
+ g
+);
+h = f(0);
+function f2(x) = (
+ h(9)
+);
+if f2(0) != 5 then (error("ret of local failed");exit());
+
+function f(x) = (
+ local g;
+ function g(y) = 19;
+ g
+);
+h = f(0);
+function f2(x) = (
+ h(9)
+);
+if f2(0) != 19 then (error("ret of local 2 failed");exit());
+
+
+function ff() = (
+ function f(x) = (
+ local g;
+ function g(y) = 999;
+ g
+ );
+ h = f(0);
+ function f2(x) = (
+ h(9)
+ );
+ if f2(0) != 999 then (error("ret of local failed");exit());
+);
+ff();
+
print("true");
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
