Bill Allombert on Sat, 20 Mar 2004 15:51:53 +0100


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

patch to speed up binomial


Hello PARI developers,

Here a patch to speed up binomial() the same way as mpfact:

Basically it use divide_conquer_prod. That should be around
60% faster on my tests.

This patch can probably be made cleaner.
I don't think it affects the accuracy of the result in any way
(when n is a real).

Cheers,
Bill.

Index: src/basemath/bibli2.c
===================================================================
RCS file: /home/cvs/pari/src/basemath/bibli2.c,v
retrieving revision 1.60
diff -u -r1.60 bibli2.c
--- src/basemath/bibli2.c	25 Jan 2004 17:50:05 -0000	1.60
+++ src/basemath/bibli2.c	20 Mar 2004 14:25:27 -0000
@@ -497,7 +497,7 @@
     if (k == 0) return gun;
     return gcopy(n);
   }
-  av = avma; y = n;
+  av = avma;
   if (typ(n) == t_INT)
   {
     if (signe(n) > 0)
@@ -512,14 +512,19 @@
         return gcopy(n);
       }
     }
-    for (i=2; i<=k; i++)
-      y = gdivgs(gmul(y,addis(n,i-1-k)), i);
+    y = cgetg(k+1,t_VEC);
+    for (i=1; i<=k; i++)
+      y[i] = lsubis(n,i-1);
+    y=divide_conquer_prod(y,mulii);
   }
   else
   {
-    for (i=2; i<=k; i++)
-      y = gdivgs(gmul(y,gaddgs(n,i-1-k)), i);
+    y = cgetg(k+1,t_VEC);
+    for (i=1; i<=k; i++)
+      y[i] = lsubgs(n,i-1);
+    y=divide_conquer_prod(y,gmul);
   }
+  y=gdiv(y,mpfact(k));
   return gerepileupto(av, y);
 }