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