1 /* mpn_mod_1s_3p (ap, n, b, cps)
2 Divide (ap,,n) by b. Return the single-limb remainder.
3 Requires that d < B / 3.
5 Contributed to the GNU project by Torbjorn Granlund.
7 THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY
8 SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
9 GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
11 Copyright 2008, 2009 Free Software Foundation, Inc.
13 This file is part of the GNU MP Library.
15 The GNU MP Library is free software; you can redistribute it and/or modify
16 it under the terms of the GNU Lesser General Public License as published by
17 the Free Software Foundation; either version 3 of the License, or (at your
18 option) any later version.
20 The GNU MP Library is distributed in the hope that it will be useful, but
21 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
22 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
23 License for more details.
25 You should have received a copy of the GNU Lesser General Public License
26 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
33 mpn_mod_1s_3p_cps (mp_limb_t cps
[6], mp_limb_t b
)
36 mp_limb_t B1modb
, B2modb
, B3modb
, B4modb
;
39 ASSERT (b
<= GMP_NUMB_MAX
/ 3);
41 count_leading_zeros (cnt
, b
);
46 B1modb
= -b
* ((bi
>> (GMP_LIMB_BITS
-cnt
)) | (CNST_LIMB(1) << cnt
));
47 ASSERT (B1modb
<= b
); /* NB: not fully reduced mod b */
48 udiv_rnd_preinv (B2modb
, B1modb
, b
, bi
);
49 udiv_rnd_preinv (B3modb
, B2modb
, b
, bi
);
50 udiv_rnd_preinv (B4modb
, B3modb
, b
, bi
);
54 cps
[2] = B1modb
>> cnt
;
55 cps
[3] = B2modb
>> cnt
;
56 cps
[4] = B3modb
>> cnt
;
57 cps
[5] = B4modb
>> cnt
;
61 mpn_mod_1s_3p (mp_srcptr ap
, mp_size_t n
, mp_limb_t b
, mp_limb_t cps
[6])
63 mp_limb_t rh
, rl
, bi
, q
, ph
, pl
, ch
, cl
, r
;
64 mp_limb_t B1modb
, B2modb
, B3modb
, B4modb
;
73 umul_ppmm (ph
, pl
, ap
[n
- 2], B1modb
);
74 add_ssaaaa (ph
, pl
, ph
, pl
, 0, ap
[n
- 3]);
75 umul_ppmm (ch
, cl
, ap
[n
- 1], B2modb
);
76 add_ssaaaa (rh
, rl
, ph
, pl
, ch
, cl
);
78 for (i
= n
- 6; i
>= 0; i
-= 3)
81 + ap[i+1] * (B mod b) <= (B-1)(b-1)
82 + ap[i+2] * (B^2 mod b) <= (B-1)(b-1)
83 + LO(rr) * (B^3 mod b) <= (B-1)(b-1)
84 + HI(rr) * (B^4 mod b) <= (B-1)(b-1)
86 umul_ppmm (ph
, pl
, ap
[i
+ 1], B1modb
);
87 add_ssaaaa (ph
, pl
, ph
, pl
, 0, ap
[i
+ 0]);
89 umul_ppmm (ch
, cl
, ap
[i
+ 2], B2modb
);
90 add_ssaaaa (ph
, pl
, ph
, pl
, ch
, cl
);
92 umul_ppmm (ch
, cl
, rl
, B3modb
);
93 add_ssaaaa (ph
, pl
, ph
, pl
, ch
, cl
);
95 umul_ppmm (rh
, rl
, rh
, B4modb
);
96 add_ssaaaa (rh
, rl
, rh
, rl
, ph
, pl
);
101 umul_ppmm (ph
, pl
, rl
, B1modb
);
102 add_ssaaaa (ph
, pl
, ph
, pl
, 0, ap
[i
+ 2]);
103 umul_ppmm (rh
, rl
, rh
, B2modb
);
104 add_ssaaaa (rh
, rl
, rh
, rl
, ph
, pl
);
107 umul_ppmm (ph
, pl
, rl
, B1modb
);
108 add_ssaaaa (ph
, pl
, ph
, pl
, 0, ap
[0]);
109 umul_ppmm (rh
, rl
, rh
, B2modb
);
110 add_ssaaaa (rh
, rl
, rh
, rl
, ph
, pl
);
118 umul_ppmm (rh
, cl
, rh
, B1modb
);
119 add_ssaaaa (rh
, rl
, rh
, rl
, 0, cl
);
120 r
= (rh
<< cnt
) | (rl
>> (GMP_LIMB_BITS
- cnt
));
122 udiv_qrnnd_preinv (q
, r
, rh
>> (GMP_LIMB_BITS
- cnt
),
123 (rh
<< cnt
) | (rl
>> (GMP_LIMB_BITS
- cnt
)), b
, bi
);
124 ASSERT (q
<= 3); /* optimize for small quotient? */
127 udiv_qrnnd_preinv (q
, r
, r
, rl
<< cnt
, b
, bi
);