1 /*****
2 * arrayop
3 * John Bowman
4 *
5 * Array operations
6 *****/
7 #ifndef ARRAYOP_H
8 #define ARRAYOP_H
9
10 #include "util.h"
11 #include "stack.h"
12 #include "array.h"
13 #include "types.h"
14 #include "fileio.h"
15 #include "callable.h"
16 #include "mathop.h"
17
18 namespace run {
19
20 using vm::pop;
21 using vm::read;
22 using vm::array;
23 using camp::tab;
24
25 vm::array *copyArray(vm::array *a);
26 vm::array *copyArray2(vm::array *a);
27
28 template<class T, class U, template <class S> class op>
arrayOp(vm::stack * s)29 void arrayOp(vm::stack *s)
30 {
31 U b=pop<U>(s);
32 array *a=pop<array*>(s);
33 size_t size=checkArray(a);
34 array *c=new array(size);
35 for(size_t i=0; i < size; i++)
36 (*c)[i]=op<T>()(read<T>(a,i),b,i);
37 s->push(c);
38 }
39
40 template<class T, class U, template <class S> class op>
opArray(vm::stack * s)41 void opArray(vm::stack *s)
42 {
43 array *a=pop<array*>(s);
44 T b=pop<T>(s);
45 size_t size=checkArray(a);
46 array *c=new array(size);
47 for(size_t i=0; i < size; i++)
48 (*c)[i]=op<U>()(b,read<U>(a,i),i);
49 s->push(c);
50 }
51
52 template<class T, template <class S> class op>
arrayArrayOp(vm::stack * s)53 void arrayArrayOp(vm::stack *s)
54 {
55 array *b=pop<array*>(s);
56 array *a=pop<array*>(s);
57 size_t size=checkArrays(a,b);
58 array *c=new array(size);
59 for(size_t i=0; i < size; i++)
60 (*c)[i]=op<T>()(read<T>(a,i),read<T>(b,i),i);
61 s->push(c);
62 }
63
64 template<class T>
sumArray(vm::stack * s)65 void sumArray(vm::stack *s)
66 {
67 array *a=pop<array*>(s);
68 size_t size=checkArray(a);
69 T sum=0;
70 for(size_t i=0; i < size; i++)
71 sum += read<T>(a,i);
72 s->push(sum);
73 }
74
75 extern const char *arrayempty;
76
77 template<class T, template <class S> class op>
binopArray(vm::stack * s)78 void binopArray(vm::stack *s)
79 {
80 array *a=pop<array*>(s);
81 size_t size=checkArray(a);
82 if(size == 0) vm::error(arrayempty);
83 T m=read<T>(a,0);
84 for(size_t i=1; i < size; i++)
85 m=op<T>()(m,read<T>(a,i));
86 s->push(m);
87 }
88
89 template<class T, template <class S> class op>
binopArray2(vm::stack * s)90 void binopArray2(vm::stack *s)
91 {
92 array *a=pop<array*>(s);
93 size_t size=checkArray(a);
94 bool empty=true;
95 T m=0;
96 for(size_t i=0; i < size; i++) {
97 array *ai=read<array*>(a,i);
98 size_t aisize=checkArray(ai);
99 if(aisize) {
100 if(empty) {
101 m=read<T>(ai,0);
102 empty=false;
103 }
104 for(size_t j=0; j < aisize; j++)
105 m=op<T>()(m,read<T>(ai,j));
106 }
107 }
108 if(empty) vm::error(arrayempty);
109 s->push(m);
110 }
111
112 template<class T, template <class S> class op>
binopArray3(vm::stack * s)113 void binopArray3(vm::stack *s)
114 {
115 array *a=pop<array*>(s);
116 size_t size=checkArray(a);
117 bool empty=true;
118 T m=0;
119 for(size_t i=0; i < size; i++) {
120 array *ai=read<array*>(a,i);
121 size_t aisize=checkArray(ai);
122 for(size_t j=0; j < aisize; j++) {
123 array *aij=read<array*>(ai,j);
124 size_t aijsize=checkArray(aij);
125 if(aijsize) {
126 if(empty) {
127 m=read<T>(aij,0);
128 empty=false;
129 }
130 for(size_t k=0; k < aijsize; k++) {
131 m=op<T>()(m,read<T>(aij,k));
132 }
133 }
134 }
135 }
136 if(empty) vm::error(arrayempty);
137 s->push(m);
138 }
139
140 template<class T, class U, template <class S> class op>
array2Op(vm::stack * s)141 void array2Op(vm::stack *s)
142 {
143 U b=pop<U>(s);
144 array *a=pop<array*>(s);
145 size_t size=checkArray(a);
146 array *c=new array(size);
147 for(size_t i=0; i < size; ++i) {
148 array *ai=read<array*>(a,i);
149 size_t aisize=checkArray(ai);
150 array *ci=new array(aisize);
151 (*c)[i]=ci;
152 for(size_t j=0; j < aisize; j++)
153 (*ci)[j]=op<T>()(read<T>(ai,j),b,0);
154 }
155 s->push(c);
156 }
157
158 template<class T, class U, template <class S> class op>
opArray2(vm::stack * s)159 void opArray2(vm::stack *s)
160 {
161 array *a=pop<array*>(s);
162 T b=pop<T>(s);
163 size_t size=checkArray(a);
164 array *c=new array(size);
165 for(size_t i=0; i < size; ++i) {
166 array *ai=read<array*>(a,i);
167 size_t aisize=checkArray(ai);
168 array *ci=new array(aisize);
169 (*c)[i]=ci;
170 for(size_t j=0; j < aisize; j++)
171 (*ci)[j]=op<U>()(read<U>(ai,j),b,0);
172 }
173 s->push(c);
174 }
175
176 template<class T, template <class S> class op>
array2Array2Op(vm::stack * s)177 void array2Array2Op(vm::stack *s)
178 {
179 array *b=pop<array*>(s);
180 array *a=pop<array*>(s);
181 size_t size=checkArrays(a,b);
182 array *c=new array(size);
183 for(size_t i=0; i < size; ++i) {
184 array *ai=read<array*>(a,i);
185 array *bi=read<array*>(b,i);
186 size_t aisize=checkArrays(ai,bi);
187 array *ci=new array(aisize);
188 (*c)[i]=ci;
189 for(size_t j=0; j < aisize; j++)
190 (*ci)[j]=op<T>()(read<T>(ai,j),read<T>(bi,j),0);
191 }
192 s->push(c);
193 }
194
195 template<class T>
Array2Equals(vm::stack * s)196 bool Array2Equals(vm::stack *s)
197 {
198 array *b=pop<array*>(s);
199 array *a=pop<array*>(s);
200 size_t n=checkArray(a);
201 if(n != checkArray(b)) return false;
202 if(n == 0) return true;
203 size_t n0=checkArray(read<array*>(a,0));
204 if(n0 != checkArray(read<array*>(b,0))) return false;
205
206 for(size_t i=0; i < n; ++i) {
207 array *ai=read<array*>(a,i);
208 array *bi=read<array*>(b,i);
209 for(size_t j=0; j < n0; ++j) {
210 if(read<T>(ai,j) != read<T>(bi,j))
211 return false;
212 }
213 }
214 return true;
215 }
216
217 template<class T>
array2Equals(vm::stack * s)218 void array2Equals(vm::stack *s)
219 {
220 s->push(Array2Equals<T>(s));
221 }
222
223 template<class T>
array2NotEquals(vm::stack * s)224 void array2NotEquals(vm::stack *s)
225 {
226 s->push(!Array2Equals<T>(s));
227 }
228
229 template<class T>
diagonal(vm::stack * s)230 void diagonal(vm::stack *s)
231 {
232 array *a=pop<array*>(s);
233 size_t n=checkArray(a);
234 array *c=new array(n);
235 for(size_t i=0; i < n; ++i) {
236 array *ci=new array(n);
237 (*c)[i]=ci;
238 for(size_t j=0; j < i; ++j)
239 (*ci)[j]=T();
240 (*ci)[i]=read<T>(a,i);
241 for(size_t j=i+1; j < n; ++j)
242 (*ci)[j]=T();
243 }
244 s->push(c);
245 }
246
247 template<class T>
248 struct compare {
operatorcompare249 bool operator() (const vm::item& a, const vm::item& b)
250 {
251 return vm::get<T>(a) < vm::get<T>(b);
252 }
253 };
254
255 template<class T>
sortArray(vm::stack * s)256 void sortArray(vm::stack *s)
257 {
258 array *c=copyArray(pop<array*>(s));
259 sort(c->begin(),c->end(),compare<T>());
260 s->push(c);
261 }
262
263 template<class T>
264 struct compare2 {
operatorcompare2265 bool operator() (const vm::item& A, const vm::item& B)
266 {
267 array *a=vm::get<array*>(A);
268 array *b=vm::get<array*>(B);
269 size_t size=a->size();
270 if(size != b->size()) return false;
271
272 for(size_t j=0; j < size; j++) {
273 if(read<T>(a,j) < read<T>(b,j)) return true;
274 if(read<T>(a,j) > read<T>(b,j)) return false;
275 }
276 return false;
277 }
278 };
279
280 // Sort the rows of a 2-dimensional array by the first column, breaking
281 // ties with successively higher columns.
282 template<class T>
sortArray2(vm::stack * s)283 void sortArray2(vm::stack *s)
284 {
285 array *c=copyArray(pop<array*>(s));
286 stable_sort(c->begin(),c->end(),compare2<T>());
287 s->push(c);
288 }
289
290 // Search a sorted ordered array a of n elements for key, returning the index i
291 // if a[i] <= key < a[i+1], -1 if key is less than all elements of a, or
292 // n-1 if key is greater than or equal to the last element of a.
293 template<class T>
searchArray(vm::stack * s)294 void searchArray(vm::stack *s)
295 {
296 T key=pop<T>(s);
297 array *a=pop<array*>(s);
298 size_t size= a->size();
299 if(size == 0 || key < read<T>(a,0)) {s->push(-1); return;}
300 size_t u=size-1;
301 if(key >= read<T>(a,u)) {s->push((Int) u); return;}
302 size_t l=0;
303
304 while (l < u) {
305 size_t i=(l+u)/2;
306 if(key < read<T>(a,i)) u=i;
307 else if(key < read<T>(a,i+1)) {s->push((Int) i); return;}
308 else l=i+1;
309 }
310 s->push(0);
311 }
312
313 extern string emptystring;
314
315 void writestring(vm::stack *s);
316
317 template<class T>
write(vm::stack * s)318 void write(vm::stack *s)
319 {
320 array *a=pop<array*>(s);
321 vm::callable *suffix=pop<vm::callable *>(s,NULL);
322 T first=pop<T>(s);
323 string S=pop<string>(s,emptystring);
324 vm::item it=pop(s);
325 bool defaultfile=isdefault(it);
326 camp::ofile *f=defaultfile ? &camp::Stdout : vm::get<camp::ofile*>(it);
327 if(!f->isOpen() || !f->enabled()) return;
328
329 size_t size=checkArray(a);
330 if(S != "") f->write(S);
331 f->write(first);
332 for(size_t i=0; i < size; ++i) {
333 f->write(tab);
334 f->write(read<T>(a,i));
335 }
336 if(f->text()) {
337 if(suffix) {
338 s->push(f);
339 suffix->call(s);
340 } else if(defaultfile) {
341 try {
342 f->writeline();
343 } catch (quit&) {
344 }
345 }
346 }
347 }
348
349 template<class T>
writeArray(vm::stack * s)350 void writeArray(vm::stack *s)
351 {
352 array *A=pop<array*>(s);
353 array *a=pop<array*>(s);
354 string S=pop<string>(s,emptystring);
355 vm::item it=pop(s);
356 bool defaultfile=isdefault(it);
357 camp::ofile *f=defaultfile ? &camp::Stdout : vm::get<camp::ofile*>(it);
358 if(!f->isOpen() || !f->enabled()) return;
359
360 size_t asize=checkArray(a);
361 size_t Asize=checkArray(A);
362 if(f->Standard()) interact::lines=0;
363 else if(!f->isOpen()) return;
364 try {
365 if(S != "") {f->write(S); f->writeline();}
366
367 size_t i=0;
368 bool cont=true;
369 while(cont) {
370 cont=false;
371 bool first=true;
372 if(i < asize) {
373 vm::item& I=(*a)[i];
374 if(defaultfile) cout << i << ":\t";
375 if(!I.empty())
376 f->write(vm::get<T>(I));
377 cont=true;
378 first=false;
379 }
380 unsigned count=0;
381 for(size_t j=0; j < Asize; ++j) {
382 array *Aj=read<array*>(A,j);
383 size_t Ajsize=checkArray(Aj);
384 if(i < Ajsize) {
385 if(f->text()) {
386 if(first && defaultfile) cout << i << ":\t";
387 for(unsigned k=0; k <= count; ++k)
388 f->write(tab);
389 count=0;
390 }
391 vm::item& I=(*Aj)[i];
392 if(!I.empty())
393 f->write(vm::get<T>(I));
394 cont=true;
395 first=false;
396 } else count++;
397 }
398 ++i;
399 if(cont && f->text()) f->writeline();
400 }
401 } catch (quit&) {
402 }
403 f->flush();
404 }
405
406 template<class T>
writeArray2(vm::stack * s)407 void writeArray2(vm::stack *s)
408 {
409 array *a=pop<array*>(s);
410 vm::item it=pop(s);
411 bool defaultfile=isdefault(it);
412 camp::ofile *f=defaultfile ? &camp::Stdout : vm::get<camp::ofile*>(it);
413 if(!f->isOpen() || !f->enabled()) return;
414
415 size_t size=checkArray(a);
416 if(f->Standard()) interact::lines=0;
417
418 try {
419 for(size_t i=0; i < size; i++) {
420 vm::item& I=(*a)[i];
421 if(!I.empty()) {
422 array *ai=vm::get<array*>(I);
423 size_t aisize=checkArray(ai);
424 for(size_t j=0; j < aisize; j++) {
425 if(j > 0 && f->text()) f->write(tab);
426 vm::item& I=(*ai)[j];
427 if(!I.empty())
428 f->write(vm::get<T>(I));
429 }
430 }
431 if(f->text()) f->writeline();
432 }
433 } catch (quit&) {
434 }
435 f->flush();
436 }
437
438 template<class T>
writeArray3(vm::stack * s)439 void writeArray3(vm::stack *s)
440 {
441 array *a=pop<array*>(s);
442 vm::item it=pop(s);
443 bool defaultfile=isdefault(it);
444 camp::ofile *f=defaultfile ? &camp::Stdout : vm::get<camp::ofile*>(it);
445 if(!f->isOpen() || !f->enabled()) return;
446
447 size_t size=checkArray(a);
448 if(f->Standard()) interact::lines=0;
449
450 try {
451 for(size_t i=0; i < size;) {
452 vm::item& I=(*a)[i];
453 if(!I.empty()) {
454 array *ai=vm::get<array*>(I);
455 size_t aisize=checkArray(ai);
456 for(size_t j=0; j < aisize; j++) {
457 vm::item& I=(*ai)[j];
458 if(!I.empty()) {
459 array *aij=vm::get<array*>(I);
460 size_t aijsize=checkArray(aij);
461 for(size_t k=0; k < aijsize; k++) {
462 if(k > 0 && f->text()) f->write(tab);
463 vm::item& I=(*aij)[k];
464 if(!I.empty())
465 f->write(vm::get<T>(I));
466 }
467 }
468 if(f->text()) f->writeline();
469 }
470 }
471 ++i;
472 if(i < size && f->text()) f->writeline();
473 }
474 } catch (quit&) {
475 }
476 f->flush();
477 }
478
479 template <class T, class S, T (*func)(S)>
arrayFunc(vm::stack * s)480 void arrayFunc(vm::stack *s)
481 {
482 array *a=pop<array*>(s);
483 size_t size=checkArray(a);
484 array *c=new array(size);
485 for(size_t i=0; i < size; i++)
486 (*c)[i]=func(read<S>(a,i));
487 s->push(c);
488 }
489
490 template <class T, class S, T (*func)(S)>
arrayFunc2(vm::stack * s)491 void arrayFunc2(vm::stack *s)
492 {
493 array *a=pop<array*>(s);
494 size_t size=checkArray(a);
495 array *c=new array(size);
496 for(size_t i=0; i < size; ++i) {
497 array *ai=read<array*>(a,i);
498 size_t aisize=checkArray(ai);
499 array *ci=new array(aisize);
500 (*c)[i]=ci;
501 for(size_t j=0; j < aisize; j++)
502 (*ci)[j]=func(read<S>(ai,j));
503 }
504 s->push(c);
505 }
506
507 vm::array *Identity(Int n);
508 camp::triple operator *(const vm::array& a, const camp::triple& v);
509 double norm(double *a, size_t n);
510 double norm(camp::triple *a, size_t n);
511
checkdimension(const vm::array * a,size_t dim)512 inline size_t checkdimension(const vm::array *a, size_t dim)
513 {
514 size_t size=checkArray(a);
515 if(dim && size != dim) {
516 ostringstream buf;
517 buf << "array of length " << dim << " expected";
518 vm::error(buf);
519 }
520 return size;
521 }
522
523 template<class T>
524 inline void copyArrayC(T* &dest, const vm::array *a, size_t dim=0,
525 GCPlacement placement=NoGC)
526 {
527 size_t size=checkdimension(a,dim);
528 dest=(placement == NoGC) ? new T[size] : new(placement) T[size];
529 for(size_t i=0; i < size; i++)
530 dest[i]=vm::read<T>(a,i);
531 }
532
533 template<class T, class A>
534 inline void copyArrayC(T* &dest, const vm::array *a, T (*cast)(A),
535 size_t dim=0, GCPlacement placement=NoGC)
536 {
537 size_t size=checkdimension(a,dim);
538 dest=(placement == NoGC) ? new T[size] : new(placement) T[size];
539 for(size_t i=0; i < size; i++)
540 dest[i]=cast(vm::read<A>(a,i));
541 }
542
543 template<typename T>
copyCArray(const size_t n,const T * p)544 inline vm::array* copyCArray(const size_t n, const T* p)
545 {
546 vm::array* a = new vm::array(n);
547 for(size_t i=0; i < n; ++i) (*a)[i] = p[i];
548 return a;
549 }
550
551 template<class T>
552 inline void copyArray2C(T* &dest, const vm::array *a, bool square=true,
553 size_t dim2=0, GCPlacement placement=NoGC)
554 {
555 size_t n=checkArray(a);
556 size_t m=(square || n == 0) ? n : checkArray(vm::read<vm::array*>(a,0));
557 if(n > 0 && dim2 && m != dim2) {
558 ostringstream buf;
559 buf << "second matrix dimension must be " << dim2;
560 vm::error(buf);
561 }
562
563 dest=(placement == NoGC) ? new T[n*m] : new(placement) T[n*m];
564 for(size_t i=0; i < n; i++) {
565 vm::array *ai=vm::read<vm::array*>(a,i);
566 size_t aisize=checkArray(ai);
567 if(aisize == m) {
568 T *desti=dest+i*m;
569 for(size_t j=0; j < m; j++)
570 desti[j]=vm::read<T>(ai,j);
571 } else
572 vm::error(square ? "matrix must be square" :
573 "matrix must be rectangular");
574 }
575 }
576
577 template<typename T>
copyCArray2(const size_t n,const size_t m,const T * p)578 inline vm::array* copyCArray2(const size_t n, const size_t m, const T* p)
579 {
580 vm::array* a=new vm::array(n);
581 for(size_t i=0; i < n; ++i) {
582 array *ai=new array(m);
583 (*a)[i]=ai;
584 for(size_t j=0; j < m; ++j)
585 (*ai)[j]=p[m*i+j];
586 }
587 return a;
588 }
589
590 } // namespace run
591
592 #endif // ARRAYOP_H
593