genius r642 - in trunk: . lib lib/equation_solving po src

[jirka svn gnome org]

Author: jirka
Date: Fri Feb 22 07:42:30 2008
New Revision: 642
URL: http://svn.gnome.org/viewvc/genius?rev=642&view=rev

Log:

Fri Feb 22 01:38:54 2008  Jiri (George) Lebl <jirka 5z com>

	* src/mpwrap.[ch]: add mpw_re_sgn and mpw_im_sgn functions and fix
	  memory leak in mpw_re and mpw_im

	* src/funclib.c: Implement QuadraticFormula internally, marginally
	  improving performance, but mainly handle special cases better
	  and avoid instability in solutions (avoid bad cancellation in most
	  cases)

	* lib/equation_solving/formulas.gel: remove QuadraticFormula from
	  here

	* src/longtest.gel: add tests for quadratic formula



Modified:
   trunk/ChangeLog
   trunk/lib/equation_solving/formulas.gel
   trunk/lib/library-strings.c
   trunk/po/ChangeLog
   trunk/po/cs.po
   trunk/src/funclib.c
   trunk/src/longtest.gel
   trunk/src/mpwrap.c
   trunk/src/mpwrap.h

Modified: trunk/lib/equation_solving/formulas.gel
==============================================================================
--- trunk/lib/equation_solving/formulas.gel	(original)
+++ trunk/lib/equation_solving/formulas.gel	Fri Feb 22 07:42:30 2008
@@ -1,13 +1,3 @@
-SetHelp ("QuadraticFormula", "equation_solving",
-         "Find roots of a quadratic polynomial (given as vector of coefficients)")
-function QuadraticFormula(p) = (
-	if not IsPoly(p) or elements(p) != 3 or p@(3) == 0 then
-		(error("QuadraticFormula: argument 1 must be a quadratic polynomial");bailout);
-	[ (- p@(2) + sqrt(p@(2)^2 - 4*p@(1)*p@(3))) / (2 * p@(3)) 
-	 (- p@(2) - sqrt(p@(2)^2 - 4*p@(1)*p@(3))) / (2 * p@(3))]
-)
-protect ("QuadraticFormula")
-
 # see http://en.wikipedia.org/wiki/Cubic_equation
 SetHelp ("CubicFormula", "equation_solving",
          "Find roots of a cubic polynomial (given as vector of coefficients)")

Modified: trunk/lib/library-strings.c
==============================================================================
--- trunk/lib/library-strings.c	(original)
+++ trunk/lib/library-strings.c	Fri Feb 22 07:42:30 2008
@@ -221,7 +221,6 @@
 char *fake = N_("Find root of a function using the Muller's method");
 char *fake = N_("Find root of a function using the secant method");
 char *fake = N_("Find roots of a polynomial (given as vector of coefficients)");
-char *fake = N_("Find roots of a quadratic polynomial (given as vector of coefficients)");
 char *fake = N_("Find roots of a quartic polynomial (given as vector of coefficients)");
 char *fake = N_("Use classical non-adaptive Runge-Kutta of fourth order method to numerically solve y'=f(x,y) for initial x0,y0 going to x1 with n increments, returns y at x1");
 char *fake = N_("Calculate average of an entire matrix");

Modified: trunk/src/funclib.c
==============================================================================
--- trunk/src/funclib.c	(original)
+++ trunk/src/funclib.c	Fri Feb 22 07:42:30 2008
@@ -4359,6 +4359,87 @@
 }
 
 static GelETree *
+QuadraticFormula_op(GelCtx *ctx, GelETree * * a, gboolean *exception)
+{
+	GelETree *n;
+	GelMatrixW *m;
+	GelMatrix *nm;
+
+	mpw_ptr aa, bb, cc;
+	mpw_t r1, r2;
+
+	if G_UNLIKELY (a[0]->type == NULL_NODE)
+		return gel_makenum_null ();
+
+	if G_UNLIKELY ( ! check_poly(a,1,"QuadraticFormula",TRUE))
+		return NULL;
+
+	m = a[0]->mat.matrix;
+
+	if G_UNLIKELY (gel_matrixw_elements (m) != 3 ||
+		       mpw_zero_p (gel_matrixw_index(m,2,0)->val.value)) {
+		gel_errorout (_("%s: argument 1 must be a quadratic polynomial"),
+			      "QuadraticFormula");
+		return NULL;
+	}
+
+	aa = gel_matrixw_index(m,2,0)->val.value;
+	bb = gel_matrixw_index(m,1,0)->val.value;
+	cc = gel_matrixw_index(m,0,0)->val.value;
+
+	mpw_init (r1);
+	mpw_init (r2);
+
+	if (mpw_zero_p (cc)) {
+		mpw_div (r1, bb, aa);
+		mpw_neg (r1, r1);
+		mpw_set_ui (r2, 0);
+	} else if (mpw_zero_p (bb)) {
+		mpw_mul (r1, aa, cc);
+		mpw_neg (r1, r1);
+		mpw_sqrt (r1, r1);
+		mpw_div (r1, r1, aa);
+		mpw_neg (r2, r1);
+	} else {
+		mpw_mul (r1, bb, bb);
+		mpw_mul (r2, aa, cc);
+		mpw_mul_ui (r2, r2, 4);
+		mpw_sub (r1, r1, r2);
+		mpw_sqrt (r1, r1);
+		/* r1 is now the sqrt of the discriminant */
+
+		/* try to avoid instability */
+		if (mpw_re_sgn (r1) != mpw_re_sgn (bb)) {
+			mpw_neg (r1, r1);
+		}
+
+		mpw_add (r1, r1, bb);
+		mpw_div_ui (r1, r1, 2);
+		mpw_neg (r1, r1);
+
+		/* r1 = (bb + s * sqrt(bb^2 - 4*aa*cc)) / (-2); */
+
+		/* set r2 first */
+		mpw_div (r2, cc, r1);
+
+		mpw_div (r1, r1, aa);
+	}
+
+	nm = gel_matrix_new ();
+	gel_matrix_set_size (nm, 1, 2, FALSE /* padding */);
+	gel_matrix_index (nm, 0, 0) = gel_makenum_use (r1);
+	gel_matrix_index (nm, 0, 1) = gel_makenum_use (r2);
+
+	GET_NEW_NODE (n);
+	n->type = MATRIX_NODE;
+	n->mat.matrix = gel_matrixw_new_with_matrix (nm);
+	n->mat.quoted = FALSE;
+
+	return n;
+}
+
+
+static GelETree *
 PolyToString_op (GelCtx *ctx, GelETree * * a, gboolean *exception)
 {
 	GelETree *n;
@@ -6073,6 +6154,8 @@
 	VFUNC (PolyToString, 2, "p,var", "polynomial", N_("Make string out of a polynomial (as vector)"));
 	FUNC (PolyToFunction, 1, "p", "polynomial", N_("Make function out of a polynomial (as vector)"));
 
+	FUNC (QuadraticFormula, 1, "p", "equation_solving", N_("Find roots of a quadratic polynomial (given as vector of coefficients)"));
+
 	FUNC (Combinations, 2, "k,n", "combinatorics", N_("Get all combinations of k numbers from 1 to n as a vector of vectors"));
 	FUNC (NextCombination, 2, "v,n", "combinatorics", N_("Get combination that would come after v in call to combinations, first combination should be [1:k]."));
 	FUNC (Permutations, 2, "k,n", "combinatorics", N_("Get all permutations of k numbers from 1 to n as a vector of vectors"));

Modified: trunk/src/longtest.gel
==============================================================================
--- trunk/src/longtest.gel	(original)
+++ trunk/src/longtest.gel	Fri Feb 22 07:42:30 2008
@@ -92,6 +92,31 @@
   true
 );
 
+function roottestquad() = (
+  for n=1 to 100 do (
+    if n <= 50 then
+    	p = randint(10,3)-5 + 1i* (randint(10,3)-5)
+    else
+        p = randint(10,3)-5;
+    if p@(3) == 0 then
+	p@(3) = 1;
+    #print (p);
+    r = QuadraticFormula (p);
+    #print (r);
+    f = PolyToFunction (p);
+    if (|f(r@(1))| >= 0.01 or
+	|f(r@(2))| >= 0.01) then (
+	#print (n);
+	#print (p);
+	#print (r);
+	#print (f(r@(1)));
+	#print (f(r@(2)));
+      return false;
+    )
+  );
+  true
+);
+
 function roottestquart() = (
   for n=1 to 100 do (
     if n <= 50 then
@@ -176,6 +201,7 @@
 	if not derivtest() then (error("error on deriv test");errors = errors + 1);
 
 	#roots
+	if not roottestquad() then (error("error on root test quadratic");errors = errors + 1);
 	if not roottestcube() then (error("error on root test cubic");errors = errors + 1);
 	if not roottestquart() then (error("error on root test quartic");errors = errors + 1);
 	if not nthroottest() then (error("error on nth root test");errors = errors + 1);

Modified: trunk/src/mpwrap.c
==============================================================================
--- trunk/src/mpwrap.c	(original)
+++ trunk/src/mpwrap.c	Fri Feb 22 07:42:30 2008
@@ -3357,6 +3357,18 @@
 	return 0;
 }
 
+int
+mpw_re_sgn(mpw_ptr op)
+{
+	return mpwl_sgn(op->r);
+}
+
+int
+mpw_im_sgn(mpw_ptr op)
+{
+	return mpwl_sgn(op->i);
+}
+
 void
 mpw_abs(mpw_ptr rop,mpw_ptr op)
 {
@@ -5363,17 +5375,28 @@
 void
 mpw_im(mpw_ptr rop, mpw_ptr op)
 {
+	if (rop == op) {
+		MAKE_IMAG(rop);
+		rop->r = rop->i;
+		rop->i = gel_zero;
+		return;
+	}
 	MAKE_REAL(rop);
-	rop->r=op->i;
-	op->i->alloc.usage++;
+	DEALLOC_MPWL (rop->r);
+	rop->r = op->i;
+	ALLOC_MPWL (rop->r);
 }
 
 void
 mpw_re(mpw_ptr rop, mpw_ptr op)
 {
 	MAKE_REAL(rop);
-	rop->r=op->r;
-	op->r->alloc.usage++;
+	if (rop == op) {
+		return;
+	}
+	DEALLOC_MPWL (rop->r);
+	rop->r = op->r;
+	ALLOC_MPWL (rop->r);
 }
 
 void

Modified: trunk/src/mpwrap.h
==============================================================================
--- trunk/src/mpwrap.h	(original)
+++ trunk/src/mpwrap.h	Fri Feb 22 07:42:30 2008
@@ -154,6 +154,11 @@
 
 int mpw_sgn(mpw_ptr op);
 
+/* sign of the real part */
+int mpw_re_sgn(mpw_ptr op);
+/* sign of the im part */
+int mpw_im_sgn(mpw_ptr op);
+
 void mpw_neg(mpw_ptr rop,mpw_ptr op);
 
 void mpw_add(mpw_ptr rop,mpw_ptr op1, mpw_ptr op2);