From 1c3d4cf54a0edd484c4405f5332d39bb17f1aee0 Mon Sep 17 00:00:00 2001 From: Mark Wooding Date: Wed, 5 Jun 2013 17:14:30 +0100 Subject: [PATCH] Overhaul `math' representation machinery. 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 | 1 + configure.ac | 64 +++++++ m4/.gitignore | 2 + m4/mdw-probe-constant.m4 | 183 ++++++++++++++++++ m4/mdw-uint-bits.m4 | 91 +++++++++ math/Makefile.am | 135 +++++++------ math/mpgen | 481 +++++++++++++++++++++++++++++++++++++++++++++++ math/typeinfo.py.in | 4 + 8 files changed, 890 insertions(+), 71 deletions(-) create mode 100644 m4/.gitignore create mode 100644 m4/mdw-probe-constant.m4 create mode 100644 m4/mdw-uint-bits.m4 create mode 100644 math/mpgen create mode 100644 math/typeinfo.py.in diff --git a/Makefile.am b/Makefile.am index 4cacfbe..e29c36e 100644 --- a/Makefile.am +++ b/Makefile.am @@ -25,6 +25,7 @@ ### MA 02111-1307, USA. include $(top_srcdir)/vars.am +ACLOCAL_AMFLAGS = -Im4 SUBDIRS = diff --git a/configure.ac b/configure.ac index aca24fe..7c66807 100644 --- a/configure.ac +++ b/configure.ac @@ -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 +#include ]) + catacomb_COMPILE_TIME_CONSTANT([hi], [$3], [$1_max], +[#include +#include ]) + 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 index 0000000..94f2b51 --- /dev/null +++ b/m4/.gitignore @@ -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 index 0000000..3916582 --- /dev/null +++ b/m4/mdw-probe-constant.m4 @@ -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 +### +### 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 . +### +### 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 . + */ + +/* 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 index 0000000..d5e678f --- /dev/null +++ b/m4/mdw-uint-bits.m4 @@ -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 +### +### 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 . +### +### 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 +#ifdef HAVE_STDINT_H +# include +#endif]) + down=0 + while :; do + mdw_PROBE_CONSTANT([answer], + [low($guess) ? -1 : low($guess - 1) ? 0 : +1], [ +#ifdef HAVE_STDINT_H +# include +#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'])]) diff --git a/math/Makefile.am b/math/Makefile.am index dfc5b58..7cabf89 100644 --- a/math/Makefile.am +++ b/math/Makefile.am @@ -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 index 0000000..4ed4f16 --- /dev/null +++ b/math/mpgen @@ -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 \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 +#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 index 0000000..3ff9b7c --- /dev/null +++ b/math/typeinfo.py.in @@ -0,0 +1,4 @@ +### -*-python-*- + +TYPEINFO = [@type_bits@] +LIMITS = [@limits@] -- 2.11.0