* sense fast. Here, we do efficient reduction modulo sparse polynomials.
* (It works for arbitrary polynomials, but isn't efficient for dense ones.)
*
- * Suppose that %$p(x) = x^n + p'(x) = \sum_{0\le i<n} p_i x^i$%, hopefully
- * with only a few other %$p_i \ne 0$%. We're going to compile %$p$% into a
- * sequence of instructions which can be used to perform reduction modulo
- * %$p$%. The important observation is that %$x^n \equiv p' \pmod p$%.
+ * Suppose that %$p = x^n + p'$% where %$p' = \sum_{0\le i<n} p_i x^i$%,
+ * hopefully with only a few %$p_i \ne 0$%. We're going to compile %$p$%
+ * into a sequence of instructions which can be used to perform reduction
+ * modulo %$p$%. The important observation is that
+ * %$x^n \equiv p' \pmod p$%.
*
* Suppose we're working with %$w$%-bit words; let %$n = N w + n'$% with
* %$0 \le n' < w$%. Let %$u(x)$% be some arbitrary polynomial. Write
- * %$u = z x^k + u'$% with %$\deg u' < k \ge n$%; then a reduction step uses
+ * %$u = z x^k + u'$% with %$\deg u' < k \ge n$%. Then a reduction step uses
* that %$u \equiv u' + z p' x^{k-n} \pmod p$%: the right hand side has
* degree %$\max \{ \deg u', k + \deg p' - n + \deg z \} < \deg u$%, so this
* makes progress towards a complete reduction.
unsigned f; /* Flags */
#define f_lsr 1u /* Overflow from previous word */
#define f_load 2u /* Outstanding @LOAD@ */
+#define f_fip 4u /* Final-pass offset is set */
instr_v iv; /* Instruction vector */
+ size_t fip; /* Offset for final-pass reduction */
size_t w; /* Currently loaded target word */
size_t wi; /* Left-shifts for current word */
+ gfreduce *r; /* Reduction context pointer */
};
#define INSTR(g_, op_, arg_) do { \
static void emit_load(struct gen *g, size_t w)
{
+ /* --- If this is not the low-order word then note final-pass start --- *
+ *
+ * Once we've eliminated the whole high-degree words, there will possibly
+ * remain a few high-degree bits. We can further reduce the subject
+ * polynomial by subtracting an appropriate multiple of %$p'$%, but if we
+ * do this naively we'll end up addressing `low-order' words beyond the
+ * bottom of our input. We solve this problem by storing an alternative
+ * start position for this final pass (which works because we scan bits
+ * right-to-left).
+ */
+
+ if (!(g->f & f_fip) && w < g->r->lim) {
+ g->fip = DA_LEN(&g->iv);
+ g->f |= f_fip;
+ }
+
+ /* --- Actually emit the instruction --- */
+
INSTR(g, GFRI_LOAD, w);
g->f |= f_load;
g->w = w;
/* --- Sort out the easy stuff --- */
+ g.r = r;
d = mp_bits(p); assert(d); d--;
r->lim = d/MPW_BITS;
dw = d%MPW_BITS;
INSTR(&g, GFRI_STORE, g.w);
}
+ /* --- Copy the instruction vector.
+ *
+ * If we've not set a final-pass offset yet then now would be an excellent
+ * time. Obviously it should be right at the end, because there's nothing
+ * for a final pass to do.
+ */
+
r->in = DA_LEN(&g.iv);
r->iv = xmalloc(r->in * sizeof(gfreduce_instr));
memcpy(r->iv, DA(&g.iv), r->in * sizeof(gfreduce_instr));
+
+ if (!(g.f & f_fip)) g.fip = DA_LEN(&g.iv);
+ r->fiv = r->iv + g.fip;
+
DA_DESTROY(&g.iv);
}
#undef f_lsr
#undef f_load
+#undef f_fip
/* --- @gfreduce_destroy@ --- *
*
/* --- @gfreduce_dump@ --- *
*
- * Arguments: @gfreduce *r@ = structure to dump
+ * Arguments: @const gfreduce *r@ = structure to dump
* @FILE *fp@ = file to dump on
*
* Returns: ---
* Use: Dumps a reduction context.
*/
-void gfreduce_dump(gfreduce *r, FILE *fp)
+void gfreduce_dump(const gfreduce *r, FILE *fp)
{
size_t i;
(unsigned long)r->lim, (unsigned long)r->mask);
for (i = 0; i < r->in; i++) {
static const char *opname[] = { "load", "lsl", "lsr", "store" };
+ if (&r->iv[i] == r->fiv)
+ fputs("final:\n", fp);
assert(r->iv[i].op < N(opname));
fprintf(fp, " %s %lu\n",
opname[r->iv[i].op],
(unsigned long)r->iv[i].arg);
}
+ if (&r->iv[i] == r->fiv)
+ fputs("final:\n", fp);
}
/* --- @gfreduce_do@ --- *
*
- * Arguments: @gfreduce *r@ = reduction context
+ * Arguments: @const gfreduce *r@ = reduction context
* @mp *d@ = destination
* @mp *x@ = source
*
}
}
-mp *gfreduce_do(gfreduce *r, mp *d, mp *x)
+mp *gfreduce_do(const gfreduce *r, mp *d, mp *x)
{
mpw *v, *vl;
const gfreduce_instr *il;
while (*vl & r->mask) {
z = *vl & r->mask;
*vl &= ~r->mask;
- run(r->iv, il, vl, z);
+ run(r->fiv, il, vl, z);
}
}
}
/* --- @gfreduce_sqrt@ --- *
*
- * Arguments: @gfreduce *r@ = pointer to reduction context
+ * Arguments: @const gfreduce *r@ = pointer to reduction context
* @mp *d@ = destination
* @mp *x@ = some polynomial
*
* Returns: The square root of @x@ modulo @r->p@, or null.
*/
-mp *gfreduce_sqrt(gfreduce *r, mp *d, mp *x)
+mp *gfreduce_sqrt(const gfreduce *r, mp *d, mp *x)
{
mp *y = MP_COPY(x);
mp *z, *spare = MP_NEW;
unsigned long m = mp_bits(r->p) - 1;
unsigned long i;
+ /* --- This is pretty easy --- *
+ *
+ * Note that %$x = x^{2^m}$%; therefore %$(x^{2^{m-1}})^2 = x^{2^m} = x$%,
+ * so %$x^{2^{m-1}}$% is the square root we seek.
+ */
+
for (i = 0; i < m - 1; i++) {
mp *t = gf_sqr(spare, y);
spare = y;
/* --- @gfreduce_trace@ --- *
*
- * Arguments: @gfreduce *r@ = pointer to reduction context
+ * Arguments: @const gfreduce *r@ = pointer to reduction context
* @mp *x@ = some polynomial
*
* Returns: The trace of @x@. (%$\Tr(x)=x + x^2 + \cdots + x^{2^{m-1}}$%
- * if %$x \in \gf{2^m}$%).
+ * if %$x \in \gf{2^m}$%). Since the trace is invariant under
+ * the Frobenius automorphism (i.e., %$\Tr(x)^2 = \Tr(x)$%), it
+ * must be an element of the base field, i.e., %$\gf{2}$%, and
+ * we only need a single bit to represent it.
*/
-int gfreduce_trace(gfreduce *r, mp *x)
+int gfreduce_trace(const gfreduce *r, mp *x)
{
mp *y = MP_COPY(x);
mp *spare = MP_NEW;
/* --- @gfreduce_halftrace@ --- *
*
- * Arguments: @gfreduce *r@ = pointer to reduction context
+ * Arguments: @const gfreduce *r@ = pointer to reduction context
* @mp *d@ = destination
* @mp *x@ = some polynomial
*
* if %$x \in \gf{2^m}$% with %$m$% odd).
*/
-mp *gfreduce_halftrace(gfreduce *r, mp *d, mp *x)
+mp *gfreduce_halftrace(const gfreduce *r, mp *d, mp *x)
{
mp *y = MP_COPY(x);
mp *spare = MP_NEW;
/* --- @gfreduce_quadsolve@ --- *
*
- * Arguments: @gfreduce *r@ = pointer to reduction context
+ * Arguments: @const gfreduce *r@ = pointer to reduction context
* @mp *d@ = destination
* @mp *x@ = some polynomial
*
- * Returns: A polynomial @y@ such that %$y^2 + y = x$%, or null.
+ * Returns: A polynomial @z@ such that %$z^2 + z = x$%, or null.
+ *
+ * Use: Solves quadratic equations in a field with characteristic 2.
+ * Suppose we have an equation %$y^2 + A y + B = 0$% where
+ * %$A \ne 0$%. (If %$A = 0$% then %$y = \sqrt{B}$% and you
+ * want @gfreduce_sqrt@ instead.) Use this function to solve
+ * %$z^2 + z = B/A^2$%; then set %$y = A z$%, since
+ * %$y^2 + y = A^2 z^2 + A^2 z = A^2 (z^2 + z) = B$% as
+ * required.
+ *
+ * The two roots are %$z$% and %$z + 1$%; this function always
+ * returns the one with zero scalar coefficient.
*/
-mp *gfreduce_quadsolve(gfreduce *r, mp *d, mp *x)
+mp *gfreduce_quadsolve(const gfreduce *r, mp *d, mp *x)
{
unsigned long m = mp_bits(r->p) - 1;
mp *t;
+ /* --- About the solutions --- *
+ *
+ * Factor %$z^2 + z = z (z + 1)$%. Therefore, if %$z^2 + z = x$% and
+ * %$z' = z + 1$% then %$z'^2 + z' = z^2 + 1 + z + 1 = z^2 + z$%, so
+ * %$z + 1$% is the other solution.
+ *
+ * A solution exists if and only if %$\Tr(x) = 0$%. To see the `only if'
+ * implication, recall that the trace function is linear, and hence
+ * $%\Tr(z^2 + z) = \Tr(z)^2 + \Tr(z) = \Tr(z) + \Tr(z) = 0$%. The `if'
+ * direction will be proven using explicit constructions captured in the
+ * code below.
+ */
+
MP_COPY(x);
- if (m & 1)
+ if (m & 1) {
+
+ /* --- A short-cut for fields with odd degree ---
+ *
+ * The method below works in all binary fields, but there's a quicker way
+ * which works whenever the degree is odd. The half-trace is
+ * %$z = \sum_{0\le i\le (m-1)/2} x^{2^{2i}}$%. Then %$z^2 + z = {}$%
+ * %$\sum_{0\le i\le (m-1)/2} (x^{2^{2i}} + x^{2^{2i+1}}) = {}$%
+ * %$\Tr(x) + x^{2^m} = \Tr(x) + x$%. This therefore gives us the
+ * solution we want whenever %$\Tr(x) = 0$%.
+ */
+
d = gfreduce_halftrace(r, d, x);
- else {
+ } else {
mp *z, *w, *rho = MP_NEW;
mp *spare = MP_NEW;
grand *fr = fibrand_create(0);
unsigned long i;
+ /* --- Unpicking the magic --- *
+ *
+ * Choose %$\rho \inr \gf{2^m}$% with %$\Tr(\rho) = 1$%. Let
+ * %$z = \sum_{0\le i<m} \rho^{2^i} \sum_{0\le j<i} x^{2^j} = {}$%
+ * %$\sum_{1\le i<m} \rho^{2^i} (x + \sum_{1\le j<i} x^{2^j} = {}$%
+ * %$\rho^2 x + \rho^4 (x + x^2) + \rho^8 (x + x^2 + x^4) + \cdots + {}$%
+ * %$\rho^{2^{m-1}} (x + x^2 + x^{2^{m-2}})$%. Then %$z^2 = {}$%
+ * %$\sum_{0\le i<m} \rho^{2^{i+1}} \sum_{0\le j<i} x^{2^{j+1}} = {}$%
+ * %$\sum_{1\le i\le m} \rho^{2^i} \sum_{1\le j<i} x^{2^j} = {}$%
+ * %$\sum_{1\le i<m} \rho^{2^i} \sum_{1\le j<i} x^{2^j} + {}$%
+ * %$\rho^{2^m} \sum_{1\le j<m} x^{2^j}$%; and, somewhat miraculously,
+ * %$z^2 + z = \sum_{1\le i<m} \rho^{2^i} x + {}$%
+ * %$\rho \sum_{1\le i<m} x^{2^i} = x (\Tr(\rho) + \rho) + {}$%
+ * %$\rho (\Tr(x) + x) = x \Tr(\rho) + \rho \Tr(x)$%. Again,
+ * this gives us the root we want whenever %$\Tr(x) = 0$%.
+ *
+ * The loop below calculates %$w = \Tr(\rho)$% and %$z$% simultaneously,
+ * since the same powers of %$\rho$% are wanted in both calculations.
+ */
+
for (;;) {
rho = mprand(rho, m, fr, 0);
z = MP_ZERO;
d = z;
}
+ /* --- Check that we calculated the right answer --- *
+ *
+ * It should be correct; if it's not then maybe the ring we're working in
+ * isn't really a field.
+ */
+
t = gf_sqr(MP_NEW, d); t = gfreduce_do(r, t, t); t = gf_add(t, t, d);
if (!MP_EQ(t, x)) {
MP_DROP(d);
}
MP_DROP(t);
MP_DROP(x);
+
+ /* --- Pick a canonical root --- *
+ *
+ * The two roots are %$z$% and %$z + 1$%; pick the one with a zero
+ * scalar coefficient just for consistency's sake.
+ */
+
if (d) d->v[0] &= ~(mpw)1;
return (d);
}
/* --- @gfreduce_exp@ --- *
*
- * Arguments: @gfreduce *gr@ = pointer to reduction context
+ * Arguments: @const gfreduce *gr@ = pointer to reduction context
* @mp *d@ = fake destination
* @mp *a@ = base
* @mp *e@ = exponent
* Returns: Result, %$a^e \bmod m$%.
*/
-mp *gfreduce_exp(gfreduce *gr, mp *d, mp *a, mp *e)
+mp *gfreduce_exp(const gfreduce *gr, mp *d, mp *a, mp *e)
{
mp *x = MP_ONE;
mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW;
#ifdef TEST_RIG
-#define MP(x) mp_readstring(MP_NEW, #x, 0, 0)
-
static int vreduce(dstr *v)
{
mp *d = *(mp **)v[0].buf;
~mdw / catacomb / blobdiff