/* -*-c-*-
*
- * $Id: mpx.c,v 1.16 2003/05/16 09:09:24 mdw Exp $
+ * $Id: mpx.c,v 1.17 2004/03/27 00:04:46 mdw Exp $
*
* Low-level multiprecision arithmetic
*
/*----- Revision history --------------------------------------------------*
*
* $Log: mpx.c,v $
+ * 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.
*
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