Updated from ../=mpn/gmp-1.910
authorroland <roland>
Mon, 4 Dec 1995 18:37:45 +0000 (18:37 +0000)
committerroland <roland>
Mon, 4 Dec 1995 18:37:45 +0000 (18:37 +0000)
stdlib/gmp-impl.h
stdlib/gmp.h
stdlib/longlong.h
sysdeps/m88k/add_n.s [new file with mode: 0644]
sysdeps/m88k/mul_1.s [new file with mode: 0644]
sysdeps/m88k/sub_n.s [new file with mode: 0644]

index 0d2a8fc..2f0956d 100644 (file)
@@ -179,24 +179,22 @@ void _mp_default_free ();
    strings in base 2..36.  */
 struct bases
 {
-  /* Number of digits in the conversion base that always fits in
-     an mp_limb.  For example, for base 10 this is 10, since
-     2**32 = 4294967296 has ten digits.  */
+  /* Number of digits in the conversion base that always fits in an mp_limb.
+     For example, for base 10 on a machine where a mp_limb has 32 bits this
+     is 9, since 10**9 is the largest number that fits into a mp_limb.  */
   int chars_per_limb;
 
   /* log(2)/log(conversion_base) */
   float chars_per_bit_exactly;
 
-  /* big_base is conversion_base**chars_per_limb, i.e. the biggest
-     number that fits a word, built by factors of conversion_base.
-     Exception: For 2, 4, 8, etc, big_base is log2(base), i.e. the
-     number of bits used to represent each digit in the base.  */
+  /* base**chars_per_limb, i.e. the biggest number that fits a word, built by
+     factors of base.  Exception: For 2, 4, 8, etc, big_base is log2(base),
+     i.e. the number of bits used to represent each digit in the base.  */
   mp_limb big_base;
 
-  /* big_base_inverted is a BITS_PER_MP_LIMB bit approximation to
-     1/big_base, represented as a fixed-point number.  Instead of
-     dividing by big_base an application can choose to multiply
-     by big_base_inverted.  */
+  /* A BITS_PER_MP_LIMB bit approximation to 1/big_base, represented as a
+     fixed-point number.  Instead of dividing by big_base an application can
+     choose to multiply by big_base_inverted.  */
   mp_limb big_base_inverted;
 };
 
index 0b2cb29..2437799 100644 (file)
@@ -23,6 +23,7 @@ the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
 #ifndef __GNU_MP__
 #define __need_size_t
 #include <stddef.h>
+#undef __need_size_t
 
 #if defined (__STDC__)
 #define __gmp_const const
@@ -40,7 +41,7 @@ the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
 typedef unsigned int           mp_limb;
 typedef int                    mp_limb_signed;
 #else
-#if _LONG_LONG_LIMB
+#ifdef _LONG_LONG_LIMB
 typedef unsigned long long int mp_limb;
 typedef long long int          mp_limb_signed;
 #else
@@ -110,11 +111,11 @@ typedef __mpq_struct mpq_t[1];
 
 typedef struct
 {
-  mp_size_t alloc;             /* Number of *limbs* allocated and pointed
-                                  to by the D field.  */
   mp_size_t prec;              /* Max precision, in number of `mp_limb's.
                                   Set by mpf_init and modified by
-                                  mpf_set_prec.  */
+                                  mpf_set_prec.  The area pointed to
+                                  by the `d' field contains `prec' + 1
+                                  limbs.  */
   mp_size_t size;              /* abs(SIZE) is the number of limbs
                                   the last field points to.  If SIZE
                                   is negative this is a negative
@@ -127,7 +128,7 @@ typedef struct
 typedef __mpf_struct mpf_t[1];
 
 /* Types for function declarations in gmp files.  */
-/* ??? Should not pollute user name space ??? */
+/* ??? Should not pollute user name space with these ??? */
 typedef __gmp_const __mpz_struct *mpz_srcptr;
 typedef __mpz_struct *mpz_ptr;
 typedef __gmp_const __mpf_struct *mpf_srcptr;
index bbb92e3..e52bf32 100644 (file)
@@ -139,6 +139,7 @@ the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
     __asm__ ("clz %0,%1"                                               \
             : "=r" ((USItype)(count))                                  \
             : "r" ((USItype)(x)))
+#define COUNT_LEADING_ZEROS_0 32
 #endif /* __a29k__ */
 
 #if defined (__alpha__) && W_TYPE_SIZE == 64
@@ -298,9 +299,9 @@ extern UDItype __udiv_qrnnd ();
           struct {USItype __h, __l;} __i;                              \
          } __xx;                                                       \
     __asm__ ("xmpyu %1,%2,%0"                                          \
-            : "=fx" (__xx.__ll)                                        \
-            : "fx" ((USItype)(u)),                                     \
-              "fx" ((USItype)(v)));                                    \
+            : "=*f" (__xx.__ll)                                        \
+            : "*f" ((USItype)(u)),                                     \
+              "*f" ((USItype)(v)));                                    \
     (wh) = __xx.__i.__h;                                               \
     (wl) = __xx.__i.__l;                                               \
   } while (0)
@@ -339,7 +340,7 @@ extern USItype __udiv_qrnnd ();
        sub             %0,%1,%0                ; Subtract it.
        " : "=r" (count), "=r" (__tmp) : "1" (x));                      \
   } while (0)
-#endif
+#endif /* hppa */
 
 #if (defined (__i370__) || defined (__mvs__)) && W_TYPE_SIZE == 32
 #define umul_ppmm(xh, xl, m0, m1) \
@@ -431,7 +432,29 @@ extern USItype __udiv_qrnnd ();
 #endif
 #endif /* 80x86 */
 
+#if defined (__i860__) && W_TYPE_SIZE == 32
+#define rshift_rhlc(r,h,l,c) \
+  __asm__ ("shr %3,r0,r0\;shrd %1,%2,%0"                               \
+          "=r" (r) : "r" (h), "r" (l), "rn" (c))
+#endif /* i860 */
+
 #if defined (__i960__) && W_TYPE_SIZE == 32
+#define add_ssaaaa(sh, sl, ah, al, bh, bl) \
+  __asm__ ("cmpo 1,0\;addc %5,%4,%1\;addc %3,%2,%0"                    \
+          : "=r" ((USItype)(sh)),                                      \
+            "=&r" ((USItype)(sl))                                      \
+          : "%dI" ((USItype)(ah)),                                     \
+            "dI" ((USItype)(bh)),                                      \
+            "%dI" ((USItype)(al)),                                     \
+            "dI" ((USItype)(bl)))
+#define sub_ddmmss(sh, sl, ah, al, bh, bl) \
+  __asm__ ("cmpo 0,0\;subc %5,%4,%1\;subc %3,%2,%0"                    \
+          : "=r" ((USItype)(sh)),                                      \
+            "=&r" ((USItype)(sl))                                      \
+          : "dI" ((USItype)(ah)),                                      \
+            "dI" ((USItype)(bh)),                                      \
+            "dI" ((USItype)(al)),                                      \
+            "dI" ((USItype)(bl)))
 #define umul_ppmm(w1, w0, u, v) \
   ({union {UDItype __ll;                                               \
           struct {USItype __l, __h;} __i;                              \
@@ -448,7 +471,39 @@ extern USItype __udiv_qrnnd ();
             : "%dI" ((USItype)(u)),                                    \
               "dI" ((USItype)(v)));                                    \
     __w; })  
-#endif /* __i960__ */
+#define udiv_qrnnd(q, r, nh, nl, d) \
+  do {                                                                 \
+    union {UDItype __ll;                                               \
+          struct {USItype __l, __h;} __i;                              \
+         } __nn;                                                       \
+    __nn.__i.__h = (nh); __nn.__i.__l = (nl);                          \
+    __asm__ ("ediv %d,%n,%0"                                           \
+          : "=d" (__rq.__ll)                                           \
+          : "dI" (__nn.__ll),                                          \
+            "dI" ((USItype)(d)));                                      \
+    (r) = __rq.__i.__l; (q) = __rq.__i.__h;                            \
+  } while (0)
+#define count_leading_zeros(count, x) \
+  do {                                                                 \
+    USItype __cbtmp;                                                   \
+    __asm__ ("scanbit %1,%0"                                           \
+            : "=r" (__cbtmp)                                           \
+            : "r" ((USItype)(x)));                                     \
+    (count) = __cbtmp ^ 31;                                            \
+  } while (0)
+#define COUNT_LEADING_ZEROS_0 (-32) /* sic */
+#if defined (__i960mx)         /* what is the proper symbol to test??? */
+#define rshift_rhlc(r,h,l,c) \
+  do {                                                                 \
+    union {UDItype __ll;                                               \
+          struct {USItype __l, __h;} __i;                              \
+         } __nn;                                                       \
+    __nn.__i.__h = (h); __nn.__i.__l = (l);                            \
+    __asm__ ("shre %2,%1,%0"                                           \
+            : "=d" (r) : "dI" (__nn.__ll), "dI" (c));                  \
+  }
+#endif /* i960mx */
+#endif /* i960 */
 
 #if (defined (__mc68000__) || defined (__mc68020__) || defined (__NeXT__) || defined(mc68020)) && W_TYPE_SIZE == 32
 #define add_ssaaaa(sh, sl, ah, al, bh, bl) \
@@ -469,7 +524,7 @@ extern USItype __udiv_qrnnd ();
             "d" ((USItype)(bh)),                                       \
             "1" ((USItype)(al)),                                       \
             "g" ((USItype)(bl)))
-#if (defined (__mc68020__) || defined (__NeXT__) || defined(mc68020)) && W_TYPE_SIZE == 32
+#if (defined (__mc68020__) || defined (__NeXT__) || defined(mc68020))
 #define umul_ppmm(w1, w0, u, v) \
   __asm__ ("mulu%.l %3,%1:%0"                                          \
           : "=d" ((USItype)(w0)),                                      \
@@ -496,8 +551,9 @@ extern USItype __udiv_qrnnd ();
   __asm__ ("bfffo %1{%b2:%b2},%0"                                      \
           : "=d" ((USItype)(count))                                    \
           : "od" ((USItype)(x)), "n" (0))
+#define COUNT_LEADING_ZEROS_0 32
 #else /* not mc68020 */
-#define umul_ppmmxx(xh, xl, a, b) \
+#define umul_ppmm(xh, xl, a, b) \
   do { USItype __umul_tmp1, __umul_tmp2;                               \
        __asm__ ("| Inlined umul_ppmm
        move%.l %5,%3
@@ -557,6 +613,7 @@ extern USItype __udiv_qrnnd ();
             : "r" ((USItype)(x)));                                     \
     (count) = __cbtmp ^ 31;                                            \
   } while (0)
+#define COUNT_LEADING_ZEROS_0 63 /* sic */
 #if defined (__m88110__)
 #define umul_ppmm(wh, wl, u, v) \
   do {                                                                 \
@@ -738,6 +795,7 @@ extern USItype __udiv_qrnnd ();
   __asm__ ("{cntlz|cntlzw} %0,%1"                                      \
           : "=r" ((USItype)(count))                                    \
           : "r" ((USItype)(x)))
+#define COUNT_LEADING_ZEROS_0 32
 #if defined (_ARCH_PPC)
 #define umul_ppmm(ph, pl, m0, m1) \
   do {                                                                 \
@@ -887,7 +945,7 @@ extern USItype __udiv_qrnnd ();
        (count) += 16;                                                  \
       }                                                                        \
   } while (0)
-#endif
+#endif /* RT/ROMP */
 
 #if defined (__sh2__) && W_TYPE_SIZE == 32
 #define umul_ppmm(w1, w0, u, v) \
@@ -1154,20 +1212,6 @@ extern USItype __udiv_qrnnd ();
     (xh) += ((((signed int) __m0 >> 15) & __m1)                                \
             + (((signed int) __m1 >> 15) & __m0));                     \
   } while (0)
-#define umul_ppmm_off(xh, xl, m0, m1) \
-  do {                                                                 \
-    union {long int __ll;                                              \
-          struct {unsigned int __h, __l;} __i;                         \
-         } __xx;                                                       \
-    __asm__ ("mult     %S0,%H3"                                        \
-            : "=r" (__xx.__i.__h),                                     \
-              "=r" (__xx.__i.__l)                                      \
-            : "%1" (m0),                                               \
-              "rQR" (m1));                                             \
-    (xh) = __xx.__i.__h + ((((signed int) m0 >> 15) & m1)              \
-                          + (((signed int) m1 >> 15) & m0));           \
-    (xl) = __xx.__i.__l;                                               \
-  } while (0)
 #endif /* __z8000__ */
 
 #endif /* __GNUC__ */
diff --git a/sysdeps/m88k/add_n.s b/sysdeps/m88k/add_n.s
new file mode 100644 (file)
index 0000000..7e4cccc
--- /dev/null
@@ -0,0 +1,103 @@
+; mc88100 __mpn_add -- Add two limb vectors of the same length > 0 and store
+; sum in a third limb vector.
+
+; Copyright (C) 1992, 1994 Free Software Foundation, Inc.
+
+; This file is part of the GNU MP Library.
+
+; The GNU MP Library is free software; you can redistribute it and/or modify
+; it under the terms of the GNU Library General Public License as published by
+; the Free Software Foundation; either version 2 of the License, or (at your
+; option) any later version.
+
+; The GNU MP Library is distributed in the hope that it will be useful, but
+; WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+; or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
+; License for more details.
+
+; You should have received a copy of the GNU Library General Public License
+; along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+; the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+
+
+; INPUT PARAMETERS
+; res_ptr      r2
+; s1_ptr       r3
+; s2_ptr       r4
+; size         r5
+
+; This code has been optimized to run one instruction per clock, avoiding
+; load stalls and writeback contention.  As a result, the instruction
+; order is not always natural.
+
+; The speed is about 4.6 clocks/limb + 18 clocks/limb-vector on an 88100,
+; but on the 88110, it seems to run much slower, 6.6 clocks/limb.
+
+       text
+       align    16
+       global   ___mpn_add_n
+___mpn_add_n:
+       ld      r6,r3,0                 ; read first limb from s1_ptr
+       extu    r10,r5,3
+       ld      r7,r4,0                 ; read first limb from s2_ptr
+
+       subu.co r5,r0,r5                ; (clear carry as side effect)
+       mak     r5,r5,3<4>
+       bcnd    eq0,r5,Lzero
+
+       or      r12,r0,lo16(Lbase)
+       or.u    r12,r12,hi16(Lbase)
+       addu    r12,r12,r5              ; r12 is address for entering in loop
+
+       extu    r5,r5,2                 ; divide by 4
+       subu    r2,r2,r5                ; adjust res_ptr
+       subu    r3,r3,r5                ; adjust s1_ptr
+       subu    r4,r4,r5                ; adjust s2_ptr
+
+       or      r8,r6,r0
+
+       jmp.n   r12
+        or     r9,r7,r0
+
+Loop:  addu    r3,r3,32
+       st      r8,r2,28
+       addu    r4,r4,32
+       ld      r6,r3,0
+       addu    r2,r2,32
+       ld      r7,r4,0
+Lzero: subu    r10,r10,1               ; add 0 + 8r limbs (adj loop cnt)
+Lbase: ld      r8,r3,4
+       addu.cio r6,r6,r7
+       ld      r9,r4,4
+       st      r6,r2,0
+       ld      r6,r3,8                 ; add 7 + 8r limbs
+       addu.cio r8,r8,r9
+       ld      r7,r4,8
+       st      r8,r2,4
+       ld      r8,r3,12                ; add 6 + 8r limbs
+       addu.cio r6,r6,r7
+       ld      r9,r4,12
+       st      r6,r2,8
+       ld      r6,r3,16                ; add 5 + 8r limbs
+       addu.cio r8,r8,r9
+       ld      r7,r4,16
+       st      r8,r2,12
+       ld      r8,r3,20                ; add 4 + 8r limbs
+       addu.cio r6,r6,r7
+       ld      r9,r4,20
+       st      r6,r2,16
+       ld      r6,r3,24                ; add 3 + 8r limbs
+       addu.cio r8,r8,r9
+       ld      r7,r4,24
+       st      r8,r2,20
+       ld      r8,r3,28                ; add 2 + 8r limbs
+       addu.cio r6,r6,r7
+       ld      r9,r4,28
+       st      r6,r2,24
+       bcnd.n  ne0,r10,Loop            ; add 1 + 8r limbs
+        addu.cio r8,r8,r9
+
+       st      r8,r2,28                ; store most significant limb
+
+       jmp.n    r1
+        addu.ci r2,r0,r0               ; return carry-out from most sign. limb
diff --git a/sysdeps/m88k/mul_1.s b/sysdeps/m88k/mul_1.s
new file mode 100644 (file)
index 0000000..35c238d
--- /dev/null
@@ -0,0 +1,128 @@
+; mc88100 __mpn_mul_1 -- Multiply a limb vector with a single limb and
+; store the product in a second limb vector.
+
+; Copyright (C) 1992, 1994 Free Software Foundation, Inc.
+
+; This file is part of the GNU MP Library.
+
+; The GNU MP Library is free software; you can redistribute it and/or modify
+; it under the terms of the GNU Library General Public License as published by
+; the Free Software Foundation; either version 2 of the License, or (at your
+; option) any later version.
+
+; The GNU MP Library is distributed in the hope that it will be useful, but
+; WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+; or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
+; License for more details.
+
+; You should have received a copy of the GNU Library General Public License
+; along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+; the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+
+
+; INPUT PARAMETERS
+; res_ptr      r2
+; s1_ptr       r3
+; size         r4
+; s2_limb      r5
+
+; Common overhead is about 11 cycles/invocation.
+
+; The speed for S2_LIMB >= 0x10000 is approximately 21 cycles/limb.  (The
+; pipeline stalls 2 cycles due to WB contention.)
+
+; The speed for S2_LIMB < 0x10000 is approximately 16 cycles/limb.  (The
+; pipeline stalls 2 cycles due to WB contention and 1 cycle due to latency.)
+
+; To enhance speed:
+; 1. Unroll main loop 4-8 times.
+; 2. Schedule code to avoid WB contention.  It might be tempting to move the
+;    ld instruction in the loops down to save 2 cycles (less WB contention),
+;    but that looses because the ultimate value will be read from outside
+;    the allocated space.  But if we handle the ultimate multiplication in
+;    the tail, we can do this.
+; 3. Make the multiplication with less instructions.  I think the code for
+;    (S2_LIMB >= 0x10000) is not minimal.
+; With these techniques the (S2_LIMB >= 0x10000) case would run in 17 or
+; less cycles/limb; the (S2_LIMB < 0x10000) case would run in 11
+; cycles/limb.  (Assuming infinite unrolling.)
+
+       text
+       align    16
+       global   ___mpn_mul_1
+___mpn_mul_1:
+
+       ; Make S1_PTR and RES_PTR point at the end of their blocks
+       ; and negate SIZE.
+       lda      r3,r3[r4]
+       lda      r6,r2[r4]              ; RES_PTR in r6 since r2 is retval
+       subu     r4,r0,r4
+
+       addu.co  r2,r0,r0               ; r2 = cy = 0
+       ld       r9,r3[r4]
+       mask     r7,r5,0xffff           ; r7 = lo(S2_LIMB)
+       extu     r8,r5,16               ; r8 = hi(S2_LIMB)
+       bcnd.n   eq0,r8,Lsmall          ; jump if (hi(S2_LIMB) == 0)
+        subu    r6,r6,4
+
+; General code for any value of S2_LIMB.
+
+       ; Make a stack frame and save r25 and r26
+       subu     r31,r31,16
+       st.d     r25,r31,8
+
+       ; Enter the loop in the middle
+       br.n    L1
+       addu     r4,r4,1
+
+Loop:
+       ld       r9,r3[r4]
+       st       r26,r6[r4]
+; bcnd ne0,r0,0                        ; bubble
+       addu     r4,r4,1
+L1:    mul      r26,r9,r5              ; low word of product   mul_1   WB ld
+       mask     r12,r9,0xffff          ; r12 = lo(s1_limb)     mask_1
+       mul      r11,r12,r7             ; r11 =  prod_0         mul_2   WB mask_1
+       mul      r10,r12,r8             ; r10 = prod_1a         mul_3
+       extu     r13,r9,16              ; r13 = hi(s1_limb)     extu_1  WB mul_1
+       mul      r12,r13,r7             ; r12 = prod_1b         mul_4   WB extu_1
+       mul      r25,r13,r8             ; r25  = prod_2         mul_5   WB mul_2
+       extu     r11,r11,16             ; r11 = hi(prod_0)      extu_2  WB mul_3
+       addu     r10,r10,r11            ;                       addu_1  WB extu_2
+; bcnd ne0,r0,0                        ; bubble                        WB addu_1
+       addu.co  r10,r10,r12            ;                               WB mul_4
+       mask.u   r10,r10,0xffff         ; move the 16 most significant bits...
+       addu.ci  r10,r10,r0             ; ...to the low half of the word...
+       rot      r10,r10,16             ; ...and put carry in pos 16.
+       addu.co  r26,r26,r2             ; add old carry limb
+       bcnd.n   ne0,r4,Loop
+        addu.ci r2,r25,r10             ; compute new carry limb
+
+       st       r26,r6[r4]
+       ld.d     r25,r31,8
+       jmp.n    r1
+        addu    r31,r31,16
+
+; Fast code for S2_LIMB < 0x10000
+Lsmall:
+       ; Enter the loop in the middle
+       br.n    SL1
+       addu     r4,r4,1
+
+SLoop:
+       ld       r9,r3[r4]              ;
+       st       r8,r6[r4]              ;
+       addu     r4,r4,1                ;
+SL1:   mul      r8,r9,r5               ; low word of product
+       mask     r12,r9,0xffff          ; r12 = lo(s1_limb)
+       extu     r13,r9,16              ; r13 = hi(s1_limb)
+       mul      r11,r12,r7             ; r11 =  prod_0
+       mul      r12,r13,r7             ; r12 = prod_1b
+       addu.cio r8,r8,r2               ; add old carry limb
+       extu     r10,r11,16             ; r11 = hi(prod_0)
+       addu     r10,r10,r12            ;
+       bcnd.n   ne0,r4,SLoop
+       extu     r2,r10,16              ; r2 = new carry limb
+
+       jmp.n    r1
+       st       r8,r6[r4]
diff --git a/sysdeps/m88k/sub_n.s b/sysdeps/m88k/sub_n.s
new file mode 100644 (file)
index 0000000..3963cd5
--- /dev/null
@@ -0,0 +1,104 @@
+; mc88100 __mpn_sub -- Subtract two limb vectors of the same length > 0 and
+; store difference in a third limb vector.
+
+; Copyright (C) 1992, 1994 Free Software Foundation, Inc.
+
+; This file is part of the GNU MP Library.
+
+; The GNU MP Library is free software; you can redistribute it and/or modify
+; it under the terms of the GNU Library General Public License as published by
+; the Free Software Foundation; either version 2 of the License, or (at your
+; option) any later version.
+
+; The GNU MP Library is distributed in the hope that it will be useful, but
+; WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+; or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
+; License for more details.
+
+; You should have received a copy of the GNU Library General Public License
+; along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
+; the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+
+
+; INPUT PARAMETERS
+; res_ptr      r2
+; s1_ptr       r3
+; s2_ptr       r4
+; size         r5
+
+; This code has been optimized to run one instruction per clock, avoiding
+; load stalls and writeback contention.  As a result, the instruction
+; order is not always natural.
+
+; The speed is about 4.6 clocks/limb + 18 clocks/limb-vector on an 88100,
+; but on the 88110, it seems to run much slower, 6.6 clocks/limb.
+
+       text
+       align    16
+       global   ___mpn_sub_n
+___mpn_sub_n:
+       ld      r6,r3,0                 ; read first limb from s1_ptr
+       extu    r10,r5,3
+       ld      r7,r4,0                 ; read first limb from s2_ptr
+
+       subu.co r5,r0,r5                ; (clear carry as side effect)
+       mak     r5,r5,3<4>
+       bcnd    eq0,r5,Lzero
+
+       or      r12,r0,lo16(Lbase)
+       or.u    r12,r12,hi16(Lbase)
+       addu    r12,r12,r5              ; r12 is address for entering in loop
+
+       extu    r5,r5,2                 ; divide by 4
+       subu    r2,r2,r5                ; adjust res_ptr
+       subu    r3,r3,r5                ; adjust s1_ptr
+       subu    r4,r4,r5                ; adjust s2_ptr
+
+       or      r8,r6,r0
+
+       jmp.n   r12
+        or     r9,r7,r0
+
+Loop:  addu    r3,r3,32
+       st      r8,r2,28
+       addu    r4,r4,32
+       ld      r6,r3,0
+       addu    r2,r2,32
+       ld      r7,r4,0
+Lzero: subu    r10,r10,1               ; subtract 0 + 8r limbs (adj loop cnt)
+Lbase: ld      r8,r3,4
+       subu.cio r6,r6,r7
+       ld      r9,r4,4
+       st      r6,r2,0
+       ld      r6,r3,8                 ; subtract 7 + 8r limbs
+       subu.cio r8,r8,r9
+       ld      r7,r4,8
+       st      r8,r2,4
+       ld      r8,r3,12                ; subtract 6 + 8r limbs
+       subu.cio r6,r6,r7
+       ld      r9,r4,12
+       st      r6,r2,8
+       ld      r6,r3,16                ; subtract 5 + 8r limbs
+       subu.cio r8,r8,r9
+       ld      r7,r4,16
+       st      r8,r2,12
+       ld      r8,r3,20                ; subtract 4 + 8r limbs
+       subu.cio r6,r6,r7
+       ld      r9,r4,20
+       st      r6,r2,16
+       ld      r6,r3,24                ; subtract 3 + 8r limbs
+       subu.cio r8,r8,r9
+       ld      r7,r4,24
+       st      r8,r2,20
+       ld      r8,r3,28                ; subtract 2 + 8r limbs
+       subu.cio r6,r6,r7
+       ld      r9,r4,28
+       st      r6,r2,24
+       bcnd.n  ne0,r10,Loop            ; subtract 1 + 8r limbs
+        subu.cio r8,r8,r9
+
+       st      r8,r2,28                ; store most significant limb
+
+       addu.ci r2,r0,r0                ; return carry-out from most sign. limb
+       jmp.n    r1
+        xor    r2,r2,1