/* -*-c-*-
*
- * $Id: mpx.c,v 1.16 2003/05/16 09:09:24 mdw Exp $
+ * $Id: mpx.c,v 1.20 2004/04/08 01:36:15 mdw Exp $
*
* Low-level multiprecision arithmetic
*
* MA 02111-1307, USA.
*/
-/*----- Revision history --------------------------------------------------*
- *
- * $Log: mpx.c,v $
- * Revision 1.16 2003/05/16 09:09:24 mdw
- * Fix @mp_lsl2c@. Turns out to be surprisingly tricky.
- *
- * Revision 1.15 2002/10/20 01:12:31 mdw
- * Two's complement I/O fixes.
- *
- * Revision 1.14 2002/10/19 18:55:08 mdw
- * Fix overflows in shift primitives.
- *
- * Revision 1.13 2002/10/19 17:56:50 mdw
- * Fix bit operations. Test them (a bit) better.
- *
- * Revision 1.12 2002/10/06 22:52:50 mdw
- * Pile of changes for supporting two's complement properly.
- *
- * Revision 1.11 2001/04/03 19:36:05 mdw
- * Add some simple bitwise operations so that Perl can use them.
- *
- * Revision 1.10 2000/10/08 12:06:12 mdw
- * Provide @mpx_ueq@ for rapidly testing equality of two integers.
- *
- * Revision 1.9 2000/06/26 07:52:50 mdw
- * Portability fix for the bug fix.
- *
- * Revision 1.8 2000/06/25 12:59:02 mdw
- * (mpx_udiv): Fix bug in quotient digit estimation.
- *
- * Revision 1.7 1999/12/22 15:49:07 mdw
- * New function for division by a small integer.
- *
- * Revision 1.6 1999/11/20 22:43:44 mdw
- * Integrate testing for MPX routines.
- *
- * Revision 1.5 1999/11/20 22:23:27 mdw
- * Add function versions of some low-level macros with wider use.
- *
- * Revision 1.4 1999/11/17 18:04:09 mdw
- * Add two's-complement functionality. Improve mpx_udiv a little by
- * performing the multiplication of the divisor by q with the subtraction
- * from r.
- *
- * Revision 1.3 1999/11/13 01:57:31 mdw
- * Remove stray debugging code.
- *
- * Revision 1.2 1999/11/13 01:50:59 mdw
- * Multiprecision routines finished and tested.
- *
- * Revision 1.1 1999/09/03 08:41:12 mdw
- * Initial import.
- *
- */
-
/*----- Header files ------------------------------------------------------*/
#include <assert.h>
/* --- 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;
}