OpenSolaris Launch

1 files changed, 2620 insertions, 0 deletions

diff --git a/usr/src/common/bignum/bignumimpl.c b/usr/src/common/bignum/bignumimpl.c
new file mode 100644
index 0000000000..cf42c030ee
--- /dev/null
+++ b/usr/src/common/bignum/bignumimpl.c
@@ -0,0 +1,2620 @@
+/*
+ * CDDL HEADER START
+ *
+ * The contents of this file are subject to the terms of the
+ * Common Development and Distribution License, Version 1.0 only
+ * (the "License").  You may not use this file except in compliance
+ * with the License.
+ *
+ * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
+ * or http://www.opensolaris.org/os/licensing.
+ * See the License for the specific language governing permissions
+ * and limitations under the License.
+ *
+ * When distributing Covered Code, include this CDDL HEADER in each
+ * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
+ * If applicable, add the following below this CDDL HEADER, with the
+ * fields enclosed by brackets "[]" replaced with your own identifying
+ * information: Portions Copyright [yyyy] [name of copyright owner]
+ *
+ * CDDL HEADER END
+ */
+/*
+ * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
+ * Use is subject to license terms.
+ */
+
+#pragma ident	"%Z%%M%	%I%	%E% SMI"
+
+#define	big_div_pos_fast big_div_pos
+
+#include "bignum.h"
+
+/*
+ * Configuration guide
+ * -------------------
+ *
+ * There are 4 preprocessor symbols used to configure the bignum
+ * implementation.  This file contains no logic to configure based on
+ * processor; we leave that to the Makefiles to specify.
+ *
+ * USE_FLOATING_POINT
+ *   Meaning: There is support for a fast floating-point implementation of
+ *   Montgomery multiply.
+ *
+ * PSR_MUL
+ *   Meaning: There are processor-specific versions of the low level
+ *   functions to implement big_mul.  Those functions are: big_mul_set_vec,
+ *   big_mul_add_vec, big_mul_vec, and big_sqr_vec.  PSR_MUL implies support
+ *   for all 4 functions.  You cannot pick and choose which subset of these
+ *   functions to support; that would lead to a rat's nest of #ifdefs.
+ *
+ * HWCAP
+ *   Meaning: Call multiply support functions through a function pointer.
+ *   On x86, there are multiple implementations for differnt hardware
+ *   capabilities, such as MMX, SSE2, etc.  Tests are made at run-time, when
+ *   a function is first used.  So, the support functions are called through
+ *   a function pointer.  There is no need for that on Sparc, because there
+ *   is only one implementation; support functions are called directly.
+ *   Later, if there were some new VIS instruction, or something, and a
+ *   run-time test were needed, rather than variant kernel modules and
+ *   libraries, then HWCAP would be defined for Sparc, as well.
+ *
+ * UMUL64
+ *   Meaning: It is safe to use generic C code that assumes the existence
+ *   of a 32 x 32 --> 64 bit unsigned multiply.  If this is not defined,
+ *   then the generic code for big_mul_add_vec() must necessarily be very slow,
+ *   because it must fall back to using 16 x 16 --> 32 bit multiplication.
+ *
+ */
+
+
+#ifdef	_KERNEL
+
+#include <sys/types.h>
+#include <sys/kmem.h>
+#include <sys/param.h>
+#include <sys/sunddi.h>
+
+#define	big_malloc(size)	kmem_alloc(size, KM_NOSLEEP)
+#define	big_free(ptr, size)	kmem_free(ptr, size)
+
+void *
+big_realloc(void *from, size_t oldsize, size_t newsize)
+{
+	void *rv;
+
+	rv = kmem_alloc(newsize, KM_NOSLEEP);
+	if (rv != NULL)
+		bcopy(from, rv, oldsize);
+	kmem_free(from, oldsize);
+	return (rv);
+}
+
+#else	/* _KERNEL */
+
+#include <stdlib.h>
+#include <stdio.h>
+
+#ifndef MALLOC_DEBUG
+
+#define	big_malloc(size)	malloc(size)
+#define	big_free(ptr, size)	free(ptr)
+
+#else
+
+void
+big_free(void *ptr, size_t size)
+{
+	printf("freed %d bytes at %p\n", size, ptr);
+	free(ptr);
+}
+
+void *
+big_malloc(size_t size)
+{
+	void *rv;
+	rv = malloc(size);
+	printf("malloced %d bytes, addr:%p\n", size, rv);
+	return (rv);
+}
+#endif /* MALLOC_DEBUG */
+
+#define	big_realloc(x, y, z) realloc((x), (z))
+
+void
+printbignum(char *aname, BIGNUM *a)
+{
+	int i;
+
+	(void) printf("\n%s\n%d\n", aname, a->sign*a->len);
+	for (i = a->len - 1; i >= 0; i--) {
+		(void) printf("%08x ", a->value[i]);
+		if ((i % 8 == 0) && (i != 0))
+		    (void) printf("\n");
+	}
+	(void) printf("\n");
+}
+
+#endif	/* _KERNEL */
+
+
+/* size in 32-bit words */
+BIG_ERR_CODE
+big_init(BIGNUM *number, int size)
+{
+	number->value = big_malloc(sizeof (uint32_t) * size);
+	if (number->value == NULL) {
+		return (BIG_NO_MEM);
+	}
+	number->size = size;
+	number->len = 0;
+	number->sign = 1;
+	number->malloced = 1;
+	return (BIG_OK);
+}
+
+/* size in 32-bit words */
+BIG_ERR_CODE
+big_init1(BIGNUM *number, int size, uint32_t *buf, int bufsize)
+{
+	if ((buf == NULL) || (size > bufsize)) {
+		number->value = big_malloc(sizeof (uint32_t) * size);
+		if (number->value == NULL) {
+			return (BIG_NO_MEM);
+		}
+		number->size = size;
+		number->malloced = 1;
+	} else {
+		number->value = buf;
+		number->size = bufsize;
+		number->malloced = 0;
+	}
+		number->len = 0;
+		number->sign = 1;
+
+	return (BIG_OK);
+}
+
+void
+big_finish(BIGNUM *number)
+{
+	if (number->malloced == 1) {
+		big_free(number->value, sizeof (uint32_t) * number->size);
+		number->malloced = 0;
+	}
+}
+
+/*
+ *  bn->size should be at least (len + 3) / 4
+ * converts from byte-big-endian format to bignum format (words in
+ * little endian order, but bytes within the words big endian)
+ */
+void
+bytestring2bignum(BIGNUM *bn, uchar_t *kn, size_t len)
+{
+	int		i, j, offs;
+	uint32_t	word;
+	uchar_t		*knwordp;
+
+#ifdef	_LP64
+	/* LINTED */
+	offs = (uint32_t)len % sizeof (uint32_t);
+	/* LINTED */
+	bn->len = (uint32_t)len / sizeof (uint32_t);
+	/* LINTED */
+	for (i = 0; i < (uint32_t)len / sizeof (uint32_t); i++) {
+#else	/* !_LP64 */
+	offs = len % sizeof (uint32_t);
+	bn->len = len / sizeof (uint32_t);
+	for (i = 0; i < len / sizeof (uint32_t); i++) {
+#endif	/* _LP64 */
+		knwordp = &(kn[len - sizeof (uint32_t) * (i + 1)]);
+		word = knwordp[0];
+		for (j = 1; j < sizeof (uint32_t); j++) {
+			word = (word << 8)+ knwordp[j];
+		}
+		bn->value[i] = word;
+	}
+	if (offs > 0) {
+		word = kn[0];
+		for (i = 1; i < offs; i++) word = (word << 8) + kn[i];
+		bn->value[bn->len++] = word;
+	}
+	while ((bn->len > 1) && (bn->value[bn->len-1] == 0)) {
+		bn->len --;
+	}
+}
+
+/*
+ * copies the least significant len bytes if
+ * len < bn->len * sizeof (uint32_t)
+ * converts from bignum format to byte-big-endian format.
+ * bignum format is words in little endian order,
+ * but bytes within words in native byte order.
+ */
+void
+bignum2bytestring(uchar_t *kn, BIGNUM *bn, size_t len)
+{
+	int		i, j, offs;
+	uint32_t	word;
+
+	if (len < sizeof (uint32_t) * bn->len) {
+#ifdef	_LP64
+		/* LINTED */
+		for (i = 0; i < (uint32_t)len / sizeof (uint32_t); i++) {
+#else	/* !_LP64 */
+		for (i = 0; i < len / sizeof (uint32_t); i++) {
+#endif	/* _LP64 */
+			word = bn->value[i];
+			for (j = 0; j < sizeof (uint32_t); j++) {
+				kn[len - sizeof (uint32_t) * i - j - 1] =
+				    word & 0xff;
+				word = word >> 8;
+			}
+		}
+#ifdef	_LP64
+		/* LINTED */
+		offs = (uint32_t)len % sizeof (uint32_t);
+#else	/* !_LP64 */
+		offs = len % sizeof (uint32_t);
+#endif	/* _LP64 */
+		if (offs > 0) {
+			word = bn->value[len / sizeof (uint32_t)];
+#ifdef	_LP64
+			    /* LINTED */
+			    for (i =  (uint32_t)len % sizeof (uint32_t);
+				i > 0; i --) {
+#else	/* !_LP64 */
+			    for (i = len % sizeof (uint32_t); i > 0; i --) {
+#endif	/* _LP64 */
+				    kn[i - 1] = word & 0xff;
+				    word = word >> 8;
+			    }
+		}
+	} else {
+		for (i = 0; i < bn->len; i++) {
+			word = bn->value[i];
+			for (j = 0; j < sizeof (uint32_t); j++) {
+				kn[len - sizeof (uint32_t) * i - j - 1] =
+				    word & 0xff;
+				word = word >> 8;
+			}
+		}
+#ifdef	_LP64
+		/* LINTED */
+		for (i = 0; i < (uint32_t)len - sizeof (uint32_t) * bn->len;
+		    i++) {
+#else	/* !_LP64 */
+		for (i = 0; i < len - sizeof (uint32_t) * bn->len; i++) {
+#endif	/* _LP64 */
+			kn[i] = 0;
+		}
+	}
+}
+
+
+int
+big_bitlength(BIGNUM *a)
+{
+	int		l, b;
+	uint32_t	c;
+
+	l = a->len - 1;
+	while ((l > 0) && (a->value[l] == 0)) {
+		l--;
+	}
+	b = sizeof (uint32_t) * 8;
+	c = a->value[l];
+	while ((b > 1) && ((c & 0x80000000) == 0)) {
+		c = c << 1;
+		b--;
+	}
+	return (l * (int)sizeof (uint32_t) * 8 + b);
+}
+
+
+BIG_ERR_CODE
+big_copy(BIGNUM *dest, BIGNUM *src)
+{
+	uint32_t *newptr;
+	int i, len;
+
+	len = src->len;
+	while ((len > 1) && (src->value[len - 1] == 0))
+		len--;
+	src->len = len;
+	if (dest->size < len) {
+		if (dest->malloced == 1) {
+			newptr = (uint32_t *)big_realloc(dest->value,
+			    sizeof (uint32_t) * dest->size,
+			    sizeof (uint32_t) * len);
+		} else {
+			newptr = (uint32_t *)
+			    big_malloc(sizeof (uint32_t) * len);
+			if (newptr != NULL) dest->malloced = 1;
+		}
+		if (newptr == NULL)
+			return (BIG_NO_MEM);
+		dest->value = newptr;
+		dest->size = len;
+	}
+	dest->len = len;
+	dest->sign = src->sign;
+	for (i = 0; i < len; i++) dest->value[i] = src->value[i];
+
+	return (BIG_OK);
+}
+
+
+BIG_ERR_CODE
+big_extend(BIGNUM *number, int size)
+{
+	uint32_t	*newptr;
+	int		i;
+
+	if (number->size >= size)
+		return (BIG_OK);
+	if (number->malloced) {
+		number->value =
+		    big_realloc(number->value,
+			sizeof (uint32_t) * number->size,
+			sizeof (uint32_t) * size);
+	} else {
+		newptr = big_malloc(sizeof (uint32_t) * size);
+		if (newptr != NULL) {
+			for (i = 0; i < number->size; i++) {
+				newptr[i] = number->value[i];
+			}
+		}
+		number->value = newptr;
+	}
+
+	if (number->value == NULL)
+		return (BIG_NO_MEM);
+
+	number->size = size;
+	number->malloced = 1;
+	return (BIG_OK);
+}
+
+
+int
+big_is_zero(BIGNUM *n)
+{
+	int i, result;
+
+	result = 1;
+	for (i = 0; i < n->len; i++)
+		if (n->value[i] != 0) result = 0;
+	return (result);
+}
+
+
+BIG_ERR_CODE
+big_add_abs(BIGNUM *result, BIGNUM *aa, BIGNUM *bb)
+{
+	int i, shorter, longer;
+	uint32_t cy, ai;
+	uint32_t *r, *a, *b, *c;
+	BIG_ERR_CODE err;
+
+	if (aa->len > bb->len) {
+		shorter = bb->len;
+		longer = aa->len;
+		c = aa->value;
+	} else {
+		shorter = aa->len;
+		longer = bb->len;
+		c = bb->value;
+	}
+	if (result->size < longer + 1) {
+		err = big_extend(result, longer + 1);
+		if (err != BIG_OK)
+			return (err);
+	}
+
+	r = result->value;
+	a = aa->value;
+	b = bb->value;
+	cy = 0;
+	for (i = 0; i < shorter; i++) {
+		ai = a[i];
+		r[i] = ai + b[i] + cy;
+		if (r[i] > ai) cy = 0;
+		else if (r[i] < ai) cy = 1;
+	}
+	for (; i < longer; i++) {
+		ai = c[i];
+		r[i] = ai + cy;
+		if (r[i] >= ai) cy = 0;
+	}
+	if (cy == 1) {
+		r[i] = cy;
+		result->len = longer + 1;
+	} else {
+		result->len = longer;
+	}
+	result->sign = 1;
+	return (BIG_OK);
+}
+
+
+/* caller must make sure that result has at least len words allocated */
+void
+big_sub_vec(uint32_t *r, uint32_t *a, uint32_t *b, int len)
+{
+	int i;
+	uint32_t cy, ai;
+
+	cy = 1;
+	for (i = 0; i < len; i++) {
+		ai = a[i];
+		r[i] = ai + (~b[i]) + cy;
+		if (r[i] > ai) cy = 0;
+		else if (r[i] < ai) cy = 1;
+	}
+}
+
+
+/* result=aa-bb  it is assumed that aa>=bb */
+BIG_ERR_CODE
+big_sub_pos(BIGNUM *result, BIGNUM *aa, BIGNUM *bb)
+{
+	int i, shorter;
+	uint32_t cy, ai;
+	uint32_t *r, *a, *b;
+	BIG_ERR_CODE err;
+
+	if (aa->len > bb->len) shorter = bb->len;
+	else shorter = aa->len;
+	if (result->size < aa->len) {
+		err = big_extend(result, aa->len);
+		if (err != BIG_OK)
+			return (err);
+	}
+
+	r = result->value;
+	a = aa->value;
+	b = bb->value;
+	result->len = aa->len;
+	cy = 1;
+	for (i = 0; i < shorter; i++) {
+		ai = a[i];
+		r[i] = ai + (~b[i]) + cy;
+		if (r[i] > ai) cy = 0;
+		else if (r[i] < ai) cy = 1;
+	}
+	for (; i < aa->len; i++) {
+		ai = a[i];
+		r[i] = ai + (~0) + cy;
+		if (r[i] < ai) cy = 1;
+	}
+	result->sign = 1;
+	if (cy == 0)
+		return (BIG_INVALID_ARGS);
+	else
+		return (BIG_OK);
+}
+
+
+/* returns -1 if |aa|<|bb|, 0 if |aa|==|bb| 1 if |aa|>|bb| */
+int
+big_cmp_abs(BIGNUM *aa, BIGNUM *bb)
+{
+	int i;
+
+	if (aa->len > bb->len) {
+		for (i = aa->len - 1; i > bb->len - 1; i--) {
+			if (aa->value[i] > 0)
+				return (1);
+		}
+	} else if (aa->len < bb->len) {
+		for (i = bb->len - 1; i > aa->len - 1; i--) {
+			if (bb->value[i] > 0)
+				return (-1);
+		}
+	} else i = aa->len-1;
+	for (; i >= 0; i--) {
+		if (aa->value[i] > bb->value[i])
+			return (1);
+		else if (aa->value[i] < bb->value[i])
+			return (-1);
+	}
+
+	return (0);
+}
+
+
+BIG_ERR_CODE
+big_sub(BIGNUM *result, BIGNUM *aa, BIGNUM *bb)
+{
+	BIG_ERR_CODE err;
+
+	if ((bb->sign == -1) && (aa->sign == 1)) {
+		if ((err = big_add_abs(result, aa, bb)) != BIG_OK)
+			return (err);
+		result->sign = 1;
+	} else if ((aa->sign == -1) && (bb->sign == 1)) {
+		if ((err = big_add_abs(result, aa, bb)) != BIG_OK)
+			return (err);
+		result->sign = -1;
+	} else if ((aa->sign == 1) && (bb->sign == 1)) {
+		if (big_cmp_abs(aa, bb) >= 0) {
+			if ((err = big_sub_pos(result, aa, bb)) != BIG_OK)
+				return (err);
+			result->sign = 1;
+		} else {
+			if ((err = big_sub_pos(result, bb, aa)) != BIG_OK)
+				return (err);
+			result->sign = -1;
+		}
+	} else {
+		if (big_cmp_abs(aa, bb) >= 0) {
+			if ((err = big_sub_pos(result, aa, bb)) != BIG_OK)
+				return (err);
+			result->sign = -1;
+		} else {
+			if ((err = big_sub_pos(result, bb, aa)) != BIG_OK)
+				return (err);
+			result->sign = 1;
+		}
+	}
+	return (BIG_OK);
+}
+
+
+
+BIG_ERR_CODE
+big_add(BIGNUM *result, BIGNUM *aa, BIGNUM *bb)
+{
+	BIG_ERR_CODE err;
+
+	if ((bb->sign == -1) && (aa->sign == -1)) {
+		if ((err = big_add_abs(result, aa, bb)) != BIG_OK)
+			return (err);
+		result->sign = -1;
+	} else if ((aa->sign == 1) && (bb->sign == 1)) {
+		if ((err = big_add_abs(result, aa, bb)) != BIG_OK)
+			return (err);
+		result->sign = 1;
+	} else if ((aa->sign == 1) && (bb->sign == -1)) {
+		if (big_cmp_abs(aa, bb) >= 0) {
+			if ((err = big_sub_pos(result, aa, bb)) != BIG_OK)
+				return (err);
+			result->sign = 1;
+		} else {
+			if ((err = big_sub_pos(result, bb, aa)) != BIG_OK)
+				return (err);
+			result->sign = -1;
+		}
+	} else {
+		if (big_cmp_abs(aa, bb) >= 0) {
+			if ((err = big_sub_pos(result, aa, bb)) != BIG_OK)
+				return (err);
+			result->sign = -1;
+		} else {
+			if ((err = big_sub_pos(result, bb, aa)) != BIG_OK)
+				return (err);
+			result->sign = 1;
+		}
+	}
+	return (BIG_OK);
+}
+
+
+/* result = aa/2 aa must be positive */
+BIG_ERR_CODE
+big_half_pos(BIGNUM *result, BIGNUM *aa)
+{
+	BIG_ERR_CODE err;
+	int i;
+	uint32_t cy, cy1;
+	uint32_t *a, *r;
+
+	if (result->size < aa->len) {
+		err = big_extend(result, aa->len);
+		if (err != BIG_OK)
+			return (err);
+	}
+
+	result->len = aa->len;
+	a = aa->value;
+	r = result->value;
+	cy = 0;
+	for (i = aa->len-1; i >= 0; i--) {
+		cy1 = a[i] << 31;
+		r[i] = (cy|(a[i] >> 1));
+		cy = cy1;
+	}
+	if (r[result->len-1] == 0) result->len--;
+	return (BIG_OK);
+}
+
+/* result  =  aa*2 aa must be positive */
+BIG_ERR_CODE
+big_double(BIGNUM *result, BIGNUM *aa)
+{
+	BIG_ERR_CODE err;
+	int i, rsize;
+	uint32_t cy, cy1;
+	uint32_t *a, *r;
+
+	if ((aa->len > 0) && ((aa->value[aa->len - 1] & 0x80000000) != 0))
+		rsize = aa->len + 1;
+	else rsize = aa->len;
+
+	if (result->size < rsize) {
+		err = big_extend(result, rsize);
+		if (err != BIG_OK)
+			return (err);
+	}
+
+	a = aa->value;
+	r = result->value;
+	if (rsize == aa->len + 1) r[rsize - 1] = 1;
+	cy = 0;
+	for (i = 0; i < aa->len; i++) {
+		cy1 = a[i] >> 31;
+		r[i] = (cy | (a[i] << 1));
+		cy = cy1;
+	}
+	result->len = rsize;
+	return (BIG_OK);
+}
+
+/* returns aa mod b, aa must be nonneg, b must be a max 16-bit integer */
+uint32_t
+big_mod16_pos(BIGNUM *aa, uint32_t b)
+{
+	int i;
+	uint32_t rem;
+
+	if (aa->len == 0)
+		return (0);
+	rem = aa->value[aa->len - 1] % b;
+	for (i = aa->len - 2; i >= 0; i--) {
+		rem = ((rem << 16) | (aa->value[i] >> 16)) % b;
+		rem = ((rem << 16) | (aa->value[i] & 0xffff)) % b;
+	}
+	return (rem);
+}
+
+
+/*
+ * result = aa - (2^32)^lendiff * bb
+ * result->size should be at least aa->len at entry
+ * aa, bb, and result should be positive
+ */
+void
+big_sub_pos_high(BIGNUM *result, BIGNUM *aa, BIGNUM *bb)
+{
+	int i, lendiff;
+	BIGNUM res1, aa1;
+
+	lendiff = aa->len - bb->len;
+	res1.size = result->size - lendiff;
+	res1.malloced = 0;
+	res1.value = result->value + lendiff;
+	aa1.size = aa->size - lendiff;
+	aa1.value = aa->value + lendiff;
+	aa1.len = bb->len;
+	aa1.sign = 1;
+	(void) big_sub_pos(&res1, &aa1, bb);
+	if (result->value != aa->value) {
+		for (i = 0; i < lendiff; i++) {
+			result->value[i] = aa->value[i];
+		}
+	}
+	result->len = aa->len;
+}
+
+
+/*
+ * returns 1, 0, or -1 depending on whether |aa| > , ==, or <
+ *							(2^32)^lendiff * |bb|
+ * aa->len should be >= bb->len
+ */
+int
+big_cmp_abs_high(BIGNUM *aa, BIGNUM *bb)
+{
+	int lendiff;
+	BIGNUM aa1;
+
+	lendiff = aa->len - bb->len;
+	aa1.len = bb->len;
+	aa1.size = aa->size - lendiff;
+	aa1.malloced = 0;
+	aa1.value = aa->value + lendiff;
+	return (big_cmp_abs(&aa1, bb));
+}
+
+
+/*
+ * result = aa * b where b is a max. 16-bit positive integer.
+ * result should have enough space allocated.
+ */
+void
+big_mul16_low(BIGNUM *result, BIGNUM *aa, uint32_t b)
+{
+	int i;
+	uint32_t t1, t2, ai, cy;
+	uint32_t *a, *r;
+
+	a = aa->value;
+	r = result->value;
+	cy = 0;
+	for (i = 0; i < aa->len; i++) {
+		ai = a[i];
+		t1 = (ai & 0xffff) * b + cy;
+		t2 = (ai >> 16) * b + (t1 >> 16);
+		r[i] = (t1 & 0xffff) | (t2 << 16);
+		cy = t2 >> 16;
+	}
+	r[i] = cy;
+	result->len = aa->len + 1;
+	result->sign = aa->sign;
+}
+
+
+/*
+ * result = aa * b * 2^16 where b is a max. 16-bit positive integer.
+ * result should have enough space allocated.
+ */
+void
+big_mul16_high(BIGNUM *result, BIGNUM *aa, uint32_t b)
+{
+	int i;
+	uint32_t t1, t2, ai, cy, ri;
+	uint32_t *a, *r;
+
+	a = aa->value;
+	r = result->value;
+	cy = 0;
+	ri = 0;
+	for (i = 0; i < aa->len; i++) {
+		ai = a[i];
+		t1 = (ai & 0xffff) * b + cy;
+		t2 = (ai >> 16) * b + (t1 >> 16);
+		r[i] = (t1 << 16) + ri;
+		ri = t2 & 0xffff;
+		cy = t2 >> 16;
+	}
+	r[i] = (cy << 16) + ri;
+	result->len = aa->len + 1;
+	result->sign = aa->sign;
+}
+
+/* it is assumed that result->size is big enough */
+void
+big_shiftleft(BIGNUM *result, BIGNUM *aa, int offs)
+{
+	int i;
+	uint32_t cy, ai;
+
+	if (offs == 0) {
+		if (result != aa) {
+			(void) big_copy(result, aa);
+		}
+		return;
+	}
+	cy = 0;
+	for (i = 0; i < aa->len; i++) {
+		ai = aa->value[i];
+		result->value[i] = (ai << offs) | cy;
+		cy = ai >> (32 - offs);
+	}
+	if (cy != 0) {
+		result->len = aa->len + 1;
+		result->value[result->len - 1] = cy;
+	} else {
+		result->len = aa->len;
+	}
+	result->sign = aa->sign;
+}
+
+/* it is assumed that result->size is big enough */
+void
+big_shiftright(BIGNUM *result, BIGNUM *aa, int offs)
+{
+	int i;
+	uint32_t cy, ai;
+
+	if (offs == 0) {
+		if (result != aa) {
+			(void) big_copy(result, aa);
+		}
+		return;
+	}
+	cy = aa->value[0] >> offs;
+	for (i = 1; i < aa->len; i++) {
+		ai = aa->value[i];
+		result->value[i-1] = (ai << (32 - offs)) | cy;
+		cy = ai >> offs;
+	}
+	result->len = aa->len;
+	result->value[result->len - 1] = cy;
+	result->sign = aa->sign;
+}
+
+
+/*
+ * result = aa/bb   remainder = aa mod bb
+ * it is assumed that aa and bb are positive
+ */
+BIG_ERR_CODE
+big_div_pos_fast(BIGNUM *result, BIGNUM *remainder, BIGNUM *aa, BIGNUM *bb)
+{
+	BIG_ERR_CODE err;
+	int i, alen, blen, tlen, rlen, offs;
+	uint32_t higha, highb, coeff;
+	uint64_t highb64;
+	uint32_t *a, *b;
+	BIGNUM bbhigh, bblow, tresult, tmp1, tmp2;
+	uint32_t tmp1value[BIGTMPSIZE];
+	uint32_t tmp2value[BIGTMPSIZE];
+	uint32_t tresultvalue[BIGTMPSIZE];
+	uint32_t bblowvalue[BIGTMPSIZE];
+	uint32_t bbhighvalue[BIGTMPSIZE];
+
+	a = aa->value;
+	b = bb->value;
+	alen = aa->len;
+	blen = bb->len;
+	while ((alen > 1) && (a[alen - 1] == 0)) alen = alen - 1;
+	aa->len = alen;
+	while ((blen > 1) && (b[blen - 1] == 0)) blen = blen - 1;
+	bb->len = blen;
+	if ((blen == 1) && (b[0] == 0))
+		return (BIG_DIV_BY_0);
+
+	if (big_cmp_abs(aa, bb) < 0) {
+		if ((remainder != NULL) &&
+		    ((err = big_copy(remainder, aa)) != BIG_OK))
+			return (err);
+		if (result != NULL) {
+			result->len = 1;
+			result->sign = 1;
+			result->value[0] = 0;
+		}
+		return (BIG_OK);
+	}
+
+	if ((err = big_init1(&bblow, blen + 1,
+	    bblowvalue, arraysize(bblowvalue))) != BIG_OK)
+		return (err);
+
+	if ((err = big_init1(&bbhigh, blen + 1,
+	    bbhighvalue, arraysize(bbhighvalue))) != BIG_OK)
+		goto ret1;
+
+	if ((err = big_init1(&tmp1, alen + 2,
+	    tmp1value, arraysize(tmp1value))) != BIG_OK)
+		goto ret2;
+
+	if ((err = big_init1(&tmp2, blen + 2,
+	    tmp2value, arraysize(tmp2value))) != BIG_OK)
+		goto ret3;
+
+	if ((err = big_init1(&tresult, alen - blen + 2,
+	    tresultvalue, arraysize(tresultvalue))) != BIG_OK)
+		goto ret4;
+
+	offs = 0;
+	if (blen > 1) {
+		highb64 = (((uint64_t)(b[blen - 1])) << 32) |
+		    ((uint64_t)(b[blen - 2]));
+	} else {
+		highb64 = (((uint64_t)(b[blen - 1])) << 32);
+	}
+	if (highb64 >= 0x1000000000000ull) {
+		highb64 = highb64 >> 16;
+		offs = 16;
+	}
+	while ((highb64 & 0x800000000000ull) == 0) {
+		highb64 = highb64 << 1;
+		offs++;
+	}
+#ifdef	_LP64
+	/* LINTED */
+	highb = (highb64 >> 32) & 0xffffffff;
+#else	/* !_LP64 */
+	highb = highb64 >> 32;
+#endif	/* _LP64 */
+
+	big_shiftleft(&bblow, bb, offs);
+	if (offs <= 15) {
+		big_shiftleft(&bbhigh, &bblow, 16);
+	} else {
+		big_shiftright(&bbhigh, &bblow, 16);
+	}
+	if (bbhigh.value[bbhigh.len - 1] == 0) {
+		bbhigh.len--;
+	} else {
+		bbhigh.value[bbhigh.len] = 0;
+	}
+
+	big_shiftleft(&tmp1, aa, offs);
+	rlen = tmp1.len - bblow.len + 1;
+	tresult.len = rlen;
+
+	tmp1.len++;
+	tlen = tmp1.len;
+	tmp1.value[tmp1.len - 1] = 0;
+	for (i = 0; i < rlen; i++) {
+		higha = (tmp1.value[tlen - 1] << 16) +
+		    (tmp1.value[tlen - 2] >> 16);
+		coeff = higha / (highb + 1);
+		big_mul16_high(&tmp2, &bblow, coeff);
+		big_sub_pos_high(&tmp1, &tmp1, &tmp2);
+		bbhigh.len++;
+		while (tmp1.value[tlen - 1] > 0) {
+			big_sub_pos_high(&tmp1, &tmp1, &bbhigh);
+			coeff++;
+		}
+		bbhigh.len--;
+		tlen--;
+		tmp1.len--;
+		while (big_cmp_abs_high(&tmp1, &bbhigh) >= 0) {
+			big_sub_pos_high(&tmp1, &tmp1, &bbhigh);
+			coeff++;
+		}
+		tresult.value[rlen - i - 1] = coeff << 16;
+		higha = tmp1.value[tlen - 1];
+		coeff = higha / (highb + 1);
+		big_mul16_low(&tmp2, &bblow, coeff);
+		tmp2.len--;
+		big_sub_pos_high(&tmp1, &tmp1, &tmp2);
+		while (big_cmp_abs_high(&tmp1, &bblow) >= 0) {
+			big_sub_pos_high(&tmp1, &tmp1, &bblow);
+			coeff++;
+		}
+		tresult.value[rlen - i - 1] =
+		    tresult.value[rlen - i - 1] + coeff;
+	}
+
+	big_shiftright(&tmp1, &tmp1, offs);
+
+	err = BIG_OK;
+
+	if ((remainder != NULL) &&
+	    ((err = big_copy(remainder, &tmp1)) != BIG_OK))
+		goto ret;
+
+	if (result != NULL)
+		err = big_copy(result, &tresult);
+
+ret:
+	big_finish(&tresult);
+ret4:
+	big_finish(&tmp1);
+ret3:
+	big_finish(&tmp2);
+ret2:
+	big_finish(&bbhigh);
+ret1:
+	big_finish(&bblow);
+	return (err);
+}
+
+/*
+ * If there is no processor-specific integer implementation of
+ * the lower level multiply functions, then this code is provided
+ * for big_mul_set_vec(), big_mul_add_vec(), big_mul_vec() and
+ * big_sqr_vec().
+ *
+ * There are two generic implementations.  One that assumes that
+ * there is hardware and C compiler support for a 32 x 32 --> 64
+ * bit unsigned multiply, but otherwise is not specific to any
+ * processor, platform, or ISA.
+ *
+ * The other makes very few assumptions about hardware capabilities.
+ * It does not even assume that there is any implementation of a
+ * 32 x 32 --> 64 bit multiply that is accessible to C code and
+ * appropriate to use.  It falls constructs 32 x 32 --> 64 bit
+ * multiplies from 16 x 16 --> 32 bit multiplies.
+ *
+ */
+
+#if !defined(PSR_MUL)
+
+#ifdef UMUL64
+
+#define	UNROLL8
+
+#define	MUL_SET_VEC_ROUND_PREFETCH(R) \
+	p = pf * d; \
+	pf = (uint64_t)a[R+1]; \
+	t = p + cy; \
+	r[R] = (uint32_t)t; \
+	cy = t >> 32
+
+#define	MUL_SET_VEC_ROUND_NOPREFETCH(R) \
+	p = pf * d; \
+	t = p + cy; \
+	r[R] = (uint32_t)t; \
+	cy = t >> 32
+
+#define	MUL_ADD_VEC_ROUND_PREFETCH(R) \
+	t = (uint64_t)r[R]; \
+	p = pf * d; \
+	pf = (uint64_t)a[R+1]; \
+	t = p + t + cy; \
+	r[R] = (uint32_t)t; \
+	cy = t >> 32
+
+#define	MUL_ADD_VEC_ROUND_NOPREFETCH(R) \
+	t = (uint64_t)r[R]; \
+	p = pf * d; \
+	t = p + t + cy; \
+	r[R] = (uint32_t)t; \
+	cy = t >> 32
+
+#ifdef UNROLL8
+
+#define	UNROLL 8
+
+/*
+ * r = a * b
+ * where r and a are vectors; b is a single 32-bit digit
+ */
+
+uint32_t
+big_mul_set_vec(uint32_t *r, uint32_t *a, int len, uint32_t b)
+{
+	uint64_t d, pf, p, t, cy;
+
+	if (len == 0)
+		return (0);
+	cy = 0;
+	d = (uint64_t)b;
+	pf = (uint64_t)a[0];
+	while (len > UNROLL) {
+		MUL_SET_VEC_ROUND_PREFETCH(0);
+		MUL_SET_VEC_ROUND_PREFETCH(1);
+		MUL_SET_VEC_ROUND_PREFETCH(2);
+		MUL_SET_VEC_ROUND_PREFETCH(3);
+		MUL_SET_VEC_ROUND_PREFETCH(4);
+		MUL_SET_VEC_ROUND_PREFETCH(5);
+		MUL_SET_VEC_ROUND_PREFETCH(6);
+		MUL_SET_VEC_ROUND_PREFETCH(7);
+		r += UNROLL;
+		a += UNROLL;
+		len -= UNROLL;
+	}
+	if (len == UNROLL) {
+		MUL_SET_VEC_ROUND_PREFETCH(0);
+		MUL_SET_VEC_ROUND_PREFETCH(1);
+		MUL_SET_VEC_ROUND_PREFETCH(2);
+		MUL_SET_VEC_ROUND_PREFETCH(3);
+		MUL_SET_VEC_ROUND_PREFETCH(4);
+		MUL_SET_VEC_ROUND_PREFETCH(5);
+		MUL_SET_VEC_ROUND_PREFETCH(6);
+		MUL_SET_VEC_ROUND_NOPREFETCH(7);
+		return ((uint32_t)cy);
+	}
+	while (len > 1) {
+		MUL_SET_VEC_ROUND_PREFETCH(0);
+		++r;
+		++a;
+		--len;
+	}
+	if (len > 0) {
+		MUL_SET_VEC_ROUND_NOPREFETCH(0);
+	}
+	return ((uint32_t)cy);
+}
+
+/*
+ * r += a * b
+ * where r and a are vectors; b is a single 32-bit digit
+ */
+
+uint32_t
+big_mul_add_vec(uint32_t *r, uint32_t *a, int len, uint32_t b)
+{
+	uint64_t d, pf, p, t, cy;
+
+	if (len == 0)
+		return (0);
+	cy = 0;
+	d = (uint64_t)b;
+	pf = (uint64_t)a[0];
+	while (len > 8) {
+		MUL_ADD_VEC_ROUND_PREFETCH(0);
+		MUL_ADD_VEC_ROUND_PREFETCH(1);
+		MUL_ADD_VEC_ROUND_PREFETCH(2);
+		MUL_ADD_VEC_ROUND_PREFETCH(3);
+		MUL_ADD_VEC_ROUND_PREFETCH(4);
+		MUL_ADD_VEC_ROUND_PREFETCH(5);
+		MUL_ADD_VEC_ROUND_PREFETCH(6);
+		MUL_ADD_VEC_ROUND_PREFETCH(7);
+		r += 8;
+		a += 8;
+		len -= 8;
+	}
+	if (len == 8) {
+		MUL_ADD_VEC_ROUND_PREFETCH(0);
+		MUL_ADD_VEC_ROUND_PREFETCH(1);
+		MUL_ADD_VEC_ROUND_PREFETCH(2);
+		MUL_ADD_VEC_ROUND_PREFETCH(3);
+		MUL_ADD_VEC_ROUND_PREFETCH(4);
+		MUL_ADD_VEC_ROUND_PREFETCH(5);
+		MUL_ADD_VEC_ROUND_PREFETCH(6);
+		MUL_ADD_VEC_ROUND_NOPREFETCH(7);
+		return ((uint32_t)cy);
+	}
+	while (len > 1) {
+		MUL_ADD_VEC_ROUND_PREFETCH(0);
+		++r;
+		++a;
+		--len;
+	}
+	if (len > 0) {
+		MUL_ADD_VEC_ROUND_NOPREFETCH(0);
+	}
+	return ((uint32_t)cy);
+}
+#endif /* UNROLL8 */
+
+void
+big_sqr_vec(uint32_t *r, uint32_t *a, int len)
+{
+	uint32_t *tr, *ta;
+	int tlen, row, col;
+	uint64_t p, s, t, t2, cy;
+	uint32_t d;
+
+	tr = r + 1;
+	ta = a;
+	tlen = len - 1;
+	tr[tlen] = big_mul_set_vec(tr, ta + 1, tlen, ta[0]);
+	while (--tlen > 0) {
+		tr += 2;
+		++ta;
+		tr[tlen] = big_mul_add_vec(tr, ta + 1, tlen, ta[0]);
+	}
+	s = (uint64_t)a[0];
+	s = s * s;
+	r[0] = (uint32_t)s;
+	cy = s >> 32;
+	p = ((uint64_t)r[1] << 1) + cy;
+	r[1] = (uint32_t)p;
+	cy = p >> 32;
+	row = 1;
+	col = 2;
+	while (row < len) {
+		s = (uint64_t)a[row];
+		s = s * s;
+		p = (uint64_t)r[col] << 1;
+		t = p + s;
+		d = (uint32_t)t;
+		t2 = (uint64_t)d + cy;
+		r[col] = (uint32_t)t2;
+		cy = (t >> 32) + (t2 >> 32);
+		if (row == len - 1)
+			break;
+		p = ((uint64_t)r[col+1] << 1) + cy;
+		r[col+1] = (uint32_t)p;
+		cy = p >> 32;
+		++row;
+		col += 2;
+	}
+	r[col+1] = (uint32_t)cy;
+}
+
+#else /* ! UMUL64 */
+
+/*
+ * r = r + a * digit, r and a are vectors of length len
+ * returns the carry digit
+ */
+uint32_t
+big_mul_add_vec(uint32_t *r, uint32_t *a, int len, uint32_t digit)
+{
+	uint32_t cy, cy1, retcy, dlow, dhigh;
+	int i;
+
+	cy1 = 0;
+	dlow = digit & 0xffff;
+	dhigh = digit >> 16;
+	for (i = 0; i < len; i++) {
+		cy = (cy1 >> 16) + dlow * (a[i] & 0xffff) + (r[i] & 0xffff);
+		cy1 = (cy >> 16) + dlow * (a[i]>>16) + (r[i] >> 16);
+		r[i] = (cy & 0xffff) | (cy1 << 16);
+	}
+	retcy = cy1 >> 16;
+
+	cy1 = r[0] & 0xffff;
+	for (i = 0; i < len - 1; i++) {
+		cy = (cy1 >> 16) + dhigh * (a[i] & 0xffff) + (r[i] >> 16);
+		r[i] = (cy1 & 0xffff) | (cy << 16);
+		cy1 = (cy >> 16) + dhigh * (a[i] >> 16) + (r[i + 1] & 0xffff);
+	}
+	cy = (cy1 >> 16) + dhigh * (a[len - 1] & 0xffff) + (r[len - 1] >> 16);
+	r[len - 1] = (cy1 & 0xffff) | (cy << 16);
+	retcy = (cy >> 16) + dhigh * (a[len - 1] >> 16) + retcy;
+
+	return (retcy);
+}
+
+/*
+ * r = a * digit, r and a are vectors of length len
+ * returns the carry digit
+ */
+uint32_t
+big_mul_set_vec(uint32_t *r, uint32_t *a, int len, uint32_t digit)
+{
+	return (big_mul_add_vec(r, a, len, digit));
+}
+
+void
+big_sqr_vec(uint32_t *r, uint32_t *a, int len)
+{
+	int i;
+
+	r[len] = big_mul_set_vec(r, a, len, a[0]);
+	for (i = 1; i < len; ++i)
+		r[len + i] = big_mul_add_vec(r+i, a, len, a[i]);
+}
+
+#endif /* UMUL64 */
+
+void
+big_mul_vec(uint32_t *r, uint32_t *a, int alen, uint32_t *b, int blen)
+{
+	int i;
+
+	r[alen] = big_mul_set_vec(r, a, alen, b[0]);
+	for (i = 1; i < blen; ++i)
+		r[alen + i] = big_mul_add_vec(r+i, a, alen, b[i]);
+}
+
+
+#endif /* ! PSR_MUL */
+
+
+/*
+ * result = aa * bb  result->value should be big enough to hold the result
+ *
+ * Implementation: Standard grammar school algorithm
+ *
+ */
+
+BIG_ERR_CODE
+big_mul(BIGNUM *result, BIGNUM *aa, BIGNUM *bb)
+{
+	BIGNUM tmp1;
+	uint32_t tmp1value[BIGTMPSIZE];
+	uint32_t *r, *t, *a, *b;
+	BIG_ERR_CODE err;
+	int i, alen, blen, rsize, sign, diff;
+
+	if (aa == bb) {
+		diff = 0;
+	} else {
+		diff = big_cmp_abs(aa, bb);
+		if (diff < 0) {
+			BIGNUM *tt;
+			tt = aa;
+			aa = bb;
+			bb = tt;
+		}
+	}
+	a = aa->value;
+	b = bb->value;
+	alen = aa->len;
+	blen = bb->len;
+	while ((alen > 1) && (a[alen - 1] == 0)) alen--;
+	aa->len = alen;
+	while ((blen > 1) && (b[blen - 1] == 0)) blen--;
+	bb->len = blen;
+
+	rsize = alen + blen;
+	if (result->size < rsize) {
+		err = big_extend(result, rsize);
+		if (err != BIG_OK)
+			return (err);
+		/* aa or bb might be an alias to result */
+		a = aa->value;
+		b = bb->value;
+	}
+	r = result->value;
+
+	if (((alen == 1) && (a[0] == 0)) || ((blen == 1) && (b[0] == 0))) {
+		result->len = 1;
+		result->sign = 1;
+		r[0] = 0;
+		return (BIG_OK);
+	}
+	sign = aa->sign * bb->sign;
+	if ((alen == 1) && (a[0] == 1)) {
+		for (i = 0; i < blen; i++) r[i] = b[i];
+		result->len = blen;
+		result->sign = sign;
+		return (BIG_OK);
+	}
+	if ((blen == 1) && (b[0] == 1)) {
+		for (i = 0; i < alen; i++) r[i] = a[i];
+		result->len = alen;
+		result->sign = sign;
+		return (BIG_OK);
+	}
+
+	err = big_init1(&tmp1, rsize, tmp1value, arraysize(tmp1value));
+	if (err != BIG_OK)
+		return (err);
+	t = tmp1.value;
+	for (i = 0; i < rsize; i++) t[i] = 0;
+
+	if (diff == 0 && alen > 2)
+		BIG_SQR_VEC(t, a, alen);
+	else if (blen > 0)
+		BIG_MUL_VEC(t, a, alen, b, blen);
+	if (t[rsize - 1] == 0)
+		--rsize;
+	tmp1.len = rsize;
+	if ((err = big_copy(result, &tmp1)) != BIG_OK)
+		return (err);
+
+	result->sign = sign;
+
+	if (tmp1.malloced) big_finish(&tmp1);
+
+	return (BIG_OK);
+}
+
+
+/*
+ * caller must ensure that  a < n,  b < n  and  ret->size >=  2 * n->len + 1
+ * and that ret is not n
+ */
+BIG_ERR_CODE
+big_mont_mul(BIGNUM *ret, BIGNUM *a, BIGNUM *b, BIGNUM *n, uint32_t n0)
+{
+	int i, j, nlen, needsubtract;
+	uint32_t *nn, *rr;
+	uint32_t digit, c;
+	BIG_ERR_CODE err;
+
+	nlen = n->len;
+	nn = n->value;
+
+	rr = ret->value;
+
+	if ((err = big_mul(ret, a, b)) != BIG_OK)
+		return (err);
+
+	rr = ret->value;
+	for (i = ret->len; i < 2 * nlen + 1; i++) rr[i] = 0;
+	for (i = 0; i < nlen; i++) {
+		digit = rr[i];
+		digit = digit * n0;
+
+		c = BIG_MUL_ADD_VEC(rr + i, nn, nlen, digit);
+		j = i + nlen;
+		rr[j] += c;
+		while (rr[j] < c) {
+			rr[j + 1] += 1;
+			j++;
+			c = 1;
+		}
+	}
+
+	needsubtract = 0;
+	if ((rr[2 * nlen]  != 0))
+		needsubtract = 1;
+	else {
+		for (i = 2 * nlen - 1; i >= nlen; i--) {
+			if (rr[i] > nn[i - nlen]) {
+				needsubtract = 1;
+				break;
+			} else if (rr[i] < nn[i - nlen]) break;
+		}
+	}
+	if (needsubtract)
+		big_sub_vec(rr, rr + nlen, nn, nlen);
+	else {
+		for (i = 0; i < nlen; i++)
+			rr[i] = rr[i + nlen];
+	}
+	for (i = nlen - 1; (i >= 0) && (rr[i] == 0); i--);
+	ret->len = i+1;
+
+	return (BIG_OK);
+}
+
+uint32_t
+big_n0(uint32_t n)
+{
+	int i;
+	uint32_t result, tmp;
+
+	result = 0;
+	tmp = 0xffffffff;
+	for (i = 0; i < 32; i++) {
+		if ((tmp & 1) == 1) {
+			result = (result >> 1) | 0x80000000;
+			tmp = tmp - n;
+		} else  result = (result>>1);
+		tmp = tmp >> 1;
+	}
+
+	return (result);
+}
+
+
+int
+big_numbits(BIGNUM *n)
+{
+	int i, j;
+	uint32_t t;
+
+	for (i = n->len - 1; i > 0; i--)
+		if (n->value[i] != 0) break;
+	t = n->value[i];
+	for (j = 32; j > 0; j--) {
+		if ((t & 0x80000000) == 0)
+			t = t << 1;
+		else
+			return (32 * i + j);
+	}
+	return (0);
+}
+
+/* caller must make sure that a < n */
+BIG_ERR_CODE
+big_mont_rr(BIGNUM *result, BIGNUM *n)
+{
+	BIGNUM rr;
+	uint32_t rrvalue[BIGTMPSIZE];
+	int len, i;
+	BIG_ERR_CODE err;
+
+	rr.malloced = 0;
+	len = n->len;
+
+	if ((err = big_init1(&rr, 2 * len + 1,
+	    rrvalue, arraysize(rrvalue))) != BIG_OK)
+		return (err);
+
+	for (i = 0; i < 2 * len; i++) rr.value[i] = 0;
+	rr.value[2 * len] = 1;
+	rr.len = 2 * len + 1;
+	if ((err = big_div_pos(NULL, &rr, &rr, n)) != BIG_OK)
+		goto ret;
+	err = big_copy(result, &rr);
+ret:
+	if (rr.malloced) big_finish(&rr);
+	return (err);
+}
+
+/* caller must make sure that a < n */
+BIG_ERR_CODE
+big_mont_conv(BIGNUM *result, BIGNUM *a, BIGNUM *n, uint32_t n0, BIGNUM *n_rr)
+{
+	BIGNUM rr;
+	uint32_t rrvalue[BIGTMPSIZE];
+	int len, i;
+	BIG_ERR_CODE err;
+
+	rr.malloced = 0;
+	len = n->len;
+
+	if ((err = big_init1(&rr, 2 * len + 1, rrvalue, arraysize(rrvalue)))
+	    != BIG_OK)
+			return (err);
+
+	if (n_rr == NULL) {
+		for (i = 0; i < 2 * len; i++) rr.value[i] = 0;
+		rr.value[2 * len] = 1;
+		rr.len = 2 * len + 1;
+		if ((err = big_div_pos(NULL, &rr, &rr, n)) != BIG_OK)
+			goto ret;
+		n_rr = &rr;
+	}
+
+	if ((err = big_mont_mul(&rr, n_rr, a, n, n0)) != BIG_OK)
+		goto ret;
+	err = big_copy(result, &rr);
+ret:
+	if (rr.malloced) big_finish(&rr);
+	return (err);
+}
+
+
+#define	MAX_EXP_BIT_GROUP_SIZE 6
+#define	APOWERS_MAX_SIZE (1 << (MAX_EXP_BIT_GROUP_SIZE - 1))
+
+#ifdef USE_FLOATING_POINT
+
+/*
+ * This version makes use of floating point for performance.
+ */
+static BIG_ERR_CODE
+_big_modexp(BIGNUM *result, BIGNUM *a, BIGNUM *e, BIGNUM *n, BIGNUM *n_rr)
+{
+	BIGNUM ma, tmp, rr;
+	uint32_t mavalue[BIGTMPSIZE];
+	uint32_t tmpvalue[BIGTMPSIZE];
+	uint32_t rrvalue[BIGTMPSIZE];
+	int i, j, k, l, m, p, bit, bitind, bitcount, nlen;
+	BIG_ERR_CODE err;
+	uint32_t n0;
+	double dn0;
+	double *dn, *dt, *d16r, *d32r;
+	uint32_t *nint, *prod;
+	double *apowers[APOWERS_MAX_SIZE];
+	int nbits, groupbits, apowerssize;
+
+	nbits = big_numbits(e);
+	if (nbits < 50) {
+		groupbits = 1;
+		apowerssize = 1;
+	} else {
+		groupbits = MAX_EXP_BIT_GROUP_SIZE;
+		apowerssize = 1 << (groupbits - 1);
+	}
+
+	if ((err = big_init1(&ma, n->len, mavalue, arraysize(mavalue)))	!=
+	    BIG_OK)
+		return (err);
+	ma.len = 1;
+	ma.value[0] = 0;
+
+	if ((err = big_init1(&tmp, 2 * n->len + 1,
+	    tmpvalue, arraysize(tmpvalue))) != BIG_OK)
+		goto ret1;
+	tmp.len = 1;
+	tmp.value[0] = 0;
+
+	rr.malloced = 0;
+	if (n_rr == NULL) {
+		if ((err = big_init1(&rr, 2 * n->len + 1,
+		    rrvalue, arraysize(rrvalue))) != BIG_OK)
+			goto ret2;
+		if (big_mont_rr(&rr, n) != BIG_OK)
+			goto ret2;
+		n_rr = &rr;
+	}
+
+	n0 = big_n0(n->value[0]);
+
+	if (big_cmp_abs(a, n) > 0) {
+		if ((err = big_div_pos(NULL, &ma, a, n)) != BIG_OK)
+			goto ret2;
+		err = big_mont_conv(&ma, &ma, n, n0, n_rr);
+	} else {
+		err = big_mont_conv(&ma, a, n, n0, n_rr);
+	}
+	if (err != BIG_OK)
+		goto ret3;
+
+	tmp.len = 1;
+	tmp.value[0] = 1;
+	if ((err = big_mont_conv(&tmp, &tmp, n, n0, n_rr)) != BIG_OK)
+		goto ret3;
+
+	nlen = n->len;
+	dn0 = (double)(n0 & 0xffff);
+
+	dn = dt = d16r = d32r = NULL;
+	nint = prod = NULL;
+	for (i = 0; i < apowerssize; i++) {
+		apowers[i] = NULL;
+	}
+
+	if ((dn = big_malloc(nlen * sizeof (double))) == NULL) {
+		err = BIG_NO_MEM;
+		goto ret;
+	}
+	if ((dt = big_malloc((4 * nlen + 2) * sizeof (double))) == NULL) {
+		err = BIG_NO_MEM;
+		goto ret;
+	}
+	if ((nint = big_malloc(nlen * sizeof (uint32_t))) == NULL) {
+		err = BIG_NO_MEM;
+		goto ret;
+	}
+	if ((prod = big_malloc((nlen + 1) * sizeof (uint32_t))) == NULL) {
+		err = BIG_NO_MEM;
+		goto ret;
+	}
+	if ((d16r = big_malloc((2 * nlen + 1) * sizeof (double))) == NULL) {
+		err = BIG_NO_MEM;
+		goto ret;
+	}
+	if ((d32r = big_malloc(nlen * sizeof (double))) == NULL) {
+		err = BIG_NO_MEM;
+		goto ret;
+	}
+	for (i = 0; i < apowerssize; i++) {
+		if ((apowers[i] = big_malloc((2 * nlen + 1) *
+		    sizeof (double))) == NULL) {
+			err = BIG_NO_MEM;
+			goto ret;
+		}
+	}
+
+	for (i = 0; i < ma.len; i++) nint[i] = ma.value[i];
+	for (; i < nlen; i++) nint[i] = 0;
+	conv_i32_to_d32_and_d16(d32r, apowers[0], nint, nlen);
+
+	for (i = 0; i < n->len; i++) nint[i] = n->value[i];
+	for (; i < nlen; i++) nint[i] = 0;
+	conv_i32_to_d32(dn, nint, nlen);
+
+	mont_mulf_noconv(prod, d32r, apowers[0], dt, dn, nint, nlen, dn0);
+	conv_i32_to_d32(d32r, prod, nlen);
+	for (i = 1; i < apowerssize; i++) {
+		mont_mulf_noconv(prod, d32r, apowers[i - 1],
+		    dt, dn, nint, nlen, dn0);
+		conv_i32_to_d16(apowers[i], prod, nlen);
+	}
+
+	for (i = 0; i < tmp.len; i++) prod[i] = tmp.value[i];
+	for (; i < nlen + 1; i++) prod[i] = 0;
+
+	bitind = nbits % 32;
+	k = 0;
+	l = 0;
+	p = 0;
+	bitcount = 0;
+	for (i = nbits / 32; i >= 0; i--) {
+		for (j = bitind - 1; j >= 0; j--) {
+			bit = (e->value[i] >> j) & 1;
+			if ((bitcount == 0) && (bit == 0)) {
+				conv_i32_to_d32_and_d16(d32r, d16r,
+				    prod, nlen);
+				mont_mulf_noconv(prod, d32r, d16r,
+				    dt, dn, nint, nlen, dn0);
+			} else {
+				bitcount++;
+				p = p * 2 + bit;
+				if (bit == 1) {
+					k = k + l + 1;
+					l = 0;
+				} else {
+					l++;
+				}
+				if (bitcount == groupbits) {
+					for (m = 0; m < k; m++) {
+						conv_i32_to_d32_and_d16(
+							d32r, d16r,
+							prod, nlen);
+						mont_mulf_noconv(prod, d32r,
+						    d16r, dt, dn, nint,
+						    nlen, dn0);
+					}
+					conv_i32_to_d32(d32r, prod, nlen);
+					mont_mulf_noconv(prod, d32r,
+					    apowers[p >> (l+1)],
+					    dt, dn, nint, nlen, dn0);
+					for (m = 0; m < l; m++) {
+						conv_i32_to_d32_and_d16(
+							d32r, d16r,
+							prod, nlen);
+						mont_mulf_noconv(prod, d32r,
+						    d16r, dt, dn, nint,
+						    nlen, dn0);
+					}
+					k = 0;
+					l = 0;
+					p = 0;
+					bitcount = 0;
+				}
+			}
+		}
+		bitind = 32;
+	}
+
+	for (m = 0; m < k; m++) {
+		conv_i32_to_d32_and_d16(d32r, d16r, prod, nlen);
+		mont_mulf_noconv(prod, d32r, d16r, dt, dn, nint, nlen, dn0);
+	}
+	if (p != 0) {
+		conv_i32_to_d32(d32r, prod, nlen);
+		mont_mulf_noconv(prod, d32r, apowers[p >> (l + 1)],
+		    dt, dn, nint, nlen, dn0);
+	}
+	for (m = 0; m < l; m++) {
+		conv_i32_to_d32_and_d16(d32r, d16r, prod, nlen);
+		mont_mulf_noconv(prod, d32r, d16r, dt, dn, nint, nlen, dn0);
+	}
+
+	ma.value[0] = 1;
+	ma.len = 1;
+	for (i = 0; i < nlen; i++) tmp.value[i] = prod[i];
+	for (i = nlen - 1; (i > 0) && (prod[i] == 0); i--);
+	tmp.len = i + 1;
+	if ((err = big_mont_mul(&tmp, &tmp, &ma, n, n0)) != BIG_OK)
+		goto ret;
+	err = big_copy(result, &tmp);
+ret:
+	for (i = apowerssize - 1; i >= 0; i--) {
+		if (apowers[i] != NULL)
+			big_free(apowers[i], (2 * nlen + 1) * sizeof (double));
+	}
+	if (d32r != NULL)
+		big_free(d32r, nlen * sizeof (double));
+	if (d16r != NULL)
+		big_free(d16r, (2 * nlen + 1) * sizeof (double));
+	if (prod != NULL)
+		big_free(prod, (nlen + 1) * sizeof (uint32_t));
+	if (nint != NULL)
+		big_free(nint, nlen * sizeof (uint32_t));
+	if (dt != NULL)
+		big_free(dt, (4 * nlen + 2) * sizeof (double));
+	if (dn != NULL)
+		big_free(dn, nlen * sizeof (double));
+
+ret3:
+	big_finish(&rr);
+ret2:
+	big_finish(&tmp);
+ret1:
+	big_finish(&ma);
+	return (err);
+
+}
+
+#ifdef _KERNEL
+
+#include <sys/sysmacros.h>
+#include <sys/regset.h>
+#include <sys/fpu/fpusystm.h>
+
+/* the alignment for block stores to save fp registers */
+#define	FPR_ALIGN	(64)
+
+extern void big_savefp(kfpu_t *);
+extern void big_restorefp(kfpu_t *);
+
+#endif /* _KERNEL */
+
+BIG_ERR_CODE
+big_modexp(BIGNUM *result, BIGNUM *a, BIGNUM *e, BIGNUM *n, BIGNUM *n_rr)
+{
+#ifdef _KERNEL
+	BIG_ERR_CODE rv;
+	uint8_t fpua[sizeof (kfpu_t) + FPR_ALIGN];
+	kfpu_t *fpu;
+
+#ifdef DEBUG
+	if (!fpu_exists)
+		return (BIG_GENERAL_ERR);
+#endif
+
+	fpu =  (kfpu_t *)P2ROUNDUP((uintptr_t)fpua, FPR_ALIGN);
+	big_savefp(fpu);
+
+	rv = _big_modexp(result, a, e, n, n_rr);
+
+	big_restorefp(fpu);
+
+	return (rv);
+#else
+	return (_big_modexp(result, a, e, n, n_rr));
+#endif	/* _KERNEL */
+}
+
+#else /* ! USE_FLOATING_POINT */
+
+/*
+ * This version uses strictly integer math and is safe in the kernel
+ * for all platforms.
+ */
+
+/*
+ * computes a^e mod n
+ * assumes a < n, n odd, result->value at least as long as n->value
+ */
+BIG_ERR_CODE
+big_modexp(BIGNUM *result, BIGNUM *a, BIGNUM *e, BIGNUM *n, BIGNUM *n_rr)
+{
+	BIGNUM ma, tmp, rr;
+	uint32_t mavalue[BIGTMPSIZE];
+	uint32_t tmpvalue[BIGTMPSIZE];
+	uint32_t rrvalue[BIGTMPSIZE];
+	BIGNUM apowers[APOWERS_MAX_SIZE];
+	int i, j, k, l, m, p,
+	    bit, bitind, bitcount, groupbits, apowerssize;
+	BIG_ERR_CODE err;
+	uint32_t n0;
+
+	int nbits;
+
+	nbits = big_numbits(e);
+	if (nbits < 50) {
+		groupbits = 1;
+		apowerssize = 1;
+	} else {
+		groupbits = MAX_EXP_BIT_GROUP_SIZE;
+		apowerssize = 1 << (groupbits - 1);
+	}
+
+	if ((err = big_init1(&ma, n->len,
+	    mavalue, arraysize(mavalue))) != BIG_OK)
+		return (err);
+	ma.len = 1;
+	ma.value[0] = 0;
+
+	if ((err = big_init1(&tmp, 2 * n->len + 1,
+	    tmpvalue, arraysize(tmpvalue))) != BIG_OK)
+		goto ret1;
+	tmp.len = 1;
+	tmp.value[0] = 1;
+
+	n0 = big_n0(n->value[0]);
+
+	rr.malloced = 0;
+	if (n_rr == NULL) {
+		if ((err = big_init1(&rr, 2 * n->len + 1,
+		    rrvalue, arraysize(rrvalue))) != BIG_OK)
+			goto ret2;
+
+		if (big_mont_rr(&rr, n) != BIG_OK)
+			goto ret3;
+		n_rr = &rr;
+	}
+
+	for (i = 0; i < apowerssize; i++) apowers[i].malloced = 0;
+	for (i = 0; i < apowerssize; i++) {
+		if ((err = big_init1(&(apowers[i]), n->len, NULL, 0)) !=
+		    BIG_OK)
+			goto ret;
+	}
+
+	if (big_cmp_abs(a, n) > 0) {
+		if ((err = big_div_pos(NULL, &ma, a, n)) != BIG_OK)
+			goto ret;
+		err = big_mont_conv(&ma, &ma, n, n0, n_rr);
+	} else {
+		err = big_mont_conv(&ma, a, n, n0, n_rr);
+	}
+	if (err != BIG_OK) goto ret;
+
+	(void) big_copy(&(apowers[0]), &ma);
+	if ((err = big_mont_mul(&tmp, &ma, &ma, n, n0)) != BIG_OK)
+		goto ret;
+	(void) big_copy(&ma, &tmp);
+
+	for (i = 1; i < apowerssize; i++) {
+		if ((err = big_mont_mul(&tmp, &ma,
+		    &(apowers[i-1]), n, n0)) != BIG_OK)
+			goto ret;
+		(void) big_copy(&apowers[i], &tmp);
+	}
+
+	tmp.len = 1;
+	tmp.value[0] = 1;
+	if ((err = big_mont_conv(&tmp, &tmp, n, n0, n_rr)) != BIG_OK)
+		goto ret;
+
+	bitind = nbits % 32;
+	k = 0;
+	l = 0;
+	p = 0;
+	bitcount = 0;
+	for (i = nbits / 32; i >= 0; i--) {
+		for (j = bitind - 1; j >= 0; j--) {
+			bit = (e->value[i] >> j) & 1;
+			if ((bitcount == 0) && (bit == 0)) {
+				if ((err = big_mont_mul(&tmp,
+				    &tmp, &tmp, n, n0)) != BIG_OK)
+					goto ret;
+			} else {
+				bitcount++;
+				p = p * 2 + bit;
+				if (bit == 1) {
+					k = k + l + 1;
+					l = 0;
+				} else {
+					l++;
+				}
+				if (bitcount == groupbits) {
+					for (m = 0; m < k; m++) {
+						if ((err = big_mont_mul(&tmp,
+						    &tmp, &tmp, n, n0)) !=
+						    BIG_OK)
+							goto ret;
+					}
+					if ((err = big_mont_mul(&tmp, &tmp,
+					    &(apowers[p >> (l + 1)]),
+					    n, n0)) != BIG_OK)
+						goto ret;
+					for (m = 0; m < l; m++) {
+						if ((err = big_mont_mul(&tmp,
+						    &tmp, &tmp, n, n0)) !=
+						    BIG_OK)
+							goto ret;
+					}
+					k = 0;
+					l = 0;
+					p = 0;
+					bitcount = 0;
+				}
+			}
+		}
+		bitind = 32;
+	}
+
+	for (m = 0; m < k; m++) {
+		if ((err = big_mont_mul(&tmp, &tmp, &tmp, n, n0)) != BIG_OK)
+			goto ret;
+	}
+	if (p != 0) {
+		if ((err = big_mont_mul(&tmp, &tmp,
+		    &(apowers[p >> (l + 1)]), n, n0)) != BIG_OK)
+			goto ret;
+	}
+	for (m = 0; m < l; m++) {
+		if ((err = big_mont_mul(&tmp, &tmp, &tmp, n, n0)) != BIG_OK)
+			goto ret;
+	}
+
+	ma.value[0] = 1;
+	ma.len = 1;
+	if ((err = big_mont_mul(&tmp, &tmp, &ma, n, n0)) != BIG_OK)
+		goto ret;
+	err = big_copy(result, &tmp);
+ret:
+	for (i = apowerssize - 1; i >= 0; i--) {
+		big_finish(&(apowers[i]));
+	}
+ret3:
+	if (rr.malloced) big_finish(&rr);
+ret2:
+	if (tmp.malloced) big_finish(&tmp);
+ret1:
+	if (ma.malloced) big_finish(&ma);
+	return (err);
+}
+
+#endif /* USE_FLOATING_POINT */
+
+
+BIG_ERR_CODE
+big_modexp_crt(BIGNUM *result, BIGNUM *a, BIGNUM *dmodpminus1,
+    BIGNUM *dmodqminus1, BIGNUM *p, BIGNUM *q, BIGNUM *pinvmodq,
+    BIGNUM *p_rr, BIGNUM *q_rr)
+{
+	BIGNUM ap, aq, tmp;
+	int alen, biglen, sign;
+	BIG_ERR_CODE err;
+
+	if (p->len > q->len) biglen = p->len;
+	else biglen = q->len;
+
+	if ((err = big_init1(&ap, p->len, NULL, 0)) != BIG_OK)
+		return (err);
+	if ((err = big_init1(&aq, q->len, NULL, 0)) != BIG_OK)
+		goto ret1;
+	if ((err = big_init1(&tmp, biglen + q->len + 1, NULL, 0)) != BIG_OK)
+		goto ret2;
+
+	/*
+	 * check whether a is too short - to avoid timing attacks
+	 */
+	alen = a->len;
+	while ((alen > p->len) && (a->value[alen - 1] == 0)) {
+		alen--;
+	}
+	if (alen < p->len + q->len) {
+		/*
+		 * a is too short, add p*q to it before
+		 * taking it modulo p and q
+		 * this will also affect timing, but this difference
+		 * does not depend on p or q, only on a
+		 * (in "normal" operation, this path will never be
+		 * taken, so it is not a performance penalty
+		 */
+		if ((err = big_mul(&tmp, p, q)) != BIG_OK)
+			goto ret;
+		if ((err = big_add(&tmp, &tmp, a)) != BIG_OK)
+			goto ret;
+		if ((err = big_div_pos(NULL, &ap, &tmp, p)) != BIG_OK)
+			goto ret;
+		if ((err = big_div_pos(NULL, &aq, &tmp, q)) != BIG_OK)
+			goto ret;
+	} else {
+		if ((err = big_div_pos(NULL, &ap, a, p)) != BIG_OK)
+			goto ret;
+		if ((err = big_div_pos(NULL, &aq, a, q)) != BIG_OK)
+			goto ret;
+	}
+
+	if ((err = big_modexp(&ap, &ap, dmodpminus1, p, p_rr)) != BIG_OK)
+		goto ret;
+	if ((err = big_modexp(&aq, &aq, dmodqminus1, q, q_rr)) != BIG_OK)
+		goto ret;
+	if ((err = big_sub(&tmp, &aq, &ap)) != BIG_OK)
+		goto ret;
+	if ((err = big_mul(&tmp, &tmp, pinvmodq)) != BIG_OK)
+		goto ret;
+	sign = tmp.sign;
+	tmp.sign = 1;
+	if ((err = big_div_pos(NULL, &aq, &tmp, q)) != BIG_OK)
+		goto ret;
+	if ((sign == -1) && (!big_is_zero(&aq))) {
+		(void) big_sub_pos(&aq, q, &aq);
+	}
+	if ((err = big_mul(&tmp, &aq, p)) != BIG_OK)
+		goto ret;
+	err = big_add_abs(result, &ap, &tmp);
+
+ret:
+	big_finish(&tmp);
+ret2:
+	big_finish(&aq);
+ret1:
+	big_finish(&ap);
+
+	return (err);
+}
+
+
+uint32_t onearr[1] = {1};
+BIGNUM One = {1, 1, 1, 0, onearr};
+
+uint32_t twoarr[1] = {2};
+BIGNUM Two = {1, 1, 1, 0, twoarr};
+
+uint32_t fourarr[1] = {4};
+BIGNUM Four = {1, 1, 1, 0, fourarr};
+
+BIG_ERR_CODE
+big_sqrt_pos(BIGNUM *result, BIGNUM *n)
+{
+	BIGNUM *high, *low, *mid, *t;
+	BIGNUM t1, t2, t3, prod;
+	uint32_t t1value[BIGTMPSIZE];
+	uint32_t t2value[BIGTMPSIZE];
+	uint32_t t3value[BIGTMPSIZE];
+	uint32_t prodvalue[BIGTMPSIZE];
+	int i, nbits, diff, nrootbits, highbits;
+	BIG_ERR_CODE err;
+
+	nbits = big_numbits(n);
+
+	if ((err = big_init1(&t1, n->len + 1,
+	    t1value, arraysize(t1value))) != BIG_OK)
+		return (err);
+	if ((err = big_init1(&t2, n->len + 1,
+	    t2value, arraysize(t2value))) != BIG_OK)
+		goto ret1;
+	if ((err = big_init1(&t3, n->len + 1,
+	    t3value, arraysize(t3value))) != BIG_OK)
+		goto ret2;
+	if ((err = big_init1(&prod, n->len + 1,
+	    prodvalue, arraysize(prodvalue))) != BIG_OK)
+		goto ret3;
+
+	nrootbits = (nbits + 1) / 2;
+	t1.len = t2.len = t3.len = (nrootbits - 1) / 32 + 1;
+	for (i = 0; i < t1.len; i++) {
+		t1.value[i] = 0;
+		t2.value[i] = 0xffffffff;
+	}
+	highbits = nrootbits - 32 * (t1.len - 1);
+	if (highbits == 32) {
+		t1.value[t1.len - 1] = 0x80000000;
+		t2.value[t2.len - 1] = 0xffffffff;
+	} else {
+		t1.value[t1.len - 1] = 1 << (highbits - 1);
+		t2.value[t2.len - 1] = 2 * t1.value[t1.len - 1] - 1;
+	}
+	high = &t2;
+	low = &t1;
+	mid = &t3;
+
+	if ((err = big_mul(&prod, high, high)) != BIG_OK)
+		goto ret;
+	diff = big_cmp_abs(&prod, n);
+	if (diff <= 0) {
+		err = big_copy(result, high);
+		goto ret;
+	}
+
+	(void) big_sub_pos(mid, high, low);
+	while (big_cmp_abs(&One, mid) != 0) {
+		(void) big_add_abs(mid, high, low);
+		(void) big_half_pos(mid, mid);
+		if ((err = big_mul(&prod, mid, mid)) != BIG_OK)
+			goto ret;
+		diff = big_cmp_abs(&prod, n);
+		if (diff > 0) {
+			t = high;
+			high = mid;
+			mid = t;
+		} else if (diff < 0) {
+			t = low;
+			low = mid;
+			mid = t;
+		} else {
+			err = big_copy(result, low);
+			goto ret;
+		}
+		(void) big_sub_pos(mid, high, low);
+	}
+
+	err = big_copy(result, low);
+ret:
+	if (prod.malloced) big_finish(&prod);
+ret3:
+	if (t3.malloced) big_finish(&t3);
+ret2:
+	if (t2.malloced) big_finish(&t2);
+ret1:
+	if (t1.malloced) big_finish(&t1);
+
+	return (err);
+}
+
+
+BIG_ERR_CODE
+big_Jacobi_pos(int *jac, BIGNUM *nn, BIGNUM *mm)
+{
+	BIGNUM *t, *tmp2, *m, *n;
+	BIGNUM t1, t2, t3;
+	uint32_t t1value[BIGTMPSIZE];
+	uint32_t t2value[BIGTMPSIZE];
+	uint32_t t3value[BIGTMPSIZE];
+	int len, err;
+
+	if (big_is_zero(nn) ||
+	    (((nn->value[0] & 1) | (mm->value[0] & 1)) == 0)) {
+		*jac = 0;
+		return (BIG_OK);
+	}
+
+	if (nn->len > mm->len) len = nn->len;
+	else len = mm->len;
+
+	if ((err = big_init1(&t1, len,
+	    t1value, arraysize(t1value))) != BIG_OK)
+		return (err);
+	if ((err = big_init1(&t2, len,
+	    t2value, arraysize(t2value))) != BIG_OK)
+		goto ret1;
+	if ((err = big_init1(&t3, len,
+	    t3value, arraysize(t3value))) != BIG_OK)
+		goto ret2;
+
+	n = &t1;
+	m = &t2;
+	tmp2 = &t3;
+
+	(void) big_copy(n, nn);
+	(void) big_copy(m, mm);
+
+	*jac = 1;
+	while (big_cmp_abs(&One, m) != 0) {
+		if (big_is_zero(n)) {
+			*jac = 0;
+			goto ret;
+		}
+		if ((m->value[0] & 1) == 0) {
+			if (((n->value[0] & 7) == 3) ||
+			    ((n->value[0] & 7) == 5)) *jac = -*jac;
+			(void) big_half_pos(m, m);
+		} else if ((n->value[0] & 1) == 0) {
+			if (((m->value[0] & 7) == 3) ||
+			    ((m->value[0] & 7) == 5)) *jac = -*jac;
+			(void) big_half_pos(n, n);
+		} else {
+			if (((m->value[0] & 3) == 3) &&
+			    ((n->value[0] & 3) == 3)) {
+				*jac = -*jac;
+			}
+			if ((err = big_div_pos(NULL, tmp2, m, n)) != BIG_OK)
+				goto ret;
+			t = tmp2;
+			tmp2 = m;
+			m = n;
+			n = t;
+		}
+	}
+	err = BIG_OK;
+
+ret:
+	if (t3.malloced) big_finish(&t3);
+ret2:
+	if (t2.malloced) big_finish(&t2);
+ret1:
+	if (t1.malloced) big_finish(&t1);
+
+	return (err);
+}
+
+
+BIG_ERR_CODE
+big_Lucas(BIGNUM *Lkminus1, BIGNUM *Lk, BIGNUM *p, BIGNUM *k, BIGNUM *n)
+{
+	int m, w, i;
+	uint32_t bit;
+	BIGNUM ki, tmp, tmp2;
+	uint32_t kivalue[BIGTMPSIZE];
+	uint32_t tmpvalue[BIGTMPSIZE];
+	uint32_t tmp2value[BIGTMPSIZE];
+	BIG_ERR_CODE err;
+
+	if (big_cmp_abs(k, &One) == 0) {
+		(void) big_copy(Lk, p);
+		(void) big_copy(Lkminus1, &Two);
+		return (BIG_OK);
+	}
+
+	if ((err = big_init1(&ki, k->len + 1,
+	    kivalue, arraysize(kivalue))) != BIG_OK)
+		return (err);
+
+	if ((err = big_init1(&tmp, 2 * n->len +1,
+	    tmpvalue, arraysize(tmpvalue))) != BIG_OK)
+		goto ret1;
+
+	if ((err = big_init1(&tmp2, n->len,
+	    tmp2value, arraysize(tmp2value))) != BIG_OK)
+		goto ret2;
+
+	m = big_numbits(k);
+	ki.len = (m - 1) / 32 + 1;
+	w = (m - 1) / 32;
+	bit = 1 << ((m - 1) % 32);
+	for (i = 0; i < ki.len; i++) ki.value[i] = 0;
+	ki.value[ki.len - 1] = bit;
+	if (big_cmp_abs(k, &ki) != 0)
+		(void) big_double(&ki, &ki);
+	(void) big_sub_pos(&ki, &ki, k);
+
+	(void) big_copy(Lk, p);
+	(void) big_copy(Lkminus1, &Two);
+
+	for (i = 0; i < m; i++) {
+		if ((err = big_mul(&tmp, Lk, Lkminus1)) != BIG_OK)
+			goto ret;
+		(void) big_add_abs(&tmp, &tmp, n);
+		(void) big_sub_pos(&tmp, &tmp, p);
+		if ((err = big_div_pos(NULL, &tmp2, &tmp, n)) != BIG_OK)
+			goto ret;
+
+		if ((ki.value[w] & bit) != 0) {
+			if ((err = big_mul(&tmp, Lkminus1, Lkminus1)) !=
+			    BIG_OK)
+				goto ret;
+			(void) big_add_abs(&tmp, &tmp, n);
+			(void) big_sub_pos(&tmp, &tmp, &Two);
+			if ((err = big_div_pos(NULL, Lkminus1, &tmp, n)) !=
+			    BIG_OK)
+				goto ret;
+			(void) big_copy(Lk, &tmp2);
+		} else {
+			if ((err = big_mul(&tmp, Lk, Lk)) != BIG_OK)
+				goto ret;
+			(void) big_add_abs(&tmp, &tmp, n);
+			(void) big_sub_pos(&tmp, &tmp, &Two);
+			if ((err = big_div_pos(NULL, Lk, &tmp, n)) != BIG_OK)
+				goto ret;
+			(void) big_copy(Lkminus1, &tmp2);
+		}
+		bit = bit >> 1;
+		if (bit == 0) {
+			bit = 0x80000000;
+			w--;
+		}
+	}
+
+	err = BIG_OK;
+
+ret:
+	if (tmp2.malloced) big_finish(&tmp2);
+ret2:
+	if (tmp.malloced) big_finish(&tmp);
+ret1:
+	if (ki.malloced) big_finish(&ki);
+
+	return (err);
+}
+
+
+BIG_ERR_CODE
+big_isprime_pos(BIGNUM *n)
+{
+	BIGNUM o, nminus1, tmp, Lkminus1, Lk;
+	uint32_t ovalue[BIGTMPSIZE];
+	uint32_t nminus1value[BIGTMPSIZE];
+	uint32_t tmpvalue[BIGTMPSIZE];
+	uint32_t Lkminus1value[BIGTMPSIZE];
+	uint32_t Lkvalue[BIGTMPSIZE];
+	BIG_ERR_CODE err;
+	int e, i, jac;
+
+	if (big_cmp_abs(n, &One) == 0)
+		return (BIG_FALSE);
+	if (big_cmp_abs(n, &Two) == 0)
+		return (BIG_TRUE);
+	if ((n->value[0] & 1) == 0)
+		return (BIG_FALSE);
+
+	if ((err = big_init1(&o, n->len, ovalue, arraysize(ovalue))) != BIG_OK)
+		return (err);
+
+	if ((err = big_init1(&nminus1, n->len,
+	    nminus1value, arraysize(nminus1value))) != BIG_OK)
+		goto ret1;
+
+	if ((err = big_init1(&tmp, 2 * n->len,
+	    tmpvalue, arraysize(tmpvalue))) != BIG_OK)
+		goto ret2;
+
+	if ((err = big_init1(&Lkminus1, n->len,
+	    Lkminus1value, arraysize(Lkminus1value))) != BIG_OK)
+		goto ret3;
+
+	if ((err = big_init1(&Lk, n->len,
+	    Lkvalue, arraysize(Lkvalue))) != BIG_OK)
+		goto ret4;
+
+	(void) big_sub_pos(&o, n, &One); 	/* cannot fail */
+	(void) big_copy(&nminus1, &o);		/* cannot fail */
+	e = 0;
+	while ((o.value[0] & 1) == 0) {
+		e++;
+		(void) big_half_pos(&o, &o);  /* cannot fail */
+	}
+	if ((err = big_modexp(&tmp, &Two, &o, n, NULL)) != BIG_OK)
+		goto ret;
+	i = 0;
+	while ((i < e) &&
+	    (big_cmp_abs(&tmp, &One) != 0) &&
+	    (big_cmp_abs(&tmp, &nminus1) != 0)) {
+		if ((err = big_modexp(&tmp, &tmp, &Two, n, NULL)) !=  BIG_OK)
+			goto ret;
+		i++;
+	}
+	if (!((big_cmp_abs(&tmp, &nminus1) == 0) ||
+	    ((i == 0) && (big_cmp_abs(&tmp, &One) == 0)))) {
+		err = BIG_FALSE;
+		goto ret;
+	}
+
+	if ((err = big_sqrt_pos(&tmp, n)) != BIG_OK)
+		goto ret;
+	if ((err = big_mul(&tmp, &tmp, &tmp)) != BIG_OK)
+		goto ret;
+	if (big_cmp_abs(&tmp, n) == 0) {
+		err = BIG_FALSE;
+		goto ret;
+	}
+
+	(void) big_copy(&o, &Two);
+	do {
+		(void) big_add_abs(&o, &o, &One);
+		if ((err = big_mul(&tmp, &o, &o)) != BIG_OK)
+			goto ret;
+		(void) big_sub_pos(&tmp, &tmp, &Four);
+		if ((err = big_Jacobi_pos(&jac, &tmp, n)) != BIG_OK)
+			goto ret;
+	} while (jac != -1);
+
+	(void) big_add_abs(&tmp, n, &One);
+	if ((err = big_Lucas(&Lkminus1, &Lk, &o, &tmp, n)) != BIG_OK)
+		goto ret;
+	if ((big_cmp_abs(&Lkminus1, &o) == 0) && (big_cmp_abs(&Lk, &Two) == 0))
+		err = BIG_TRUE;
+	else err = BIG_FALSE;
+
+ret:
+	if (Lk.malloced) big_finish(&Lk);
+ret4:
+	if (Lkminus1.malloced) big_finish(&Lkminus1);
+ret3:
+	if (tmp.malloced) big_finish(&tmp);
+ret2:
+	if (nminus1.malloced) big_finish(&nminus1);
+ret1:
+	if (o.malloced) big_finish(&o);
+
+	return (err);
+}
+
+
+#define	SIEVESIZE 1000
+
+uint32_t smallprimes[] =
+{
+3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
+51, 53, 59, 61, 67, 71, 73, 79, 83, 89, 91, 97
+};
+
+
+BIG_ERR_CODE
+big_nextprime_pos(BIGNUM *result, BIGNUM *n)
+{
+	BIG_ERR_CODE err;
+	int sieve[SIEVESIZE];
+	int i;
+	uint32_t off, p;
+
+	if ((err = big_copy(result, n)) != BIG_OK)
+		return (err);
+	result->value[0] |= 1;
+	/* LINTED */
+	while (1) {
+		for (i = 0; i < SIEVESIZE; i++) sieve[i] = 0;
+		for (i = 0;
+		    i < sizeof (smallprimes) / sizeof (uint32_t); i++) {
+			p = smallprimes[i];
+			off = big_mod16_pos(result, p);
+			off = p - off;
+			if ((off % 2) == 1) off = off + p;
+			off = off/2;
+			while (off < SIEVESIZE) {
+				sieve[off] = 1;
+				off = off + p;
+			}
+		}
+
+		for (i = 0; i < SIEVESIZE; i++) {
+			if (sieve[i] == 0) {
+				err = big_isprime_pos(result);
+				if (err != BIG_FALSE) {
+					if (err != BIG_TRUE)
+						return (err);
+					else
+						return (BIG_OK);
+				}
+
+			}
+			if ((err = big_add_abs(result, result, &Two)) !=
+			    BIG_OK)
+				return (err);
+		}
+	}
+	return (BIG_OK);
+}
+
+
+BIG_ERR_CODE
+big_nextprime_pos_slow(BIGNUM *result, BIGNUM *n)
+{
+	BIG_ERR_CODE err;
+
+
+	if ((err = big_copy(result, n)) != BIG_OK)
+		return (err);
+	result->value[0] |= 1;
+	while ((err = big_isprime_pos(result)) != BIG_TRUE) {
+		if (err != BIG_FALSE)
+			return (err);
+		if ((err = big_add_abs(result, result, &Two)) != BIG_OK)
+			return (err);
+	}
+	return (BIG_OK);
+}
+
+
+/*
+ * given m and e, computes the rest in the equation
+ * gcd(m, e) = cm * m + ce * e
+ */
+BIG_ERR_CODE
+big_ext_gcd_pos(BIGNUM *gcd, BIGNUM *cm, BIGNUM *ce, BIGNUM *m, BIGNUM *e)
+{
+	BIGNUM *xi, *ri, *riminus1, *riminus2, *t,
+	    *vmi, *vei, *vmiminus1, *veiminus1;
+	BIGNUM t1, t2, t3, t4, t5, t6, t7, t8, tmp;
+	uint32_t t1value[BIGTMPSIZE];
+	uint32_t t2value[BIGTMPSIZE];
+	uint32_t t3value[BIGTMPSIZE];
+	uint32_t t4value[BIGTMPSIZE];
+	uint32_t t5value[BIGTMPSIZE];
+	uint32_t t6value[BIGTMPSIZE];
+	uint32_t t7value[BIGTMPSIZE];
+	uint32_t t8value[BIGTMPSIZE];
+	uint32_t tmpvalue[BIGTMPSIZE];
+	BIG_ERR_CODE err;
+	int len;
+
+	if (big_cmp_abs(m, e) >= 0) len = m->len;
+	else len = e->len;
+
+	if ((err = big_init1(&t1, len,
+	    t1value, arraysize(t1value))) != BIG_OK)
+		return (err);
+	if ((err = big_init1(&t2, len,
+	    t2value, arraysize(t2value))) != BIG_OK)
+			goto ret1;
+	if ((err = big_init1(&t3, len,
+	    t3value, arraysize(t3value))) != BIG_OK)
+			goto ret2;
+	if ((err = big_init1(&t4, len,
+	    t4value, arraysize(t3value))) != BIG_OK)
+			goto ret3;
+	if ((err = big_init1(&t5, len,
+	    t5value, arraysize(t5value))) != BIG_OK)
+			goto ret4;
+	if ((err = big_init1(&t6, len,
+	    t6value, arraysize(t6value))) != BIG_OK)
+			goto ret5;
+	if ((err = big_init1(&t7, len,
+	    t7value, arraysize(t7value))) != BIG_OK)
+			goto ret6;
+	if ((err = big_init1(&t8, len,
+	    t8value, arraysize(t8value))) != BIG_OK)
+			goto ret7;
+
+	if ((err = big_init1(&tmp, 2 * len,
+	    tmpvalue, arraysize(tmpvalue))) != BIG_OK)
+		goto ret8;
+
+	ri = &t1;
+	ri->value[0] = 1;
+	ri->len = 1;
+	xi = &t2;
+	riminus1 = &t3;
+	riminus2 = &t4;
+	vmi = &t5;
+	vei = &t6;
+	vmiminus1 = &t7;
+	veiminus1 = &t8;
+
+	(void) big_copy(vmiminus1, &One);
+	(void) big_copy(vmi, &One);
+	(void) big_copy(veiminus1, &One);
+	(void) big_copy(xi, &One);
+	vei->len = 1;
+	vei->value[0] = 0;
+
+	(void) big_copy(riminus1, m);
+	(void) big_copy(ri, e);
+
+	while (!big_is_zero(ri)) {
+		t = riminus2;
+		riminus2 = riminus1;
+		riminus1 = ri;
+		ri = t;
+		if ((err = big_mul(&tmp, vmi, xi)) != BIG_OK)
+			goto ret;
+		if ((err = big_sub(vmiminus1, vmiminus1, &tmp)) != BIG_OK)
+			goto ret;
+		t = vmiminus1;
+		vmiminus1 = vmi;
+		vmi = t;
+		if ((err = big_mul(&tmp, vei, xi)) != BIG_OK)
+			goto ret;
+		if ((err = big_sub(veiminus1, veiminus1, &tmp)) != BIG_OK)
+			goto ret;
+		t = veiminus1;
+		veiminus1 = vei;
+		vei = t;
+		if ((err = big_div_pos(xi, ri, riminus2, riminus1)) != BIG_OK)
+			goto ret;
+	}
+	if ((gcd != NULL) && ((err = big_copy(gcd, riminus1)) != BIG_OK))
+		goto ret;
+	if ((cm != NULL) && ((err = big_copy(cm, vmi)) != BIG_OK))
+		goto ret;
+	if (ce != NULL)
+		err = big_copy(ce, vei);
+ret:
+	if (tmp.malloced) big_finish(&tmp);
+ret8:
+	if (t8.malloced) big_finish(&t8);
+ret7:
+	if (t7.malloced) big_finish(&t7);
+ret6:
+	if (t6.malloced) big_finish(&t6);
+ret5:
+	if (t5.malloced) big_finish(&t5);
+ret4:
+	if (t4.malloced) big_finish(&t4);
+ret3:
+	if (t3.malloced) big_finish(&t3);
+ret2:
+	if (t2.malloced) big_finish(&t2);
+ret1:
+	if (t1.malloced) big_finish(&t1);
+
+	return (err);
+}