GNOME.org
 

[genius] Mon Dec 17 13:45:21 2012 Jiri (George) Lebl <jirka 5z com>

commit 4dea9c7044955d4f23b5b3072e5588cc86cd74ec
Author: Jiri (George) Lebl <jirka 5z com>
Date:   Mon Dec 17 13:45:31 2012 -0600

    Mon Dec 17 13:45:21 2012  Jiri (George) Lebl <jirka 5z com>
    
    	* src/eval.c: when forloop (sum/prod/for) is with floats then check
    	  for overrun, and if within 2^-20 times the step size, then still
    	  execute the last body (assume things were because of roundoff
    	  error), this makes things like for n=0 to 1 by 0.1 actually work
    	  as expected.
    
    	* src/mpwrap.[ch], src/eval.c, src/funclib.c, src/matop.c: minor
    	  optimization in make_float and new macro to check real part for
    	  being floating point
    
    	* src/geniustests.txt, src/testfourier.gel: update testsuite

 ChangeLog           |   14 ++++++++++
 NEWS                |    3 ++
 src/eval.c          |   68 ++++++++++++++++++++++++++++++++++++++++++++------
 src/funclib.c       |    2 +-
 src/geniustests.txt |    8 ++++++
 src/matop.c         |    2 +-
 src/mpwrap.c        |   10 +++---
 src/mpwrap.h        |    4 ++-
 src/testfourier.gel |    4 +-
 9 files changed, 96 insertions(+), 19 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 1749356..620cd37 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,17 @@
+Mon Dec 17 13:45:21 2012  Jiri (George) Lebl <jirka 5z com>
+
+	* src/eval.c: when forloop (sum/prod/for) is with floats then check
+	  for overrun, and if within 2^-20 times the step size, then still
+	  execute the last body (assume things were because of roundoff
+	  error), this makes things like for n=0 to 1 by 0.1 actually work
+	  as expected.
+
+	* src/mpwrap.[ch], src/eval.c, src/funclib.c, src/matop.c: minor
+	  optimization in make_float and new macro to check real part for
+	  being floating point
+
+	* src/geniustests.txt, src/testfourier.gel: update testsuite
+
 Tue Dec 11 00:10:02 2012  Jiri (George) Lebl <jirka 5z com>
 
 	* configure.in, src/Makefile.am: exorcise last vestiges of the
diff --git a/NEWS b/NEWS
index a9d39b0..541f29b 100644
--- a/NEWS
+++ b/NEWS
@@ -17,6 +17,9 @@ Changes to 1.0.16
   color is more useful when zoomed out.
 * Simpler output when typing "help foo" when foo is neither defined nor
   documented.
+* When for/sum/prod loops are in terms of floating point numbers and
+  the final number is within 2^-20 times the step size of the goal,
+  assume there were roundoff errors and still execute the body
 * Handle wider matrices than 2^15 columns in expansion
 * Fix flicker when plotting surfaces to allow animations with 3d plots
 * Fix possible uninitialized crash when reading badly formed standard library
diff --git a/src/eval.c b/src/eval.c
index 6a3d1a3..81384ec 100644
--- a/src/eval.c
+++ b/src/eval.c
@@ -4099,12 +4099,47 @@ iter_pop_stack(GelCtx *ctx)
 		case GE_FOR:
 			{
 				GelEvalFor *evf = data;
-				if(evf->by)
-					mpw_add(evf->x,evf->x,evf->by);
+				gboolean done = FALSE;
+				if (evf->by)
+					mpw_add (evf->x, evf->x, evf->by);
 				else
-					mpw_add_ui(evf->x,evf->x,1);
-				/*if done*/
-				if(mpw_cmp(evf->x,evf->to) == -evf->init_cmp) {
+					mpw_add_ui (evf->x, evf->x, 1);
+				/* we know we aren't dealing with complexes */
+				if (mpw_is_real_part_float (evf->x)) {
+					int thecmp = mpw_cmp (evf->x, evf->to);
+					if (mpw_cmp (evf->x, evf->to) == -evf->init_cmp) {
+						/* maybe we just missed it, let's look back within 2^-20 of the by and see */
+						mpw_t tmp;
+						if (evf->by != NULL) {
+							mpfr_ptr f;
+							/* by is definitely mpfr */
+							mpw_init_set (tmp, evf->by);
+							f = mpw_peek_real_mpf (tmp);
+							mpfr_mul_2si (f, f, -20, GMP_RNDN);
+						} else {
+							mpw_init (tmp);
+							mpw_set_d (tmp, 1.0/1048576.0 /* 2^-20 */);
+						}
+
+						mpw_sub (tmp, evf->x, tmp);
+
+						done = (mpw_cmp(tmp,evf->to) == -evf->init_cmp);
+
+						/* don't use x, but use the to, x might be too far */
+						if ( ! done) {
+							mpw_set (evf->x, evf->to);
+						}
+
+						mpw_clear (tmp);
+					} else {
+						done = FALSE;
+					}
+				} else {
+					/*if done*/
+					done = (mpw_cmp(evf->x,evf->to) == -evf->init_cmp);
+				}
+
+				if (done) {
 					GelETree *res;
 					GE_POP_STACK(ctx,data,flag);
 					g_assert ((flag & GE_MASK) == GE_POST);
@@ -4990,8 +5025,8 @@ iter_forloop (GelCtx *ctx, GelETree *n, gboolean *repushed)
 	
 	init_cmp = mpw_cmp(from->val.value,to->val.value);
 	
-	/*if no iterations*/
 	if(!by) {
+		/*if no iterations*/
 		if(init_cmp>0) {
 			d_addfunc(d_makevfunc(ident->id.id,gel_copynode(from)));
 			freetree_full(n,TRUE,FALSE);
@@ -5007,10 +5042,17 @@ iter_forloop (GelCtx *ctx, GelETree *n, gboolean *repushed)
 		} else if(init_cmp==0) {
 			init_cmp = -1;
 		}
+		if (mpw_is_real_part_float (from->val.value) ||
+		    mpw_is_real_part_float (to->val.value)) {
+			/* ensure all float */
+			mpw_make_float (to->val.value);
+			mpw_make_float (from->val.value);
+		}
 		evf = evf_new(type, from->val.value,to->val.value,NULL,init_cmp,
 			      gel_copynode(body),body,ident->id.id);
 	} else {
 		int sgn = mpw_sgn(by->val.value);
+		/*if no iterations*/
 		if((sgn>0 && init_cmp>0) || (sgn<0 && init_cmp<0)) {
 			d_addfunc(d_makevfunc(ident->id.id,gel_copynode(from)));
 			freetree_full(n,TRUE,FALSE);
@@ -5026,6 +5068,14 @@ iter_forloop (GelCtx *ctx, GelETree *n, gboolean *repushed)
 		}
 		if(init_cmp == 0)
 			init_cmp = -sgn;
+		if (mpw_is_real_part_float (from->val.value) ||
+		    mpw_is_real_part_float (to->val.value) ||
+		    mpw_is_real_part_float (by->val.value)) {
+			/* ensure all float */
+			mpw_make_float (to->val.value);
+			mpw_make_float (from->val.value);
+			mpw_make_float (by->val.value);
+		}
 		evf = evf_new(type, from->val.value,to->val.value,by->val.value,
 			      init_cmp,gel_copynode(body),body,ident->id.id);
 	}
@@ -8213,7 +8263,7 @@ op_precalc_1 (GelETree *n,
 	if (l->type != GEL_VALUE_NODE ||
 	    (respect_type &&
 	     (mpw_is_complex (l->val.value) ||
-	      mpw_is_float (l->val.value))))
+	      mpw_is_real_part_float (l->val.value))))
 		return;
 	mpw_init(res);
 	(*func)(res,l->val.value);
@@ -8239,8 +8289,8 @@ op_precalc_2 (GelETree *n,
 	    (respect_type &&
 	     (mpw_is_complex (l->val.value) ||
 	      mpw_is_complex (r->val.value) ||
-	      mpw_is_float (l->val.value) ||
-	      mpw_is_float (r->val.value))))
+	      mpw_is_real_part_float (l->val.value) ||
+	      mpw_is_real_part_float (r->val.value))))
 		return;
 	mpw_init(res);
 	(*func)(res,l->val.value,r->val.value);
diff --git a/src/funclib.c b/src/funclib.c
index 845ff7a..6bfc0f9 100644
--- a/src/funclib.c
+++ b/src/funclib.c
@@ -2006,7 +2006,7 @@ IsFloat_op(GelCtx *ctx, GelETree * * a, gboolean *exception)
 	if(a[0]->type!=GEL_VALUE_NODE ||
 	   mpw_is_complex(a[0]->val.value))
 		return gel_makenum_bool (0);
-	else if(mpw_is_float(a[0]->val.value))
+	else if(mpw_is_real_part_float(a[0]->val.value))
 		return gel_makenum_bool (1);
 	else
 		return gel_makenum_bool (0);
diff --git a/src/geniustests.txt b/src/geniustests.txt
index 6fe0e47..ec00e7f 100644
--- a/src/geniustests.txt
+++ b/src/geniustests.txt
@@ -1112,6 +1112,14 @@ sinc(5)==sin(5)/5						true
 A=[1;2];B=[3;4;5;6;7];[A,B,0,null,4]+""				"[1,3,0,4;2,4,0,4;1,5,0,4;2,6,0,4;1,7,0,4]"
 A=[1,2;3,4];B=[5;6;7];[A,B]+""					"[1,2,5;3,4,6;1,2,7]"
 [[1;[1,2]],[3,[7;4]]]+""					"[1,1,3,7;1,2,3,4]"
+sum n=0 to 1 by 0.100000001 do n				5.500000045
+sum n=0 to 0.95 by 0.1 do n					4.5
+sum n=0 to 0.9 by 0.1 do n					4.5
+sum n=0 to 1 by 0.1 do n					5.5
+for n=0 to 1 by 0.1 do n					1.0
+for n=0 to 1 by 0.100000001 do n				1.0
+for n=0 to 1 by 0.10000001 do n					0.90000009
+for n=1 to 9.99999999999 do n					9.99999999999
 load "nullspacetest.gel"					true
 load "longtest.gel"						true
 load "testprec.gel"						true
diff --git a/src/matop.c b/src/matop.c
index 5d99471..f33a702 100644
--- a/src/matop.c
+++ b/src/matop.c
@@ -151,7 +151,7 @@ gel_is_matrix_value_only_rational (GelMatrixW *m)
 			if (n != NULL &&
 			    (n->type != GEL_VALUE_NODE ||
 			     mpw_is_complex (n->val.value) ||
-			     mpw_is_float (n->val.value))) {
+			     mpw_is_real_part_float (n->val.value))) {
 				m->cached_value_only_rational = 1;
 				m->value_only_rational = 0;
 				return FALSE;
diff --git a/src/mpwrap.c b/src/mpwrap.c
index c1c5563..8e234e4 100644
--- a/src/mpwrap.c
+++ b/src/mpwrap.c
@@ -1,5 +1,5 @@
 /* GENIUS Calculator
- * Copyright (C) 1997-2011 Jiri (George) Lebl
+ * Copyright (C) 1997-2012 Jiri (George) Lebl
  *
  * Author: Jiri (George) Lebl
  *
@@ -5120,13 +5120,13 @@ mpw_make_int(mpw_ptr rop)
 void
 mpw_make_float(mpw_ptr rop)
 {
-	if (MPW_IS_REAL (rop)) {
+	if (rop->r->type != MPW_FLOAT) {
 		MAKE_COPY(rop->r);
 		mpwl_make_float(rop->r);
-	} else {
-		MAKE_COPY(rop->r);
+	}
+	if ( ! MPW_IS_REAL (rop) &&
+	    rop->i->type != MPW_FLOAT) {
 		MAKE_COPY(rop->i);
-		mpwl_make_float(rop->r);
 		mpwl_make_float(rop->i);
 	}
 }
diff --git a/src/mpwrap.h b/src/mpwrap.h
index 9ebbf09..7afb57b 100644
--- a/src/mpwrap.h
+++ b/src/mpwrap.h
@@ -1,5 +1,5 @@
 /* GENIUS Calculator
- * Copyright (C) 1997-2008 Jiri (George) Lebl
+ * Copyright (C) 1997-2012 Jiri (George) Lebl
  *
  * Author: Jiri (George) Lebl
  *
@@ -303,6 +303,8 @@ gboolean mpw_is_integer(mpw_ptr op);
 gboolean mpw_is_rational(mpw_ptr op);
 gboolean mpw_is_float(mpw_ptr op);
 
+#define mpw_is_real_part_float(op) ((op)->r->type == MPW_FLOAT)
+
 #define mpw_is_complex_float(op) \
 	 ( ((op)->r->type == MPW_FLOAT) || \
 	   (MPW_IS_COMPLEX (op) && ((op)->i->type == MPW_FLOAT)) )
diff --git a/src/testfourier.gel b/src/testfourier.gel
index 2dc567b..2654c65 100644
--- a/src/testfourier.gel
+++ b/src/testfourier.gel
@@ -14,10 +14,10 @@ for x=-0.7 to 0.7 by 0.1 do (
 function ff(x) = x^3+x^2;
 f = NumericalFourierSineSeriesFunction(ff,1,30);
 for x=0.3 to 0.7 by 0.1 do (
-	if |f(x)-ff(x)| >= 0.03 then (error("Fourier test 7 fail at x=" + x);exit());
+	if |f(x)-ff(x)| >= x/15 then (error("Fourier test 7 fail at x=" + x);exit());
 );
 for x=0.3 to 0.7 by 0.1 do (
-	if |-f(-x)-ff(x)| >= 0.03 then (error("Fourier test 8 fail at x=" + x);exit());
+	if |-f(-x)-ff(x)| >= x/15 then (error("Fourier test 8 fail at x=" + x);exit());
 );
 
 function ff(x) = x^3+x^2;
[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]