1 /* ---------------------------------------------------------------------
2 *
3 * -- Integer Matrix Library (IML)
4 * (C) Copyright 2004, 2006 All Rights Reserved
5 *
6 * -- IML routines -- Version 1.0.1 -- November, 2006
7 *
8 * Author : Zhuliang Chen
9 * Contributor(s) : Arne Storjohann
10 * University of Waterloo -- School of Computer Science
11 * Waterloo, Ontario, N2L3G1 Canada
12 *
13 * ---------------------------------------------------------------------
14 *
15 * -- Copyright notice and Licensing terms:
16 *
17 * Redistribution and use in source and binary forms, with or without
18 * modification, are permitted provided that the following conditions
19 * are met:
20 *
21 * 1. Redistributions of source code must retain the above copyright
22 * notice, this list of conditions and the following disclaimer.
23 * 2. Redistributions in binary form must reproduce the above copyright
24 * notice, this list of conditions, and the following disclaimer in
25 * the documentation and/or other materials provided with the distri-
26 * bution.
27 * 3. The name of the University, the IML group, or the names of its
28 * contributors may not be used to endorse or promote products deri-
29 * ved from this software without specific written permission.
30 *
31 * -- Disclaimer:
32 *
33 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
34 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
35 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
36 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
37 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPE-
38 * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
39 * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
40 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
41 * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (IN-
42 * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
43 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
44 *
45 */
46
47
48
49 #include "basisop.h"
50
51 /*
52 * Calling Sequence:
53 * Dmod(p, A, n, m, lda)
54 *
55 * Summary:
56 * Perform mod p operation inplace over a double matrix/submatrix
57 *
58 * Description:
59 * Given a n x m matrix or submatrix A and an integer p, this function uses
60 * fmod function to do operation A mod p inplace. Although the type of input
61 * modulus p in this function is Double, it is a cast of an integer.
62 *
63 * Input:
64 * p: Double, module
65 * A: 1-dim Double array length n*m, representation of a n x m double
66 * matrix/submatrix
67 * n: long, row dimension of A
68 * m: long, column dimension of A
69 * lda: long, number of entries between two continuous rows of A, help track
70 * the pointer if A is a submatrix
71 *
72 */
73
74 void
Dmod(const Double p,Double * A,const long n,const long m,const long lda)75 Dmod (const Double p, Double *A, const long n, const long m, const long lda)
76 {
77 long i, j;
78 Double temp;
79
80 for (i = 0; i < n; i++)
81 for (j = 0; j < m; j++)
82 *(A+i*lda+j) = (temp = fmod(*(A+i*lda+j), p)) >= 0 ? temp : p+temp;
83 }
84
85
86 /*
87 * Calling Sequence:
88 * DCopy(n, m, A, lda, B, ldb)
89 *
90 * Summary:
91 * Copy Double matrix/submatrix A to another matrix/submatrix B
92 *
93 * Input:
94 * n: long, row dimension of matrix/submatrix A and B
95 * m: long, column dimension of matrix/submatrix A and B
96 * A: 1-dim Double array length n*m, representation of a n x m Double
97 * matrix/submatrix
98 * lda: long, stride of two continuous rows of A
99 * B: 1-dim Double array length n*m, representation of a n x m Double
100 * matrix/submatrix
101 * ldb: long, stride of two continuous rows of B
102 *
103 */
104
105 void
DCopy(const long n,const long m,const Double * A,const long lda,Double * B,const long ldb)106 DCopy (const long n, const long m, const Double* A, const long lda, \
107 Double* B, const long ldb)
108 {
109 long i, j;
110
111 for (i = 0; i < n; i++)
112 for (j = 0; j < m; j++)
113 *(B+i*ldb+j) = *(A+i*lda+j);
114 }
115
116
117
118 /*
119 * Calling Sequence:
120 * randomDb(n, m, bd)
121 *
122 * Summary:
123 * Generate a random Double dense matrix
124 *
125 * Description:
126 * Given the bound bd, this function generates a random Double dense matrix
127 * such that any entry e satisfies 0 <= e <= 2^bd-1
128 *
129 * Input:
130 * n: long, row dimension of the generated matrix
131 * m: long, column dimension of the generated matrix
132 * bd: long, magnitude bound of entries
133 *
134 * Return:
135 * M: 1-dim Double array length n*m, representation of a n x m random matrix
136 *
137 */
138
139 Double *
randomDb(const long n,const long m,const long bd)140 randomDb (const long n, const long m, const long bd)
141 {
142 long i, j;
143 mpz_t mp_rand;
144 gmp_randstate_t state;
145 unsigned long seed;
146 FILE* devrandom;
147 Double* M;
148 static unsigned long inc = 0;
149
150 #if HAVE_TIME_H
151 time_t t1;
152 #endif
153
154 M = XMALLOC(Double, n*m);
155 mpz_init(mp_rand);
156 gmp_randinit_default(state);
157 mpz_set_ui(mp_rand, 5);
158
159 seed = 3828173;
160
161 /* generate random seed using /dev/random */
162 if ((devrandom = fopen("/dev/urandom", "r")) != NULL)
163 {
164 fread(&seed, sizeof(seed), 1, devrandom);
165 fclose(devrandom);
166 }
167 else
168 {
169 #if HAVE_TIME_H
170 time(&t1);
171 seed = (unsigned long)t1;
172 #endif
173 }
174
175 seed += inc;
176 inc += 1;
177 gmp_randseed_ui(state, seed);
178 for (i = 0; i < n*m; i++)
179 {
180 mpz_urandomb(mp_rand, state, bd);
181 M[i] = mpz_get_d(mp_rand);
182 }
183 mpz_clear(mp_rand);
184 gmp_randclear(state);
185 return M;
186 }
187
188
189
190 /*
191 * Calling Sequence:
192 * RandomPrime(lb, hb)
193 *
194 * Summary:
195 * Generate a random prime p satisfying 2^lb <= p <= 2^hb-1
196 *
197 * Input:
198 * lb: FiniteField, lower bound of the prime
199 * hb: FiniteField, upper bound of the prime
200 *
201 * Return:
202 * p: FiniteField, randomly generated prime
203 *
204 * Note:
205 * lb and hb can not exceed the bit-length of unsigned long.
206 *
207 */
208
209 FiniteField
RandPrime(const FiniteField lb,const FiniteField hb)210 RandPrime (const FiniteField lb, const FiniteField hb)
211 {
212 mpz_t mp_rand, mp_temp, mp_lb, mp_hb;
213 gmp_randstate_t state;
214 FiniteField p;
215 unsigned long seed;
216 FILE* devrandom;
217 static unsigned long inc = 0;
218
219 #if HAVE_TIME_H
220 time_t t1;
221 #endif
222
223 { mpz_init(mp_rand); mpz_init(mp_temp); }
224 { mpz_init(mp_lb); mpz_init(mp_hb); }
225 mpz_ui_pow_ui(mp_lb, 2, lb);
226 mpz_ui_pow_ui(mp_hb, 2, hb);
227 mpz_sub(mp_temp, mp_hb, mp_lb);
228 gmp_randinit_default(state);
229
230 seed = 3828173;
231
232 /* generate random seed using /dev/urandom */
233 if ((devrandom = fopen("/dev/urandom", "r")) != NULL)
234 {
235 fread(&seed, sizeof(seed), 1, devrandom);
236 fclose(devrandom);
237 }
238 else
239 {
240 #if HAVE_TIME_H
241 time(&t1);
242 seed = (unsigned long)t1;
243 #endif
244 }
245
246 seed += inc;
247 inc += 1;
248 gmp_randseed_ui(state, seed);
249 mpz_urandomm(mp_rand, state, mp_temp); /* 0 <= mp_rand <= 2^hb-2^lb-1 */
250 mpz_add(mp_rand, mp_rand, mp_lb); /* 2^lb <= mp_rand <= 2^hb-1 */
251 while (mpz_probab_prime_p(mp_rand, 10) == 0)
252 mpz_sub_ui(mp_rand, mp_rand, 1);
253 p = mpz_get_ui(mp_rand);
254
255 { mpz_clear(mp_rand); mpz_clear(mp_temp); }
256 { mpz_clear(mp_lb); mpz_clear(mp_hb); }
257 gmp_randclear(state);
258 return p;
259 }
260
261
262
263 /*
264 * Calling Sequence:
265 * max <-- maxMagnLong(A, n, m, lda)
266 *
267 * Summary:
268 * Compute maximum magnitude of a n x m signed long matrix/submatrix A
269 *
270 * Input:
271 * A: 1-dim long array length n*m, representation of a n x m
272 * matrix/submatrix
273 * n: long, row dimension of A
274 * m: long, column dimension of A
275 * lda: long, number of entries of two continuous rows of A (normally m)
276 *
277 * Return:
278 * max: long, maximum magnitude of A
279 *
280 */
281
282 long
maxMagnLong(const long * A,const long n,const long m,const long lda)283 maxMagnLong (const long *A, const long n, const long m, const long lda)
284 {
285 long i, j, temp, max=0;
286
287 for (i = 0; i < n; i++)
288 for (j = 0; j < m; j++)
289 if ((temp = labs(A[i*lda+j])) > max) { max = temp; }
290 return max;
291 }
292
293
294
295 /*
296 *
297 * Calling Sequence:
298 * maxMagnMP(mp_A, n, m, lda, mp_max)
299 *
300 * Summary:
301 * Compute maximum magnitude of a n x m mpz_t matrix/submatrix mp_A
302 *
303 * Input:
304 * mp_A: 1-dim mpz_t array length n*m, representation of a n x m
305 * matrix/submatrix
306 * n: long, row dimension of mp_A
307 * m: long, column dimension of mp_A
308 * lda: long, number of entries of two continuous rows of mp_A (normally m)
309 *
310 * Output:
311 * mp_max: mpz_t, maximum magnitude of mp_A
312 *
313 */
314
315 void
maxMagnMP(mpz_t * mp_A,const long n,const long m,const long lda,mpz_t mp_max)316 maxMagnMP (mpz_t *mp_A, const long n, const long m, const long lda, \
317 mpz_t mp_max)
318 {
319 long i, j;
320
321 mpz_set_ui(mp_max, 0);
322 for (i = 0; i < n; i++)
323 for (j = 0; j < m; j++)
324 if (mpz_cmpabs(mp_A[i*lda+j], mp_max) > 0)
325 mpz_abs(mp_max, mp_A[i*lda+j]);
326 return;
327 }
328
329
330 /*
331 * Calling Sequence:
332 * scalCpMP(n, m, lda, ldm, mp_a, mp_A, mp_M)
333 *
334 * Summary:
335 * Copy a submatrix with scale from one matrix to the other one
336 *
337 * Description:
338 * Copy with scalar n x m submatrix from A to M. The submatrix starts
339 * from position pointed by mp_A with stride lda. Then the submatrix is
340 * copied to position pointed by mp_M with stride ldm, meanwhile, the
341 * entries of the submatrix are scaled by factor mp_a.
342 *
343 * Input:
344 * n: long, row dimension of the submatrix
345 * m: long, column dimension of the submatrix
346 * lda: long, stride of two continuous rows of the submatrix in A
347 * ldm: long, stride of two continuous rows of the submatrix in M
348 * mp_a: mpz_t, scalar factor
349 * mp_A: mpz_t pointer, start point of the submatrix in A
350 * mp_M: mpz_t pointer, start point of the submatrix in M
351 *
352 */
353
354 void
scalCpMP(const long n,const long m,const long lda,const long ldm,const mpz_t mp_a,mpz_t * mp_A,mpz_t * mp_M)355 scalCpMP (const long n, const long m, const long lda, const long ldm, \
356 const mpz_t mp_a, mpz_t *mp_A, mpz_t *mp_M)
357 {
358 long i, j;
359
360 if (mpz_cmp_ui(mp_a, 1) != 0)
361 for (i = 0; i < n; i++)
362 for (j = 0; j < m; j++)
363 {
364 mpz_set(mp_M[i*ldm+j], mp_A[i*lda+j]);
365 mpz_mul(mp_M[i*ldm+j], mp_M[i*ldm+j], mp_a);
366 }
367 else
368 for (i = 0; i < n; i++)
369 for (j = 0; j < m; j++)
370 mpz_set(mp_M[i*ldm+j], mp_A[i*lda+j]);
371 }
372
373