Annotation of OpenXM_contrib/gmp/mpn/generic/mode1o.c, Revision 1.1.1.1
1.1 ohara 1: /* mpn_modexact_1c_odd -- mpn by limb exact division style remainder.
2:
3: THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY. THEY'RE ALMOST
4: CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
5: FUTURE GNU MP RELEASES.
6:
7: Copyright 2000, 2001, 2002 Free Software Foundation, Inc.
8:
9: This file is part of the GNU MP Library.
10:
11: The GNU MP Library is free software; you can redistribute it and/or modify
12: it under the terms of the GNU Lesser General Public License as published by
13: the Free Software Foundation; either version 2.1 of the License, or (at your
14: option) any later version.
15:
16: The GNU MP Library is distributed in the hope that it will be useful, but
17: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19: License for more details.
20:
21: You should have received a copy of the GNU Lesser General Public License
22: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
23: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
24: MA 02111-1307, USA.
25: */
26:
27: #include "gmp.h"
28: #include "gmp-impl.h"
29: #include "longlong.h"
30:
31:
32: /* Calculate an r satisfying
33:
34: r*b^k + a - c == q*d
35:
36: where b=2^BITS_PER_MP_LIMB, a is {src,size}, k is either size or size-1
37: (the caller won't know which), and q is the quotient (discarded). d must
38: be odd, c can be any limb value.
39:
40: If c<d then r will be in the range 0<=r<d, or if c>=d then 0<=r<=d.
41:
42: This slightly strange function suits the initial Nx1 reduction for GCDs
43: or Jacobi symbols since the factors of 2 in b^k can be ignored, leaving
44: -r == a mod d (by passing c=0). For a GCD the factor of -1 on r can be
45: ignored, or for the Jacobi symbol it can be accounted for. The function
46: also suits divisibility and congruence testing since if r=0 (or r=d) is
47: obtained then a==c mod d.
48:
49:
50: r is a bit like the remainder returned by mpn_divexact_by3c, and is the
51: sort of remainder a hypothetical mpn_divexact_1 might return. Like
52: mpn_divexact_by3c, r represents a borrow, since effectively quotient
53: limbs are chosen so that subtracting that multiple of d from src at each
54: step will produce a zero limb.
55:
56: A long calculation can be done piece by piece from low to high by passing
57: the return value from one part as the carry parameter to the next part.
58: The effective final k becomes anything between size and size-n, if n
59: pieces are used.
60:
61:
62: A similar sort of routine could be constructed based on adding multiples
63: of d at each limb, much like redc in mpz_powm does. Subtracting however
64: has a small advantage that when subtracting to cancel out l there's never
65: a borrow into h, whereas using an addition would put a carry into h
66: depending whether l==0 or l!=0.
67:
68:
69: In terms of efficiency, this function is similar to a mul-by-inverse
70: mpn_mod_1. Both are essentially two multiplies and are best suited to
71: CPUs with low latency multipliers (in comparison to a divide instruction
72: at least.) But modexact has a few less supplementary operations, only
73: needs low part and high part multiplies, and has fewer working quantities
74: (helping CPUs with few registers).
75:
76:
77: In the main loop it will be noted that the new carry (call it r) is the
78: sum of the high product h and any borrow from l=s-c. If c<d then we will
79: have r<d too, for the following reasons. Let q=l*inverse be the quotient
80: limb, so that q*d = b*h + l. Now if h=d-1 then
81:
82: l = q*d - b*(d-1) <= (b-1)*d - b*(d-1) = b-d
83:
84: But if l=s-c produces a borrow when c<d, then l>=b-d+1 and hence will
85: never have h=d-1 and so r=h+borrow <= d-1.
86:
87: When c>=d, on the other hand, h=d-1 can certainly occur together with a
88: borrow, thereby giving only r<=d, as per the function definition above.
89:
90: As a design decision it's left to the caller to check for r=d if it might
91: be passing c>=d. Several applications have c<d initially so the extra
92: test is often unnecessary, for example the GCDs or a plain divisibility
93: d|a test will pass c=0.
94:
95:
96: The special case for size==1 is so that it can be assumed c<=d in the
97: high<=divisor test at the end. c<=d is only guaranteed after at least
98: one iteration of the main loop. There's also a decent chance one % is
99: faster than a modlimb_invert, though that will depend on the processor.
100:
101: A CPU specific implementation might want to omit the size==1 code or the
102: high<divisor test. mpn/x86/k6/mode1o.asm for instance finds neither
103: useful. */
104:
105:
106: mp_limb_t
107: mpn_modexact_1c_odd (mp_srcptr src, mp_size_t size, mp_limb_t d, mp_limb_t c)
108: {
109: mp_limb_t s, h, l, inverse, dummy, dmul;
110: mp_size_t i;
111:
112: ASSERT (size >= 1);
113: ASSERT (d & 1);
114: ASSERT_MPN (src, size);
115: ASSERT_LIMB (d);
116: ASSERT_LIMB (c);
117:
118: if (size == 1)
119: {
120: s = src[0];
121: if (s > c)
122: {
123: l = s-c;
124: h = l % d;
125: if (h != 0)
126: h = d - h;
127: }
128: else
129: {
130: l = c-s;
131: h = l % d;
132: }
133: return h;
134: }
135:
136:
137: modlimb_invert (inverse, d);
138: dmul = d << GMP_NAIL_BITS;
139:
140: i = 0;
141: do
142: {
143: s = src[i];
144: SUBC_LIMB (c, l, s, c);
145: l = (l * inverse) & GMP_NUMB_MASK;
146: umul_ppmm (h, dummy, l, dmul);
147: c += h;
148: }
149: while (++i < size-1);
150:
151:
152: s = src[i];
153: if (s <= d)
154: {
155: /* With high<=d the final step can be a subtract and addback. If c==0
156: then the addback will restore to l>=0. If c==d then will get l==d
157: if s==0, but that's ok per the function definition. */
158:
159: l = c - s;
160: if (l > c)
161: l += d;
162:
163: ASSERT (l <= d);
164: return l;
165: }
166: else
167: {
168: /* Can't skip a divide, just do the loop code once more. */
169:
170: SUBC_LIMB (c, l, s, c);
171: l = (l * inverse) & GMP_NUMB_MASK;
172: umul_ppmm (h, dummy, l, dmul);
173: c += h;
174:
175: ASSERT (c <= d);
176: return c;
177: }
178: }
179:
180:
181:
182: #if 0
183:
184: /* The following is an alternate form that might shave one cycle on a
185: superscalar processor since it takes c+=h off the dependent chain,
186: leaving just a low product, high product, and a subtract.
187:
188: This is for CPU specific implementations to consider. A special case for
189: high<divisor and/or size==1 can be added if desired.
190:
191: Notice that c is only ever 0 or 1, since if s-c produces a borrow then
192: x=0xFF..FF and x-h cannot produce a borrow. The c=(x>s) could become
193: c=(x==0xFF..FF) too, if that helped. */
194:
195: mp_limb_t
196: mpn_modexact_1c_odd (mp_srcptr src, mp_size_t size, mp_limb_t d, mp_limb_t h)
197: {
198: mp_limb_t s, x, y, inverse, dummy, dmul, c1, c2;
199: mp_limb_t c = 0;
200: mp_size_t i;
201:
202: ASSERT (size >= 1);
203: ASSERT (d & 1);
204:
205: modlimb_invert (inverse, d);
206: dmul = d << GMP_NAIL_BITS;
207:
208: for (i = 0; i < size; i++)
209: {
210: ASSERT (c==0 || c==1);
211:
212: s = src[i];
213: SUBC_LIMB (c1, x, s, c);
214:
215: SUBC_LIMB (c2, y, x, h);
216: c = c1 + c2;
217:
218: y = (y * inverse) & GMP_NUMB_MASK;
219: umul_ppmm (h, dummy, y, dmul);
220: }
221:
222: h += c;
223: return h;
224: }
225:
226: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to mode1o.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpn / generic