xref: /dflybsd-src/contrib/gmp/mpn/generic/redc_2.c (revision 86d7f5d305c6adaa56ff4582ece9859d73106103)
186d7f5d3SJohn Marino /* mpn_redc_2.  Set cp[] <- up[]/R^n mod mp[].  Clobber up[].
286d7f5d3SJohn Marino    mp[] is n limbs; up[] is 2n limbs.
386d7f5d3SJohn Marino 
486d7f5d3SJohn Marino    THIS IS AN INTERNAL FUNCTION WITH A MUTABLE INTERFACE.  IT IS ONLY
586d7f5d3SJohn Marino    SAFE TO REACH THIS FUNCTION THROUGH DOCUMENTED INTERFACES.
686d7f5d3SJohn Marino 
786d7f5d3SJohn Marino Copyright (C) 2000, 2001, 2002, 2004, 2008 Free Software Foundation, Inc.
886d7f5d3SJohn Marino 
986d7f5d3SJohn Marino This file is part of the GNU MP Library.
1086d7f5d3SJohn Marino 
1186d7f5d3SJohn Marino The GNU MP Library is free software; you can redistribute it and/or modify
1286d7f5d3SJohn Marino it under the terms of the GNU Lesser General Public License as published by
1386d7f5d3SJohn Marino the Free Software Foundation; either version 3 of the License, or (at your
1486d7f5d3SJohn Marino option) any later version.
1586d7f5d3SJohn Marino 
1686d7f5d3SJohn Marino The GNU MP Library is distributed in the hope that it will be useful, but
1786d7f5d3SJohn Marino WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1886d7f5d3SJohn Marino or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
1986d7f5d3SJohn Marino License for more details.
2086d7f5d3SJohn Marino 
2186d7f5d3SJohn Marino You should have received a copy of the GNU Lesser General Public License
2286d7f5d3SJohn Marino along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
2386d7f5d3SJohn Marino 
2486d7f5d3SJohn Marino #include "gmp.h"
2586d7f5d3SJohn Marino #include "gmp-impl.h"
2686d7f5d3SJohn Marino #include "longlong.h"
2786d7f5d3SJohn Marino 
2886d7f5d3SJohn Marino 
2986d7f5d3SJohn Marino #if GMP_NAIL_BITS != 0
3086d7f5d3SJohn Marino you lose
3186d7f5d3SJohn Marino #endif
3286d7f5d3SJohn Marino 
3386d7f5d3SJohn Marino /* For testing purposes, define our own mpn_addmul_2 if there is none already
3486d7f5d3SJohn Marino    available.  */
3586d7f5d3SJohn Marino #ifndef HAVE_NATIVE_mpn_addmul_2
3686d7f5d3SJohn Marino mp_limb_t
mpn_addmul_2(mp_ptr rp,mp_srcptr up,mp_size_t n,mp_srcptr vp)3786d7f5d3SJohn Marino mpn_addmul_2 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_srcptr vp)
3886d7f5d3SJohn Marino {
3986d7f5d3SJohn Marino   rp[n] = mpn_addmul_1 (rp, up, n, vp[0]);
4086d7f5d3SJohn Marino   return mpn_addmul_1 (rp + 1, up, n, vp[1]);
4186d7f5d3SJohn Marino }
4286d7f5d3SJohn Marino #endif
4386d7f5d3SJohn Marino 
4486d7f5d3SJohn Marino #if defined (__GNUC__) && defined (__ia64) && W_TYPE_SIZE == 64
4586d7f5d3SJohn Marino #define umul2low(ph, pl, uh, ul, vh, vl) \
4686d7f5d3SJohn Marino   do {									\
4786d7f5d3SJohn Marino     mp_limb_t _ph, _pl;							\
4886d7f5d3SJohn Marino     __asm__ ("xma.hu %0 = %3, %5, f0\n\t"				\
4986d7f5d3SJohn Marino 	     "xma.l %1 = %3, %5, f0\n\t"				\
5086d7f5d3SJohn Marino 	     ";;\n\t"							\
5186d7f5d3SJohn Marino 	     "xma.l %0 = %3, %4, %0\n\t"				\
5286d7f5d3SJohn Marino 	     ";;\n\t"							\
5386d7f5d3SJohn Marino 	     "xma.l %0 = %2, %5, %0"					\
5486d7f5d3SJohn Marino 	     : "=&f" (ph), "=&f" (pl)					\
5586d7f5d3SJohn Marino 	     : "f" (uh), "f" (ul), "f" (vh), "f" (vl));			\
5686d7f5d3SJohn Marino   } while (0)
5786d7f5d3SJohn Marino #endif
5886d7f5d3SJohn Marino 
5986d7f5d3SJohn Marino #ifndef umul2low
6086d7f5d3SJohn Marino #define umul2low(ph, pl, uh, ul, vh, vl) \
6186d7f5d3SJohn Marino   do {									\
6286d7f5d3SJohn Marino     mp_limb_t _ph, _pl;							\
6386d7f5d3SJohn Marino     umul_ppmm (_ph, _pl, ul, vl);					\
6486d7f5d3SJohn Marino     (ph) = _ph + (ul) * (vh) + (uh) * (vl);				\
6586d7f5d3SJohn Marino     (pl) = _pl;								\
6686d7f5d3SJohn Marino   } while (0)
6786d7f5d3SJohn Marino #endif
6886d7f5d3SJohn Marino 
6986d7f5d3SJohn Marino void
mpn_redc_2(mp_ptr rp,mp_ptr up,mp_srcptr mp,mp_size_t n,mp_srcptr mip)7086d7f5d3SJohn Marino mpn_redc_2 (mp_ptr rp, mp_ptr up, mp_srcptr mp, mp_size_t n, mp_srcptr mip)
7186d7f5d3SJohn Marino {
7286d7f5d3SJohn Marino   mp_limb_t q[2];
7386d7f5d3SJohn Marino   mp_size_t j;
7486d7f5d3SJohn Marino   mp_limb_t upn;
7586d7f5d3SJohn Marino   mp_limb_t cy;
7686d7f5d3SJohn Marino 
7786d7f5d3SJohn Marino   ASSERT (n > 0);
7886d7f5d3SJohn Marino   ASSERT_MPN (up, 2*n);
7986d7f5d3SJohn Marino 
8086d7f5d3SJohn Marino   if ((n & 1) != 0)
8186d7f5d3SJohn Marino     {
8286d7f5d3SJohn Marino       up[0] = mpn_addmul_1 (up, mp, n, (up[0] * mip[0]) & GMP_NUMB_MASK);
8386d7f5d3SJohn Marino       up++;
8486d7f5d3SJohn Marino     }
8586d7f5d3SJohn Marino 
8686d7f5d3SJohn Marino   for (j = n - 2; j >= 0; j -= 2)
8786d7f5d3SJohn Marino     {
8886d7f5d3SJohn Marino       umul2low (q[1], q[0], mip[1], mip[0], up[1], up[0]);
8986d7f5d3SJohn Marino       upn = up[n];		/* mpn_addmul_2 overwrites this */
9086d7f5d3SJohn Marino       up[1] = mpn_addmul_2 (up, mp, n, q);
9186d7f5d3SJohn Marino       up[0] = up[n];
9286d7f5d3SJohn Marino       up[n] = upn;
9386d7f5d3SJohn Marino       up += 2;
9486d7f5d3SJohn Marino     }
9586d7f5d3SJohn Marino   cy = mpn_add_n (rp, up, up - n, n);
9686d7f5d3SJohn Marino   if (cy != 0)
9786d7f5d3SJohn Marino     mpn_sub_n (rp, rp, mp, n);
9886d7f5d3SJohn Marino }
99