xref: /netbsd-src/external/mit/isl/dist/isl_seq.c (revision 5971e316fdea024efff6be8f03536623db06833e)
1*5971e316Smrg /*
2*5971e316Smrg  * Copyright 2008-2009 Katholieke Universiteit Leuven
3*5971e316Smrg  * Copyright 2011      INRIA Saclay
4*5971e316Smrg  *
5*5971e316Smrg  * Use of this software is governed by the MIT license
6*5971e316Smrg  *
7*5971e316Smrg  * Written by Sven Verdoolaege, K.U.Leuven, Departement
8*5971e316Smrg  * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
9*5971e316Smrg  */
10*5971e316Smrg 
11*5971e316Smrg #include <isl_ctx_private.h>
12*5971e316Smrg #include <isl_seq.h>
13*5971e316Smrg 
isl_seq_clr(isl_int * p,unsigned len)14*5971e316Smrg void isl_seq_clr(isl_int *p, unsigned len)
15*5971e316Smrg {
16*5971e316Smrg 	int i;
17*5971e316Smrg 	for (i = 0; i < len; ++i)
18*5971e316Smrg 		isl_int_set_si(p[i], 0);
19*5971e316Smrg }
20*5971e316Smrg 
isl_seq_set_si(isl_int * p,int v,unsigned len)21*5971e316Smrg void isl_seq_set_si(isl_int *p, int v, unsigned len)
22*5971e316Smrg {
23*5971e316Smrg 	int i;
24*5971e316Smrg 	for (i = 0; i < len; ++i)
25*5971e316Smrg 		isl_int_set_si(p[i], v);
26*5971e316Smrg }
27*5971e316Smrg 
isl_seq_set(isl_int * p,isl_int v,unsigned len)28*5971e316Smrg void isl_seq_set(isl_int *p, isl_int v, unsigned len)
29*5971e316Smrg {
30*5971e316Smrg 	int i;
31*5971e316Smrg 	for (i = 0; i < len; ++i)
32*5971e316Smrg 		isl_int_set(p[i], v);
33*5971e316Smrg }
34*5971e316Smrg 
isl_seq_neg(isl_int * dst,isl_int * src,unsigned len)35*5971e316Smrg void isl_seq_neg(isl_int *dst, isl_int *src, unsigned len)
36*5971e316Smrg {
37*5971e316Smrg 	int i;
38*5971e316Smrg 	for (i = 0; i < len; ++i)
39*5971e316Smrg 		isl_int_neg(dst[i], src[i]);
40*5971e316Smrg }
41*5971e316Smrg 
isl_seq_cpy(isl_int * dst,isl_int * src,unsigned len)42*5971e316Smrg void isl_seq_cpy(isl_int *dst, isl_int *src, unsigned len)
43*5971e316Smrg {
44*5971e316Smrg 	int i;
45*5971e316Smrg 	for (i = 0; i < len; ++i)
46*5971e316Smrg 		isl_int_set(dst[i], src[i]);
47*5971e316Smrg }
48*5971e316Smrg 
isl_seq_submul(isl_int * dst,isl_int f,isl_int * src,unsigned len)49*5971e316Smrg void isl_seq_submul(isl_int *dst, isl_int f, isl_int *src, unsigned len)
50*5971e316Smrg {
51*5971e316Smrg 	int i;
52*5971e316Smrg 	for (i = 0; i < len; ++i)
53*5971e316Smrg 		isl_int_submul(dst[i], f, src[i]);
54*5971e316Smrg }
55*5971e316Smrg 
isl_seq_addmul(isl_int * dst,isl_int f,isl_int * src,unsigned len)56*5971e316Smrg void isl_seq_addmul(isl_int *dst, isl_int f, isl_int *src, unsigned len)
57*5971e316Smrg {
58*5971e316Smrg 	int i;
59*5971e316Smrg 	for (i = 0; i < len; ++i)
60*5971e316Smrg 		isl_int_addmul(dst[i], f, src[i]);
61*5971e316Smrg }
62*5971e316Smrg 
isl_seq_swp_or_cpy(isl_int * dst,isl_int * src,unsigned len)63*5971e316Smrg void isl_seq_swp_or_cpy(isl_int *dst, isl_int *src, unsigned len)
64*5971e316Smrg {
65*5971e316Smrg 	int i;
66*5971e316Smrg 	for (i = 0; i < len; ++i)
67*5971e316Smrg 		isl_int_swap_or_set(dst[i], src[i]);
68*5971e316Smrg }
69*5971e316Smrg 
isl_seq_scale(isl_int * dst,isl_int * src,isl_int m,unsigned len)70*5971e316Smrg void isl_seq_scale(isl_int *dst, isl_int *src, isl_int m, unsigned len)
71*5971e316Smrg {
72*5971e316Smrg 	int i;
73*5971e316Smrg 	for (i = 0; i < len; ++i)
74*5971e316Smrg 		isl_int_mul(dst[i], src[i], m);
75*5971e316Smrg }
76*5971e316Smrg 
isl_seq_scale_down(isl_int * dst,isl_int * src,isl_int m,unsigned len)77*5971e316Smrg void isl_seq_scale_down(isl_int *dst, isl_int *src, isl_int m, unsigned len)
78*5971e316Smrg {
79*5971e316Smrg 	int i;
80*5971e316Smrg 	for (i = 0; i < len; ++i)
81*5971e316Smrg 		isl_int_divexact(dst[i], src[i], m);
82*5971e316Smrg }
83*5971e316Smrg 
isl_seq_cdiv_q(isl_int * dst,isl_int * src,isl_int m,unsigned len)84*5971e316Smrg void isl_seq_cdiv_q(isl_int *dst, isl_int *src, isl_int m, unsigned len)
85*5971e316Smrg {
86*5971e316Smrg 	int i;
87*5971e316Smrg 	for (i = 0; i < len; ++i)
88*5971e316Smrg 		isl_int_cdiv_q(dst[i], src[i], m);
89*5971e316Smrg }
90*5971e316Smrg 
isl_seq_fdiv_q(isl_int * dst,isl_int * src,isl_int m,unsigned len)91*5971e316Smrg void isl_seq_fdiv_q(isl_int *dst, isl_int *src, isl_int m, unsigned len)
92*5971e316Smrg {
93*5971e316Smrg 	int i;
94*5971e316Smrg 	for (i = 0; i < len; ++i)
95*5971e316Smrg 		isl_int_fdiv_q(dst[i], src[i], m);
96*5971e316Smrg }
97*5971e316Smrg 
isl_seq_fdiv_r(isl_int * dst,isl_int * src,isl_int m,unsigned len)98*5971e316Smrg void isl_seq_fdiv_r(isl_int *dst, isl_int *src, isl_int m, unsigned len)
99*5971e316Smrg {
100*5971e316Smrg 	int i;
101*5971e316Smrg 	for (i = 0; i < len; ++i)
102*5971e316Smrg 		isl_int_fdiv_r(dst[i], src[i], m);
103*5971e316Smrg }
104*5971e316Smrg 
isl_seq_combine(isl_int * dst,isl_int m1,isl_int * src1,isl_int m2,isl_int * src2,unsigned len)105*5971e316Smrg void isl_seq_combine(isl_int *dst, isl_int m1, isl_int *src1,
106*5971e316Smrg 			isl_int m2, isl_int *src2, unsigned len)
107*5971e316Smrg {
108*5971e316Smrg 	int i;
109*5971e316Smrg 	isl_int tmp;
110*5971e316Smrg 
111*5971e316Smrg 	if (dst == src1 && isl_int_is_one(m1)) {
112*5971e316Smrg 		if (isl_int_is_zero(m2))
113*5971e316Smrg 			return;
114*5971e316Smrg 		for (i = 0; i < len; ++i)
115*5971e316Smrg 			isl_int_addmul(src1[i], m2, src2[i]);
116*5971e316Smrg 		return;
117*5971e316Smrg 	}
118*5971e316Smrg 
119*5971e316Smrg 	isl_int_init(tmp);
120*5971e316Smrg 	for (i = 0; i < len; ++i) {
121*5971e316Smrg 		isl_int_mul(tmp, m1, src1[i]);
122*5971e316Smrg 		isl_int_addmul(tmp, m2, src2[i]);
123*5971e316Smrg 		isl_int_set(dst[i], tmp);
124*5971e316Smrg 	}
125*5971e316Smrg 	isl_int_clear(tmp);
126*5971e316Smrg }
127*5971e316Smrg 
128*5971e316Smrg /* Eliminate element "pos" from "dst" using "src".
129*5971e316Smrg  * In particular, let d = dst[pos] and s = src[pos], then
130*5971e316Smrg  * dst is replaced by (|s| dst - sgn(s)d src)/gcd(s,d),
131*5971e316Smrg  * such that dst[pos] is zero after the elimination.
132*5971e316Smrg  * If "m" is not NULL, then *m is multiplied by |s|/gcd(s,d).
133*5971e316Smrg  * That is, it is multiplied by the same factor as "dst".
134*5971e316Smrg  */
isl_seq_elim(isl_int * dst,isl_int * src,unsigned pos,unsigned len,isl_int * m)135*5971e316Smrg void isl_seq_elim(isl_int *dst, isl_int *src, unsigned pos, unsigned len,
136*5971e316Smrg 		  isl_int *m)
137*5971e316Smrg {
138*5971e316Smrg 	isl_int a;
139*5971e316Smrg 	isl_int b;
140*5971e316Smrg 
141*5971e316Smrg 	if (isl_int_is_zero(dst[pos]))
142*5971e316Smrg 		return;
143*5971e316Smrg 
144*5971e316Smrg 	isl_int_init(a);
145*5971e316Smrg 	isl_int_init(b);
146*5971e316Smrg 
147*5971e316Smrg 	isl_int_gcd(a, src[pos], dst[pos]);
148*5971e316Smrg 	isl_int_divexact(b, dst[pos], a);
149*5971e316Smrg 	if (isl_int_is_pos(src[pos]))
150*5971e316Smrg 		isl_int_neg(b, b);
151*5971e316Smrg 	isl_int_divexact(a, src[pos], a);
152*5971e316Smrg 	isl_int_abs(a, a);
153*5971e316Smrg 	isl_seq_combine(dst, a, dst, b, src, len);
154*5971e316Smrg 
155*5971e316Smrg 	if (m)
156*5971e316Smrg 		isl_int_mul(*m, *m, a);
157*5971e316Smrg 
158*5971e316Smrg 	isl_int_clear(a);
159*5971e316Smrg 	isl_int_clear(b);
160*5971e316Smrg }
161*5971e316Smrg 
isl_seq_eq(isl_int * p1,isl_int * p2,unsigned len)162*5971e316Smrg int isl_seq_eq(isl_int *p1, isl_int *p2, unsigned len)
163*5971e316Smrg {
164*5971e316Smrg 	int i;
165*5971e316Smrg 	for (i = 0; i < len; ++i)
166*5971e316Smrg 		if (isl_int_ne(p1[i], p2[i]))
167*5971e316Smrg 			return 0;
168*5971e316Smrg 	return 1;
169*5971e316Smrg }
170*5971e316Smrg 
isl_seq_cmp(isl_int * p1,isl_int * p2,unsigned len)171*5971e316Smrg int isl_seq_cmp(isl_int *p1, isl_int *p2, unsigned len)
172*5971e316Smrg {
173*5971e316Smrg 	int i;
174*5971e316Smrg 	int cmp;
175*5971e316Smrg 	for (i = 0; i < len; ++i)
176*5971e316Smrg 		if ((cmp = isl_int_cmp(p1[i], p2[i])) != 0)
177*5971e316Smrg 			return cmp;
178*5971e316Smrg 	return 0;
179*5971e316Smrg }
180*5971e316Smrg 
isl_seq_is_neg(isl_int * p1,isl_int * p2,unsigned len)181*5971e316Smrg int isl_seq_is_neg(isl_int *p1, isl_int *p2, unsigned len)
182*5971e316Smrg {
183*5971e316Smrg 	int i;
184*5971e316Smrg 
185*5971e316Smrg 	for (i = 0; i < len; ++i) {
186*5971e316Smrg 		if (isl_int_abs_ne(p1[i], p2[i]))
187*5971e316Smrg 			return 0;
188*5971e316Smrg 		if (isl_int_is_zero(p1[i]))
189*5971e316Smrg 			continue;
190*5971e316Smrg 		if (isl_int_eq(p1[i], p2[i]))
191*5971e316Smrg 			return 0;
192*5971e316Smrg 	}
193*5971e316Smrg 	return 1;
194*5971e316Smrg }
195*5971e316Smrg 
isl_seq_first_non_zero(isl_int * p,unsigned len)196*5971e316Smrg int isl_seq_first_non_zero(isl_int *p, unsigned len)
197*5971e316Smrg {
198*5971e316Smrg 	int i;
199*5971e316Smrg 
200*5971e316Smrg 	for (i = 0; i < len; ++i)
201*5971e316Smrg 		if (!isl_int_is_zero(p[i]))
202*5971e316Smrg 			return i;
203*5971e316Smrg 	return -1;
204*5971e316Smrg }
205*5971e316Smrg 
isl_seq_last_non_zero(isl_int * p,unsigned len)206*5971e316Smrg int isl_seq_last_non_zero(isl_int *p, unsigned len)
207*5971e316Smrg {
208*5971e316Smrg 	int i;
209*5971e316Smrg 
210*5971e316Smrg 	for (i = len - 1; i >= 0; --i)
211*5971e316Smrg 		if (!isl_int_is_zero(p[i]))
212*5971e316Smrg 			return i;
213*5971e316Smrg 	return -1;
214*5971e316Smrg }
215*5971e316Smrg 
isl_seq_abs_max(isl_int * p,unsigned len,isl_int * max)216*5971e316Smrg void isl_seq_abs_max(isl_int *p, unsigned len, isl_int *max)
217*5971e316Smrg {
218*5971e316Smrg 	int i;
219*5971e316Smrg 
220*5971e316Smrg 	isl_int_set_si(*max, 0);
221*5971e316Smrg 
222*5971e316Smrg 	for (i = 0; i < len; ++i)
223*5971e316Smrg 		if (isl_int_abs_gt(p[i], *max))
224*5971e316Smrg 			isl_int_abs(*max, p[i]);
225*5971e316Smrg }
226*5971e316Smrg 
isl_seq_abs_min_non_zero(isl_int * p,unsigned len)227*5971e316Smrg int isl_seq_abs_min_non_zero(isl_int *p, unsigned len)
228*5971e316Smrg {
229*5971e316Smrg 	int i, min = isl_seq_first_non_zero(p, len);
230*5971e316Smrg 	if (min < 0)
231*5971e316Smrg 		return -1;
232*5971e316Smrg 	for (i = min + 1; i < len; ++i) {
233*5971e316Smrg 		if (isl_int_is_zero(p[i]))
234*5971e316Smrg 			continue;
235*5971e316Smrg 		if (isl_int_abs_lt(p[i], p[min]))
236*5971e316Smrg 			min = i;
237*5971e316Smrg 	}
238*5971e316Smrg 	return min;
239*5971e316Smrg }
240*5971e316Smrg 
isl_seq_gcd(isl_int * p,unsigned len,isl_int * gcd)241*5971e316Smrg void isl_seq_gcd(isl_int *p, unsigned len, isl_int *gcd)
242*5971e316Smrg {
243*5971e316Smrg 	int i, min = isl_seq_abs_min_non_zero(p, len);
244*5971e316Smrg 
245*5971e316Smrg 	if (min < 0) {
246*5971e316Smrg 		isl_int_set_si(*gcd, 0);
247*5971e316Smrg 		return;
248*5971e316Smrg 	}
249*5971e316Smrg 	isl_int_abs(*gcd, p[min]);
250*5971e316Smrg 	for (i = 0; isl_int_cmp_si(*gcd, 1) > 0 && i < len; ++i) {
251*5971e316Smrg 		if (i == min)
252*5971e316Smrg 			continue;
253*5971e316Smrg 		if (isl_int_is_zero(p[i]))
254*5971e316Smrg 			continue;
255*5971e316Smrg 		isl_int_gcd(*gcd, *gcd, p[i]);
256*5971e316Smrg 	}
257*5971e316Smrg }
258*5971e316Smrg 
isl_seq_normalize(struct isl_ctx * ctx,isl_int * p,unsigned len)259*5971e316Smrg void isl_seq_normalize(struct isl_ctx *ctx, isl_int *p, unsigned len)
260*5971e316Smrg {
261*5971e316Smrg 	if (len == 0)
262*5971e316Smrg 		return;
263*5971e316Smrg 	isl_seq_gcd(p, len, &ctx->normalize_gcd);
264*5971e316Smrg 	if (!isl_int_is_zero(ctx->normalize_gcd) &&
265*5971e316Smrg 	    !isl_int_is_one(ctx->normalize_gcd))
266*5971e316Smrg 		isl_seq_scale_down(p, p, ctx->normalize_gcd, len);
267*5971e316Smrg }
268*5971e316Smrg 
isl_seq_lcm(isl_int * p,unsigned len,isl_int * lcm)269*5971e316Smrg void isl_seq_lcm(isl_int *p, unsigned len, isl_int *lcm)
270*5971e316Smrg {
271*5971e316Smrg 	int i;
272*5971e316Smrg 
273*5971e316Smrg 	if (len == 0) {
274*5971e316Smrg 		isl_int_set_si(*lcm, 1);
275*5971e316Smrg 		return;
276*5971e316Smrg 	}
277*5971e316Smrg 	isl_int_set(*lcm, p[0]);
278*5971e316Smrg 	for (i = 1; i < len; ++i)
279*5971e316Smrg 		isl_int_lcm(*lcm, *lcm, p[i]);
280*5971e316Smrg }
281*5971e316Smrg 
isl_seq_inner_product(isl_int * p1,isl_int * p2,unsigned len,isl_int * prod)282*5971e316Smrg void isl_seq_inner_product(isl_int *p1, isl_int *p2, unsigned len,
283*5971e316Smrg 			   isl_int *prod)
284*5971e316Smrg {
285*5971e316Smrg 	int i;
286*5971e316Smrg 	if (len == 0) {
287*5971e316Smrg 		isl_int_set_si(*prod, 0);
288*5971e316Smrg 		return;
289*5971e316Smrg 	}
290*5971e316Smrg 	isl_int_mul(*prod, p1[0], p2[0]);
291*5971e316Smrg 	for (i = 1; i < len; ++i)
292*5971e316Smrg 		isl_int_addmul(*prod, p1[i], p2[i]);
293*5971e316Smrg }
294*5971e316Smrg 
isl_seq_hash(isl_int * p,unsigned len,uint32_t hash)295*5971e316Smrg uint32_t isl_seq_hash(isl_int *p, unsigned len, uint32_t hash)
296*5971e316Smrg {
297*5971e316Smrg 	int i;
298*5971e316Smrg 	for (i = 0; i < len; ++i) {
299*5971e316Smrg 		if (isl_int_is_zero(p[i]))
300*5971e316Smrg 			continue;
301*5971e316Smrg 		hash *= 16777619;
302*5971e316Smrg 		hash ^= (i & 0xFF);
303*5971e316Smrg 		hash = isl_int_hash(p[i], hash);
304*5971e316Smrg 	}
305*5971e316Smrg 	return hash;
306*5971e316Smrg }
307*5971e316Smrg 
308*5971e316Smrg /* Given two affine expressions "p" of length p_len (including the
309*5971e316Smrg  * denominator and the constant term) and "subs" of length subs_len,
310*5971e316Smrg  * plug in "subs" for the variable at position "pos".
311*5971e316Smrg  * The variables of "subs" and "p" are assumed to match up to subs_len,
312*5971e316Smrg  * but "p" may have additional variables.
313*5971e316Smrg  * "v" is an initialized isl_int that can be used internally.
314*5971e316Smrg  *
315*5971e316Smrg  * In particular, if "p" represents the expression
316*5971e316Smrg  *
317*5971e316Smrg  *	(a i + g)/m
318*5971e316Smrg  *
319*5971e316Smrg  * with i the variable at position "pos" and "subs" represents the expression
320*5971e316Smrg  *
321*5971e316Smrg  *	f/d
322*5971e316Smrg  *
323*5971e316Smrg  * then the result represents the expression
324*5971e316Smrg  *
325*5971e316Smrg  *	(a f + d g)/(m d)
326*5971e316Smrg  *
327*5971e316Smrg  */
isl_seq_substitute(isl_int * p,int pos,isl_int * subs,int p_len,int subs_len,isl_int v)328*5971e316Smrg void isl_seq_substitute(isl_int *p, int pos, isl_int *subs,
329*5971e316Smrg 	int p_len, int subs_len, isl_int v)
330*5971e316Smrg {
331*5971e316Smrg 	isl_int_set(v, p[1 + pos]);
332*5971e316Smrg 	isl_int_set_si(p[1 + pos], 0);
333*5971e316Smrg 	isl_seq_combine(p + 1, subs[0], p + 1, v, subs + 1, subs_len - 1);
334*5971e316Smrg 	isl_seq_scale(p + subs_len, p + subs_len, subs[0], p_len - subs_len);
335*5971e316Smrg 	isl_int_mul(p[0], p[0], subs[0]);
336*5971e316Smrg }
337*5971e316Smrg 
isl_seq_get_hash(isl_int * p,unsigned len)338*5971e316Smrg uint32_t isl_seq_get_hash(isl_int *p, unsigned len)
339*5971e316Smrg {
340*5971e316Smrg 	uint32_t hash = isl_hash_init();
341*5971e316Smrg 
342*5971e316Smrg 	return isl_seq_hash(p, len, hash);
343*5971e316Smrg }
344*5971e316Smrg 
isl_seq_get_hash_bits(isl_int * p,unsigned len,unsigned bits)345*5971e316Smrg uint32_t isl_seq_get_hash_bits(isl_int *p, unsigned len, unsigned bits)
346*5971e316Smrg {
347*5971e316Smrg 	uint32_t hash;
348*5971e316Smrg 
349*5971e316Smrg 	hash = isl_seq_get_hash(p, len);
350*5971e316Smrg 	return isl_hash_bits(hash, bits);
351*5971e316Smrg }
352*5971e316Smrg 
isl_seq_dump(isl_int * p,unsigned len)353*5971e316Smrg void isl_seq_dump(isl_int *p, unsigned len)
354*5971e316Smrg {
355*5971e316Smrg 	int i;
356*5971e316Smrg 
357*5971e316Smrg 	for (i = 0; i < len; ++i) {
358*5971e316Smrg 		if (i)
359*5971e316Smrg 			fprintf(stderr, " ");
360*5971e316Smrg 		isl_int_print(stderr, p[i], 0);
361*5971e316Smrg 	}
362*5971e316Smrg 	fprintf(stderr, "\n");
363*5971e316Smrg }
364