From 51f5bbe0bfe8a344c3cc1f5a9fddf885fa75d057 Mon Sep 17 00:00:00 2001 From: Mark Wooding Date: Fri, 9 Aug 2013 10:30:48 +0100 Subject: [PATCH] math/gfreduce.c: Refactor and document. There was an unfortunate amount of code duplication in the old code, particularly around the issuing of LSR instruction sequences. Gather the relevant state into a structure, and split the repeated code into separate functions. Also add a bunch of commentary explaining the underlying theory. --- math/gfreduce.c | 239 ++++++++++++++++++++++++++++++++++++-------------------- math/t/gfreduce | 6 ++ 2 files changed, 160 insertions(+), 85 deletions(-) diff --git a/math/gfreduce.c b/math/gfreduce.c index aed3f22..96e0ced 100644 --- a/math/gfreduce.c +++ b/math/gfreduce.c @@ -47,26 +47,22 @@ DA_DECL(instr_v, gfreduce_instr); * * Let's face it, @gfx_div@ sucks. It works (I hope), but it's not in any * sense fast. Here, we do efficient reduction modulo sparse polynomials. - * - * Suppose we have a polynomial @X@ we're trying to reduce mod @P@. If we - * take the topmost nonzero word of @X@, call it @w@, then we can eliminate - * it by subtracting off %$w P x^{k}$% for an appropriate value of @k@. The - * trick is in observing that if @P@ is sparse we can do this multiplication - * and subtraction efficiently, just by XORing appropriate shifts of @w@ into - * @X@. - * - * The first tricky bit is in working out when to stop. I'll use eight-bit - * words to demonstrate what I'm talking about. - * - * xxxxxx xxxxxxxx xxxxxxxx xxxxxxxx xxxxxxxx xxxxxxxx - * 001ppppp pppppppp pppppppp pppppppp - * || - * |<------------ bp ------------->| - * |<------------ nw --------------->| - * - * The trick of taking whole words off of @X@ stops working when there are - * only @nw@ words left. Then we have to mask off the bottom bits of @w@ - * before continuing. + * (It works for arbitrary polynomials, but isn't efficient for dense ones.) + * + * Suppose that %$p(x) = x^n + p'(x) = \sum_{0\le iiv; \ + size_t _i = DA_LEN(_iv); \ + \ + DA_ENSURE(_iv, 1); \ + DA(_iv)[_i].op = (op_); \ + DA(_iv)[_i].arg = (arg_); \ + DA_EXTEND(_iv, 1); \ +} while (0) + +static void emit_load(struct gen *g, size_t w) +{ + INSTR(g, GFRI_LOAD, w); + g->f |= f_load; + g->w = w; +} + +static void emit_right_shifts(struct gen *g) +{ + gfreduce_instr *ip; + size_t i, wl; + + /* --- Close off the current word --- * + * + * If we shifted into this current word with a nonzero bit offset, then + * we'll also need to arrange to perform a sequence of right shifts into + * the following word, which we might as well do by scanning the + * instruction sequence (which starts at @wi@). + * + * Either way, we leave a @LOAD@ unmatched if there was one before, in the + * hope that callers have an easier time; @g->w@ is updated to reflect the + * currently open word. + */ + + if (!(g->f & f_lsr)) + return; + + wl = DA_LEN(&g->iv); + INSTR(g, GFRI_STORE, g->w); + emit_load(g, g->w - 1); + for (i = g->wi; i < wl; i++) { + ip = &DA(&g->iv)[i]; + assert(ip->op == GFRI_LSL); + if (ip->arg) + INSTR(g, GFRI_LSR, MPW_BITS - ip->arg); + } + g->f &= ~f_lsr; +} + +static void ensure_loaded(struct gen *g, size_t w) +{ + if (!(g->f & f_load)) { + emit_load(g, w); + g->wi = DA_LEN(&g->iv); + } else if (w != g->w) { + emit_right_shifts(g); + if (w != g->w) { + INSTR(g, GFRI_STORE, g->w); + emit_load(g, w); + } + g->wi = DA_LEN(&g->iv); + } +} + void gfreduce_create(gfreduce *r, mp *p) { - instr_v iv = DA_INIT; + struct gen g = { 0, DA_INIT }; unsigned long d; unsigned dw; mpscan sc; unsigned long i; - gfreduce_instr *ip; - unsigned f = 0; - size_t w, ww, wi, wl, ll, bb; + size_t w, bb; /* --- Sort out the easy stuff --- */ @@ -103,79 +171,80 @@ void gfreduce_create(gfreduce *r, mp *p) } r->p = mp_copy(p); - /* --- Stash a new instruction --- */ - -#define INSTR(op_, arg_) do { \ - DA_ENSURE(&iv, 1); \ - DA(&iv)[DA_LEN(&iv)].op = (op_); \ - DA(&iv)[DA_LEN(&iv)].arg = (arg_); \ - DA_EXTEND(&iv, 1); \ -} while (0) - -#define f_lsr 1u + /* --- How this works --- * + * + * The instruction sequence is run with two ambient parameters: a pointer + * (usually) just past the most significant word of the polynomial to be + * reduced; and a word %$z$% which is the multiple of %$p'$% we are meant + * to add. + * + * The sequence visits each word of the polynomial at most once. Suppose + * %$u = z x^{w N} + u'$%; our pointer points just past the end of %$u'$%. + * Word %$I$% of %$u'$% will be affected by modulus bits %$p_i$% where + * %$(N - I - 1) w + 1 \le i \le (N - I + 1) w - 1$%, so %$p_i$% affects + * word %$I = \lceil (n - i + 1)/w \rceil$% and (if %$i$% is not a multiple + * of %$w$%) also word %$I - 1$%. + * + * We have four instructions: @LOAD@ reads a specified word of %$u$% into an + * accumulator, and @STORE@ stores it back (we'll always store back to the + * same word we most recently read, but this isn't a requirement); and + * @LSL@ and @LSR@, which XOR in appropriately shifted copies of %$z$% into + * the accumulator. So a typical program will contain sequences of @LSR@ + * and @LSL@ instructions sandwiched between @LOAD@/@STORE@ pairs. + * + * We do a single right-to-left pass across %$p$%. + */ - w = (d + MPW_BITS - 1)/MPW_BITS; - INSTR(GFRI_LOAD, w); - wi = DA_LEN(&iv); - f = 0; - ll = 0; bb = MPW_BITS - dw; + for (i = 0, mp_scan(&sc, p); mp_step(&sc) && i < d; i++) { if (!mp_bit(&sc)) continue; - ww = (d - i + MPW_BITS - 1)/MPW_BITS; - if (ww != w) { - wl = DA_LEN(&iv); - INSTR(GFRI_STORE, w); - if (!ll) - ll = DA_LEN(&iv); - if (!(f & f_lsr)) - INSTR(GFRI_LOAD, ww); - else { - INSTR(GFRI_LOAD, w - 1); - for (; wi < wl; wi++) { - ip = &DA(&iv)[wi]; - assert(ip->op == GFRI_LSL); - if (ip->arg) - INSTR(GFRI_LSR, MPW_BITS - ip->arg); - } - if (w - 1 != ww) { - INSTR(GFRI_STORE, w - 1); - INSTR(GFRI_LOAD, ww); - } - f &= ~f_lsr; - } - w = ww; - wi = DA_LEN(&iv); - } - INSTR(GFRI_LSL, (bb + i)%MPW_BITS); + + /* --- We've found a set bit, so work out which word it affects --- * + * + * In general, a bit affects two words: it needs to be shifted left into + * one, and shifted right into the next. We find the former here. + */ + + w = (d - i + MPW_BITS - 1)/MPW_BITS; + + /* --- Concentrate on the appropriate word --- */ + + ensure_loaded(&g, w); + + /* --- Accumulate a new @LSL@ instruction --- * + * + * If this was a nonzero shift, then we'll need to arrange to do right + * shifts into the following word. + */ + + INSTR(&g, GFRI_LSL, (bb + i)%MPW_BITS); if ((bb + i)%MPW_BITS) - f |= f_lsr; - } - wl = DA_LEN(&iv); - INSTR(GFRI_STORE, w); - if (!ll) - ll = DA_LEN(&iv); - if (f & f_lsr) { - INSTR(GFRI_LOAD, w - 1); - for (; wi < wl; wi++) { - ip = &DA(&iv)[wi]; - assert(ip->op == GFRI_LSL); - if (ip->arg) - INSTR(GFRI_LSR, MPW_BITS - ip->arg); - } - INSTR(GFRI_STORE, w - 1); + g.f |= f_lsr; } -#undef INSTR + /* --- Wrapping up --- * + * + * We probably need a final @STORE@, and maybe a sequence of right shifts. + */ - r->in = DA_LEN(&iv); + if (g.f & f_load) { + emit_right_shifts(&g); + INSTR(&g, GFRI_STORE, g.w); + } + + r->in = DA_LEN(&g.iv); r->iv = xmalloc(r->in * sizeof(gfreduce_instr)); - r->liv = r->iv + ll; - memcpy(r->iv, DA(&iv), r->in * sizeof(gfreduce_instr)); - DA_DESTROY(&iv); + memcpy(r->iv, DA(&g.iv), r->in * sizeof(gfreduce_instr)); + DA_DESTROY(&g.iv); } +#undef INSTR + +#undef f_lsr +#undef f_load + /* --- @gfreduce_destroy@ --- * * * Arguments: @gfreduce *r@ = structure to free diff --git a/math/t/gfreduce b/math/t/gfreduce index f712316..ecece92 100644 --- a/math/t/gfreduce +++ b/math/t/gfreduce @@ -13,6 +13,12 @@ reduce { 0xb2ca471b0867d5fae2e4f27a2d2706da 0xf254423fef93d5d7a76ecf22c656c1352c53257875945d33 0x582f783fc210f72814780e69b0bd29ff; + + # --- Bugs discovered --- + + 0x20000000000000000000000000000000000000004000000000000000001 + 0x110414154054140445011511541540514401111414505115044104145451505001151441450000541004550554154040500000411050400055041 + 0x1bf878a39fbee9cf20a2f6f41eadda756518f1669b2c7d8f9234965b6b3; } modexp { -- 2.11.0