[Nickle] nickle: Branch 'master'
Keith Packard
keithp at keithp.com
Mon Jul 5 21:30:26 PDT 2010
Makefile.am | 3 +
builtin-math.c | 10 ++++++
builtin.5c | 2 +
compile.c | 18 ++++++++++-
factorial.5c | 85 +++++++++++++++++++++++++++++++++++++++++++++++++++++
gram.y | 6 +++
math.5c | 3 +
nickle.h | 1
prime_sieve.5c | 90 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
9 files changed, 215 insertions(+), 3 deletions(-)
New commits:
commit b249d719eeb86e8834ccfa2c91c4ead07d497e09
Author: Keith Packard <keithp at keithp.com>
Date: Tue Jul 6 00:28:41 2010 -0400
Replace naïve factorial algorithm with prime sieve from Peter Luschny
diff --git a/Makefile.am b/Makefile.am
index 7611173..a66d6b6 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -12,7 +12,8 @@ BUILD_DATE := $(shell sh $(top_builddir)/date-sh)
NICKLEFILES = builtin.5c math.5c scanf.5c mutex.5c \
arc4.5c prng.5c command.5c abort.5c \
printf.5c history.5c ctype.5c string.5c socket.5c \
- file.5c parse-args.5c svg.5c process.5c
+ file.5c parse-args.5c svg.5c process.5c \
+ prime_sieve.5c factorial.5c
DEBIAN = debian/nickle.install debian/changelog debian/compat \
debian/control debian/copyright debian/rules debian/lintian.override
diff --git a/builtin-math.c b/builtin-math.c
index 5066474..6b2f61c 100644
--- a/builtin-math.c
+++ b/builtin-math.c
@@ -28,6 +28,10 @@ import_Math_namespace()
" int popcount (int i)\n"
"\n"
" Return the number of '1' bits in 'i'.\n" },
+ { do_Math_factorial, "factorial", "i", "i", "\n"
+ " int factorial (int i)\n"
+ "\n"
+ " Return the factorial of 'i'.\n" },
{ 0 }
};
static const struct fbuiltin_2 funcs_2[] = {
@@ -146,3 +150,9 @@ do_Math_popcount(Value v) {
ENTER ();
RETURN (Popcount (v));
}
+
+Value
+do_Math_factorial(Value v) {
+ ENTER ();
+ RETURN (Factorial (v));
+}
diff --git a/builtin.5c b/builtin.5c
index ae36829..af07fce 100644
--- a/builtin.5c
+++ b/builtin.5c
@@ -117,6 +117,8 @@ library "command.5c";
# Note that some libraries extend namespaces, and
# thus aren't autoload/autoimport candidates.
library "math.5c";
+library "prime_sieve.5c"
+library "factorial.5c"
import Math;
library "scanf.5c";
library "socket.5c";
diff --git a/compile.c b/compile.c
index 82d7add..c0c6c65 100644
--- a/compile.c
+++ b/compile.c
@@ -2528,7 +2528,23 @@ _CompileExpr (ObjPtr obj, ExprPtr expr, Bool evaluate, ExprPtr stat, CodePtr cod
case UMINUS: obj = CompileUnOp (obj, expr, NegateOp, stat, code); break;
case LNOT: obj = CompileUnFunc (obj, expr, Lnot, stat, code,"~"); break;
case BANG: obj = CompileUnFunc (obj, expr, Not, stat, code,"!"); break;
- case FACT: obj = CompileUnFunc (obj, expr, Factorial, stat, code,"!"); break;
+ case FACT:
+ obj = _CompileExpr (obj, expr->tree.left, True, stat, code);
+ SetPush (obj);
+ obj = _CompileExpr (obj, expr->tree.right, True, stat, code);
+ expr->base.type = TypeCombineUnary (expr->tree.right->base.type,
+ expr->base.tag);
+ if (!expr->base.type)
+ {
+ CompileError (obj, stat, "Incompatible type, value '%T', for ! operation",
+ expr->tree.right->base.type);
+ expr->base.type = typePoly;
+ break;
+ }
+ BuildInst (obj, OpCall, inst, stat);
+ inst->ints.value = 1;
+ BuildInst (obj, OpNoop, inst, stat);
+ break;
case INC:
if (expr->tree.left)
{
diff --git a/factorial.5c b/factorial.5c
new file mode 100644
index 0000000..7b5b8be
--- /dev/null
+++ b/factorial.5c
@@ -0,0 +1,85 @@
+/*
+ * Copyright © 2010 Keith Packard and Bart Massey.
+ * All Rights Reserved. See the file COPYING in this directory
+ * for licensing information.
+ */
+
+extend namespace Math {
+ /*
+ * PrimeSwing factorial algorithm by Peter Luschny
+ *
+ * Translated from the Java version as found here:
+ *
+ * http://www.luschny.de/math/factorial/FastFactorialFunctions.htm
+ *
+ * That code was released under the MIT license using this reference:
+ *
+ * http://www.opensource.org/licenses/mit-license.php
+ *
+ * This code is functionally equivalent, although is not a direct copy
+ */
+ public int factorial(int n)
+ /*
+ * Return the factorial of 'n'
+ */
+ {
+
+ /* For small values, just do it the naïve way */
+ if (n < 20) {
+ int r = 1;
+ while (n > 0) {
+ r *= n;
+ n--;
+ }
+ return r;
+ }
+
+ int plen = 2 * floor(sqrt(n) + n / (log2(n) - 1));
+ int[plen] prime_list;
+ int[] primes = PrimeSieve::primes(n);
+
+ static int[] smallOddSwing = {
+ 1, 1, 1, 3, 3, 15, 5, 35, 35, 315, 63, 693, 231, 3003, 429,
+ 6435, 6435, 109395, 12155, 230945, 46189, 969969, 88179,
+ 2028117, 676039, 16900975, 1300075, 35102025, 5014575,
+ 145422675, 9694845, 300540195, 300540195
+ };
+
+ /* Compute the swing factorial of 'n' */
+ int swing(int n) {
+ if (n < dim(smallOddSwing))
+ return smallOddSwing[n];
+ int sqrt_n = floor(sqrt(n));
+ PrimeSieve::range_t r1 = PrimeSieve::primes_range(&primes, 3, sqrt_n);
+ PrimeSieve::range_t r2 = PrimeSieve::primes_range(&primes, sqrt_n + 1, n // 3);
+ int count = 0;
+ for (int i_p = r1.start; i_p < r1.end; i_p++) {
+ int prime = primes[i_p];
+ int q = n, p = 1;
+ while ((q //= prime) > 0)
+ if ((q & 1) == 1)
+ p *= prime;
+ if (p > 1)
+ prime_list[count++] = p;
+ }
+ for (int i_p = r2.start; i_p < r2.end; i_p++) {
+ int prime = primes[i_p];
+ if (((n // prime) & 1) == 1)
+ prime_list[count++] = prime;
+ }
+ int prod = PrimeSieve::primorial(&primes, n // 2 + 1, n);
+ for (int i = 0; i < count; i++)
+ prod *= prime_list[i];
+ return prod;
+ }
+
+ /* Final factorial computation */
+ int recursive_factorial(int n) {
+ if (n < 2)
+ return 1;
+ return recursive_factorial(n>>1) ** 2 * swing(n);
+ }
+
+ return recursive_factorial(n) << (n - popcount(n));
+ }
+}
diff --git a/gram.y b/gram.y
index 4f358f7..d99bfc5 100644
--- a/gram.y
+++ b/gram.y
@@ -1305,7 +1305,11 @@ simpleexpr : simpleexpr assignop simpleexpr %prec ASSIGN
| BANG simpleexpr
{ $$ = NewExprTree(BANG, $2, (Expr *) 0); }
| simpleexpr BANG %prec FACT
- { $$ = NewExprTree(FACT, $1, (Expr *) 0); }
+ {
+ $$ = NewExprTree(FACT,
+ BuildName ("Math", "factorial"),
+ $1);
+ }
| INC simpleexpr
{ $$ = NewExprTree(INC, $2, (Expr *) 0); }
| simpleexpr INC
diff --git a/math.5c b/math.5c
index a6cd170..b25028e 100644
--- a/math.5c
+++ b/math.5c
@@ -648,8 +648,10 @@ extend namespace Math {
a4 = a**4;
iter = prec + 8;
v = 0;
+ int niter = 0;
while (iter-- > 0)
{
+ niter++;
term = ai/i! - aj/j!;
v += term;
if (exponent (v) > exponent (term) + iprec)
@@ -659,6 +661,7 @@ extend namespace Math {
i += 4;
j += 4;
}
+ printf ("cos: %d iters\n", niter);
return imprecise (v + term);
}
diff --git a/nickle.h b/nickle.h
index 610019d..25c3ade 100644
--- a/nickle.h
+++ b/nickle.h
@@ -940,6 +940,7 @@ Value do_xor (Value, Value);
Value do_Math_pow (Value, Value);
Value do_Math_assignpow (Value, Value);
Value do_Math_popcount (Value);
+Value do_Math_factorial (Value);
Value do_File_putb (Value, Value);
Value do_File_putc (Value, Value);
Value do_File_ungetb (Value, Value);
diff --git a/prime_sieve.5c b/prime_sieve.5c
new file mode 100644
index 0000000..751455b
--- /dev/null
+++ b/prime_sieve.5c
@@ -0,0 +1,90 @@
+/*
+ * Copyright © 2010 Keith Packard and Bart Massey.
+ * All Rights Reserved. See the file COPYING in this directory
+ * for licensing information.
+ */
+
+public namespace PrimeSieve {
+
+ int num_to_pos(int n) = (n - 3) // 2;
+ int pos_to_num(int p) = (p * 2) + 3;
+
+ /* Generate a list of primes from a boolean array of composite
+ * test results
+ */
+ int[] sieve_to_array(&bool[] composite, int n_prime) {
+ int[n_prime] primes;
+ int p = 0;
+ for (int i = 0; i < dim(composite); i++)
+ if (!composite[i])
+ primes[p++] = pos_to_num(i);
+ return primes;
+ }
+
+ /* Use the sieve of Eratosthenes to compute the set of
+ * primes <= max
+ */
+ public int[] primes (int max)
+ {
+ int n = num_to_pos(max+1);
+ bool[n] composite = { false ... };
+ int n_prime = 0;
+
+ for (int p = 0; p < n; p++) {
+ if (!composite[p]) {
+ int prime = pos_to_num(p);
+ int step = prime;
+ int init = p + prime;
+ for (int v = init; v < n; v += step)
+ composite[v] = true;
+ n_prime++;
+ }
+ }
+ return sieve_to_array(&composite, n_prime);
+ }
+
+ public typedef struct {
+ int start, end;
+ } range_t;
+
+ /* Binary search in a list of numbers for the value <= 'v' */
+ public int primes_search(&int[] primes, int v) {
+ int min = 0, max = dim(primes);
+ while (min < max) {
+ int mid = (max + min) >> 1;
+ int p = primes[mid];
+ if (p == v)
+ return mid;
+ if (p > v)
+ max = mid;
+ else {
+ if (min == mid)
+ return mid;
+ min = mid;
+ }
+ }
+ return min;
+ }
+
+ /* Return the indices within a list of primes such that
+ * primes[start] <= min && max < primes[end]
+ */
+ public range_t primes_range(&int[] primes, int min, int max) {
+ int start = primes_search(&primes, min);
+ if (primes[start] < min) start++;
+ int end = primes_search(&primes, max) + 1;
+ return (range_t) { start = start, end = end };
+ }
+
+ /*
+ * Multiply all primes between min and max within the
+ * supplied list.
+ */
+ public int primorial(&int[] primes, int min, int max) {
+ range_t r = primes_range(&primes, min, max);
+ int v = 1;
+ for (int i = r.start; i < r.end; i++)
+ v *= primes[i];
+ return v;
+ }
+}
More information about the Nickle
mailing list