commit c0bc7b6e2d090554c9d940bc3614e089a688503a
parent d133f10c9a98c10f60df827e134f45f56e6ad0a0
Author: Mattias Andrée <maandree@kth.se>
Date: Thu, 3 Mar 2016 10:33:29 +0100
Add zabs, zadd, zdiv, zmod, zmodmul, zmodpow, zneg, zpow, zsub, and the newly introduced zmodsqr
Signed-off-by: Mattias Andrée <maandree@kth.se>
Diffstat:
15 files changed, 392 insertions(+), 4 deletions(-)
diff --git a/man/zmodmul.3 b/man/zmodmul.3
@@ -1,6 +1,6 @@
.TH ZMODMUL 3 libzahl
.SH NAME
-zmodmul - Calculate the modular product of two big integer
+zmodmul - Calculate a modular product of two big integer
.SH SYNOPSIS
.nf
#include <zahl.h>
@@ -34,6 +34,7 @@ It is possible to calculate the modular product
with a faster algorithm than calculating the
product and than the modulus of that product.
.SH SEE ALSO
+.BR zmodsqr (3),
.BR zmodpow (3),
.BR zstr (3),
.BR zadd (3),
diff --git a/man/zmodpow.3 b/man/zmodpow.3
@@ -34,6 +34,7 @@ It is possible to calculate the modular power
with a faster algorithm than calculating the
power and than the modulus of that power.
.SH SEE ALSO
+.BR zmodsqr (3),
.BR zmodmul (3),
.BR zsqr (3),
.BR zstr (3),
diff --git a/man/zmodsqr.3 b/man/zmodsqr.3
@@ -0,0 +1,45 @@
+.TH ZMODSQR 3 libzahl
+.SH NAME
+zsqr - Calculate a modular square of a big integer
+.SH SYNOPSIS
+.nf
+#include <zahl.h>
+
+void zmodsqr(z_t \fIsquare\fP, z_t \fIinteger\fP, z_t \fImodulator\fP);
+.fi
+.SH DESCRIPTION
+.B zmodsqr
+calculates the square of an
+.IR integer ,
+modulus a
+.IR modulator ,
+and stores the result in
+.IR square .
+That is,
+.I square
+gets
+.IR integer ².
+Mod
+.IR modulator .
+.P
+It is safe to call
+.B zmodsqr
+with non-unique parameters.
+.SH RATIONALE
+See rationle for
+.BR zmodmul (3),
+and
+.BR zsqr (3).
+.SH SEE ALSO
+.BR zmodmul (3),
+.BR zmodpow (3),
+.BR zsqr (3),
+.BR zstr (3),
+.BR zadd (3),
+.BR zsub (3),
+.BR zmul (3),
+.BR zdiv (3),
+.BR zmod (3),
+.BR zneg (3),
+.BR zabs (3),
+.BR zpow (3)
diff --git a/src/internals.h b/src/internals.h
@@ -3,6 +3,7 @@
#define BITS_PER_CHAR 32
#define LB_BITS_PER_CHAR 5
+#define ZAHL_CHAR_MAX UINT32_MAX
#define FLOOR_BITS_TO_CHARS(bits) ((bits) >> LB_BITS_PER_CHAR)
#define CEILING_BITS_TO_CHARS(bits) (((bits) + (BITS_PER_CHAR - 1)) >> LB_BITS_PER_CHAR)
@@ -15,7 +16,14 @@
X(libzahl_tmp_str_div)\
X(libzahl_tmp_str_rem)\
X(libzahl_tmp_gcd_u)\
- X(libzahl_tmp_gcd_v)
+ X(libzahl_tmp_gcd_v)\
+ X(libzahl_tmp_modmul)\
+ X(libzahl_tmp_div)\
+ X(libzahl_tmp_mod)\
+ X(libzahl_tmp_pow_b)\
+ X(libzahl_tmp_pow_c)\
+ X(libzahl_tmp_pow_d)\
+ X(libzahl_tmp_modsqr)
#define LIST_CONSTS\
X(libzahl_const_1e19, zsetu, 10000000000000000000ULL) /* The largest power of 10 < 2⁶⁴. */\
diff --git a/src/zabs.c b/src/zabs.c
@@ -0,0 +1,11 @@
+/* See LICENSE file for copyright and license details. */
+#include "internals"
+
+
+void
+zabs(z_t a, z_t b)
+{
+ if (a != b)
+ zset(a, b);
+ SET_SIGNUM(a, !zzero(a));
+}
diff --git a/src/zadd.c b/src/zadd.c
@@ -0,0 +1,94 @@
+/* See LICENSE file for copyright and license details. */
+#include "internals"
+
+#include <stdlib.h>
+#include <string.h>
+
+
+void
+zadd_unsigned(z_t a, z_t b, z_t c)
+{
+ size_t i, size, n;
+ uint32_t carry[] = {0, 0};
+ zahl_char_t *addend;
+
+ if (a == c) {
+ zadd_unsigned(a, c, b);
+ return;
+ }
+
+ size = b->used > c->used ? b->used : c->used;
+ if (a->alloced < size + 1) {
+ a->alloced = size + 1;
+ a->chars = realloc(a->chars, a->alloced * sizeof(*(a->chars)));
+ }
+ a->chars[size] = 0;
+
+ n = b->used + c->used - size;
+ if (a == b) {
+ if (a->used < c->used) {
+ n = c->used;
+ memset(a->chars + a->used, 0, (n - a->used) * sizeof(*(a->chars)));
+ }
+ addend = c->chars;
+ } else if (a == c) {
+ if (a->used < b->used) {
+ n = b->used;
+ memset(a->chars + a->used, 0, (n - a->used) * sizeof(*(a->chars)));
+ }
+ addend = b->chars;
+ } else if (b->used > c->used) {
+ memcpy(a->chars, b->chars, b->used * sizeof(*(a->chars)));
+ a->used = b->used;
+ addend = c->chars;
+ } else {
+ memcpy(a->chars, c->chars, c->used * sizeof(*(a->chars)));
+ a->used = c->used;
+ addend = b->chars;
+ }
+
+ for (i = 0; i < n; i++) {
+ if (carry[i & 1])
+ carry[~i & 1] = (ZAHL_CHAR_MAX - a->chars[i] <= addend[i]);
+ else
+ carry[~i & 1] = (ZAHL_CHAR_MAX - a->chars[i] < addend[i]);
+ a->chars[i] += addend[i] + carry[i & 1];
+ }
+
+ while (carry[~i & 1]) {
+ carry[i & 1] = a->chars[i] == ZAHL_CHAR_MAX;
+ a->chars[i++] += 1;
+ }
+
+ if (a->used < i)
+ a->used = i;
+ SET_SIGNUM(a, !zzero(b) | !zzero(c));
+}
+
+
+void
+zadd(z_t a, z_t b, z_t c)
+{
+ if (zzero(b)) {
+ if (a != b)
+ zset(a, c);
+ } else if (zzero(c)) {
+ if (a != b)
+ zset(a, b);
+ } else if (b == c) {
+ zlsh(a, b, 1);
+ } else if ((zsignum(b) | zsignum(c)) < 0) {
+ if (zsignum(b) < 0) {
+ if (zsignum(c) < 0) {
+ zadd_unsigned(a, b, c);
+ SET_SIGNUM(a, -zsignum(a));
+ } else {
+ zsub_unsigned(a, c, b);
+ }
+ } else {
+ zsub_unsigned(a, b, c);
+ }
+ } else {
+ zadd_unsigned(a, b, c);
+ }
+}
diff --git a/src/zdiv.c b/src/zdiv.c
@@ -0,0 +1,9 @@
+/* See LICENSE file for copyright and license details. */
+#include "internals"
+
+
+void
+zdiv(z_t a, z_t b, z_t c)
+{
+ zdivmod(a, libzahl_tmp_div, b, c);
+}
diff --git a/src/zmod.c b/src/zmod.c
@@ -0,0 +1,9 @@
+/* See LICENSE file for copyright and license details. */
+#include "internals"
+
+
+void
+zmod(z_t a, z_t b, z_t c)
+{
+ zdivmod(libzahl_tmp_mod, a, b, c);
+}
diff --git a/src/zmodmul.c b/src/zmodmul.c
@@ -0,0 +1,17 @@
+/* See LICENSE file for copyright and license details. */
+#include "internals"
+
+
+void
+zmodmul(z_t a, z_t b, z_t c)
+{
+ /* TODO Montgomery modular multiplication */
+ if (a == d) {
+ zset(libzahl_tmp_modmul, d);
+ zmul(a, b, c);
+ zmod(a, a, libzahl_tmp_modmul);
+ } else {
+ zmul(a, b, c);
+ zmod(a, a, d);
+ }
+}
diff --git a/src/zmodpow.c b/src/zmodpow.c
@@ -0,0 +1,50 @@
+/* See LICENSE file for copyright and license details. */
+#include "internals"
+
+#include <errno.h>
+
+#define tb libzahl_tmp_pow_b
+#define tc libzahl_tmp_pow_c
+#define td libzahl_tmp_pow_d
+
+
+void
+zmodpow(z_t a, z_t b, z_t c, z_t d)
+{
+ size_t i, n;
+
+ if (zsignum(c) <= 0) {
+ if (zzero(c)) {
+ if (zzero(b)) {
+ errno = EDOM; /* Indeterminate form: 0:th power of 0 */
+ FAILURE_JUMP();
+ }
+ zsetu(a, 1);
+ } else if (zzero(b) || zzero(d)) {
+ errno = EDOM; /* Undefined form: Division by 0 */
+ FAILURE_JUMP();
+ } else {
+ SET_SIGNUM(a, 0);
+ }
+ return;
+ } else if (zzero(d)) {
+ errno = EDOM; /* Undefined form: Division by 0 */
+ FAILURE_JUMP();
+ } else if (zzero(b)) {
+ SET_SIGNUM(a, 0);
+ return;
+ }
+
+ n = zbits(c);
+
+ zmod(tb, b, d);
+ zset(tc, c);
+ zset(td, d);
+ zsetu(a, 1);
+
+ for (i = 0; i < n; i++) {
+ if (zbtest(tc, i))
+ zmodmul(a, a, tb, td);
+ zmodsqr(tb, tb, td);
+ }
+}
diff --git a/src/zmodsqr.c b/src/zmodsqr.c
@@ -0,0 +1,17 @@
+/* See LICENSE file for copyright and license details. */
+#include "internals"
+
+
+void
+zmodsqr(z_t a, z_t b, z_t c)
+{
+ /* TODO What is the fastest way to do zmodsqr? */
+ if (a == c) {
+ zset(libzahl_tmp_modsqr, c);
+ zsqr(a, b);
+ zmod(a, a, libzahl_tmp_modsqr);
+ } else {
+ zsqr(a, b);
+ zmod(a, a, c);
+ }
+}
diff --git a/src/zneg.c b/src/zneg.c
@@ -0,0 +1,11 @@
+/* See LICENSE file for copyright and license details. */
+#include "internals"
+
+
+void
+zneg(z_t a, z_t b)
+{
+ if (a != b)
+ zset(a, b);
+ SET_SIGNUM(a, -zsignum(a));
+}
diff --git a/src/zpow.c b/src/zpow.c
@@ -0,0 +1,45 @@
+/* See LICENSE file for copyright and license details. */
+#include "internals"
+
+#include <errno.h>
+
+#define tb libzahl_tmp_pow_b
+#define tc libzahl_tmp_pow_c
+
+
+void
+zpow(z_t a, z_t b, z_t c)
+{
+ size_t i, n;
+
+ if (zsignum(c) <= 0) {
+ if (zzero(c)) {
+ if (zzero(b)) {
+ errno = EDOM; /* Indeterminate form: 0:th power of 0 */
+ FAILURE_JUMP();
+ }
+ zsetu(a, 1);
+ } else if (zzero(b)) {
+ errno = EDOM; /* Undefined form: Division by 0 */
+ FAILURE_JUMP();
+ } else {
+ SET_SIGNUM(a, 0);
+ }
+ return;
+ } else if (zzero(b)) {
+ SET_SIGNUM(a, 0);
+ return;
+ }
+
+ n = zbits(c);
+
+ zset(tb, b);
+ zset(tc, c);
+ zsetu(a, 1);
+
+ for (i = 0; i < n; i++) {
+ if (zbtest(tc, i))
+ zmul(a, a, tb);
+ zsqr(tb, tb);
+ }
+}
diff --git a/src/zsub.c b/src/zsub.c
@@ -0,0 +1,70 @@
+/* See LICENSE file for copyright and license details. */
+#include "internals"
+
+#include <stdlib.h>
+#include <string.h>
+
+
+void
+zsub_unsigned(z_t a, z_t b, z_t c)
+{
+ int magcmp = zcmpmag(b, c);
+ z_t s;
+ size_t i, n;
+ zahl_char_t carry[] = {0, 0};
+
+ if (magcmp <= 0) {
+ if (magcmp == 0) {
+ SET_SIGNUM(a, 0);
+ return;
+ }
+ if (a != c)
+ zset(a, c);
+ *subtrahend = *b;
+ } else {
+ if (a != b)
+ zset(a, b);
+ *subtrahend = *c;
+ }
+
+ n = a->used < s->used ? a->used : s->used;
+
+ for (i = 0; i < n; i++) {
+ carry[~i & 1] = carry[i & 1] ? (a->chars[i] <= s->chars[i]) : (a->chars[i] < s->chars[i]);
+ a->chars[i] -= s->chars[i];
+ a->chars[i] -= carry[i & 1];
+ }
+
+ if (carry[i] & 1) {
+ while (!a->chars[i])
+ a->chars[i++] = ZAHL_CHAR_MAX;
+ a->chars[i] -= 1;
+ }
+ SET_SIGNUM(a, magcmp);
+}
+
+void
+zsub(z_t a, z_t b, z_t c)
+{
+ if (b == c) {
+ SET_SIGNUM(a, 0);
+ } else if (zzero(b)) {
+ zneg(a, c);
+ } else if (zzero(c)) {
+ if (a != b)
+ zset(a, b);
+ } else if ((zsignum(b) | zsignum(c)) < 0) {
+ if (zsignum(b) < 0) {
+ if (zsignum(c) < 0) {
+ zsub_unsigned(a, c, b);
+ } else {
+ zadd_unsigned(a, b, c);
+ SET_SIGNUM(a, -zsignum(a));
+ }
+ } else {
+ zadd_unsigned(a, b, c);
+ }
+ } else {
+ zsub_unsigned(a, b, c);
+ }
+}
diff --git a/zahl.h b/zahl.h
@@ -73,6 +73,7 @@ void zdiv(z_t, z_t, z_t); /* a := b / c */
void zdivmod(z_t, z_t, z_t, z_t); /* a := c / d, b = c % d */
void zmod(z_t, z_t, z_t); /* a := b % c */
void zsqr(z_t, z_t); /* a := b² */
+void zmodsqr(z_t, z_t, z_t); /* a := b² % c */
void zneg(z_t, z_t); /* a := -b */
void zabs(z_t, z_t); /* a := |b| */
void zpow(z_t, z_t, z_t); /* a := b ↑ c */
@@ -80,8 +81,7 @@ void zmodpow(z_t, z_t, z_t, z_t); /* a := (b ↑ c) % d */
/* These are used internally and may be removed in a future version. */
void zadd_unsigned(z_t, z_t, z_t); /* a := |b| + |c|, b and c must not be the same reference. */
-void zsub_unsigned(z_t, z_t, z_t); /* a := |b| - |c|, b and c must not be the same reference. */
-void zsub_positive(z_t, z_t, z_t); /* a := |b| - |c|, assumes b ≥ c and that, and b is not c. */
+void zsub_unsigned(z_t, z_t, z_t); /* a := |b| - |c| */
/* Bitwise operations. */
