- if (MP_LEN(a) > MP_LEN(b)) {
- const mp *t = a; a = b; b = t;
- }
- n = MP_LEN(mm->m);
-
- MP_MODIFY(d, 2 * n + 1);
- dv = d->v; dvl = d->vl;
- MPX_ZERO(dv, dvl);
- av = a->v; avl = a->vl;
- bv = b->v; bvl = b->vl;
- mv = mm->m->v; mvl = mm->m->vl;
- y = *bv;
+ a = MP_COPY(a);
+ b = MP_COPY(b);
+ MP_MODIFY(d, 2 * n + 1);
+ dv = d->v; dvl = d->vl;
+ MPX_ZERO(dv, dvl);
+ av = a->v; avl = a->vl;
+ bv = b->v; bvl = b->vl;
+ mv = mm->m->v; mvl = mm->m->vl;
+ y = *bv;
+
+ /* --- Montgomery multiplication phase --- */
+
+ i = 0;
+ while (i < n && av < avl) {
+ mpw x = *av++;
+ mpw u = MPW((*dv + x * y) * mm->mi);
+ MPX_UMLAN(dv, dvl, bv, bvl, x);
+ MPX_UMLAN(dv, dvl, mv, mvl, u);
+ dv++;
+ i++;
+ }