Commit 4dea9c70 authored by Jiri (George) Lebl's avatar Jiri (George) Lebl

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
parent 1b252cfd
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
......
......@@ -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
......
......@@ -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);
......
......@@ -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);
......
......@@ -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
......
......@@ -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;
......
/* 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);
}
}
......
/* 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)) )
......
......@@ -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;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment