| Bill Allombert on Fri, 27 Jan 2012 21:49:21 +0100 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
| for and sum |
Hello PARI developers, I found out why for() is slower than sum(): (20:54) gp > for(i=1,10^8,0); (20:54) gp > ## *** last result computed in 7,965 ms. (20:54) gp > sum(i=1,10^8,0); (20:55) gp > ## *** last result computed in 6,584 ms. This is due to for() allowing the lower bound not to be a t_INT. So of course I made a branch bill-forparii which adds a function forparii that handle the case where a is an integer faster. (PAtch in attachement, made with git format-patch). (21:33) gp > for(i=1,10^8,0); (21:33) gp > ## *** last result computed in 5,981 ms. Cheers, Bill.
>From 668b9477d4903054f2e1fcc330199c52e0a78a6c Mon Sep 17 00:00:00 2001
From: Bill Allombert <Bill.Allombert@math.u-bordeaux1.fr>
Date: Fri, 27 Jan 2012 21:45:25 +0100
Subject: [PATCH] Add forparii when a is a integer.
---
src/language/sumiter.c | 29 ++++++++++++++++++++++++++++-
1 files changed, 28 insertions(+), 1 deletions(-)
diff --git a/src/language/sumiter.c b/src/language/sumiter.c
index 4cc80d5..e52164c 100644
--- a/src/language/sumiter.c
+++ b/src/language/sumiter.c
@@ -59,17 +59,44 @@ iferrnamepari(const char *err, GEN a, GEN b)
/** **/
/********************************************************************/
+static void
+forparii(GEN a, GEN b, GEN code)
+{
+ pari_sp av, av0 = avma, lim;
+ GEN aa;
+ if (gcmp(b,a) < 0) return;
+ b = gfloor(b);
+ aa = a = setloop(a);
+ av=avma; lim = stack_lim(av,1);
+ push_lex(a,code);
+ for(;;)
+ {
+ closure_evalvoid(code); if (loop_break()) break;
+ if (cmpii(a,b) >= 0) break;
+ a = get_lex(-1);
+ a = a==aa? incloop(a): gaddgs(a,1);
+ if (low_stack(lim, stack_lim(av,1)))
+ {
+ if (DEBUGMEM>1) pari_warn(warnmem,"forparii");
+ a = gerepileupto(av,a);
+ }
+ set_lex(-1,a);
+ }
+ pop_lex(1); avma = av0;
+}
+
void
forpari(GEN a, GEN b, GEN code)
{
pari_sp ltop=avma, av, lim;
+ if (typ(a) == t_INT) { forparii(a,b,code); return; }
b = gcopy(b); /* Kludge to work-around the a+(a=2) bug */
av=avma; lim = stack_lim(av,1);
push_lex(a,code);
while (gcmp(a,b) <= 0)
{
closure_evalvoid(code); if (loop_break()) break;
- a = get_lex(-1); a = typ(a) == t_INT? addis(a, 1): gaddgs(a,1);
+ a = get_lex(-1); a = gaddgs(a,1);
if (low_stack(lim, stack_lim(av,1)))
{
if (DEBUGMEM>1) pari_warn(warnmem,"forpari");
--
1.7.2.5