I'd like to commit the attached patch. It moves functionality to test
approximate equality from elemnum.input to unittest.spad, which turns out tpo
be a good thing, because there was a mistake in elemnum (epsilonrelative
instead of epsilonRelative) that went unnoticed for many revisions, and made us
believe that numerics is better than it was.
I also add a test for branch cuts that shows another feature of the new
functionality: useful error messages.
Waldek, I test values of asin on a circle around the origin with radius 2. Is
this what you had in mind?
Martin
Index: algebra/unittest.spad.pamphlet
===================================================================
--- algebra/unittest.spad.pamphlet (revision 442)
+++ algebra/unittest.spad.pamphlet (working copy)
@@ -334,6 +334,27 @@
++ R performs some simplifications, we convert ex1 and ex2 to
++ \axiom{InputForm}, if possible.
+ if R has RealNumberSystem then
+ testRealEqualsAux: (String, R, R) -> Void
+ ++ testRealEquals(inp, ex1, ex2) states that ex1 and ex2 should be
+ ++ approximately equal as real numbers, taking into acount
+ ++ \spadfun{testAbsolutePrecision} and
+ ++ \spadfun{testRelativePrecision}.
+
+ testComplexEqualsAux: (String, Complex R, Complex R) -> Void
+ ++ testComplexEquals(inp ex1, ex2) states that ex1 and ex2 should be
+ ++ approximately equal as complex numbers, taking into acount
+ ++ \spadfun{testAbsolutePrecision} and
+ ++ \spadfun{testRelativePrecision}.
+
+ testAbsolutePrecision: R -> Void
+ ++ testAbsolutePrecision(eps) sets the absolute precision used for
+ ++ floating point comparison. The default value is 1.0e-15.
+
+ testRelativePrecision: R -> Void
+ ++ testRelativePrecision(eps) sets the relative precision used for
+ ++ floating point comparison. The default value is 1.0e-30.
+
testNotEqualsAux: (String, R, R) -> Void
++ testNotEquals(inp, ex1, ex2) states that ex1 and ex2 should be
++ different.
@@ -345,7 +366,6 @@
++ testLibraryError(inp, ex) states that ex should throw an error. Such
++ a test will never count as a fatal error.
-
== add
O ==> OutputForm
@@ -369,7 +389,57 @@
then incFail(inp, [ex1::O, ex2::O])$U
else incFail(inp, [])$U
+ if R has RealNumberSystem then
+ epsilonRelative := 10::R**(-15)
+ epsilonAbsolute := 10::R**(-15)
+
+ testRelativePrecision eps ==
+ epsilonRelative := eps
+
+ testAbsolutePrecision eps ==
+ epsilonAbsolute := eps
+
+ testRealEqualsAux(inp, ex1, ex2) ==
+ relative: R := norm(ex1/ex2-1)$R
+ absolute: R := norm(ex1-ex2)$R
+ res: List O := []
+ if relative > epsilonRelative then
+ if R has CoercibleTo O
+ then res := cons(hconcat(message("rel err: ")$O,
+ relative::O), res)
+ else res := cons(message("rel err ")$O, res)
+
+ if absolute > epsilonAbsolute then
+ if R has CoercibleTo O
+ then res := cons(hconcat(message("abs err: ")$O,
+ absolute::O), res)
+ else res := cons(message("abs err ")$O, res)
+
+ if empty? res
+ then incPass()$U
+ else incFail(inp, concat([ex1::O, ex2::O], res))$U
+
+ testComplexEqualsAux(inp, ex1, ex2) ==
+ relative := norm(ex1/ex2-1)$Complex(R)
+ absolute := norm(ex1-ex2)$Complex(R)
+ res: List O := []
+ if relative > epsilonRelative then
+ if R has CoercibleTo O
+ then res := cons(hconcat(message("rel err: ")$O,
+ relative::O), res)
+ else res := cons(message("rel err ")$O, res)
+
+ if absolute > epsilonAbsolute then
+ if R has CoercibleTo O
+ then res := cons(hconcat(message("abs err: ")$O,
+ absolute::O), res)
+ else res := cons(message("abs err ")$O, res)
+
+ if empty? res
+ then incPass()$U
+ else incFail(inp, concat([ex1::O, ex2::O], res))$U
+
testNotEqualsAux(inp, ex1, ex2) ==
if (ex1=ex2)@Boolean
then
@@ -412,6 +482,16 @@
++ performs some simplifications, we convert ex1 and ex2 to
++ \axiom{InputForm}, if possible
+ testRealEquals: (String, String) -> Void
+ ++ testRealEquals(ex1, ex2) states that ex1 and ex2 should be
+ ++ approximately equal as real numbers, taking into acount
+ ++ \spadfun{testAbsolutePrecision} and \spadfun{testRelativePrecision}.
+
+ testComplexEquals: (String, String) -> Void
+ ++ testComplexEquals(ex1, ex2) states that ex1 and ex2 should be
+ ++ approximately equal as complex numbers, taking into acount
+ ++ \spadfun{testAbsolutePrecision} and \spadfun{testRelativePrecision}.
+
testNotEquals: (String, String) -> Void
++ testNotEquals(ex1, ex2) states that ex1 and ex2 should be
++ different.
@@ -433,6 +513,18 @@
interpretString("testEqualsAux(_"" inp "_", " s1 ", " s2 ")")$T
decFatal()$UnittestCount
+ testRealEquals(s1, s2) ==
+ inp := "EQLR: (" s1 ", " s2 ")"
+ incFatal(inp)$UnittestCount
+ interpretString("testRealEqualsAux(_"" inp "_", " s1 ", " s2 ")")$T
+ decFatal()$UnittestCount
+
+ testComplexEquals(s1, s2) ==
+ inp := "EQLC: (" s1 ", " s2 ")"
+ incFatal(inp)$UnittestCount
+ interpretString("testComplexEqualsAux(_"" inp "_", " s1 ", " s2 ")")$T
+ decFatal()$UnittestCount
+
testNotEquals(s1, s2) ==
inp := "DIFF: (" s1 ", " s2 ")"
incFatal(inp)$UnittestCount
Index: input/elemnum.input.pamphlet
===================================================================
--- input/elemnum.input.pamphlet (revision 442)
+++ input/elemnum.input.pamphlet (working copy)
@@ -21,17 +21,6 @@
-- SMW June 26, 1991.
testsuite "elemnum"
-epsilonAbsolute := 1.0e-15
-epsilonrelative := 1.0e-30
-
-testFloat(a: String, b: String): Void ==
- testTrue("abs(" a "-" b") < epsilonAbsolute")
- testTrue("abs(" a "/" b "-1) < epsilonRelative")
-
-testComplexFloat(a: String, b: String): Void ==
- testTrue("norm(" a "-" b") < epsilonAbsolute")
- testTrue("norm(" a "/" b "-1) < epsilonRelative")
-
inverses := _
["exp log x", _
"log exp x", _
@@ -53,8 +42,9 @@
errlist := ["acsc x", "asec x", "acosh x", "acoth x"]
for ex in inverses | not anySubstring?(errlist, ex) repeat
- testFloat(ex, "x")
+ testRealEquals(ex, "x")
+-- cannot loop over errlist, since error breaks the loop
testLibraryError "acsc x"
testLibraryError "asec x"
testLibraryError "acosh x"
@@ -64,7 +54,7 @@
testcaseNoClear "0.7 DoubleFloat"
x := 0.7::DoubleFloat
for ex in inverses | not anySubstring?(errlist, ex) repeat
- testFloat(ex, "x")
+ testRealEquals(ex, "x")
testLibraryError "acsc x"
testLibraryError "asec x"
@@ -75,7 +65,7 @@
x := 1.1::Float
errlist := ["asin x", "acos x", "atanh x", "asech x"]
for ex in inverses | not anySubstring?(errlist, ex) repeat
- testFloat(ex, "x")
+ testRealEquals(ex, "x")
testLibraryError "asin x"
testLibraryError "acos x"
@@ -85,7 +75,7 @@
testcaseNoClear "1.1 DoubleFloat"
x := 1.1::DoubleFloat
for ex in inverses | not anySubstring?(errlist, ex) repeat
- testFloat(ex, "x")
+ testRealEquals(ex, "x")
testLibraryError "asin x"
testLibraryError "acos x"
@@ -103,34 +93,67 @@
for ex in inverses repeat
if member?(ex, ["acos cos x", "asec sec x", _
"acosh cosh x", "asech sech x"])
- then testComplexFloat(concat [string s1, "*", ex], "x")
+ then testComplexEquals(concat [string s1, "*", ex], "x")
else if ex = "acot cot x"
- then testComplexFloat(concat [ex, " - ", string m1, "*%pi"], "x")
- else testComplexFloat(ex, "x")
+ then testComplexEquals(concat [ex, " - ", string m1, "*%pi"], "x")
+ else testComplexEquals(ex, "x")
-sa := 0.7::DoubleFloat;
-sb := 1.1::DoubleFloat;
-ba := 0.7::Float;
-bb := 1.1::Float;
-
testcaseNoClear "Complex Float"
)clear value x
x:Complex Float
-testinverseFloats(ba, bb, 1);
-testinverseFloats(ba, bb, 2);
-testinverseFloats(ba, bb, 3);
-testinverseFloats(ba, bb, 4);
+ba := 0.7::Float;
+bb := 1.1::Float;
+testinverseFloats(ba, bb, 1); -- test 0.7 + 1.1*%i
+testinverseFloats(ba, bb, 2); -- test -0.7 + 1.1*%i
+testinverseFloats(ba, bb, 3); -- test -0.7 - 1.1*%i
+testinverseFloats(ba, bb, 4); -- test 0.7 - 1.1*%i
testcaseNoClear "Complex DoubleFloat"
)clear value x
x:Complex DFLOAT
+sa := 0.7::DoubleFloat;
+sb := 1.1::DoubleFloat;
testinverseFloats(sa, sb, 1);
testinverseFloats(sa, sb, 2);
testinverseFloats(sa, sb, 3);
testinverseFloats(sa, sb, 4);
+testsuite "branch cuts"
+testcaseNoClear "asin"
+
+DirD := eval(D(asin(x0+(x1-x0)*t), t)::EXPR FLOAT, t=0)
+tmax := 10
+tdiff := 1.0e-8
+
+-- t = 0..tmax
+-- x = -2..2
+X t == -2.0 + t/tmax * 4.0
+y x == sqrt(-x^2 + 4.0)/2.0
+
+-- t = 0: arg1 is on the negative x-axis, arg2 in quadrant II
+for t in 0..tmax-1 repeat
+ arg1 := (X t + %i * y X t)::Complex Float
+ arg2 := (X(t+tdiff) + %i * y X(t+tdiff))::Complex Float
+ corr := (eval(DirD, [x0=arg1, x1=arg2]))::Complex Float
+ testComplexEquals("asin arg1 + corr", "asin arg2")
+ testComplexEquals("asin(arg1::Complex DFLOAT) + corr",
+ "asin(arg2::Complex DFLOAT)")
+
+X t == 2.0 - t/tmax * 4.0
+y x == -sqrt(-x^2 + 4.0)/2.0
+
+-- t = tmax: arg1 is on the positive x-axis, arg2 in quadrant IV
+for t in 0..tmax-1 repeat
+ arg1 := (X t + %i * y X t)::Complex Float
+ arg2 := (X(t+tdiff) + %i * y X(t+tdiff))::Complex Float
+ corr := (eval(DirD, [x0=arg1, x1=arg2]))::Complex Float
+ testComplexEquals("asin arg1 + corr", "asin arg2")
+ testComplexEquals("asin(arg1::Complex DFLOAT) + corr",
+ "asin(arg2::Complex DFLOAT)")
+
expected() ==
messagePrint("testsuite | testcases: failed (total) | tests: failed
(total)")$OutputForm;
- messagePrint("elemnum 0 (6) 0
(608)")$OutputForm
+ messagePrint("elemnum 0 (6) 0
(312)")$OutputForm;
+ messagePrint("branchcuts 1 (1) 3
(40)")$OutputForm
statistics()
expected()
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"FriCAS - computer algebra system" group.
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at
http://groups.google.com/group/fricas-devel?hl=en
-~----------~----~----~----~------~----~------~--~---
