/* -*-c-*-
*
- * $Id: mpx.c,v 1.16 2003/05/16 09:09:24 mdw Exp $
+ * $Id: mpx.c,v 1.19 2004/04/03 03:29:40 mdw Exp $
*
* Low-level multiprecision arithmetic
*
/*----- Revision history --------------------------------------------------*
*
* $Log: mpx.c,v $
+ * Revision 1.19 2004/04/03 03:29:40 mdw
+ * Fix overrun in @mpx_lsr@.
+ *
+ * Revision 1.18 2004/04/01 12:50:09 mdw
+ * Add cyclic group abstraction, with test code. Separate off exponentation
+ * functions for better static linking. Fix a buttload of bugs on the way.
+ * Generally ensure that negative exponents do inversion correctly. Add
+ * table of standard prime-field subgroups. (Binary field subgroups are
+ * currently unimplemented but easy to add if anyone ever finds a good one.)
+ *
+ * Revision 1.17 2004/03/27 00:04:46 mdw
+ * Implement efficient reduction for pleasant-looking primes.
+ *
* Revision 1.16 2003/05/16 09:09:24 mdw
* Fix @mp_lsl2c@. Turns out to be surprisingly tricky.
*
/* --- Single bit shifting --- */
else if (n == 1) {
- mpw w = *av++ >> 1;
+ mpw w = av < avl ? *av++ >> 1 : 0;
while (av < avl) {
mpw t;
if (dv >= dvl)
void mpx_uaddn(mpw *dv, mpw *dvl, mpw n) { MPX_UADDN(dv, dvl, n); }
+/* --- @mpx_uaddnlsl@ --- *
+ *
+ * Arguments: @mpw *dv, *dvl@ = destination and first argument vector
+ * @mpw a@ = second argument
+ * @unsigned o@ = offset in bits
+ *
+ * Returns: ---
+ *
+ * Use: Computes %$d + 2^o a$%. If the result overflows then
+ * high-order bits are discarded, as usual. We must have
+ * @0 < o < MPW_BITS@.
+ */
+
+void mpx_uaddnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o)
+{
+ mpd x = (mpd)a << o;
+
+ while (x && dv < dvl) {
+ x += *dv;
+ *dv++ = MPW(x);
+ x >>= MPW_BITS;
+ }
+}
+
/* --- @mpx_usub@ --- *
*
* Arguments: @mpw *dv, *dvl@ = destination vector base and limit
void mpx_usubn(mpw *dv, mpw *dvl, mpw n) { MPX_USUBN(dv, dvl, n); }
+/* --- @mpx_uaddnlsl@ --- *
+ *
+ * Arguments: @mpw *dv, *dvl@ = destination and first argument vector
+ * @mpw a@ = second argument
+ * @unsigned o@ = offset in bits
+ *
+ * Returns: ---
+ *
+ * Use: Computes %$d + 2^o a$%. If the result overflows then
+ * high-order bits are discarded, as usual. We must have
+ * @0 < o < MPW_BITS@.
+ */
+
+void mpx_usubnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o)
+{
+ mpw b = a >> (MPW_BITS - o);
+ a <<= o;
+
+ if (dv < dvl) {
+ mpd x = (mpd)*dv - (mpd)a;
+ *dv++ = MPW(x);
+ if (x >> MPW_BITS)
+ b++;
+ MPX_USUBN(dv, dvl, b);
+ }
+}
+
/* --- @mpx_umul@ --- *
*
* Arguments: @mpw *dv, *dvl@ = destination vector base and limit
d = dvl[-1];
for (b = MPW_BITS / 2; b; b >>= 1) {
- if (d < (MPW_MAX >> b)) {
+ if (d <= (MPW_MAX >> b)) {
d <<= b;
norm += b;
}