1 /* Steven Andrews, 10/17/01.
2 Simple routines for manipulating vectors of integers.
3 See documentation called Zn_doc.doc.
4 Copyright 2003-2007 by Steven Andrews. This work is distributed under the terms
5 of the Gnu Lesser General Public License (LGPL). */
6
7 #include <string.h>
8 #include "Zn.h"
9 #include "math2.h"
10 #include "random2.h"
11
setstdZV(int * c,int n,int k)12 int *setstdZV(int *c,int n,int k) {
13 int i;
14
15 if(k==0) for(i=0;i<n;i++) c[i]=0;
16 else if(k==1) for(i=0;i<n;i++) c[i]=1;
17 else if(k<0) {
18 for(i=0;i<n;i++) c[i]=0;
19 c[-k]=1; }
20 else if(k==2) for(i=0;i<n;i++) c[i]=i;
21 else if(k==3) for(i=0;i<n;i++) c[i]=intrand(2);
22 return c; }
23
24
printZV(int * c,int n)25 int *printZV(int *c,int n) {
26 int i;
27
28 for(i=0;i<n;i++) printf("%i ",c[i]);
29 printf("\n");
30 return c; }
31
32
Zn_vect2csvstring(int * c,int n,char * string)33 char *Zn_vect2csvstring(int *c,int n,char *string) {
34 int i;
35
36 if(n>0) snprintf(string, n*sizeof(string),"%i",c[0]);
37 for(i=1;i<n;i++) snprintf(string+strlen(string),sizeof(string)-strlen(string),",%i",c[i]);
38 return string; }
39
40
fprintZV(FILE * stream,int * c,int n)41 int *fprintZV(FILE *stream,int *c,int n) {
42 int i;
43
44 if(n) fprintf(stream,"%i",c[0]);
45 for(i=1;i<n;i++) fprintf(stream," %i",c[i]);
46 fprintf(stream,"\n");
47 return c; }
48
49
minZV(int * a,int n)50 int minZV(int *a,int n) {
51 int min,i;
52
53 min=a[0];
54 for(i=1;i<n;i++) if(a[i]<min) min=a[i];
55 return min; }
56
57
maxZV(int * a,int n)58 int maxZV(int *a,int n) {
59 int max,i;
60
61 max=a[0];
62 for(i=1;i<n;i++) if(a[i]>max) max=a[i];
63 return max; }
64
65
copyZV(int * a,int * c,int n)66 int *copyZV(int *a,int *c,int n) {
67 for(n--;n>=0;n--) c[n]=a[n];
68 return c; }
69
70
sumZV(float ax,int * a,float bx,int * b,int * c,int n)71 int *sumZV(float ax,int *a,float bx,int *b,int *c,int n) {
72 int i;
73
74 for(i=0;i<n;i++) c[i]=(int)(ax*a[i]+bx*b[i]);
75 return c; }
76
77
deriv1ZV(int * a,int * c,int n)78 int *deriv1ZV(int *a,int *c,int n) {
79 int i;
80
81 c[0]=(-3*a[0]+4*a[1]-a[2])/2;
82 for(i=1;i<n-1;i++) c[i]=(a[i+1]-a[i-1])/2;
83 c[n-1]=(a[n-3]-4*a[n-2]+3*a[n-1])/2;
84 return c; }
85
86
deriv2ZV(int * a,int * c,int n)87 int *deriv2ZV(int *a,int *c,int n) {
88 int i;
89
90 c[0]=a[0]+a[2]-2*a[1];
91 for(i=1;i<n-1;i++) c[i]=a[i-1]+a[i+1]-2*a[i];
92 c[n-1]=a[n-3]+a[n-1]-2*a[n-2];
93 return c; }
94
95
productZV(int * a,int n)96 int productZV(int *a,int n) {
97 int i,p;
98
99 p=1;
100 for(i=0;i<n;i++) p*=a[i];
101 return p; }
102
103
intfindZV(int * a,int n,int i)104 int intfindZV(int *a,int n,int i) {
105 int j;
106
107 for(j=0;j<n&&a[j]!=i;j++);
108 return j<n?j:-1; }
109
110
indx2addZV(int * indx,int * dim,int rank)111 int indx2addZV(int *indx,int *dim,int rank) {
112 int i,add;
113
114 add=indx[0];
115 for(i=1;i<rank;i++) add=dim[i]*add+indx[i];
116 return add; }
117
118
add2indxZV(int add,int * indx,int * dim,int rank)119 int *add2indxZV(int add,int *indx,int *dim,int rank) {
120 int i;
121
122 for(i=rank-1;i>0;i--) {
123 indx[i]=add%dim[i];
124 add/=dim[i]; }
125 indx[0]=add;
126 return indx; }
127
128
indx2add3ZV(int * indx,int rank)129 int indx2add3ZV(int *indx,int rank) {
130 int i,add;
131
132 add=indx[0];
133 for(i=1;i<rank;i++) add=3*add+indx[i];
134 return add; }
135
136
nextaddZV(int add,int * indx1,int * indx2,int * dim,int rank)137 int nextaddZV(int add,int *indx1,int *indx2,int *dim,int rank) {
138 int i,done;
139
140 done=0;
141 for(i=rank-1;i>=0&&!done;i--) {
142 if(add%dim[i]<indx2[i]) done=1;
143 else add/=dim[i]; }
144 if(!done) return 1+indx2addZV(indx2,dim,rank);
145 add++;
146 for(i+=2;i<rank;i++) add=dim[i]*add+indx1[i];
147 return add; }
148
149
add2indx3ZV(int add,int * indx,int rank)150 int *add2indx3ZV(int add,int *indx,int rank) {
151 int i;
152
153 for(i=rank-1;i>0;i--) {
154 indx[i]=add%3;
155 add/=3; }
156 indx[0]=add;
157 return indx; }
158
159
neighborZV(int * indx,int * c,int * dim,int rank,int type,int * wrap,int * mid)160 int neighborZV(int *indx,int *c,int *dim,int rank,int type,int *wrap,int *mid) {
161 int i,n,p,ad,k,j;
162 static int *ofst=NULL,*wp=NULL,rankold=0;
163
164 n=0;
165 if(type==-1||rank!=rankold) {
166 if(ofst) freeZV(ofst);
167 if(wp) freeZV(wp);
168 ofst=wp=NULL;
169 rankold=rank; }
170
171 if(type==0) {
172 p=productZV(dim,rank);
173 ad=indx2addZV(indx,dim,rank);
174 for(i=0;i<rank;i++) {
175 p/=dim[i];
176 if(indx[i]>0) c[n++]=ad-p; }
177 if(mid) *mid=n;
178 for(i=rank-1;i>=0;i--) {
179 if(indx[i]<dim[i]-1) c[n++]=ad+p;
180 p*=dim[i]; }}
181 else if(type==1) {
182 p=productZV(dim,rank);
183 ad=indx2addZV(indx,dim,rank);
184 for(i=0;i<rank;i++) {
185 p/=dim[i];
186 if(indx[i]>0) c[n++]=ad-p;
187 else c[n++]=ad+(dim[i]-1)*p;
188 for(k=0;k<n-1;k++) if(c[k]==c[n-1]) n--; }
189 if(mid) *mid=n;
190 for(i=rank-1;i>=0;i--) {
191 if(indx[i]<dim[i]-1) c[n++]=ad+p;
192 else c[n++]=ad-(dim[i]-1)*p;
193 p*=dim[i];
194 for(k=0;k<n-1;k++) if(c[k]==c[n-1]) n--; }}
195 else if(type==2) {
196 if(!ofst) ofst=allocZV(rank);
197 if(!ofst) return -1;
198 p=intpower(3,rank);
199 for(i=0;i<p;i++) {
200 if(i==p/2) {if(mid) *mid=n;i++;}
201 add2indx3ZV(i,ofst,rank);
202 for(j=0;j<rank;j++) {
203 ofst[j]+=indx[j]-1;
204 if(ofst[j]<0||ofst[j]>=dim[j]) j=rank+1; }
205 if(j==rank) c[n++]=indx2addZV(ofst,dim,rank); }}
206 else if(type==3) {
207 if(!ofst) ofst=allocZV(rank);
208 if(!ofst) return -1;
209 p=intpower(3,rank);
210 for(i=0;i<p;i++) {
211 if(i==p/2) {if(mid) *mid=n;i++;}
212 add2indx3ZV(i,ofst,rank);
213 for(j=0;j<rank;j++) {
214 ofst[j]+=indx[j]-1;
215 if(ofst[j]<0) ofst[j]=dim[j]-1;
216 else if(ofst[j]>=dim[j]) ofst[j]=0; }
217 c[n++]=indx2addZV(ofst,dim,rank);
218 for(k=0;k<n-1;k++) if(c[k]==c[n-1]) n--; }}
219 else if(type==4) {
220 p=productZV(dim,rank);
221 ad=indx2addZV(indx,dim,rank);
222 for(i=0;i<rank;i++) {
223 p/=dim[i];
224 if(indx[i]>0) c[n++]=ad-p;
225 else if(wrap[i]) c[n++]=ad+(dim[i]-1)*p;
226 for(k=0;k<n-1;k++) if(c[k]==c[n-1]) n--; }
227 if(mid) *mid=n;
228 for(i=rank-1;i>=0;i--) {
229 if(indx[i]<dim[i]-1) c[n++]=ad+p;
230 else if(wrap[i]) c[n++]=ad-(dim[i]-1)*p;
231 p*=dim[i];
232 for(k=0;k<n-1;k++) if(c[k]==c[n-1]) n--; }}
233 else if(type==5) {
234 if(!ofst) ofst=allocZV(rank);
235 if(!ofst) return -1;
236 p=intpower(3,rank);
237 for(i=0;i<p;i++) {
238 if(i==p/2) {if(mid) *mid=n;i++;}
239 add2indx3ZV(i,ofst,rank);
240 for(j=0;j<rank;j++) {
241 ofst[j]+=indx[j]-1;
242 if(ofst[j]<0&&wrap[j]) ofst[j]=dim[j]-1;
243 else if(ofst[j]>=dim[j]&&wrap[j]) ofst[j]=0;
244 else if(ofst[j]<0||ofst[j]>=dim[j]) j=rank+1; }
245 if(j==rank) c[n++]=indx2addZV(ofst,dim,rank);
246 for(k=0;k<n-1;k++) if(c[k]==c[n-1]) n--; }}
247 else if(type==6) {
248 p=productZV(dim,rank);
249 ad=indx2addZV(indx,dim,rank);
250 if(!wp) wp=allocZV(rank);
251 if(!wp) return -1;
252 copyZV(wrap,wp,rank);
253 setstdZV(wrap,2*rank,0);
254 for(i=0;i<rank;i++) {
255 p/=dim[i];
256 if(indx[i]>0) c[n++]=ad-p;
257 else if(wp[i]) {
258 wrap[n]=1<<2*i;
259 c[n++]=ad+(dim[i]-1)*p; }}
260 if(mid) *mid=n;
261 for(i=rank-1;i>=0;i--) {
262 if(indx[i]<dim[i]-1) c[n++]=ad+p;
263 else if(wp[i]) {
264 wrap[n]=2<<2*i;
265 c[n++]=ad-(dim[i]-1)*p; }
266 p*=dim[i]; }}
267 else if(type==7) {
268 if(!wp) wp=allocZV(rank);
269 if(!ofst) ofst=allocZV(rank);
270 if(!ofst||!wp) return -1;
271 p=intpower(3,rank);
272 copyZV(wrap,wp,rank);
273 for(i=0;i<p;i++) {
274 wrap[n]=0;
275 if(i==p/2) {if(mid) *mid=n;i++;}
276 add2indx3ZV(i,ofst,rank);
277 for(j=0;j<rank;j++) {
278 ofst[j]+=indx[j]-1;
279 if(ofst[j]<0&&wp[j]) {
280 ofst[j]=dim[j]-1;wrap[n]+=1<<2*j; }
281 else if(ofst[j]>=dim[j]&&wp[j]) {
282 ofst[j]=0;wrap[n]+=2<<2*j; }
283 else if(ofst[j]<0||ofst[j]>=dim[j]) j=rank+1; }
284 if(j==rank) c[n++]=indx2addZV(ofst,dim,rank); }}
285 return n; }
286
287
Zn_sameset(const int * a,const int * b,int * work,int n)288 int Zn_sameset(const int *a,const int *b,int *work,int n) {
289 int i,j,sum;
290
291 for(j=0;j<n;j++) work[j]=0;
292 sum=0;
293 for(i=0;i<n;i++)
294 for(j=0;j<n;j++)
295 if(b[j]==a[i] && work[j]==0) {
296 work[j]=1;
297 sum++;
298 j=n; }
299 return (sum==n); }
300
301
Zn_sort(int * a,int * b,int n)302 void Zn_sort(int *a,int *b,int n) {
303 int i,j,ir,l;
304 int az,bz;
305
306 if(!b) b=a;
307 for(i=0;i<n-1&&a[i]<a[i+1];i++)
308 ;
309 if(i==n-1) return;
310 for(i=0;i<n-1&&a[i]>a[i+1];i++)
311 ;
312 if(i==n-1) {
313 for(i=0;i<n/2;i++) {
314 az=a[i];
315 bz=b[i];
316 a[i]=a[j=n-i-1];
317 b[i]=b[j];
318 a[j]=az;
319 b[j]=bz; }
320 return; }
321 l=(n>>1)+1;
322 ir=n;
323 for(;;) {
324 if(l>1) {
325 az=a[--l-1];
326 bz=b[l-1]; }
327 else {
328 az=a[ir-1];
329 bz=b[ir-1];
330 a[ir-1]=a[0];
331 b[ir-1]=b[0];
332 if(--ir==1) {
333 a[0]=az;
334 b[0]=bz;
335 return; }}
336 i=l;
337 j=l<<1;
338 while(j<=ir) {
339 if(j<ir&&a[j-1]<a[j]) ++j;
340 if(az<a[j-1]) {
341 a[i-1]=a[j-1];
342 b[i-1]=b[j-1];
343 j+=(i=j); }
344 else
345 j=ir+1; }
346 a[i-1]=az;
347 b[i-1]=bz; }}
348
349
350
Zn_issort(int * a,int n)351 int Zn_issort(int *a,int n) {
352 int i;
353
354 if(n<=1) return 1;
355 for(i=1;i<n&&a[i-1]==a[i];i++);
356 if(i==n) return 1;
357 for(i=1;i<n&&a[i-1]<a[i];i++);
358 if(i==n) return 3;
359 for(;i<n&&a[i-1]<=a[i];i++);
360 if(i==n) return 2;
361 for(i=1;i<n&&a[i-1]>a[i];i++);
362 if(i==n) return -3;
363 for(;i<n&&a[i-1]>=a[i];i++);
364 if(i==n) return -2;
365 return 0; }
366
367
Zn_incrementcounter(int * a,int digits,int base)368 int Zn_incrementcounter(int *a,int digits,int base) {
369 int i;
370
371 i=0;
372 a[0]++;
373 while(a[i]==base) {
374 a[i]=0;
375 i++;
376 if(i==digits) return 1;
377 a[i]++; }
378 return 0; }
379
380
Zn_permute(int * a,int * b,int n,int k)381 int Zn_permute(int *a,int *b,int n,int k) {
382 int ans;
383
384 if(n==0) ans=0;
385
386 else if(n==1) {b[0]=a[0];ans=0; }
387
388 else if(n==2) {
389 if(k==0) {b[0]=a[0];b[1]=a[1];ans=1; }
390 else {b[0]=a[1];b[1]=a[0];ans=0; }
391 if(a[0]==a[1]) ans=0; }
392
393 else if(n==3) {
394 if(k==0) {b[0]=a[0];b[1]=a[1];b[2]=a[2];ans=1;}
395 else if(k==1) {b[0]=a[0];b[1]=a[2];b[2]=a[1];ans=2;}
396 else if(k==2) {b[0]=a[1];b[1]=a[0];b[2]=a[2];ans=3;}
397 else if(k==3) {b[0]=a[1];b[1]=a[2];b[2]=a[0];ans=4;}
398 else if(k==4) {b[0]=a[2];b[1]=a[0];b[2]=a[1];ans=5;}
399 else {b[0]=a[2];b[1]=a[1];b[2]=a[0];ans=0;}
400 if(a[1]==a[2]&&ans==1) ans=2;
401 else if(a[1]==a[2]&&(ans==4||ans==5)) ans=0;
402 if(a[0]==a[1]&&(ans==2||ans==3)) ans=4;
403 else if(a[0]==a[1]&&ans==5) ans=0;
404 if(a[0]==a[2]&&(ans==3||ans==4||ans==5)) ans=0; }
405
406 else ans=-1;
407
408 return ans; }
409
410
411 /******** combinatorics *********/
412
413
Zn_permutelex(int * seq,int n)414 int Zn_permutelex(int *seq,int n) {
415 int i,j;
416 int temp;
417
418 i=n-1;
419 while(i>0 && seq[i-1]>=seq[i]) i--;
420 if(i==0) { // input was final sequence
421 i=0;
422 j=n-1;
423 while(i<j) {
424 temp=seq[i]; // swap values at positions i and j
425 seq[i]=seq[j];
426 seq[j]=temp;
427 i++;
428 j--; }
429 return 2; }
430
431 j=n;
432 while(seq[j-1]<=seq[i-1]) j--;
433
434 temp=seq[i-1]; // swap values at positions (i-1) and (j-1)
435 seq[i-1]=seq[j-1];
436 seq[j-1]=temp;
437
438 i++;
439 j=n;
440 while(i<j) {
441 temp=seq[i-1]; // swap values at positions (i-1) and (j-1)
442 seq[i-1]=seq[j-1];
443 seq[j-1]=temp;
444 i++;
445 j--; }
446
447 i=n-1;
448 while(i>0 && seq[i-1]>=seq[i]) i--;
449 if(i==0) return 1; // at final sequence
450
451 return 0; }
452
453
454
455
456
457
458
459
460
461
462
463
464
465