Overhaul `math' representation machinery.
authorMark Wooding <mdw@distorted.org.uk>
Wed, 5 Jun 2013 16:14:30 +0000 (17:14 +0100)
committerMark Wooding <mdw@distorted.org.uk>
Tue, 25 Jun 2013 23:37:42 +0000 (00:37 +0100)
Collect type information from the C compiler at configuration
time (using a rather complicated hack so that it works with
cross-compilers).  Read this from a Python script `mpgen' which is now
responsible for knowing all of the `mp' representation details.

Since `mpgen' generates all of the constant tables directly, we no
longer have any need for the programs `genlimits' or `mpdump' -- or the
random collection of `awk' scripts for turning `mumbletab.in' files into
`mumbletab.c' files.  And this means that we can kill `libmpbase.la'.

With this change, Catacomb is finally safe for cross-compilation.

Makefile.am
configure.ac
m4/.gitignore [new file with mode: 0644]
m4/mdw-probe-constant.m4 [new file with mode: 0644]
m4/mdw-uint-bits.m4 [new file with mode: 0644]
math/Makefile.am
math/mpgen [new file with mode: 0644]
math/typeinfo.py.in [new file with mode: 0644]

index 4cacfbe..e29c36e 100644 (file)
@@ -25,6 +25,7 @@
 ### MA 02111-1307, USA.
 
 include $(top_srcdir)/vars.am
+ACLOCAL_AMFLAGS                 = -Im4
 
 SUBDIRS                         =
 
index aca24fe..7c66807 100644 (file)
@@ -61,6 +61,70 @@ AC_SEARCH_LIBS([sqrt], [m])
 AC_SUBST([CATACOMB_LIBS], [$LIBS])
 LIBS=$mdw_ORIG_LIBS
 
+dnl Find out whether very long integer types are available.
+AC_CHECK_HEADERS([stdint.h])
+AC_SUBST([have_stdint_h])
+AC_C_LONG_LONG
+
+dnl Find the bit lengths of the obvious integer types.  This will be useful
+dnl when deciding on a representation for multiprecision integers.
+type_bits="" type_bits_sep=""
+AC_DEFUN([catacomb_UINT_BITS],
+  [mdw_UINT_BITS([$2], [$1])
+   type_bits="$type_bits$type_bits_sep('$1', $[]$1_bits)"
+   type_bits_sep=", "])
+catacomb_UINT_BITS([uchar], [unsigned char])
+catacomb_UINT_BITS([ushort], [unsigned short])
+catacomb_UINT_BITS([uint], [unsigned int])
+catacomb_UINT_BITS([ulong], [unsigned long])
+if test "$ac_cv_c_long_long" = "yes"; then
+  catacomb_UINT_BITS([ullong], [unsigned long long])
+fi
+if test "$ac_cv_header_stdint_h" = "yes"; then
+  catacomb_UINT_BITS([uintmax], [uintmax_t])
+fi
+AC_SUBST([type_bits])
+
+dnl Determine the limits of common C integer types.
+limits="" limits_sep=""
+AC_DEFUN([catacomb_COMPILE_TIME_CONSTANT],
+  [case "$2" in
+     =*)
+       $1="$2"; $1=${$1#=}
+       ;;
+     *)
+       AC_CACHE_CHECK([compile-time value of $2], [mdw_cv_constant_$3],
+                     [mdw_PROBE_CONSTANT([mdw_cv_constant_$3], [$2], [$4])])
+       $1=$mdw_cv_constant_$3
+       ;;
+   esac])
+AC_DEFUN([catacomb_LIMIT],
+[catacomb_COMPILE_TIME_CONSTANT([lo], [$2], [$1_min],
+[#include <limits.h>
+#include <stddef.h>])
+ catacomb_COMPILE_TIME_CONSTANT([hi], [$3], [$1_max],
+[#include <limits.h>
+#include <stddef.h>])
+ limits="$limits$limits_sep('$1', $lo, $hi)" limits_sep=", "])
+catacomb_LIMIT([SCHAR],                [SCHAR_MIN],    [SCHAR_MAX])
+catacomb_LIMIT([CHAR],         [CHAR_MIN],     [CHAR_MAX])
+catacomb_LIMIT([UCHAR],                [=0],           [UCHAR_MAX])
+catacomb_LIMIT([UINT8],                [=0],           [=0xff])
+catacomb_LIMIT([SHRT],         [SHRT_MIN],     [SHRT_MAX])
+catacomb_LIMIT([USHRT],                [=0],           [USHRT_MAX])
+catacomb_LIMIT([UINT16],       [=0],           [=0xffff])
+catacomb_LIMIT([INT],          [INT_MIN],      [INT_MAX])
+catacomb_LIMIT([UINT],         [=0],           [UINT_MAX])
+catacomb_LIMIT([LONG],         [LONG_MIN],     [LONG_MAX])
+catacomb_LIMIT([ULONG],                [=0],           [ULONG_MAX])
+catacomb_LIMIT([UINT32],       [=0],           [=0xffffffff])
+if test "$ac_cv_c_long_long" = "yes"; then
+  catacomb_LIMIT([LLONG],      [LLONG_MIN],    [LLONG_MAX])
+  catacomb_LIMIT([ULLONG],     [=0],           [ULLONG_MAX])
+fi
+catacomb_LIMIT([SIZET],                [=0],           [~(size_t)0])
+AC_SUBST([limits])
+
 dnl Functions used for noise-gathering.
 AC_CHECK_FUNCS([setgroups])
 AC_CACHE_CHECK([whether the freewheel noise generator will work],
diff --git a/m4/.gitignore b/m4/.gitignore
new file mode 100644 (file)
index 0000000..94f2b51
--- /dev/null
@@ -0,0 +1,2 @@
+libtool.m4
+lt*.m4
diff --git a/m4/mdw-probe-constant.m4 b/m4/mdw-probe-constant.m4
new file mode 100644 (file)
index 0000000..3916582
--- /dev/null
@@ -0,0 +1,183 @@
+dnl -*-autoconf-*-
+
+### SYNOPSIS
+###
+###   mdw_PROBE_CONSTANT(VAR, EXPR, [PREAMBLE], [IF-FAILED])
+###
+### DESCRIPTION
+###
+###   Extracts the value of a a constant integer expression from the
+###   compiler.  This works even if the compiler in question doesn't target
+###   the current architecture.  The value must be in the range -10^244 < x <
+###   10^244; this is probably fair enough.  In the extraordinarily unliklely
+###   event that the constant value is outside these bounds, the macro will
+###   fail rather than silently giving a wrong answer.
+###
+###   The result of the macro is that the shell variable VAR has the value of
+###   the expression EXPR, in decimal.  The PREAMBLE, if given, is inserted
+###   before EXPR is evaluated; it should contain #include and #define
+###   directives which are used to compute the value of the expression.
+###
+###   The idea for this macro came from the AC_C_COMPILE_VALUE macro by
+###   Ilguiz Latypov; this implementation has a number of advantages:
+###
+###     * it has an immense range of representable values, notably including
+###       negative numbers; and
+###
+###     * it returns the value directly in a shell variable rather than
+###       inventing an AC_DEFINE for it.
+###
+### LICENSE
+###
+###   Copyright (c) 2013 Mark Wooding <mdw@distorted.org.uk>
+###
+###   This program is free software: you can redistribute it and/or modify it
+###   under the terms of the GNU General Public License as published by the
+###   Free Software Foundation, either version 2 of the License, or (at your
+###   option) any later version.
+###
+###   This program is distributed in the hope that it will be useful, but
+###   WITHOUT ANY WARRANTY; without even the implied warranty of
+###   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+###   General Public License for more details.
+###
+###   You should have received a copy of the GNU General Public License along
+###   with this program. If not, see <http://www.gnu.org/licenses/>.
+###
+###   In particular, no exception to the GPL is granted regarding generated
+###   `configure' scripts which are the output of Autoconf.
+
+# Serial 1
+AC_COPYRIGHT([
+Portions copyright (c) 2013 Mark Wooding.
+
+This configure script is free software: you can redistribute it and/or
+modify it under he terms of the GNU General Public License as published
+by the Free Software Foundation, either version 2 of the License, or
+(at your option) any later version.])
+
+AC_DEFUN([mdw__PROBE_CONSTANT_SETUP],
+  [mdw__probe_constant_body="[
+
+/* The following program is copyright (c) 2013 Mark Wooding.  It is free
+ * software: you can redistribute it and/or modify it under the terms of the
+ * GNU General Public License as published by the Free Software Foundation,
+ * either version 2 of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ * for more details.
+ *
+ * You should have received a copy of the GNU General Public License along
+ * with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+/* The constant: 1 billion.  We'll pull digits out in groups of nine, since
+ * we can work with constants of at least the size of a C \`long'.
+ */
+#define MDW__G 1000000000
+
+/* An empty macro, used as an argument sometimes. */
+#define MDW__E
+
+/* A cheesy compile-time assertion.  If X is zero, then we try to declare
+ * an array with a negative number of elements.  Wrap this up in an anonymous
+ * struct so that we don't have to worry about naming things if we make
+ * more than one assertion.
+ */
+#define MDW__ASSERT(x) struct { int v[1 - 2*!(x)]; }
+
+/* Return the value of X/DIV, with further divisions D applied, truncating
+ * towards zero.  DIV must be greater than one.  This works even if X is
+ * negative, never tries to divide negative numbers, and doesn't try to
+ * negate the most-negative value.  There are three cases: if X <= -DIV then
+ * X/DIV = -(X + DIV)/DIV - 1, and X + DIV is less negative than X so this is
+ * a safe negation; if -DIV < X < 0 then the result is zero; otherwise, X
+ * is nonnegative so the straightforward division is safe.  Because DIV > 1,
+ * X/DIV is safe to negate, and we can apply the remaining divisions to it.
+ */
+#define MDW__SHIFT(x, div, d) \\
+  ((x) >= 0 ? ((x)/div d) : \\
+   (x) <= -(div) ? -((-((x) + (div))/(div) + 1) d) : 0)
+
+/* Extract the bottommost digit of X, as an integer: i.e., the value of
+ * abs(X) mod 10.  This works even if X is negative, never tries to divide
+ * negative numbers, and doesn't try to divide the most-negative value.
+ */
+#define MDW__RAW_DIGIT(x) (((x) < 0 ? -((x) + 1) % 10 + 1 : (x)) % 10)
+
+/* Extract the bottommost digit of X, as a character; if X is zero, then
+ * produce a space instead.  This avoids leading zeroes which can be
+ * misinterpreted by callers.
+ */
+#define MDW__TEXT_DIGIT(x) ((x) ? '0' + MDW__RAW_DIGIT(x) : ' ')
+
+/* Extract the bottommost digit of the probe value, after dividing by DIV
+ * and applying the divisons D.
+ */
+#define MDW__DIGIT(div, d) \\
+  MDW__TEXT_DIGIT(MDW__SHIFT(MDW__PROBE_EXPR, div, d))
+
+/* Extract the bottommost six digits of the probe value after dividing by 10
+ * and then applying the divisions D.
+ */
+#define MDW__9DIGITS(d) \\
+  MDW__DIGIT(1000000000, d), \\
+  MDW__DIGIT( 100000000, d), \\
+  MDW__DIGIT(  10000000, d), \\
+  MDW__DIGIT(   1000000, d), \\
+  MDW__DIGIT(    100000, d), \\
+  MDW__DIGIT(     10000, d), \\
+  MDW__DIGIT(      1000, d), \\
+  MDW__DIGIT(       100, d), \\
+  MDW__DIGIT(        10, d)
+
+/* Increasingly huge divisions.  PN divides by 10^(9*2^N). */
+#define MDW__P0 /MDW__G
+#define MDW__P1 MDW__P0 MDW__P0
+#define MDW__P2 MDW__P1 MDW__P1
+#define MDW__P3 MDW__P2 MDW__P2
+#define MDW__P4 MDW__P3 MDW__P3
+#define MDW__P5 MDW__P4 MDW__P4
+
+/* Increasingly long sequences of digits.  DN(P) produces the 9 * 2^N
+ * digits after applying divisions P.
+ */
+#define MDW__D0(p) MDW__9DIGITS(p MDW__P0), MDW__9DIGITS(p MDW__E)
+#define MDW__D1(p) MDW__D0(p MDW__P1), MDW__D0(p)
+#define MDW__D2(p) MDW__D1(p MDW__P2), MDW__D1(p)
+#define MDW__D3(p) MDW__D2(p MDW__P3), MDW__D2(p)
+#define MDW__D4(p) MDW__D3(p MDW__P4), MDW__D3(p)
+
+/* Ensure that our exponential cascade is sufficient to represent the
+ * expression.
+ */
+MDW__ASSERT(MDW__SHIFT(MDW__PROBE_EXPR, 10, MDW__P5) == 0);
+
+/* Format the output.  Everything is taken care of except the bottommost
+ * digit, which we handle seaprately because we actually want a \`leading'
+ * zero here if the constant value is actually zero.  Format it so that
+ * we can extract it from the object file.
+ */
+const char mdw__probe_output[] = {
+  '\\n',
+  'm', 'd', 'w', '-',
+  'p', 'r', 'o', 'b', 'e', '-',
+  'v', 'a', 'l', 'u', 'e', '=', '\"',
+  (MDW__PROBE_EXPR < 0 ? '-' : ' '),
+  MDW__D4(MDW__E),
+  '0' + MDW__RAW_DIGIT(MDW__PROBE_EXPR),
+  '\"', '\\n'
+};]"])
+
+AC_DEFUN([mdw_PROBE_CONSTANT],
+  [AC_REQUIRE([mdw__PROBE_CONSTANT_SETUP])
+   AC_COMPILE_IFELSE([AC_LANG_SOURCE([[$3
+#define MDW__PROBE_EXPR ($2)
+$mdw__probe_constant_body]])],
+     [$1=$(sed -n \
+         's:^mdw-probe-value="\(-\|\) *\([[0-9]]*\)"$:\1\2:p' conftest.o)],
+     [m4_if([$4], [],
+           [AC_MSG_FAILURE([failed to evaluate expression])],
+           [$4])])])
diff --git a/m4/mdw-uint-bits.m4 b/m4/mdw-uint-bits.m4
new file mode 100644 (file)
index 0000000..d5e678f
--- /dev/null
@@ -0,0 +1,91 @@
+dnl -*-autoconf-*-
+
+### SYNOPSIS
+###
+###   dnl mdw_UINT_BITS(TYPE, ABBREV)
+###
+### DESCRIPTION
+###
+###   Defines ABBREV_BITS to be the number of bits representable by the
+###   unsigned type named TYPE.  This works even if the compiler in question
+###   is a cross-compiler, and correctly handles systems without 8-bit bytes,
+###   or which have hidden bits in their integer representations.
+###
+### LICENSE
+###
+###   Copyright (c) 2013 Mark Wooding <mdw@distorted.org.uk>
+###
+###   This program is free software: you can redistribute it and/or modify it
+###   under the terms of the GNU General Public License as published by the
+###   Free Software Foundation, either version 2 of the License, or (at your
+###   option) any later version.
+###
+###   This program is distributed in the hope that it will be useful, but
+###   WITHOUT ANY WARRANTY; without even the implied warranty of
+###   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+###   General Public License for more details.
+###
+###   You should have received a copy of the GNU General Public License along
+###   with this program. If not, see <http://www.gnu.org/licenses/>.
+###
+###   In particular, no exception to the GPL is granted regarding generated
+###   `configure' scripts which are the output of Autoconf.
+
+dnl Algorithm:
+dnl
+dnl We first find an upper bound on the number of bits in the type as
+dnl sizeof(TYPE) * CHAR_BIT.  This is usually a good guess at the actual
+dnl value, since most modern systems don't have hidden bits, so we test to
+dnl see if it's correct.  If not, we start a binary search to find the true
+dnl value.
+dnl
+dnl To test a guess n, we form X = (TYPE)~(TYPE)0, which is the largest value
+dnl representable as a TYPE, and attempt to shift it right by n places.  This
+dnl is a little subtle, since shifting by more than the word length is
+dnl undefined behaviour.  The implementation isn't completely perfect.  It
+dnl assumes that a shift of 15 bits must be allowed (which is reasonable,
+dnl since an int must be at least 16 bits, and the standard integer
+dnl promotions will map shorter integer types to an int).  If n > 15, we
+dnl shift down in two stages, once by dividing by 2^15 = 32768, and once by
+dnl shifting by n - 15; otherwise we just shift.  (This strange chicanery is
+dnl to deal with Microsoft compilers which incorrectly `optimize' adjacent
+dnl correct shifts into a single invalid one.)  If X >> n is nonzero then n
+dnl is too small; if X >> (n - 1) is zero then n is too large; if neither
+dnl condition holds then n is correct.
+dnl
+dnl This is, unfortunately, logarithmic if the initial guess is wrong, but
+dnl that will happen rarely on interesting platforms.  Sites wanting to speed
+dnl up the configuration process can pre-seed the configuration cache.  If
+dnl anyone really cares, we can detect a native compiler (as opposed to a
+dnl cross-compiler) and do the binary-search in C.
+
+# Serial 1
+AC_DEFUN([mdw_UINT_BITS],
+  [AC_CACHE_CHECK([size of type `$1' in bits], [mdw_cv_$2_bits],
+    [mdw_PROBE_CONSTANT([guess], [sizeof($1) * CHAR_BIT], [
+#include <limits.h>
+#ifdef HAVE_STDINT_H
+#  include <stdint.h>
+#endif])
+     down=0
+     while :; do
+       mdw_PROBE_CONSTANT([answer],
+                         [low($guess) ? -1 : low($guess - 1) ? 0 : +1], [
+#ifdef HAVE_STDINT_H
+#  include <stdint.h>
+#endif
+#define srs(x, n) ((n) <= 15 ? (x) >> (n) : ((x)/32768) >> ((n) - 15))
+#define max (($1)~($1)0 | 0u)
+#define low(b) (srs(max, b) != 0)
+])
+       case "$answer" in
+        0) break;;
+        1) up=$guess;;
+        -1) down=$guess;;
+       esac
+       guess=$(expr \( $up + $down \) / 2)
+     done
+     mdw_cv_$2_bits=$guess])
+   AS_TR_SH($2_bits)=$mdw_cv_$2_bits
+   AC_DEFINE_UNQUOTED(AS_TR_CPP([$2_BITS]), [$mdw_cv_$2_bits],
+                     [number of bits in an object of type `$1'])])
index dfc5b58..7cabf89 100644 (file)
@@ -34,19 +34,44 @@ libmath_la_LIBADD    =
 TEST_LIBS               = libmath.la
 
 ###--------------------------------------------------------------------------
-### Main multiprecision integer library.
+### Representation of multiprecision integers.
+
+## The `mpgen' tool for dealing with these things.
+mpgen                   = $(srcdir)/mpgen
+EXTRA_DIST             += mpgen
+AM_V_MPGEN              = $(AM_V_MPGEN_$(V))
+AM_V_MPGEN_             = $(AM_V_MPGEN_$(AM_DEFAULT_VERBOSITY))
+AM_V_MPGEN_0            = @echo "  MPGEN  $@";
+MPGEN                   = $(AM_V_MPGEN)$(PYTHON) $(mpgen)
+
+## The type information collected by `configure'.
+CLEANFILES             += typeinfo.py
+EXTRA_DIST             += typeinfo.py.in
+typeinfo.py: $(srcdir)/typeinfo.py.in Makefile
+       $(SUBST) $(srcdir)/typeinfo.py.in >typeinfo.py.new \
+               type_bits="@type_bits@" \
+               limits="@limits@" && \
+       mv typeinfo.py.new typeinfo.py
+
+## The header file containing our representation choices.
+BUILT_SOURCES          += mptypes.h
+CLEANFILES             += mptypes.h
+nodist_archinclude_HEADERS += mptypes.h
+mptypes.h: $(mpgen) typeinfo.py
+       $(MPGEN) mptypes >mptypes.h.in && mv mptypes.h.in mptypes.h
+
+## Limits of C types as multiprecision integers.
+BUILT_SOURCES          += mplimits.h mplimits.c
+CLEANFILES             += mplimits.h mplimits.c
+nodist_archinclude_HEADERS += mplimits.h
+nodist_libmath_la_SOURCES += mplimits.c
+mplimits.h: $(mpgen) typeinfo.py
+       $(MPGEN) mplimits_h >mplimits.h.in && mv mplimits.h.in mplimits.h
+mplimits.c: $(mpgen) typeinfo.py
+       $(MPGEN) mplimits_c >mplimits.c.in && mv mplimits.c.in mplimits.c
 
-## This library is unfortunately intertwined with some of the code generation
-## programs, so we must be rather careful.  The important bits of the maths
-## library needed by these programs is separated out into `libmpbase'.  There
-## is work going on to fix this unpleasant situation by generating the
-## relevant files from Python scripts rather than C programs, using
-## information gathered by `configure'.
-noinst_LTLIBRARIES     += libmpbase.la
-libmath_la_LIBADD      += libmpbase.la
-libmpbase_la_LIBADD     = $(mLib_LIBS)
-libmpbase_la_SOURCES    =
-$(libmpbase_la_OBJECTS): mptypes.h
+###--------------------------------------------------------------------------
+### Main multiprecision integer library.
 
 ## Additional buffer I/O functions for mathematical objects.
 pkginclude_HEADERS     += buf.h
@@ -58,17 +83,17 @@ libmath_la_SOURCES  += exp.c
 
 ## Main user-visible multiprecision arithmetic.
 pkginclude_HEADERS     += mp.h
-libmpbase_la_SOURCES   += mp-arith.c
+libmath_la_SOURCES     += mp-arith.c
 TESTS                  += mp-arith.$t
-libmpbase_la_SOURCES   += mp-const.c
+libmath_la_SOURCES     += mp-const.c
 libmath_la_SOURCES     += mp-exp.c mp-exp.h
 libmath_la_SOURCES     += mp-gcd.c
 TESTS                  += mp-gcd.$t
-libmpbase_la_SOURCES   += mp-io.c
+libmath_la_SOURCES     += mp-io.c
 libmath_la_SOURCES     += mp-jacobi.c
 TESTS                  += mp-jacobi.$t
-libmpbase_la_SOURCES   += mp-mem.c
-libmpbase_la_SOURCES   += mp-misc.c
+libmath_la_SOURCES     += mp-mem.c
+libmath_la_SOURCES     += mp-misc.c
 libmath_la_SOURCES     += mp-modexp.c
 TESTS                  += mp-modexp.$t
 libmath_la_SOURCES     += mp-modsqrt.c
@@ -85,7 +110,7 @@ TESTS                        += mp-fibonacci.$t
 
 ## Special memory allocation for multiprecision integers.
 pkginclude_HEADERS     += mparena.h
-libmpbase_la_SOURCES   += mparena.c
+libmath_la_SOURCES     += mparena.c
 
 ## Barrett reduction, an efficient method for modular reduction.
 pkginclude_HEADERS     += mpbarrett.h
@@ -108,22 +133,6 @@ libmath_la_SOURCES += mpint.c
 TESTS                  += mpint.$t
 EXTRA_DIST             += t/mpint
 
-## Table of upper and lower limits of various types of machine integers, as
-## multiprecision integers.
-nodist_archinclude_HEADERS += mplimits.h
-nodist_libmath_la_SOURCES += mplimits.c
-CLEANFILES             += mplimits.h mplimits.c
-noinst_PROGRAMS                += genlimits
-genlimits_LDADD                 = libmpbase.la
-mplimits.c: genlimits$e
-       $(AM_V_GEN)./genlimits c >mplimits.c.new && \
-               mv mplimits.c.new mplimits.c
-mplimits.h: genlimits$e
-       $(AM_V_GEN)./genlimits h >mplimits.h.new && \
-               mv mplimits.h.new mplimits.h
-$(genlimits_OBJECTS): mptypes.h
-mplimits.lo: mplimits.h
-
 ## Montgomery reduction, a clever method for modular arithmetic.
 pkginclude_HEADERS     += mpmont.h
 libmath_la_SOURCES     += mpmont.c
@@ -150,32 +159,24 @@ EXTRA_DIST                += t/mpreduce
 
 ## Iteratiion over the bianry representation of multiprecision integers.
 pkginclude_HEADERS     += mpscan.h
-libmpbase_la_SOURCES   += mpscan.c
+libmath_la_SOURCES     += mpscan.c
 
 ## Conversion between multiprecision integers and their textual
 ## representations.
 pkginclude_HEADERS     += mptext.h
-libmpbase_la_SOURCES   += mptext.c
+libmath_la_SOURCES     += mptext.c
 TESTS                  += mptext.$t
 libmath_la_SOURCES     += mptext-dstr.c
 libmath_la_SOURCES     += mptext-file.c
 libmath_la_SOURCES     += mptext-len.c
-libmpbase_la_SOURCES   += mptext-string.c
+libmath_la_SOURCES     += mptext-string.c
 EXTRA_DIST             += t/mptext
 
-## Basic types used in the representation of multiprecision integers.
-nodist_archinclude_HEADERS += mptypes.h
-BUILT_SOURCES          += mptypes.h
-CLEANFILES             += mptypes.h
-noinst_PROGRAMS                += mptypes
-mptypes.h: mptypes$e
-       $(AM_V_GEN)./mptypes >mptypes.h.new && mv mptypes.h.new mptypes.h
-
 ## Low-level multiprecision arithmetic.
 pkginclude_HEADERS     += mpx.h bitops.h mpw.h
-libmpbase_la_SOURCES   += mpx.c
+libmath_la_SOURCES     += mpx.c
 TESTS                  += mpx.$t
-libmpbase_la_SOURCES   += karatsuba.h mpx-kmul.c mpx-ksqr.c
+libmath_la_SOURCES     += karatsuba.h mpx-kmul.c mpx-ksqr.c
 TESTS                  += mpx-kmul.$t mpx-ksqr.$t
 noinst_PROGRAMS                += bittest
 TESTS                  += bittest
@@ -338,28 +339,21 @@ libmath_la_SOURCES        += f-prime.c
 
 ## Table of built-in binary fields.
 pkginclude_HEADERS     += bintab.h
-libmath_la_SOURCES     += bintab.c
+nodist_libmath_la_SOURCES += bintab.c
 CLEANFILES             += bintab.c
-EXTRA_DIST             += bintab.in bin-gentab.awk
-bintab.c: bintab.in bin-gentab.awk mpdump$e
-       $(AM_V_GEN)awk -f $(srcdir)/bin-gentab.awk \
-               <$(srcdir)/bintab.in >bintab.c.new && \
-       mv bintab.c.new bintab.c
+EXTRA_DIST             += bintab.in
+bintab.c: $(mpgen) typeinfo.py bintab.in
+       $(MPGEN) bintab $(srcdir)/bintab.in >bintab.c.new && \
+               mv bintab.c.new bintab.c
 
 ## Table of built-in prime fields.
 pkginclude_HEADERS     += ptab.h
-libmath_la_SOURCES     += ptab.c
+nodist_libmath_la_SOURCES += ptab.c
 CLEANFILES             += ptab.c
-EXTRA_DIST             += ptab.in p-gentab.awk
-ptab.c: ptab.in p-gentab.awk mpdump$e
-       $(AM_V_GEN)awk -f $(srcdir)/p-gentab.awk \
-               <$(srcdir)/ptab.in >ptab.c.new && \
-       mv ptab.c.new ptab.c
-
-## A utility for building multiprecision integer constants.
-noinst_PROGRAMS                += mpdump
-mpdump_LDADD            = libmpbase.la
-$(mpdump_OBJECTS): mptypes.h
+EXTRA_DIST             += ptab.in
+ptab.c: $(mpgen) typeinfo.py ptab.in
+       $(MPGEN) ptab $(srcdir)/ptab.in >ptab.c.new && \
+               mv ptab.c.new ptab.c
 
 ###--------------------------------------------------------------------------
 ### Elliptic curve arithmetic.
@@ -389,14 +383,13 @@ pkginclude_HEADERS        += ec-test.h
 libmath_la_SOURCES     += ec-test.c
 TESTS                  += ec-test.$t
 
-## A table of built-in elliptic curves.
+## Table of built-in elliptic-curve groups.
 pkginclude_HEADERS     += ectab.h
-libmath_la_SOURCES     += ectab.c
+nodist_libmath_la_SOURCES += ectab.c
 CLEANFILES             += ectab.c
-EXTRA_DIST             += ectab.in ec-gentab.awk
-ectab.c: ectab.in ec-gentab.awk mpdump$e
-       $(AM_V_GEN)awk -f $(srcdir)/ec-gentab.awk \
-               <$(srcdir)/ectab.in >ectab.c.new && \
-       mv ectab.c.new ectab.c
+EXTRA_DIST             += ectab.in
+ectab.c: $(mpgen) typeinfo.py ectab.in
+       $(MPGEN) ectab $(srcdir)/ectab.in >ectab.c.new && \
+               mv ectab.c.new ectab.c
 
 ###----- That's all, folks --------------------------------------------------
diff --git a/math/mpgen b/math/mpgen
new file mode 100644 (file)
index 0000000..4ed4f16
--- /dev/null
@@ -0,0 +1,481 @@
+#! @PYTHON@
+###
+### Generate multiprecision integer representations
+###
+### (c) 2013 Straylight/Edgeware
+###
+
+###----- Licensing notice ---------------------------------------------------
+###
+### This file is part of Catacomb.
+###
+### Catacomb is free software; you can redistribute it and/or modify
+### it under the terms of the GNU Library General Public License as
+### published by the Free Software Foundation; either version 2 of the
+### License, or (at your option) any later version.
+###
+### Catacomb is distributed in the hope that it will be useful,
+### but WITHOUT ANY WARRANTY; without even the implied warranty of
+### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+### GNU Library General Public License for more details.
+###
+### You should have received a copy of the GNU Library General Public
+### License along with Catacomb; if not, write to the Free
+### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+### MA 02111-1307, USA.
+
+from __future__ import with_statement
+
+import re as RX
+import optparse as OP
+import types as TY
+
+from sys import stdout
+
+###--------------------------------------------------------------------------
+### Random utilities.
+
+def write_header(mode, name):
+  stdout.write("""\
+/* -*-c-*- GENERATED by mpgen (%s)
+ *
+ * %s
+ */
+
+""" % (mode, name))
+
+def write_banner(text):
+  stdout.write("/*----- %s %s*/\n" % (text, '-' * (66 - len(text))))
+
+class struct (object): pass
+
+R_IDBAD = RX.compile('[^0-9A-Za-z]')
+def fix_name(name): return R_IDBAD.sub('_', name)
+
+###--------------------------------------------------------------------------
+### Determining the appropriate types.
+
+TYPEMAP = {}
+
+class IntClass (type):
+  def __new__(cls, name, supers, dict):
+    c = type.__new__(cls, name, supers, dict)
+    try: TYPEMAP[c.tag] = c
+    except AttributeError: pass
+    return c
+
+class BasicIntType (object):
+  __metaclass__ = IntClass
+  preamble = ''
+  typedef_prefix = ''
+  literalfmt = '%su'
+  def __init__(me, bits, rank):
+    me.bits = bits
+    me.rank = rank
+    me.litwd = len(me.literal(0))
+  def literal(me, value, fmt = None):
+    if fmt is None: fmt = '0x%0' + str((me.bits + 3)//4) + 'x'
+    return me.literalfmt % (fmt % value)
+
+class UnsignedCharType (BasicIntType):
+  tag = 'uchar'
+  name = 'unsigned char'
+
+class UnsignedShortType (BasicIntType):
+  tag = 'ushort'
+  name = 'unsigned short'
+
+class UnsignedIntType (BasicIntType):
+  tag = 'uint'
+  name = 'unsigned int'
+
+class UnsignedLongType (BasicIntType):
+  tag = 'ulong'
+  name = 'unsigned long'
+  literalfmt = '%sul'
+
+class UnsignedLongLongType (BasicIntType):
+  tag = 'ullong'
+  name = 'unsigned long long'
+  preamble = """
+#if __GNUC__ > 2 || (__GNUC__ == 2 && __GNUC_MINOR__ >= 91)
+#  define CATACOMB_GCC_EXTENSION __extension__
+#else
+#  define CATACOMB_GCC_EXTENSION
+#endif
+"""
+  typedef_prefix = 'CATACOMB_GCC_EXTENSION '
+  literalfmt = 'CATACOMB_GCC_EXTENSION %sull'
+
+class UIntMaxType (BasicIntType):
+  tag = 'uintmax'
+  name = 'uintmax_t'
+  preamble = "\n#include <stdint.h>\n"
+
+class TypeChoice (object):
+  def __init__(me, tifile):
+
+    ## Load the captured type information.
+    me.ti = TY.ModuleType('typeinfo')
+    execfile(opts.typeinfo, me.ti.__dict__)
+
+    ## Build a map of the available types.
+    tymap = {}
+    byrank = []
+    for tag, bits in me.ti.TYPEINFO:
+      rank = len(byrank)
+      tymap[tag] = rank
+      byrank.append(TYPEMAP[tag](bits, rank))
+
+    ## First pass: determine a suitable word size.  The criteria are (a)
+    ## there exists another type at least twice as long (so that we can do a
+    ## single x single -> double multiplication), and (b) operations on a
+    ## word are efficient (so we'd prefer a plain machine word).  We'll start
+    ## at `int' and work down.  Maybe this won't work: there's a plan B.
+    mpwbits = 0
+    i = tymap['uint']
+    while not mpwbits and i >= 0:
+      ibits = byrank[i].bits
+      for j in xrange(i + 1, len(byrank)):
+        if byrank[j].bits >= 2*ibits:
+          mpwbits = ibits
+          break
+
+    ## If that didn't work, then we'll start with the largest type available
+    ## and go with half its size.
+    if not mpwbits:
+      mpwbits = byrank[-1].bits//2
+
+    ## Make sure we've not ended up somewhere really silly.
+    if mpwbits < 16:
+      raise Exception, "`mpw' type is too small: your C environment is weird"
+
+    ## Now figure out suitable types for `mpw' and `mpd'.
+    def find_type(bits, what):
+      for ty in byrank:
+        if ty.bits >= bits: return ty
+      raise Exception, \
+          "failed to find suitable %d-bit type, for %s" % (bits, what)
+
+    ## Store our decisions.
+    me.mpwbits = mpwbits
+    me.mpw = find_type(mpwbits, 'mpw')
+    me.mpd = find_type(mpwbits*2, 'mpd')
+
+###--------------------------------------------------------------------------
+### Outputting constant multiprecision integers.
+
+MARGIN = 72
+
+def write_preamble():
+  stdout.write("""
+#include <mLib/macros.h>
+#define MP_(name, flags) \\
+  { (/*unconst*/ mpw *)name##__mpw, \\
+    (/*unconst*/ mpw *)name##__mpw + N(name##__mpw), \\
+    N(name##__mpw), 0, MP_CONST | flags, 0 }
+#define ZERO_MP { 0, 0, 0, 0, MP_CONST, 0 }
+#define POS_MP(name) MP_(name, 0)
+#define NEG_MP(name) MP_(name, MP_NEG)
+""")
+
+def write_limbs(name, x):
+  if not x: return
+  stdout.write("\nstatic const mpw %s__mpw[] = {" % name)
+  sep = ''
+  pos = MARGIN
+  if x < 0: x = -x
+  mask = (1 << TC.mpwbits) - 1
+
+  while x > 0:
+    w, x = x & mask, x >> TC.mpwbits
+    f = TC.mpw.literal(w)
+    if pos + 2 + len(f) <= MARGIN:
+      stdout.write(sep + ' ' + f)
+    else:
+      pos = 2
+      stdout.write(sep + '\n  ' + f)
+    pos += len(f) + 2
+    sep = ','
+
+  stdout.write("\n};\n")
+
+def mp_body(name, x):
+  return "%s_MP(%s)" % (x >= 0 and "POS" or "NEG", name)
+
+###--------------------------------------------------------------------------
+### Mode definition machinery.
+
+MODEMAP = {}
+
+def defmode(func):
+  name = func.func_name
+  if name.startswith('m_'): name = name[2:]
+  MODEMAP[name] = func
+  return func
+
+###--------------------------------------------------------------------------
+### The basic types header.
+
+@defmode
+def m_mptypes():
+  write_header("mptypes", "mptypes.h")
+  stdout.write("""\
+#ifndef CATACOMB_MPTYPES_H
+#define CATACOMB_MPTYPES_H
+""")
+
+  have = set([TC.mpw, TC.mpd])
+  for t in have:
+    stdout.write(t.preamble)
+
+  for label, t, bits in [('mpw', TC.mpw, TC.mpwbits),
+                         ('mpd', TC.mpd, TC.mpwbits*2)]:
+    LABEL = label.upper()
+    stdout.write("\n%stypedef %s %s;\n" % (t.typedef_prefix, t.name, label))
+    stdout.write("#define %s_BITS %d\n" % (LABEL, bits))
+    i = 1
+    while 2*i < bits: i *= 2
+    stdout.write("#define %s_P2 %d\n" % (LABEL, i))
+    stdout.write("#define %s_MAX %s\n" % (LABEL,
+                                          t.literal((1 << bits) - 1, "%d")))
+
+  stdout.write("\n#endif\n")
+
+###--------------------------------------------------------------------------
+### Constant tables.
+
+@defmode
+def m_mplimits_c():
+  write_header("mplimits_c", "mplimits.c")
+  stdout.write('#include "mplimits.h"\n')
+  write_preamble()
+  seen = {}
+  v = []
+  def write(x):
+    if not x or x in seen: return
+    seen[x] = 1
+    write_limbs('limits_%d' % len(v), x)
+    v.append(x)
+  for tag, lo, hi in TC.ti.LIMITS:
+    write(lo)
+    write(hi)
+
+  stdout.write("\nmp mp_limits[] = {")
+  i = 0
+  sep = "\n  "
+  for x in v:
+    stdout.write("%s%s_MP(limits_%d)\n" % (sep, x < 0 and "NEG" or "POS", i))
+    i += 1
+    sep = ",\n  "
+  stdout.write("\n};\n");
+
+@defmode
+def m_mplimits_h():
+  write_header("mplimits_h", "mplimits.h")
+  stdout.write("""\
+#ifndef CATACOMB_MPLIMITS_H
+#define CATACOMB_MPLIMITS_H
+
+#ifndef CATACOMB_MP_H
+#  include "mp.h"
+#endif
+
+extern mp mp_limits[];
+
+""")
+
+  seen = { 0: "MP_ZERO" }
+  slot = [0]
+  def find(x):
+    try:
+      r = seen[x]
+    except KeyError:
+      r = seen[x] = '(&mp_limits[%d])' % slot[0]
+      slot[0] += 1
+    return r
+  for tag, lo, hi in TC.ti.LIMITS:
+    stdout.write("#define MP_%s_MIN %s\n" % (tag, find(lo)))
+    stdout.write("#define MP_%s_MAX %s\n" % (tag, find(hi)))
+
+  stdout.write("\n#endif\n")
+
+###--------------------------------------------------------------------------
+### Group tables.
+
+class GroupTableClass (type):
+  def __new__(cls, name, supers, dict):
+    c = type.__new__(cls, name, supers, dict)
+    try: mode = c.mode
+    except AttributeError: pass
+    else: MODEMAP[c.mode] = c.run
+    return c
+
+class GroupTable (object):
+  __metaclass__ = GroupTableClass
+  keyword = 'group'
+  slots = []
+  def __init__(me):
+    me.st = st = struct()
+    st.nextmp = 0
+    st.mpmap = { None: 'NO_MP', 0: 'ZERO_MP' }
+    st.d = {}
+    me.st.name = None
+    me._names = []
+    me._defs = set()
+    me._slotmap = dict([(s.name, s) for s in me.slots])
+    me._headslots = [s for s in me.slots if s.headline]
+  def _flush(me):
+    if me.st.name is None: return
+    stdout.write("/* --- %s --- */\n" % me.st.name)
+    for s in me.slots: s.setup(me.st)
+    stdout.write("\nstatic %s c_%s = {" % (me.data_t, fix_name(me.st.name)))
+    sep = "\n  "
+    for s in me.slots:
+      stdout.write(sep)
+      s.write(me.st)
+      sep = ",\n  "
+    stdout.write("\n};\n\n")
+    me.st.d = {}
+    me.st.name = None
+  @classmethod
+  def run(cls, input):
+    me = cls()
+    write_header(me.mode, me.filename)
+    stdout.write('#include "%s"\n' % me.header)
+    write_preamble()
+    stdout.write("#define NO_MP { 0, 0, 0, 0, 0, 0 }\n\n")
+    write_banner("Group data")
+    stdout.write('\n')
+    with open(input) as file:
+      for line in file:
+        ff = line.split()
+        if not ff or ff[0].startswith('#'): continue
+        if ff[0] == 'alias':
+          if len(ff) != 3: raise Exception, "wrong number of alias arguments"
+          me._flush()
+          me._names.append((ff[1], ff[2]))
+        elif ff[0] == me.keyword:
+          if len(ff) < 2 or len(ff) > 2 + len(me._headslots):
+            raise Exception, "bad number of headline arguments"
+          me._flush()
+          me.st.name = name = ff[1]
+          me._defs.add(name)
+          me._names.append((name, name))
+          for f, s in zip(ff[2:], me._headslots): s.set(me.st, f)
+        elif ff[0] in me._slotmap:
+          if len(ff) != 2:
+            raise Exception, "bad number of values for slot `%s'" % ff[0]
+          me._slotmap[ff[0]].set(me.st, ff[1])
+        else:
+          raise Exception, "unknown keyword `%s'" % ff[0]
+    me._flush()
+    write_banner("Main table")
+    stdout.write("\nconst %s %s[] = {\n" % (me.entry_t, me.tabname))
+    for a, n in me._names:
+      if n not in me._defs:
+        raise Exception, "alias `%s' refers to unknown group `%s'" % (a, n)
+      stdout.write('  { "%s", &c_%s },\n' % (a, fix_name(n)))
+    stdout.write("  { 0, 0 }\n};\n\n")
+    write_banner("That's all, folks")
+
+class BaseSlot (object):
+  def __init__(me, name, headline = False, omitp = None, allowp = None):
+    me.name = name
+    me.headline = headline
+    me._omitp = None
+    me._allowp = None
+  def set(me, st, value):
+    if me._allowp and not me._allowp(st, value):
+      raise Exception, "slot `%s' not allowed here" % me.name
+    st.d[me] = value
+  def setup(me, st):
+    if me not in st.d and (not me._omitp or not me._omitp(st)):
+      raise Exception, "missing slot `%s'" % me.name
+
+class EnumSlot (BaseSlot):
+  def __init__(me, name, prefix, values, **kw):
+    super(EnumSlot, me).__init__(name, **kw)
+    me._values = set(values)
+    me._prefix = prefix
+  def set(me, st, value):
+    if value not in me._values:
+      raise Exception, "invalid %s value `%s'" % (me.name, value)
+    super(EnumSlot, me).set(st, value)
+  def write(me, st):
+    try: stdout.write('%s_%s' % (me._prefix, st.d[me].upper()))
+    except KeyError: stdout.write('0')
+
+class MPSlot (BaseSlot):
+  def setup(me, st):
+    v = st.d.get(me)
+    if v not in st.mpmap:
+      write_limbs('v%d' % st.nextmp, v)
+      st.mpmap[v] = mp_body('v%d' % st.nextmp, v)
+      st.nextmp += 1
+  def set(me, st, value):
+    super(MPSlot, me).set(st, long(value, 0))
+  def write(me, st):
+    stdout.write(st.mpmap[st.d.get(me)])
+
+class BinaryGroupTable (GroupTable):
+  mode = 'bintab'
+  filename = 'bintab.c'
+  header = 'bintab.h'
+  data_t = 'bindata'
+  entry_t = 'binentry'
+  tabname = 'bintab'
+  slots = [MPSlot('p'), MPSlot('q'), MPSlot('g')]
+
+class EllipticCurveTable (GroupTable):
+  mode = 'ectab'
+  filename = 'ectab.c'
+  header = 'ectab.h'
+  keyword = 'curve'
+  data_t = 'ecdata'
+  entry_t = 'ecentry'
+  tabname = 'ectab'
+  slots = [EnumSlot('type', 'FTAG',
+                    ['prime', 'niceprime', 'binpoly', 'binnorm'],
+                    headline = True),
+           MPSlot('p'),
+           MPSlot('beta',
+                  allowp = lambda st, _: st.d['type'] == 'binnorm',
+                  omitp = lambda st: st.d['type'] != 'binnorm'),
+           MPSlot('a'), MPSlot('b'), MPSlot('r'), MPSlot('h'),
+           MPSlot('gx'), MPSlot('gy')]
+
+class PrimeGroupTable (GroupTable):
+  mode = 'ptab'
+  filename = 'ptab.c'
+  header = 'ptab.h'
+  data_t = 'pdata'
+  entry_t = 'pentry'
+  tabname = 'ptab'
+  slots = [MPSlot('p'), MPSlot('q'), MPSlot('g')]
+
+###--------------------------------------------------------------------------
+### Main program.
+
+op = OP.OptionParser(
+  description = 'Generate multiprecision integer representations',
+  usage = 'usage: %prog [-t TYPEINFO] MODE [ARGS ...]',
+  version = 'Catacomb, version @VERSION@')
+for shortopt, longopt, kw in [
+  ('-t', '--typeinfo', dict(
+      action = 'store', metavar = 'PATH', dest = 'typeinfo',
+      help = 'alternative typeinfo file'))]:
+  op.add_option(shortopt, longopt, **kw)
+op.set_defaults(typeinfo = './typeinfo.py')
+opts, args = op.parse_args()
+
+if len(args) < 1: op.error('missing MODE')
+mode = args[0]
+
+TC = TypeChoice(opts.typeinfo)
+
+try: modefunc = MODEMAP[mode]
+except KeyError: op.error("unknown mode `%s'" % mode)
+modefunc(*args[1:])
+
+###----- That's all, folks --------------------------------------------------
diff --git a/math/typeinfo.py.in b/math/typeinfo.py.in
new file mode 100644 (file)
index 0000000..3ff9b7c
--- /dev/null
@@ -0,0 +1,4 @@
+### -*-python-*-
+
+TYPEINFO = [@type_bits@]
+LIMITS = [@limits@]