lbn16.c

包含哈希,对称以及非对称的经典算法 包含经典事例

/*
 * lbn16.c - Low-level bignum routines, 16-bit version.
 *
 * Copyright (c) 1995  Colin Plumb.  All rights reserved.
 * For licensing and other legal details, see the file legal.c.
 *
 * NOTE: the magic constants "16" and "32" appear in many places in this
 * file, including inside identifiers.  Because it is not possible to
 * ask "#ifdef" of a macro expansion, it is not possible to use the
 * preprocessor to conditionalize these properly.  Thus, this file is
 * intended to be edited with textual search and replace to produce
 * alternate word size versions.  Any reference to the number of bits
 * in a word must be the string "16", and that string must not appear
 * otherwise.  Any reference to twice this number must appear as "32",
 * which likewise must not appear otherwise.  Is that clear?
 *
 * Remember, when doubling the bit size replace the larger number (32)
 * first, then the smaller (16).  When halving the bit size, do the
 * opposite.  Otherwise, things will get wierd.  Also, be sure to replace
 * every instance that appears.  (:%s/foo/bar/g in vi)
 *
 * These routines work with a pointer to the least-significant end of
 * an array of WORD16s.  The BIG(x), LITTLE(y) and BIGLTTLE(x,y) macros
 * defined in lbn.h (which expand to x on a big-edian machine and y on a
 * little-endian machine) are used to conditionalize the code to work
 * either way.  If you have no assembly primitives, it doesn't matter.
 * Note that on a big-endian machine, the least-significant-end pointer
 * is ONE PAST THE END.  The bytes are ptr[-1] through ptr[-len].
 * On little-endian, they are ptr[0] through ptr[len-1].  This makes
 * perfect sense if you consider pointers to point *between* bytes rather
 * than at them.
 *
 * Because the array index values are unsigned integers, ptr[-i]
 * may not work properly, since the index -i is evaluated as an unsigned,
 * and if pointers are wider, zero-extension will produce a positive
 * number rahter than the needed negative.  The expression used in this
 * code, *(ptr-i) will, however, work.  (The array syntax is equivalent
 * to *(ptr+-i), which is a pretty subtle difference.)
 *
 * Many of these routines will get very unhappy if fed zero-length inputs.
 * They use assert() to enforce this.  An higher layer of code must make
 * sure that these aren't called with zero-length inputs.
 *
 * Any of these routines can be replaced with more efficient versions
 * elsewhere, by just #defining their names.  If one of the names
 * is #defined, the C code is not compiled in and no declaration is
 * made.  Use the BNINCLUDE file to do that.  Typically, you compile
 * asm subroutines with the same name and just, e.g.
 * #define lbnMulAdd1_16 lbnMulAdd1_16
 *
 * If you want to write asm routines, start with lbnMulAdd1_16().
 * This is the workhorse of modular exponentiation.  lbnMulN1_16() is
 * also used a fair bit, although not as much and it's defined in terms
 * of lbnMulAdd1_16 if that has a custom version.  lbnMulSub1_16 and
 * lbnDiv21_16 are used in the usual division and remainder finding.
 * (Not the Montgomery reduction used in modular exponentiation, though.)
 * Once you have lbnMulAdd1_16 defined, writing the other two should
 * be pretty easy.  (Just make sure you get the sign of the subtraction
 * in lbnMulSub1_16 right - it's dest = dest - source * k.)
 *
 * The only definitions that absolutely need a double-word (BNWORD32)
 * type are lbnMulAdd1_16 and lbnMulSub1_16; if those are provided,
 * the rest follows.  lbnDiv21_16, however, is a lot slower unless you
 * have them, and lbnModQ_16 takes after it.  That one is used quite a
 * bit for prime sieving.
 */

#if HAVE_CONFIG_H
#include "config.h"
#endif

#if !NO_ASSERT_H
#include <assert.h>
#else
#define assert(x) (void)0
#endif

#if !NO_STRING_H
#include <string.h>	/* For memcpy */
#elif HAVE_STRINGS_H
#include <strings.h>
#endif
#if NEED_MEMORY_H
#include <memory.h>
#endif

#include "lbn.h"
#include "lbn16.h"
#include "lbnmem.h"
#include "legal.h"

#include "kludge.h"

#ifndef BNWORD16
#error 16-bit bignum library requires a 16-bit data type
#endif

/* Make sure the copyright notice gets included */
volatile const char * volatile const lbnCopyright = bnCopyright;

/*
 * Most of the multiply (and Montgomery reduce) routines use an outer
 * loop that iterates over one of the operands - a so-called operand
 * scanning approach.  One big advantage of this is that the assembly
 * support routines are simpler.  The loops can be rearranged to have
 * an outer loop that iterates over the product, a so-called product
 * scanning approach.  This has the advantage of writing less data
 * and doing fewer adds to memory, so is supposedly faster.  Some
 * code has been written using a product-scanning approach, but
 * it appears to be slower, so it is turned off by default.  Some
 * experimentation would be appreciated.
 *
 * (The code is also annoying to get right and not very well commented,
 * one of my pet peeves about math libraries.  I'm sorry.)
 */
#ifndef PRODUCT_SCAN
#define PRODUCT_SCAN 0
#endif

/*
 * Copy an array of words.  <Marvin mode on>  Thrilling, isn't it? </Marvin>
 * This is a good example of how the byte offsets and BIGLITTLE() macros work.
 * Another alternative would have been
 * memcpy(dest BIG(-len), src BIG(-len), len*sizeof(BNWORD16)), but I find that
 * putting operators into conditional macros is confusing.
 */
#ifndef lbnCopy_16
void
lbnCopy_16(BNWORD16 *dest, BNWORD16 const *src, unsigned len)
{
	memcpy(BIGLITTLE(dest-len,dest), BIGLITTLE(src-len,src),
	       len * sizeof(*src));
}
#endif /* !lbnCopy_16 */

/*
 * Fill n words with zero.  This does it manually rather than calling
 * memset because it can assume alignment to make things faster while
 * memset can't.  Note how big-endian numbers are naturally addressed
 * using predecrement, while little-endian is postincrement.
 */
#ifndef lbnZero_16
void
lbnZero_16(BNWORD16 *num, unsigned len)
{
	while (len--)
		BIGLITTLE(*--num,*num++) = 0;
}
#endif /* !lbnZero_16 */

/*
 * Negate an array of words.
 * Negation is subtraction from zero.  Negating low-order words
 * entails doing nothing until a non-zero word is hit.  Once that
 * is negated, a borrow is generated and never dies until the end
 * of the number is hit.  Negation with borrow, -x-1, is the same as ~x.
 * Repeat that until the end of the number.
 *
 * Doesn't return borrow out because that's pretty useless - it's
 * always set unless the input is 0, which is easy to notice in
 * normalized form.
 */
#ifndef lbnNeg_16
void
lbnNeg_16(BNWORD16 *num, unsigned len)
{
	assert(len);

	/* Skip low-order zero words */
	while (BIGLITTLE(*--num,*num) == 0) {
		if (!--len)
			return;
		LITTLE(num++;)
	}
	/* Negate the lowest-order non-zero word */
	*num = -*num;
	/* Complement all the higher-order words */
	while (--len) {
		BIGLITTLE(--num,++num);
		*num = ~*num;
	}
}
#endif /* !lbnNeg_16 */


/*
 * lbnAdd1_16: add the single-word "carry" to the given number.
 * Used for minor increments and propagating the carry after
 * adding in a shorter bignum.
 *
 * Technique: If we have a double-width word, presumably the compiler
 * can add using its carry in inline code, so we just use a larger
 * accumulator to compute the carry from the first addition.
 * If not, it's more complex.  After adding the first carry, which may
 * be > 1, compare the sum and the carry.  If the sum wraps (causing a
 * carry out from the addition), the result will be less than each of the
 * inputs, since the wrap subtracts a number (2^16) which is larger than
 * the other input can possibly be.  If the sum is >= the carry input,
 * return success immediately.
 * In either case, if there is a carry, enter a loop incrementing words
 * until one does not wrap.  Since we are adding 1 each time, the wrap
 * will be to 0 and we can test for equality.
 */
#ifndef lbnAdd1_16	/* If defined, it's provided as an asm subroutine */
#ifdef BNWORD32
BNWORD16
lbnAdd1_16(BNWORD16 *num, unsigned len, BNWORD16 carry)
{
	BNWORD32 t;
	assert(len > 0);	/* Alternative: if (!len) return carry */

	t = (BNWORD32)BIGLITTLE(*--num,*num) + carry;
	BIGLITTLE(*num,*num++) = (BNWORD16)t;
	if ((t >> 16) == 0)
		return 0;
	while (--len) {
		if (++BIGLITTLE(*--num,*num++) != 0)
			return 0;
	}
	return 1;
}
#else /* no BNWORD32 */
BNWORD16
lbnAdd1_16(BNWORD16 *num, unsigned len, BNWORD16 carry)
{
	assert(len > 0);	/* Alternative: if (!len) return carry */

	if ((BIGLITTLE(*--num,*num++) += carry) >= carry)
		return 0;
	while (--len) {
		if (++BIGLITTLE(*--num,*num++) != 0)
			return 0;
	}
	return 1;
}
#endif
#endif/* !lbnAdd1_16 */

/*
 * lbnSub1_16: subtract the single-word "borrow" from the given number.
 * Used for minor decrements and propagating the borrow after
 * subtracting a shorter bignum.
 *
 * Technique: Similar to the add, above.  If there is a double-length type,
 * use that to generate the first borrow.
 * If not, after subtracting the first borrow, which may be > 1, compare
 * the difference and the *negative* of the carry.  If the subtract wraps
 * (causing a borrow out from the subtraction), the result will be at least
 * as large as -borrow.  If the result < -borrow, then no borrow out has
 * appeared and we may return immediately, except when borrow == 0.  To
 * deal with that case, use the identity that -x = ~x+1, and instead of
 * comparing < -borrow, compare for <= ~borrow.
 * Either way, if there is a borrow out, enter a loop decrementing words
 * until a non-zero word is reached.
 *
 * Note the cast of ~borrow to (BNWORD16).  If the size of an int is larger
 * than BNWORD16, C rules say the number is expanded for the arithmetic, so
 * the inversion will be done on an int and the valie won't be quite what
 * is expected.
 */
#ifndef lbnSub1_16	/* If defined, it's provided as an asm subroutine */
#ifdef BNWORD32
BNWORD16
lbnSub1_16(BNWORD16 *num, unsigned len, BNWORD16 borrow)
{
	BNWORD32 t;
	assert(len > 0);	/* Alternative: if (!len) return borrow */

	t = (BNWORD32)BIGLITTLE(*--num,*num) - borrow;
	BIGLITTLE(*num,*num++) = (BNWORD16)t;
	if ((t >> 16) == 0)
		return 0;
	while (--len) {
		if ((BIGLITTLE(*--num,*num++))-- != 0)
			return 0;
	}
	return 1;
}
#else /* no BNWORD32 */
BNWORD16
lbnSub1_16(BNWORD16 *num, unsigned len, BNWORD16 borrow)
{
	assert(len > 0);	/* Alternative: if (!len) return borrow */

	if ((BIGLITTLE(*--num,*num++) -= borrow) <= (BNWORD16)~borrow)
		return 0;
	while (--len) {
		if ((BIGLITTLE(*--num,*num++))-- != 0)
			return 0;
	}
	return 1;
}
#endif
#endif /* !lbnSub1_16 */

/*
 * lbnAddN_16: add two bignums of the same length, returning the carry (0 or 1).
 * One of the building blocks, along with addn1, of adding two bignums of
 * differing lengths.
 *
 * Technique: Maintain a word of carry.  If there is no double-width type,
 * use the same technique as in addn1, above, to maintain the carry by
 * comparing the inputs.  Adding the carry sources is used as an OR operator;
 * at most one of the two comparisons can possibly be true.  The first can
 * only be true if carry == 1 and x, the result, is 0.  In that case the
 * second can't possibly be true.
 */
#ifndef lbnAddN_16
#ifdef BNWORD32
BNWORD16
lbnAddN_16(BNWORD16 *num1, BNWORD16 const *num2, unsigned len)
{
	BNWORD32 t;

	assert(len > 0);	/* Alternative: change loop to test at start */

	t = (BNWORD32)BIGLITTLE(*--num1,*num1) + BIGLITTLE(*--num2,*num2++);
	BIGLITTLE(*num1,*num1++) = (BNWORD16)t;
	while (--len) {
		t = (BNWORD32)BIGLITTLE(*--num1,*num1) +
		    (BNWORD32)BIGLITTLE(*--num2,*num2++) + (t >> 16);
		BIGLITTLE(*num1,*num1++) = (BNWORD16)t;
	}

	return (BNWORD16)(t>>16);
}
#else /* no BNWORD32 */
BNWORD16
lbnAddN_16(BNWORD16 *num1, BNWORD16 const *num2, unsigned len)
{
	BNWORD16 x, carry = 0;

	assert(len > 0);	/* Alternative: change loop to test at start */

	do {
		x = BIGLITTLE(*--num2,*num2++);
		carry = (x += carry) < carry;
		carry += (BIGLITTLE(*--num1,*num1++) += x) < x;
	} while (--len);

	return carry;
}
#endif
#endif /* !lbnAddN_16 */

/*
 * lbnSubN_16: add two bignums of the same length, returning the carry (0 or 1).
 * One of the building blocks, along with subn1, of subtracting two bignums of
 * differing lengths.
 *
 * Technique: If no double-width type is availble, maintain a word of borrow.
 * First, add the borrow to the subtrahend (did you have to learn all those
 * awful words in elementary school, too?), and if it overflows, set the
 * borrow again.  Then subtract the modified subtrahend from the next word
 * of input, using the same technique as in subn1, above.
 * Adding the borrows is used as an OR operator; at most one of the two
 * comparisons can possibly be true.  The first can only be true if
 * borrow == 1 and x, the result, is 0.  In that case the second can't
 * possibly be true.
 *
 * In the double-word case, (BNWORD16)-(t>>16) is subtracted, rather than
 * adding t>>16, because the shift would need to sign-extend and that's
 * not guaranteed to happen in ANSI C, even with signed types.
 */
#ifndef lbnSubN_16
#ifdef BNWORD32
BNWORD16
lbnSubN_16(BNWORD16 *num1, BNWORD16 const *num2, unsigned len)
{
	BNWORD32 t;

	assert(len > 0);	/* Alternative: change loop to test at start */

	t = (BNWORD32)BIGLITTLE(*--num1,*num1) - BIGLITTLE(*--num2,*num2++);
	BIGLITTLE(*num1,*num1++) = (BNWORD16)t;

	while (--len) {
		t = (BNWORD32)BIGLITTLE(*--num1,*num1) -
		    (BNWORD32)BIGLITTLE(*--num2,*num2++) - (BNWORD16)-(t >> 16);
		BIGLITTLE(*num1,*num1++) = (BNWORD16)t;
	}

	return -(BNWORD16)(t>>16);
}
#else
BNWORD16
lbnSubN_16(BNWORD16 *num1, BNWORD16 const *num2, unsigned len)
{
	BNWORD16 x, borrow = 0;

	assert(len > 0);	/* Alternative: change loop to test at start */

	do {
		x = BIGLITTLE(*--num2,*num2++);
		borrow = (x += borrow) < borrow;
		borrow += (BIGLITTLE(*--num1,*num1++) -= x) > (BNWORD16)~x;
	} while (--len);

	return borrow;
}
#endif
#endif /* !lbnSubN_16 */

#ifndef lbnCmp_16
/*
 * lbnCmp_16: compare two bignums of equal length, returning the sign of
 * num1 - num2. (-1, 0 or +1).
 * 
 * Technique: Change the little-endian pointers to big-endian pointers
 * and compare from the most-significant end until a difference if found.
 * When it is, figure out the sign of the difference and return it.
 */
int
lbnCmp_16(BNWORD16 const *num1, BNWORD16 const *num2, unsigned len)
{
	BIGLITTLE(num1 -= len, num1 += len);
	BIGLITTLE(num2 -= len, num2 += len);

	while (len--) {
		if (BIGLITTLE(*num1++ != *num2++, *--num1 != *--num2)) {
			if (BIGLITTLE(num1[-1] < num2[-1], *num1 < *num2))
				return -1;
			else
				return 1;
		}
	}
	return 0;
}
#endif /* !lbnCmp_16 */

/*
 * mul16_ppmmaa(ph,pl,x,y,a,b) is an optional routine that
 * computes (ph,pl) = x * y + a + b.  mul16_ppmma and mul16_ppmm
 * are simpler versions.  If you want to be lazy, all of these
 * can be defined in terms of the others, so here we create any
 * that have not been defined in terms of the ones that have been.
 */

/* Define ones with fewer a's in terms of ones with more a's */
#if !defined(mul16_ppmma) && defined(mul16_ppmmaa)
#define mul16_ppmma(ph,pl,x,y,a) mul16_ppmmaa(ph,pl,x,y,a,0)
#endif