/* -*-c-*-
*
- * $Id: mpx.c,v 1.2 1999/11/13 01:50:59 mdw Exp $
+ * $Id: mpx.c,v 1.4 1999/11/17 18:04:09 mdw Exp $
*
* Low-level multiprecision arithmetic
*
/*----- Revision history --------------------------------------------------*
*
* $Log: mpx.c,v $
+ * 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.
*
/*----- Unsigned arithmetic -----------------------------------------------*/
+/* --- @mpx_2c@ --- *
+ *
+ * Arguments: @mpw *dv, *dvl@ = destination vector
+ * @const mpw *v, *vl@ = source vector
+ *
+ * Returns: ---
+ *
+ * Use: Calculates the two's complement of @v@.
+ */
+
+void mpx_2c(mpw *dv, mpw *dvl, const mpw *v, const mpw *vl)
+{
+ mpw c = 0;
+ while (dv < dvl && v < vl)
+ *dv++ = c = MPW(~*v++);
+ if (dv < dvl) {
+ if (c > MPW_MAX / 2)
+ c = MPW(~0);
+ while (dv < dvl)
+ *dv++ = c;
+ }
+ MPX_UADDN(dv, dvl, 1);
+}
+
/* --- @mpx_ucmp@ --- *
*
* Arguments: @const mpw *av, *avl@ = first argument vector base and limit
/* --- Work out the square at this point in the proceedings --- */
{
- mpw d = *dvv;
mpd x = (mpd)a * (mpd)a + *dvv;
*dvv++ = MPW(x);
c = MPW(x >> MPW_BITS);
avv++;
while (dvv < dvl && avv < avl) {
- mpw aa = *avv;
mpd x = (mpd)a * (mpd)*avv++;
mpd y = ((x << 1) & MPW_MAX) + c + *dvv;
c = (x >> (MPW_BITS - 1)) + (y >> MPW_BITS);
* requiring the dividend to be in the result position but it
* does make some sense really. The remainder must have
* headroom for at least two extra words. The scratch space
- * must be at least two words larger than twice the size of the
- * divisor.
+ * must be at least one word larger than the divisor.
*/
void mpx_udiv(mpw *qv, mpw *qvl, mpw *rv, mpw *rvl,
/* --- Normalize the dividend/remainder to match --- */
if (norm) {
- mpw *svvl = sv + (dvl - dv) + 1;
mpx_lsl(rv, rvl, rv, rvl, norm);
- mpx_lsl(sv, svvl, dv, dvl, norm);
+ mpx_lsl(sv, svl, dv, dvl, norm);
dv = sv;
- sv = svvl;
- dvl = svvl;
+ dvl = svl;
MPX_SHRINK(dv, dvl);
}
{
mpw *svv;
const mpw *dvv;
- mpw c = 0;
+ mpw mc = 0, sc = 0;
- /* --- Calculate the size of the chunk --- */
+ /* --- Calculate the size of the chunk --- *
+ *
+ * This does the whole job of calculating @r >> scale - qd@.
+ */
- for (svv = sv, dvv = dv; dvv < dvl; svv++, dvv++) {
- mpd x = (mpd)*dvv * (mpd)q + c;
+ for (svv = rv + scale, dvv = dv;
+ dvv < dvl && svv < rvl;
+ svv++, dvv++) {
+ mpd x = (mpd)*dvv * (mpd)q + mc;
+ mc = x >> MPW_BITS;
+ x = (mpd)*svv - MPW(x) - sc;
*svv = MPW(x);
- c = x >> MPW_BITS;
+ if (x >> MPW_BITS)
+ sc = 1;
+ else
+ sc = 0;
+ }
+
+ if (svv < rvl) {
+ mpd x = (mpd)*svv - mc - sc;
+ *svv++ = MPW(x);
+ if (x >> MPW_BITS)
+ sc = MPW_MAX;
+ else
+ sc = 0;
+ while (svv < rvl)
+ *svv++ = sc;
}
- if (c)
- *svv++ = c;
- /* --- Now make sure that we can cope with the difference --- *
+ /* --- Fix if the quotient was too large --- *
*
- * Take advantage of the fact that subtraction works two's-
- * complement.
+ * This doesn't seem to happen very often.
*/
- mpx_usub(rv + scale, rvl, rv + scale, rvl, sv, svv);
if (rvl[-1] > MPW_MAX / 2) {
mpx_uadd(rv + scale, rvl, rv + scale, rvl, dv, dvl);
q--;