genius r653 - in trunk: . src

From: jirka svn gnome org
To: svn-commits-list gnome org
Subject: genius r653 - in trunk: . src
Date: Tue, 20 May 2008 05:25:02 +0000 (UTC)
Author: jirka
Date: Tue May 20 05:25:01 2008
New Revision: 653
URL: http://svn.gnome.org/viewvc/genius?rev=653&view=rev

Log:

Tue May 20 00:25:10 2008  Jiri (George) Lebl <jirka 5z com>

	* src/matrixw.[ch], matop.c, src/eval.c, src/funclib.c:  Optimize
	  matrix manipulation a bit.  Don't be overly conservative with
	  making things private.  And OMG we used Gauss-Jordan instead of
	  backaddition.  Also when the matrix is not rational, do pivotting
	  (use largest entry, not first nonzero one)

	* src/mpwrap.c: fix mpw_abs if called with same arguments for in and
	  out and add mpw_abs_sq for getting the absolute value squared,
	  which doesn't involve a sqrt

	* src/genius.c, src/gnome-genius.c, src/eval.[ch], src/matrixw.c: 
	  init the_zero during the startup

	* src/mpzextra.c: very minor optimizations

	* src/Makefile.am: fix BUILDDIR setup



Modified:
   trunk/ChangeLog
   trunk/NEWS
   trunk/src/Makefile.am
   trunk/src/eval.c
   trunk/src/eval.h
   trunk/src/funclib.c
   trunk/src/genius.c
   trunk/src/gnome-genius.c
   trunk/src/matop.c
   trunk/src/matrixw.c
   trunk/src/matrixw.h
   trunk/src/mpwrap.c
   trunk/src/mpwrap.h
   trunk/src/mpzextra.c

Modified: trunk/NEWS
==============================================================================
--- trunk/NEWS	(original)
+++ trunk/NEWS	Tue May 20 05:25:01 2008
@@ -6,6 +6,9 @@
 * QuadraticFormula built in and more numerically stable
 * CHANGE: It's Fibonacci in correct spelling, short name is still fib
 * Calling internal functions is now slightly faster
+* Gaussian elimination is now faster, and more stable when nonrational
+  matrices are involved
+* Other minor optimizations
 * Fix crash related to returning custom functions from functions
 * Fix some memory leaks
 * Documentation updates

Modified: trunk/src/Makefile.am
==============================================================================
--- trunk/src/Makefile.am	(original)
+++ trunk/src/Makefile.am	Tue May 20 05:25:01 2008
@@ -13,7 +13,7 @@
 	-DG_LOG_DOMAIN=\"Genius\"				\
 	-DDATADIR=\""$(datadir)"\"				\
 	-DLIBEXECDIR=\""$(libexecdir)"\"			\
-	-DBUILDDIR=\""$(abs_top_builddir)"\"			\
+	-DBUILDDIR=\""@abs_top_builddir@"\"			\
 	-I$(srcdir)						\
 	-I$(top_srcdir)						\
 	-I$(top_srcdir)/ve					\

Modified: trunk/src/eval.c
==============================================================================
--- trunk/src/eval.c	(original)
+++ trunk/src/eval.c	Tue May 20 05:25:01 2008
@@ -259,6 +259,13 @@
 	return 0;
 }
 
+void
+gel_init (void)
+{
+	if (the_zero == NULL)
+		the_zero = gel_makenum_ui (0);
+}
+
 mpw_ptr
 gel_find_pre_function_modulo (GelCtx *ctx)
 {
@@ -987,7 +994,7 @@
 		} else if(gel_matrixw_width(et->mat.matrix) == 1) {
 			int xx;
 			int h = gel_matrixw_height(et->mat.matrix);
-			gel_matrixw_make_private(et->mat.matrix);
+			gel_matrixw_make_private (et->mat.matrix, FALSE /* kill_type_caches */);
 			for(x=0;x<h;x++) {
 				gel_matrix_index(dest,i,di+x) =
 					gel_matrixw_get_index(et->mat.matrix,0,x);
@@ -1007,7 +1014,8 @@
 			int h = gel_matrixw_height(et->mat.matrix);
 			int w = gel_matrixw_width(et->mat.matrix);
 
-			gel_matrixw_make_private(et->mat.matrix);
+			gel_matrixw_make_private(et->mat.matrix,
+						 FALSE /* kill_type_caches */);
 
 			for(x=0;x<h;x++) {
 				GelETree *n;
@@ -1198,7 +1206,7 @@
 		/* never should be reached */
 	}
 
-	gel_matrixw_make_private (nm);
+	gel_matrixw_make_private (nm, FALSE /* kill_type_caches */);
 
 	m = gel_matrix_new();
 	gel_matrix_set_size(m, w, h, TRUE /* padding */);
@@ -1767,7 +1775,7 @@
 		node = l;
 	}
 
-	gel_matrixw_make_private(m);
+	gel_matrixw_make_private(m, TRUE /* kill_type_caches */);
 
 	for(j=0;j<gel_matrixw_height(m);j++) {
 		for(i=0;i<gel_matrixw_width(m);i++) {
@@ -1829,7 +1837,7 @@
 		return TRUE;
 	}
 
-	gel_matrixw_make_private(m);
+	gel_matrixw_make_private(m, TRUE /* kill_type_caches */);
 
 	for (i = 0; i < gel_matrixw_width (m); i++) {
 		GelETree *t = gel_matrixw_get_index(m,i,i);
@@ -1860,7 +1868,7 @@
 	int i,j;
 	GelMatrixW *m = l->mat.matrix;
 
-	gel_matrixw_make_private(m);
+	gel_matrixw_make_private(m, TRUE /* kill_type_caches */);
 
 	for(j=0;j<gel_matrixw_height(m);j++) {
 		for(i=0;i<gel_matrixw_width(m);i++) {
@@ -1927,7 +1935,7 @@
 		return TRUE;
 	}
 	l->mat.quoted = l->mat.quoted || r->mat.quoted;
-	gel_matrixw_make_private(m1);
+	gel_matrixw_make_private(m1, TRUE /* kill_type_caches */);
 	for(j=0;j<gel_matrixw_height(m1);j++) {
 		for(i=0;i<gel_matrixw_width(m1);i++) {
 			GelETree *t = gel_matrixw_get_index (m1, i, j);
@@ -2315,7 +2323,7 @@
 	int w,h;
 
 	/*make us a private copy!*/
-	gel_matrixw_make_private(m);
+	gel_matrixw_make_private(m, TRUE /* kill_type_caches */);
 
 	w = gel_matrixw_width (m);
 	h = gel_matrixw_height (m);
@@ -2394,11 +2402,9 @@
 static gboolean
 conjugate_transpose_matrix (GelCtx *ctx, GelETree *n, GelETree *l)
 {
-	if (gel_is_matrix_value_only_real (l->mat.matrix)) {
-		l->mat.matrix->tr = !(l->mat.matrix->tr);
-	} else {
-		gel_matrix_conjugate_transpose (l->mat.matrix);
-	}
+	/* handles real case nicely */
+	gel_matrix_conjugate_transpose (l->mat.matrix);
+
 	/*remove from arglist*/
 	n->op.args = NULL;
 	replacenode(n,l);
@@ -4408,7 +4414,7 @@
 			    t->type != USERTYPE_NODE) {
 				if ( ! pushed) {
 					/*make us a private copy!*/
-					gel_matrixw_make_private (m);
+					gel_matrixw_make_private (m, TRUE /* kill_type_caches */);
 
 					/* it will be a copy */
 					t = gel_matrixw_get_index (m, x, y);
@@ -6651,12 +6657,13 @@
 			}
 		}
 	} else if (n->type == MATRIX_NODE &&
-		   n->mat.matrix != NULL) {
+		   n->mat.matrix != NULL &&
+		   ! gel_is_matrix_value_only (n->mat.matrix)) {
 		int i,j;
 		int w,h;
 		w = gel_matrixw_width (n->mat.matrix);
 		h = gel_matrixw_height (n->mat.matrix);
-		gel_matrixw_make_private (n->mat.matrix);
+		gel_matrixw_make_private (n->mat.matrix, TRUE /* kill_type_caches */);
 		for (i = 0; i < w; i++) {
 			for(j = 0; j < h; j++) {
 				GelETree *t = gel_matrixw_get_index
@@ -7028,11 +7035,13 @@
 	} else if(n->type==MATRIX_NODE) {
 		int i,j;
 		int w,h;
-		if(!n->mat.matrix)
+		if(!n->mat.matrix ||
+		   gel_is_matrix_value_only (n->mat.matrix)) {
 			goto gather_comparisons_end;
+		}
 		w = gel_matrixw_width(n->mat.matrix);
 		h = gel_matrixw_height(n->mat.matrix);
-		gel_matrixw_make_private(n->mat.matrix);
+		gel_matrixw_make_private(n->mat.matrix, TRUE /* kill_type_caches */);
 		for(j=0;j<h;j++) {
 			for(i=0;i<w;i++) {
 				GelETree *t = gel_matrixw_get_index(n->mat.matrix,i,j);
@@ -7102,12 +7111,13 @@
 			}
 		}
 	} else if (n->type == MATRIX_NODE &&
-		   n->mat.matrix != NULL) {
+		   n->mat.matrix != NULL &&
+		   ! gel_is_matrix_value_only (n->mat.matrix)) {
 		int i,j;
 		int w,h;
 		w = gel_matrixw_width (n->mat.matrix);
 		h = gel_matrixw_height (n->mat.matrix);
-		gel_matrixw_make_private (n->mat.matrix);
+		gel_matrixw_make_private (n->mat.matrix, TRUE /* kill_type_caches */);
 		for(j = 0; j < h; j++) {
 			for (i = 0; i < w; i++) {
 				GelETree *t = gel_matrixw_get_index
@@ -7153,12 +7163,13 @@
 			args = args->any.next;
 		}
 	} else if (n->type == MATRIX_NODE &&
-		   n->mat.matrix != NULL) {
+		   n->mat.matrix != NULL &&
+		   ! gel_is_matrix_value_only (n->mat.matrix)) {
 		int i,j;
 		int w,h;
 		w = gel_matrixw_width (n->mat.matrix);
 		h = gel_matrixw_height (n->mat.matrix);
-		gel_matrixw_make_private (n->mat.matrix);
+		gel_matrixw_make_private (n->mat.matrix, TRUE /* kill_type_caches */);
 		for(j = 0; j < h; j++) {
 			for (i = 0; i < w; i++) {
 				GelETree *t = gel_matrixw_get_index
@@ -7212,12 +7223,13 @@
 			}
 		}
 	} else if (n->type == MATRIX_NODE &&
-		   n->mat.matrix != NULL) {
+		   n->mat.matrix != NULL &&
+		   ! gel_is_matrix_value_only (n->mat.matrix)) {
 		int i,j;
 		int w,h;
 		w = gel_matrixw_width (n->mat.matrix);
 		h = gel_matrixw_height (n->mat.matrix);
-		gel_matrixw_make_private (n->mat.matrix);
+		gel_matrixw_make_private (n->mat.matrix, TRUE /* kill_type_caches */);
 		for(j = 0; j < h; j++) {
 			for (i = 0; i < w; i++) {
 				GelETree *t = gel_matrixw_get_index
@@ -7460,10 +7472,12 @@
 	} else if(n->type==MATRIX_NODE) {
 		int i,j;
 		int w,h;
-		if(!n->mat.matrix) return;
+		if (n->mat.matrix == NULL ||
+		    gel_is_matrix_value_only (n->mat.matrix))
+			return;
 		w = gel_matrixw_width(n->mat.matrix);
 		h = gel_matrixw_height(n->mat.matrix);
-		gel_matrixw_make_private(n->mat.matrix);
+		gel_matrixw_make_private (n->mat.matrix, TRUE /* kill_type_caches */);
 		for(j=0;j<h;j++) {
 			for(i=0;i<w;i++) {
 				GelETree *t = gel_matrixw_get_index(n->mat.matrix,i,j);
@@ -7944,10 +7958,12 @@
 	} else if(n->type==MATRIX_NODE) {
 		int i,j;
 		int w,h;
-		if(!n->mat.matrix) return;
+		if (n->mat.matrix == NULL ||
+		    gel_is_matrix_value_only (n->mat.matrix))
+			return;
 		w = gel_matrixw_width(n->mat.matrix);
 		h = gel_matrixw_height(n->mat.matrix);
-		gel_matrixw_make_private(n->mat.matrix);
+		gel_matrixw_make_private (n->mat.matrix, TRUE /* kill_type_caches */);
 		for(j=0;j<h;j++) {
 			for(i=0;i<w;i++) {
 				GelETree *t = gel_matrixw_get_index(n->mat.matrix,i,j);

Modified: trunk/src/eval.h
==============================================================================
--- trunk/src/eval.h	(original)
+++ trunk/src/eval.h	Tue May 20 05:25:01 2008
@@ -1,5 +1,5 @@
 /* GENIUS Calculator
- * Copyright (C) 1997-2007 Jiri (George) Lebl
+ * Copyright (C) 1997-2008 Jiri (George) Lebl
  *
  * Author: Jiri (George) Lebl
  *
@@ -128,6 +128,8 @@
 	GelOperPrim prim[OP_TABLE_LEN];
 };
 
+void gel_init (void);
+
 /*functions for manipulating a tree*/
 GelETree * gel_makenum(mpw_t num);
 GelETree * gel_makenum_use(mpw_t num); /*don't create a new number*/

Modified: trunk/src/funclib.c
==============================================================================
--- trunk/src/funclib.c	(original)
+++ trunk/src/funclib.c	Tue May 20 05:25:01 2008
@@ -200,7 +200,7 @@
 
 	GET_NEW_NODE (n);
 	n->type = MATRIX_NODE;
-	n->mat.matrix = gel_matrixw_new_with_matrix (m);
+	n->mat.matrix = gel_matrixw_new_with_matrix_value_only_integer (m);
 	n->mat.quoted = FALSE;
 
 	return n;
@@ -560,7 +560,7 @@
 
 		GET_NEW_NODE (n);
 		n->type = MATRIX_NODE;
-		n->mat.matrix = gel_matrixw_new_with_matrix (m);
+		n->mat.matrix = gel_matrixw_new_with_matrix_value_only_real_nonrational (m);
 		n->mat.quoted = FALSE;
 
 		return n;
@@ -594,7 +594,7 @@
 
 		GET_NEW_NODE (n);
 		n->type = MATRIX_NODE;
-		n->mat.matrix = gel_matrixw_new_with_matrix (m);
+		n->mat.matrix = gel_matrixw_new_with_matrix_value_only_real_nonrational (m);
 		n->mat.quoted = FALSE;
 
 		return n;
@@ -666,7 +666,7 @@
 
 		GET_NEW_NODE (n);
 		n->type = MATRIX_NODE;
-		n->mat.matrix = gel_matrixw_new_with_matrix (m);
+		n->mat.matrix = gel_matrixw_new_with_matrix_value_only_integer (m);
 		n->mat.quoted = FALSE;
 
 		return n;
@@ -711,7 +711,7 @@
 
 		GET_NEW_NODE (n);
 		n->type = MATRIX_NODE;
-		n->mat.matrix = gel_matrixw_new_with_matrix (m);
+		n->mat.matrix = gel_matrixw_new_with_matrix_value_only_integer (m);
 		n->mat.quoted = FALSE;
 
 		return n;
@@ -3015,7 +3015,7 @@
 
 	GET_NEW_NODE (n);
 	n->type = MATRIX_NODE;
-	n->mat.matrix = gel_matrixw_new_with_matrix (nm);
+	n->mat.matrix = gel_matrixw_new_with_matrix_value_only_integer (nm);
 	if (a[0]->type == MATRIX_NODE)
 		n->mat.quoted = a[0]->mat.quoted;
 	else
@@ -3393,7 +3393,7 @@
 
 	GET_NEW_NODE (n);
 	n->type = MATRIX_NODE;
-	n->mat.matrix = gel_matrixw_new_with_matrix (nm);
+	n->mat.matrix = gel_matrixw_new_with_matrix_value_only_integer (nm);
 	n->mat.quoted = FALSE;
 
 	g_free (cols);
@@ -3508,7 +3508,7 @@
 
 	GET_NEW_NODE (n);
 	n->type = MATRIX_NODE;
-	n->mat.matrix = gel_matrixw_new_with_matrix (nm);
+	n->mat.matrix = gel_matrixw_new_with_matrix_value_only (nm);
 	n->mat.quoted = FALSE;
 
 	return n;
@@ -4264,7 +4264,7 @@
 	gel_matrixw_set_size (mn, size, 1);
 
 	rem = gel_matrixw_copy (pm);
-	gel_matrixw_make_private (rem);
+	gel_matrixw_make_private (rem, TRUE /* kill_type_caches */);
 
 	/* we know ql can't be zero */
 	ql = gel_matrixw_index (qm, sizeq-1, 0);
@@ -4487,7 +4487,7 @@
 
 	GET_NEW_NODE (n);
 	n->type = MATRIX_NODE;
-	n->mat.matrix = gel_matrixw_new_with_matrix (nm);
+	n->mat.matrix = gel_matrixw_new_with_matrix_value_only (nm);
 	n->mat.quoted = FALSE;
 
 	return n;
@@ -4916,7 +4916,7 @@
 
 	GET_NEW_NODE (n);
 	n->type = MATRIX_NODE;
-	n->mat.matrix = gel_matrixw_new_with_matrix (mm);
+	n->mat.matrix = gel_matrixw_new_with_matrix_value_only_integer (mm);
 	n->mat.quoted = FALSE;
 
 	return n;

Modified: trunk/src/genius.c
==============================================================================
--- trunk/src/genius.c	(original)
+++ trunk/src/genius.c	Tue May 20 05:25:01 2008
@@ -1,5 +1,5 @@
 /* GENIUS Calculator
- * Copyright (C) 1997-2007 Jiri (George) Lebl
+ * Copyright (C) 1997-2008 Jiri (George) Lebl
  *
  * Author: Jiri (George) Lebl
  *
@@ -440,6 +440,8 @@
 	  dictionary*/
 	d_singlecontext ();
 
+	gel_init ();
+
 	if ( ! (do_compile || do_gettext)) {
 		/*
 		 * Read main library

Modified: trunk/src/gnome-genius.c
==============================================================================
--- trunk/src/gnome-genius.c	(original)
+++ trunk/src/gnome-genius.c	Tue May 20 05:25:01 2008
@@ -1,5 +1,5 @@
 /* GENIUS Calculator
- * Copyright (C) 1997-2007 Jiri (George) Lebl
+ * Copyright (C) 1997-2008 Jiri (George) Lebl
  *
  * Author: Jiri (George) Lebl
  *
@@ -4085,6 +4085,8 @@
 	  dictionary*/
 	d_singlecontext ();
 
+	gel_init ();
+
 	gel_add_graph_functions ();
 
 	/*

Modified: trunk/src/matop.c
==============================================================================
--- trunk/src/matop.c	(original)
+++ trunk/src/matop.c	Tue May 20 05:25:01 2008
@@ -1,5 +1,5 @@
 /* GENIUS Calculator
- * Copyright (C) 1997-2007 Jiri (George) Lebl
+ * Copyright (C) 1997-2008 Jiri (George) Lebl
  *
  * Author: Jiri (George) Lebl
  *
@@ -201,7 +201,13 @@
 gel_matrix_conjugate_transpose (GelMatrixW *m)
 {
 	int i, j, w, h;
-	gel_matrixw_make_private (m);
+
+	if (gel_is_matrix_value_only_real (m)) {
+		m->tr = !m->tr;
+		return;
+	}
+
+	gel_matrixw_make_private (m, FALSE /* kill_type_caches */);
 	m->tr = !m->tr;
 	w = gel_matrixw_width (m);
 	h = gel_matrixw_height (m);
@@ -238,7 +244,7 @@
 	int i, j, k, w, h, m1w;
 	mpw_t tmp;
 	mpw_init(tmp);
-	gel_matrixw_make_private(res);
+	gel_matrixw_make_private(res, TRUE /* kill_type_caches */);
 
 	w = gel_matrixw_width (res);
 	h = gel_matrixw_height (res);
@@ -283,13 +289,15 @@
 	mpw_clear(tmp);
 }
 
+/* m must be made private before */
 static void
 swap_rows(GelMatrixW *m, int row, int row2)
 {
 	int i, w;
 	if(row==row2) return;
 	
-	gel_matrixw_make_private(m);
+	/* no need for this, only used within gauss and matrix is already private
+	gel_matrixw_make_private(m);*/
 	
 	w = gel_matrixw_width (m);
 	for (i = 0; i < w; i++) {
@@ -299,6 +307,7 @@
 	}
 }
 
+/* m must be made private before */
 static gboolean
 div_row (GelCtx *ctx, GelMatrixW *m, int row, mpw_t div)
 {
@@ -308,7 +317,8 @@
 	if (mpw_eql_ui (div, 1))
 		return TRUE;
 
-	gel_matrixw_make_private(m);
+	/* no need for this, only used within gauss and matrix is already private
+	gel_matrixw_make_private(m);*/
 	
 	w = gel_matrixw_width (m);
 	for (i = 0; i < w; i++) {
@@ -326,6 +336,7 @@
 	return ret;
 }
 
+/* m must be made private before */
 static gboolean
 mul_sub_row (GelCtx *ctx, GelMatrixW *m, int row, mpw_t mul, int row2)
 {
@@ -333,7 +344,8 @@
 	mpw_t tmp;
 	gboolean ret = TRUE;
 	
-	gel_matrixw_make_private(m);
+	/* no need for this, only used within gauss and matrix is already private
+	gel_matrixw_make_private(m);*/
 	
 	mpw_init(tmp);
 	w = gel_matrixw_width(m);
@@ -378,34 +390,84 @@
 	mpw_t tmp;
 	gboolean ret = TRUE;
 	gboolean made_private = FALSE;
+	gboolean matrix_rational = FALSE;
+	int *pivots = NULL;
+	int pivots_max = -1;
 	
 	if(detop)
 		mpw_set_ui(detop,1);
 
+	matrix_rational = gel_is_matrix_value_only_rational (m);
+
 	mpw_init(tmp);
 	d = 0;
 	w = gel_matrixw_width (m);
         h = gel_matrixw_height (m);
+
+	if (reduce) {
+		pivots = g_new0 (int, h);
+	}
+
 	for (i = 0; i < w && d < h; i++) {
-		for (j = d; j < h; j++) {
-			GelETree *t = gel_matrixw_get_index(m,i,j);
-			if (t != NULL &&
-			    ! mpw_zero_p (t->val.value))
-				break;
+		if (matrix_rational) {
+			for (j = d; j < h; j++) {
+				GelETree *t = gel_matrixw_get_index(m,i,j);
+				if (t != NULL &&
+				    ! mpw_zero_p (t->val.value))
+					break;
+			}
+		} else {
+			/* kind of a hack */
+			int bestj = h;
+			mpw_t best_abs_sq;
+			mpw_init (best_abs_sq);
+			GelETree *bestpiv = NULL;
+			for (j = d; j < h; j++) {
+				GelETree *t = gel_matrixw_get_index(m,i,j);
+				if (t != NULL &&
+				    ! mpw_zero_p (t->val.value)) {
+					if (bestpiv == NULL) {
+						bestpiv = t;
+						bestj = j;
+					} else {
+						mpw_abs_sq (tmp, t->val.value);
+						if (mpw_cmp (tmp, best_abs_sq) > 0) {
+							bestpiv = t;
+							bestj = j;
+						}
+					}
+				}
+			}
+			mpw_clear (best_abs_sq);
+
+			j = bestj;
 		}
 
 		if (j == h) {
 			ret = FALSE;
 			if(stopsing) {
 				mpw_clear(tmp);
+				g_free (pivots);
 				return FALSE;
 			}
 			continue;
 		}
 
 		if ( ! made_private) {
-			gel_matrixw_make_private(m);
+			gel_matrixw_make_private (m, TRUE /* kill_type_caches */);
+			if (simul)
+				gel_matrixw_make_private (simul, TRUE /* kill_type_caches */);
 			made_private = TRUE;
+
+			/* the matrix will be value only */
+			m->cached_value_only = 1;
+			m->value_only = 1;
+
+			if (matrix_rational) {
+				/* the matrix will be value only rational */
+				m->cached_value_only_rational = 1;
+				m->value_only_rational = 1;
+			}
 		}
 		
 		if (j > d) {
@@ -431,9 +493,10 @@
 					return FALSE;
 				}
 				if(simul) {
-					if ( !  mul_sub_row (ctx, simul, d, tmp, j) &&
+					if ( ! mul_sub_row (ctx, simul, d, tmp, j) &&
 					     stopsing) {
 						mpw_clear(tmp);
+						g_free (pivots);
 						return FALSE;
 					}
 				}
@@ -448,27 +511,11 @@
 		}
 
 		if(reduce) {
-			for(j=0;j<d;j++) {
-				GelETree *t = gel_matrixw_get_index(m,i,j);
-				if (t != NULL &&
-				    ! mpw_zero_p (t->val.value)) {
-					mpw_div(tmp,t->val.value,piv->val.value);
-					if ( ! mul_sub_row (ctx, m, d, tmp, j) &&
-					     stopsing) {
-						mpw_clear(tmp);
-						return FALSE;
-					}
-					if(simul) {
-						if ( !  mul_sub_row (ctx, simul, d, tmp, j) &&
-						     stopsing) {
-							mpw_clear(tmp);
-							return FALSE;
-						}
-					}
-				}
-			}
+			pivots[d] = i;
+			pivots_max = d;
 		}
 
+		/* make pivot 1 */
 		for (ii = i+1; ii < w; ii++) {
 			GelETree *t = gel_matrixw_get_index(m,ii,d);
 			if(t) {
@@ -481,6 +528,7 @@
 					    t != NULL &&
 					    t->type != VALUE_NODE) {
 						mpw_clear(tmp);
+						g_free (pivots);
 						return FALSE;
 					}
 				}
@@ -492,20 +540,50 @@
 			if ( ! div_row (ctx, simul, d, piv->val.value) &&
 			    stopsing) {
 				mpw_clear(tmp);
+				g_free (pivots);
 				return FALSE;
 			}
 		}
-
 		mpw_set_ui(piv->val.value,1);
 		d++;
 	}
 
+	if(reduce) {
+		for(d = pivots_max; d >= 0; d--) {
+			i = pivots[d];
+			for(j=0;j<d;j++) {
+				GelETree *t = gel_matrixw_get_index(m,i,j);
+				if (t != NULL &&
+				    ! mpw_zero_p (t->val.value)) {
+					/* subtle: must copy t->val.value,
+					 * else we wipe it out */
+					mpw_set (tmp, t->val.value);
+					if ( ! mul_sub_row (ctx, m, d, tmp, j) &&
+					     stopsing) {
+						mpw_clear(tmp);
+						g_free (pivots);
+						return FALSE;
+					}
+					if(simul) {
+						if ( ! mul_sub_row (ctx, simul, d, tmp, j) &&
+						     stopsing) {
+							mpw_clear(tmp);
+							g_free (pivots);
+							return FALSE;
+						}
+					}
+				}
+			}
+		}
+	}
+
 	if (detop && ctx->modulo != NULL) {
 		/* FIXME: this may fail!!! */
 		gel_mod_integer_rational (detop, ctx->modulo);
 	}
 	
 	mpw_clear(tmp);
+	g_free (pivots);
 	return ret;
 }
 

Modified: trunk/src/matrixw.c
==============================================================================
--- trunk/src/matrixw.c	(original)
+++ trunk/src/matrixw.c	Tue May 20 05:25:01 2008
@@ -1,5 +1,5 @@
 /* GENIUS Calculator
- * Copyright (C) 1997-2007 Jiri (George) Lebl
+ * Copyright (C) 1997-2008 Jiri (George) Lebl
  *
  * Author: George Lebl
  *
@@ -45,7 +45,8 @@
 	GelMatrixWFreeList *next;
 };
 
-static void internal_make_private (GelMatrixW *m, int w, int h);
+static void internal_make_private (GelMatrixW *m, int w, int h,
+				   gboolean kill_type_caches);
 
 static GelMatrixWFreeList *free_matrices = NULL;
 
@@ -142,9 +143,6 @@
 	m->regw = m->m->width;
 	m->regh = m->m->height;
 	
-	if (the_zero == NULL)
-		the_zero = gel_makenum_ui (0);
-	
 	return m;
 }
 GelMatrixW *
@@ -159,7 +157,6 @@
 	m->m = mat;
 	m->m->use++;
 
-	/* clear caches as we're gonna mess with this matrix */
 	m->cached_value_only = 0;
 	m->cached_value_only_real = 0;
 	m->cached_value_only_rational = 0;
@@ -173,8 +170,100 @@
 	m->regw = m->m->width;
 	m->regh = m->m->height;
 	
-	if (the_zero == NULL)
-		the_zero = gel_makenum_ui (0);
+	return m;
+}
+
+GelMatrixW *
+gel_matrixw_new_with_matrix_value_only (GelMatrix *mat)
+{
+	GelMatrixW *m;
+	NEW_MATRIX (m);
+#ifdef MATRIX_DEBUG
+	/*debug*/printf ("%s\n", G_GNUC_PRETTY_FUNCTION);
+#endif
+	
+	m->m = mat;
+	m->m->use++;
+
+	m->cached_value_only = 1;
+	m->value_only = 1;
+	m->cached_value_only_real = 0;
+	m->cached_value_only_rational = 0;
+	m->cached_value_only_integer = 0;
+	m->cached_value_or_bool_only = 0;
+	m->rref = 0;
+	
+	m->tr = 0;
+	m->regx = NULL;
+	m->regy = NULL;
+	m->regw = m->m->width;
+	m->regh = m->m->height;
+	
+	return m;
+}
+
+GelMatrixW *
+gel_matrixw_new_with_matrix_value_only_integer (GelMatrix *mat)
+{
+	GelMatrixW *m;
+	NEW_MATRIX (m);
+#ifdef MATRIX_DEBUG
+	/*debug*/printf ("%s\n", G_GNUC_PRETTY_FUNCTION);
+#endif
+	
+	m->m = mat;
+	m->m->use++;
+
+	m->cached_value_only = 1;
+	m->value_only = 1;
+	m->cached_value_only_real = 1;
+	m->value_only_real = 1;
+	m->cached_value_only_rational = 1;
+	m->value_only_rational = 1;
+	m->cached_value_only_integer = 1;
+	m->value_only_integer = 1;
+	m->cached_value_or_bool_only = 1;
+	m->value_or_bool_only = 1;
+	m->rref = 0;
+	
+	m->tr = 0;
+	m->regx = NULL;
+	m->regy = NULL;
+	m->regw = m->m->width;
+	m->regh = m->m->height;
+	
+	return m;
+}
+
+GelMatrixW *
+gel_matrixw_new_with_matrix_value_only_real_nonrational (GelMatrix *mat)
+{
+	GelMatrixW *m;
+	NEW_MATRIX (m);
+#ifdef MATRIX_DEBUG
+	/*debug*/printf ("%s\n", G_GNUC_PRETTY_FUNCTION);
+#endif
+	
+	m->m = mat;
+	m->m->use++;
+
+	m->cached_value_only = 1;
+	m->value_only = 1;
+	m->cached_value_only_real = 1;
+	m->value_only_real = 1;
+	m->cached_value_only_rational = 1;
+	m->value_only_rational = 0;
+	m->cached_value_only_integer = 1;
+	m->value_only_integer = 0;
+	m->cached_value_or_bool_only = 1;
+	m->value_or_bool_only = 1;
+	m->rref = 0;
+	
+	m->tr = 0;
+	m->regx = NULL;
+	m->regy = NULL;
+	m->regw = m->m->width;
+	m->regh = m->m->height;
 	
 	return m;
 }
@@ -436,7 +525,9 @@
 	} else if (m->m->use > 1) {
 		/* since the use is greater than 1, we WILL get a copy of
 		 * this matrix at the right size */
-		internal_make_private (m, width, height);
+		/* it may seem we could leave caches alone, but changing size
+		 * could make those values different */
+		internal_make_private (m, width, height, TRUE /* kill_type_caches */);
 		g_assert (m->m->use == 1);
 	} else /* m->m->use == 1 */{
 		ensure_at_least_size (m, width, height);
@@ -463,7 +554,7 @@
 
 	if (width > m->regw || height > m->regh) {
 		/* FIXME: this may be a bit inefficent */
-		gel_matrixw_make_private (m);
+		gel_matrixw_make_private (m, TRUE /* kill_type_caches */);
 		make_us_a_copy (m, MAX (width, m->regw),MAX (height, m->regh));
 		ensure_at_least_size (m, width, height);
 	}
@@ -483,9 +574,11 @@
 #endif
 
 	if (m->tr)
-		internal_make_private (m, MAX(m->regh, y+1), MAX(m->regw, x+1));
+		internal_make_private (m, MAX(m->regh, y+1), MAX(m->regw, x+1),
+				       TRUE /* kill_type_caches */);
 	else
-		internal_make_private (m, MAX(m->regw, x+1), MAX(m->regh, y+1));
+		internal_make_private (m, MAX(m->regw, x+1), MAX(m->regh, y+1),
+				       TRUE /* kill_type_caches */);
 	gel_matrixw_set_at_least_size (m, x+1, y+1);
 	t = gel_matrixw_get_index (m, x, y);
 	if (t != NULL)
@@ -511,7 +604,8 @@
 			x = i % m->regh;
 			y = i / m->regh;
 		}
-		internal_make_private (m, MAX(m->regh, y+1), MAX(m->regw, x+1));
+		internal_make_private (m, MAX(m->regh, y+1), MAX(m->regw, x+1),
+				       TRUE /* kill_type_caches */);
 	} else {
 		if (m->regh == 1) {
 			x = i;
@@ -520,7 +614,8 @@
 			x = i % m->regw;
 			y = i / m->regw;
 		}
-		internal_make_private (m, MAX(m->regw, x+1), MAX(m->regh, y+1));
+		internal_make_private (m, MAX(m->regw, x+1), MAX(m->regh, y+1),
+				       TRUE /* kill_type_caches */);
 	}
 	gel_matrixw_set_at_least_size (m, x+1, y+1);
 	t = gel_matrixw_get_index (m, x, y);
@@ -663,7 +758,8 @@
 
 	/* FIXME: we will copy some nodes we don't need to copy
 	 * as we will free them just below */
-	internal_make_private (m, MAX (xmax+1, m->regw), MAX (ymax+1, m->regh));
+	internal_make_private (m, MAX (xmax+1, m->regw), MAX (ymax+1, m->regh),
+			       TRUE /* kill_type_caches */);
 	ensure_at_least_size (m, xmax+1, ymax+1);
 	/* assume that's what ensure/make_us_a_copy does */
 	g_assert (m->regx == NULL && m->regy == NULL);
@@ -726,7 +822,8 @@
 	xmax = getmax (destx, w);
 	ymax = getmax (desty, h);
 
-	internal_make_private (m, MAX (xmax+1, m->regw), MAX (ymax+1, m->regh));
+	internal_make_private (m, MAX (xmax+1, m->regw), MAX (ymax+1, m->regh),
+			       TRUE /* kill_type_caches */);
 	ensure_at_least_size (m, xmax+1, ymax+1);
 	/* assume that's what ensure/make_us_a_copy does */
 	g_assert (m->regx == NULL && m->regy == NULL);
@@ -899,18 +996,21 @@
 
 /*make private copy of the GelMatrix*/
 static void
-internal_make_private (GelMatrixW *m, int w, int h)
+internal_make_private (GelMatrixW *m, int w, int h, gboolean kill_type_caches)
 {
 #ifdef MATRIX_DEBUG
 	/*debug*/printf ("%s\n", G_GNUC_PRETTY_FUNCTION);
 #endif
 
-	/* clear caches as we're gonna mess with this matrix */
-	m->cached_value_only = 0;
-	m->cached_value_only_real = 0;
-	m->cached_value_only_rational = 0;
-	m->cached_value_only_integer = 0;
-	m->cached_value_or_bool_only = 0;
+	if (kill_type_caches) {
+		/* clear caches as we're gonna mess with this matrix */
+		m->cached_value_only = 0;
+		m->cached_value_only_real = 0;
+		m->cached_value_only_rational = 0;
+		m->cached_value_only_integer = 0;
+		m->cached_value_or_bool_only = 0;
+	}
+
 	m->rref = 0;
 
 #ifdef MATRIX_DEBUG
@@ -934,14 +1034,14 @@
 
 /*make private copy of the GelMatrix*/
 void
-gel_matrixw_make_private (GelMatrixW *m)
+gel_matrixw_make_private (GelMatrixW *m, gboolean kill_type_caches)
 {
 	g_return_if_fail(m != NULL);
 #ifdef MATRIX_DEBUG
 	/*debug*/printf ("%s\n", G_GNUC_PRETTY_FUNCTION);
 #endif
 
-	internal_make_private (m, m->regw, m->regh);
+	internal_make_private (m, m->regw, m->regh, kill_type_caches);
 }
 
 /*free a matrix*/
@@ -995,11 +1095,13 @@
 	if (m->tr) {
 		int i;
 		if (m->regw == 1) {
-			internal_make_private (m, m->regw, MAX (max+1, m->regh));
+			internal_make_private (m, m->regw, MAX (max+1, m->regh),
+					       TRUE /* kill_type_caches */);
 			ensure_at_least_size (m, m->regw, max+1);
 		} else {
 			int minw = (max / m->regh) + 1;
-			internal_make_private (m, MAX (minw, m->regw), m->regh);
+			internal_make_private (m, MAX (minw, m->regw), m->regh,
+					       TRUE /* kill_type_caches */);
 			ensure_at_least_size (m, minw, m->regh);
 		}
 		/* assume that's what ensure/make_us_a_copy does */
@@ -1032,11 +1134,13 @@
 	} else {
 		int i;
 		if (m->regh == 1) {
-			internal_make_private (m, MAX (max+1, m->regw), m->regh);
+			internal_make_private (m, MAX (max+1, m->regw), m->regh,
+					       TRUE /* kill_type_caches */);
 			ensure_at_least_size (m, max+1, m->regh);
 		} else {
 			int minh = (max / m->regw) + 1;
-			internal_make_private (m, m->regw, MAX (minh, m->regh));
+			internal_make_private (m, m->regw, MAX (minh, m->regh),
+					       TRUE /* kill_type_caches */);
 			ensure_at_least_size (m, m->regw, minh);
 		}
 		/* assume that's what ensure/make_us_a_copy does */
@@ -1086,7 +1190,8 @@
 		int minw = (max / m->regh) + 1;
 		int i;
 
-		internal_make_private (m, MAX (minw, m->regw), m->regh);
+		internal_make_private (m, MAX (minw, m->regw), m->regh,
+				       TRUE /* kill_type_caches */);
 		ensure_at_least_size (m, minw, m->regh);
 		/* assume that's what ensure/make_us_a_copy does */
 		g_assert (m->regx == NULL && m->regy == NULL);
@@ -1104,7 +1209,8 @@
 		int minh = (max / m->regw) + 1;
 		int i;
 
-		internal_make_private (m, m->regw, MAX (minh, m->regh));
+		internal_make_private (m, m->regw, MAX (minh, m->regh),
+				       TRUE /* kill_type_caches */);
 		ensure_at_least_size (m, m->regw, minh);
 		/* assume that's what ensure/make_us_a_copy does */
 		g_assert (m->regx == NULL && m->regy == NULL);

Modified: trunk/src/matrixw.h
==============================================================================
--- trunk/src/matrixw.h	(original)
+++ trunk/src/matrixw.h	Tue May 20 05:25:01 2008
@@ -1,5 +1,5 @@
 /* GENIUS Calculator
- * Copyright (C) 1997-2007 Jiri (George) Lebl
+ * Copyright (C) 1997-2008 Jiri (George) Lebl
  *
  * Author: George Lebl
  *
@@ -54,6 +54,9 @@
 /*new matrix*/
 GelMatrixW *gel_matrixw_new(void);
 GelMatrixW *gel_matrixw_new_with_matrix(GelMatrix *mat);
+GelMatrixW *gel_matrixw_new_with_matrix_value_only (GelMatrix *mat);
+GelMatrixW *gel_matrixw_new_with_matrix_value_only_integer (GelMatrix *mat);
+GelMatrixW *gel_matrixw_new_with_matrix_value_only_real_nonrational (GelMatrix *mat);
 
 /*set size of a matrix*/
 void gel_matrixw_set_size(GelMatrixW *m, int width, int height);
@@ -93,7 +96,7 @@
 void gel_matrixw_free(GelMatrixW *m);
 
 /*make private copy of the GelMatrix*/
-void gel_matrixw_make_private(GelMatrixW *m);
+void gel_matrixw_make_private(GelMatrixW *m, gboolean kill_type_caches);
 
 extern GelETree *the_zero;
 

Modified: trunk/src/mpwrap.c
==============================================================================
--- trunk/src/mpwrap.c	(original)
+++ trunk/src/mpwrap.c	Tue May 20 05:25:01 2008
@@ -3380,7 +3380,6 @@
 	} else {
 		MpwRealNum t = {{NULL}};
 
-		MAKE_REAL(rop);
 		if (rop != op) {
 			MAKE_EMPTY(rop->r, op->r->type);
 		} else {
@@ -3396,6 +3395,40 @@
 		mpwl_free(&t,TRUE);
 
 		mpwl_sqrt(rop->r,rop->r);
+
+		MAKE_REAL (rop);
+	}
+}
+
+void
+mpw_abs_sq (mpw_ptr rop,mpw_ptr op)
+{
+	if (MPW_IS_REAL (op)) {
+		if(mpwl_sgn(op->r)<0)
+			mpw_neg(rop,op);
+		else
+			mpw_set(rop,op);
+
+		/* have to actually square now */
+		mpw_mul (rop, rop, rop);
+	} else {
+		MpwRealNum t = {{NULL}};
+
+		if (rop != op) {
+			MAKE_EMPTY(rop->r, op->r->type);
+		} else {
+			MAKE_COPY (rop->r);
+		}
+		
+		mpwl_init_type (&t, op->i->type);
+		
+		mpwl_mul(rop->r,op->r,op->r);
+		mpwl_mul(&t,op->i,op->i);
+		mpwl_add(rop->r,rop->r,&t);
+
+		mpwl_free(&t,TRUE);
+
+		MAKE_REAL (rop);
 	}
 }
 

Modified: trunk/src/mpwrap.h
==============================================================================
--- trunk/src/mpwrap.h	(original)
+++ trunk/src/mpwrap.h	Tue May 20 05:25:01 2008
@@ -151,6 +151,7 @@
 
 
 void mpw_abs(mpw_ptr rop,mpw_ptr op);
+void mpw_abs_sq(mpw_ptr rop,mpw_ptr op);
 
 int mpw_sgn(mpw_ptr op);
 

Modified: trunk/src/mpzextra.c
==============================================================================
--- trunk/src/mpzextra.c	(original)
+++ trunk/src/mpzextra.c	Tue May 20 05:25:01 2008
@@ -1,11 +1,13 @@
 /* GENIUS Calculator
- * Copyright (C) 1997-2003 George Lebl
+ * Copyright (C) 1997-2008 Jiri (George) Lebl
  *
- * Author: George Lebl
+ * Author: Jiri (George) Lebl
  *
- * This program is free software; you can redistribute it and/or modify
+ * This file is part of Genius.
+ *
+ * Genius is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 2 of the License, or
+ * the Free Software Foundation, either version 3 of the License, or
  * (at your option) any later version.
  *
  * This program is distributed in the hope that it will be useful,
@@ -14,9 +16,7 @@
  * GNU General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the  Free Software
- * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
- * USA.
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
 #include "config.h"
@@ -165,7 +165,7 @@
 	  if (evalnode_hook != NULL) {
 		  (*evalnode_hook)();
 	  }
-	  if (interrupted) {
+	  if G_UNLIKELY (interrupted) {
 		  is_prime = 0;
 		  break;
 	  }
@@ -243,7 +243,7 @@
 
 	if (evalnode_hook != NULL)
 		(*evalnode_hook)();
-	if (interrupted)
+	if G_UNLIKELY (interrupted)
 		return 0;
 
 	return mpz_millerrabin (n, miller_rabin_reps-1);
@@ -387,7 +387,7 @@
 		      i = 0;
 	      }
       }
-      if (interrupted) {
+      if G_UNLIKELY (interrupted) {
 	      mpz_set_ui (n, 1);
 	      continue;
       }
@@ -523,7 +523,7 @@
 
 	mpz_clear (n);
 
-	if (interrupted) {
+	if G_UNLIKELY (interrupted) {
 		mympz_factorization_free (fact);
 		return NULL;
 	}
[Date Prev][Date Next] [Thread Prev][Thread Next] [Thread Index] [Date Index] [Author Index]