source: git/ntl/src/vec_lzz_p.c @ 287cc8

spielwiese
Last change on this file since 287cc8 was 287cc8, checked in by Hans Schönemann <hannes@…>, 14 years ago
ntl 5.5.2 git-svn-id: file:///usr/local/Singular/svn/trunk@12402 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 4.4 KB
Line 
1
2#include <NTL/vec_lzz_p.h>
3
4#include <NTL/new.h>
5
6NTL_START_IMPL
7
8NTL_vector_impl(zz_p,vec_zz_p)
9
10NTL_eq_vector_impl(zz_p,vec_zz_p)
11
12void conv(vec_zz_p& x, const vec_ZZ& a)
13{
14   long i, n;
15
16   n = a.length();
17   x.SetLength(n);
18
19   zz_p* xp = x.elts();
20   const ZZ* ap = a.elts();
21
22   for (i = 0; i < n; i++)
23      conv(xp[i], ap[i]);
24}
25
26void conv(vec_ZZ& x, const vec_zz_p& a)
27{
28   long n = a.length();
29   x.SetLength(n);
30   long i;
31   for (i = 0; i < n; i++)
32      x[i] = rep(a[i]);
33}
34
35
36
37
38void InnerProduct(zz_p& x, const vec_zz_p& a, const vec_zz_p& b)
39{
40   long n = min(a.length(), b.length());
41   long i;
42   zz_p accum, t;
43
44   clear(accum);
45   for (i = 0; i < n; i++) {
46      mul(t, a[i], b[i]);
47      add(accum, accum, t);
48   }
49
50   x = accum;
51}
52
53void InnerProduct(zz_p& x, const vec_zz_p& a, const vec_zz_p& b,
54                  long offset)
55{
56   if (offset < 0) Error("InnerProduct: negative offset");
57   if (NTL_OVERFLOW(offset, 1, 0)) Error("InnerProduct: offset too big");
58
59   long n = min(a.length(), b.length()+offset);
60   long i;
61   zz_p accum, t;
62
63   clear(accum);
64   for (i = offset; i < n; i++) {
65      mul(t, a[i], b[i-offset]);
66      add(accum, accum, t);
67   }
68
69   x = accum;
70}
71
72long CRT(vec_ZZ& gg, ZZ& a, const vec_zz_p& G)
73{
74   long n = gg.length();
75   if (G.length() != n) Error("CRT: vector length mismatch");
76
77   long p = zz_p::modulus();
78
79   ZZ new_a;
80   mul(new_a, a, p);
81
82   long a_inv;
83   a_inv = rem(a, p);
84   a_inv = InvMod(a_inv, p);
85
86   long p1;
87   p1 = p >> 1;
88
89   ZZ a1;
90   RightShift(a1, a, 1);
91
92   long p_odd = (p & 1);
93
94   long modified = 0;
95
96   long h;
97
98   ZZ g;
99   long i;
100   for (i = 0; i < n; i++) {
101      if (!CRTInRange(gg[i], a)) {
102         modified = 1;
103         rem(g, gg[i], a);
104         if (g > a1) sub(g, g, a);
105      }
106      else
107         g = gg[i];
108
109      h = rem(g, p);
110      h = SubMod(rep(G[i]), h, p);
111      h = MulMod(h, a_inv, p);
112      if (h > p1)
113         h = h - p;
114
115      if (h != 0) {
116         modified = 1;
117
118         if (!p_odd && g > 0 && (h == p1))
119            MulSubFrom(g, a, h);
120         else
121            MulAddTo(g, a, h);
122      }
123
124      gg[i] = g;
125   }
126
127   a = new_a;
128
129   return modified;
130}
131
132
133
134void mul(vec_zz_p& x, const vec_zz_p& a, zz_p b)
135{
136   long n = a.length();
137   x.SetLength(n);
138
139   long i;
140
141   if (n <= 1) {
142
143      for (i = 0; i < n; i++)
144         mul(x[i], a[i], b);
145
146   }
147   else {
148
149      long p = zz_p::modulus();
150      double pinv = zz_p::ModulusInverse();
151      long bb = rep(b);
152      mulmod_precon_t bpinv = PrepMulModPrecon(bb, p, pinv);
153
154
155      const zz_p *ap = a.elts();
156      zz_p *xp = x.elts();
157
158      for (i = 0; i < n; i++)
159         xp[i].LoopHole() = MulModPrecon(rep(ap[i]), bb, p, bpinv);
160
161   }
162}
163
164void mul(vec_zz_p& x, const vec_zz_p& a, long b_in)
165{
166   zz_p b;
167   b = b_in;
168   mul(x, a, b);
169}
170
171
172
173void add(vec_zz_p& x, const vec_zz_p& a, const vec_zz_p& b)
174{
175   long n = a.length();
176   if (b.length() != n) Error("vector add: dimension mismatch");
177
178   x.SetLength(n);
179   long i;
180   for (i = 0; i < n; i++)
181      add(x[i], a[i], b[i]);
182}
183
184void sub(vec_zz_p& x, const vec_zz_p& a, const vec_zz_p& b)
185{
186   long n = a.length();
187   if (b.length() != n) Error("vector sub: dimension mismatch");
188   x.SetLength(n);
189   long i;
190   for (i = 0; i < n; i++)
191      sub(x[i], a[i], b[i]);
192}
193
194void clear(vec_zz_p& x)
195{
196   long n = x.length();
197   long i;
198   for (i = 0; i < n; i++)
199      clear(x[i]);
200}
201
202void negate(vec_zz_p& x, const vec_zz_p& a)
203{
204   long n = a.length();
205   x.SetLength(n);
206   long i;
207   for (i = 0; i < n; i++)
208      negate(x[i], a[i]);
209}
210
211
212long IsZero(const vec_zz_p& a)
213{
214   long n = a.length();
215   long i;
216
217   for (i = 0; i < n; i++)
218      if (!IsZero(a[i]))
219         return 0;
220
221   return 1;
222}
223
224vec_zz_p operator+(const vec_zz_p& a, const vec_zz_p& b)
225{
226   vec_zz_p res;
227   add(res, a, b);
228   NTL_OPT_RETURN(vec_zz_p, res);
229}
230
231vec_zz_p operator-(const vec_zz_p& a, const vec_zz_p& b)
232{
233   vec_zz_p res;
234   sub(res, a, b);
235   NTL_OPT_RETURN(vec_zz_p, res);
236}
237
238
239vec_zz_p operator-(const vec_zz_p& a)
240{
241   vec_zz_p res;
242   negate(res, a);
243   NTL_OPT_RETURN(vec_zz_p, res);
244}
245
246
247zz_p operator*(const vec_zz_p& a, const vec_zz_p& b)
248{
249   zz_p res;
250   InnerProduct(res, a, b);
251   return res;
252}
253
254
255void VectorCopy(vec_zz_p& x, const vec_zz_p& a, long n)
256{
257   if (n < 0) Error("VectorCopy: negative length");
258   if (NTL_OVERFLOW(n, 1, 0)) Error("overflow in VectorCopy");
259
260   long m = min(n, a.length());
261
262   x.SetLength(n);
263
264   long i;
265
266   for (i = 0; i < m; i++)
267      x[i] = a[i];
268
269   for (i = m; i < n; i++)
270      clear(x[i]);
271}
272
273NTL_END_IMPL
Note: See TracBrowser for help on using the repository browser.