YAP 7.1.0
newm.c
1/*************************************************************************
2 * *
3 * YAP Prolog *
4 * *
5 * Yap Prolog was developed at NCCUP - Universidade do Porto *
6 * *
7 * Copyright L.Damas, V.S.Costa and Universidade do Porto 1985-1997 *
8 * *
9 **************************************************************************
10 * *
11 * File: matrix.c * Last rev:
12 ** mods: * comments: numerical arrays *
13 * *
14 *************************************************************************/
15
16#include "YapConfig.h"
17#include "YapInterface.h"
18#include <math.h>
19#if defined(__MINGW32__) || _MSC_VER
20#include <windows.h>
21#endif
22#if HAVE_STRING_H
23#include <string.h>
24#endif
25
26/*
27 A matrix is something of the form
28
29 TYPE = {int,double}
30 BASE = integer
31 #DIMS = an int
32 DIM1
33 ...
34 DIMn
35 DATA in C format.
36
37 floating point matrixes may need to be aligned, so we always have an
38 extra element at the end.
39*/
40
41/* maximal number of dimensions, 1024 should be enough */
42#define MAX_DIMS 1024
43
44typedef enum { INT_MATRIX, FLOAT_MATRIX } mat_data_type;
45
46typedef enum {
47 MAT_TYPE = 0,
48 MAT_BASE = 1,
49 MAT_NDIMS = 2,
50 MAT_SIZE = 3,
51 MAT_ALIGN = 4,
52 MAT_DIMS = 5,
53} mat_type;
54
55typedef enum {
56 MAT_PLUS = 0,
57 MAT_SUB = 1,
58 MAT_TIMES = 2,
59 MAT_DIV = 3,
60 MAT_IDIV = 4,
61 MAT_ZDIV = 5,
62 MAT_LOG = 6,
63 MAT_EXP = 7
64} op_type;
65
66
67typedef struct m {
68 bool floats, c_ord;
69 union {
70 double *data;
71 intptr_t *ls;
72 };
73 int ndims;
74 size_t sz;
75 size_t *dims, base;
76} M;
77
78
79
80YAP_Functor FunctorM, FunctorFloats;
81YAP_Atom AtomC;
82
83static YAP_Int *matrix_long_data(int *mat, int ndims) {
84 return (YAP_Int *)(mat + (MAT_DIMS + ndims));
85}
86
87static double *matrix_double_data(int *mat, int ndims) {
88 return (double *)(mat + (MAT_DIMS + ndims));
89}
90
91static unsigned int matrix_get_offset(int *mat, int *indx) {
92 unsigned int i, pos = mat[MAT_SIZE], off = 0;
93
94 /* find where we are */
95 for (i = 0; i < mat[MAT_NDIMS]; i++) {
96 int v;
97 pos /= mat[MAT_DIMS + i];
98 v = indx[i] - mat[MAT_BASE];
99 if (v >= mat[MAT_DIMS + i]) {
100 return off;
101 }
102 off += pos * v;
103 }
104 return off;
105}
106
107static void matrix_get_index(int *mat, unsigned int offset, int *indx) {
108 unsigned int i, pos = mat[MAT_SIZE];
109
110 /* find where we are */
111 for (i = 0; i < mat[MAT_NDIMS]; i++) {
112 pos /= mat[MAT_DIMS + i];
113 indx[i] = offset / pos;
114 offset = offset % pos;
115 }
116}
117
118static size_t GET_OFFSET(YAP_Term t, M *mat) {
119 unsigned int i, off = 0;
120 if (YAP_IsIntTerm(YAP_ARG1))
121 return YAP_IntOfTerm(YAP_ARG1);
122 /* find where we are */
123 for (i = 0; i < mat->ndims; i++) {
124 int pos, v;
125 if (!YAP_IsPairTerm(t))
126 YAP_Error(TYPE_ERROR, t, NULL);
127 pos /= mat->dims[i];
128 v = YAP_IntOfTerm(YAP_HeadOfTerm(t)) - mat->base;
129 if (v >= mat->dims[ i]) {
130 return off;
131 }
132 off += pos * v;
133 }
134 return off;
135}
136
137static void matrix_next_index(int *dims, int ndims, int *indx) {
138 unsigned int i;
139
140 /* find where we are */
141 for (i = ndims; i > 0;) {
142 i--;
143 indx[i]++;
144 if (indx[i] != dims[i])
145 return;
146 indx[i] = 0;
147 }
148}
149
150
151static bool GET_MATRIX(YAP_Term inp, M *o)
152 {
153 int *mat;
154 if ((mat= YAP_BlobOfTerm(inp))) {
155 o->floats = mat[MAT_TYPE]==FLOAT_MATRIX;
156 o->c_ord = true;
157 o->sz = mat[MAT_SIZE];
158 o->ndims = mat[MAT_NDIMS];
159 o->data = (double *)(mat + (MAT_DIMS + o->ndims));
160 o->dims=mat+MAT_DIMS;
161 return true;
162 } else if (YAP_IsAtomTerm(inp) && (
163 o->data = YAP_FetchArray(inp, &o->sz, &o->floats) )) {
164 o->ndims = 1;
165 o->dims = &o->sz;
166 o->c_ord = true;
167 return true;
168 } else {
169 o->sz = YAP_IntOfTerm(YAP_ArgOfTerm(1,inp));
170 o->c_ord = true;
171 o->ndims = 1;
172 o->dims = &o->sz;
173 o->floats = true;
174 o->data = (double *)YAP_IntOfTerm(YAP_ArgOfTerm(2,inp));
175 return true;
176 }
177 return false;
178 }
179
180
181
182static YAP_Term new_int_matrix(int ndims, int dims[], YAP_Int data[]) {
183 unsigned int sz;
184 unsigned int i, nelems = 1;
185 YAP_Term blob;
186 int *mat;
187 YAP_Int *bdata;
188 int idims[MAX_DIMS];
189
190 /* in case we don't have enough room and need to shift the stack, we can't
191 really afford to keep a pointer to the global */
192 for (i = 0; i < ndims; i++) {
193 idims[i] = dims[i];
194 nelems *= dims[i];
195 }
196 sz = ((MAT_DIMS + 1) * sizeof(int) + ndims * sizeof(int) +
197 nelems * sizeof(YAP_Int)) /
198 sizeof(YAP_CELL);
199
200 blob = YAP_MkBlobTerm(sz);
201 if (blob == YAP_TermNil()) {
202 return blob;
203 }
204 mat = (int *)YAP_BlobOfTerm(blob);
205 mat[MAT_TYPE] = INT_MATRIX;
206 mat[MAT_BASE] = 0;
207 mat[MAT_NDIMS] = ndims;
208 mat[MAT_SIZE] = nelems;
209 for (i = 0; i < ndims; i++) {
210 mat[MAT_DIMS + i] = idims[i];
211 }
212 bdata = matrix_long_data(mat, ndims);
213 if (data)
214 memmove((void *)bdata, (void *)data, sizeof(double) * nelems);
215 return blob;
216}
217
218static YAP_Term new_float_matrix(int ndims, int dims[], double data[]) {
219 unsigned int sz;
220 unsigned int i, nelems = 1;
221 YAP_Term blob;
222 int *mat;
223 double *bdata;
224 int idims[MAX_DIMS];
225
226 /* in case we don't have enough room and need to shift the stack, we can't
227 really afford to keep a pointer to the global */
228 for (i = 0; i < ndims; i++) {
229 idims[i] = dims[i];
230 nelems *= dims[i];
231 }
232 sz = ((MAT_DIMS + 1) * sizeof(int) + ndims * sizeof(int) +
233 (nelems + 1) * sizeof(double) + (sizeof(YAP_CELL) - 1)) /
234 sizeof(YAP_CELL);
235 blob = YAP_MkBlobTerm(sz);
236 if (blob == YAP_TermNil())
237 return blob;
238 mat = YAP_BlobOfTerm(blob);
239 mat[MAT_TYPE] = FLOAT_MATRIX;
240 mat[MAT_BASE] = 0;
241 mat[MAT_NDIMS] = ndims;
242 mat[MAT_SIZE] = nelems;
243 for (i = 0; i < ndims; i++) {
244 mat[MAT_DIMS + i] = idims[i];
245 }
246 bdata = matrix_double_data(mat, ndims);
247 if (data)
248 memmove((void *)bdata, (void *)data, sizeof(double) * nelems);
249 return blob;
250}
251
252static YAP_Bool scan_dims(int ndims, YAP_Term tl, int dims[MAX_DIMS]) {
253 int i;
254
255 for (i = 0; i < ndims; i++) {
256 YAP_Term th;
257 int d;
258
259 if (!YAP_IsPairTerm(tl)) {
260 return FALSE;
261 }
262 th = YAP_HeadOfTerm(tl);
263 if (!YAP_IsIntTerm(th)) {
264 /* ERROR */
265 return FALSE;
266 }
267 d = YAP_IntOfTerm(th);
268 if (d < 0) {
269 /* ERROR */
270 return FALSE;
271 }
272 dims[i] = d;
273 tl = YAP_TailOfTerm(tl);
274 }
275 if (tl != YAP_TermNil()) {
276 /* ERROR */
277 return FALSE;
278 }
279 return TRUE;
280}
281
282static YAP_Bool scan_dims_args(int ndims, YAP_Term tl, int dims[MAX_DIMS]) {
283 int i;
284
285 for (i = 0; i < ndims; i++) {
286 YAP_Term th;
287 int d;
288
289 th = YAP_ArgOfTerm(2 + i, tl);
290 if (!YAP_IsIntTerm(th)) {
291 /* ERROR */
292 return FALSE;
293 }
294 d = YAP_IntOfTerm(th);
295 if (d < 0) {
296 /* ERROR */
297 return FALSE;
298 }
299 dims[i] = d;
300 }
301 return TRUE;
302}
303
304static YAP_Bool cp_int_matrix(YAP_Term tl, YAP_Term matrix) {
305 int *mat = (int *)YAP_BlobOfTerm(matrix);
306 int i, nelems = mat[MAT_SIZE];
307 YAP_Int *j = matrix_long_data(mat, mat[MAT_NDIMS]);
308
309 for (i = 0; i < nelems; i++) {
310 YAP_Term th;
311 int d;
312
313 if (!YAP_IsPairTerm(tl)) {
314 return FALSE;
315 }
316 th = YAP_HeadOfTerm(tl);
317 if (!YAP_IsIntTerm(th)) {
318 /* ERROR */
319 return FALSE;
320 }
321 d = YAP_IntOfTerm(th);
322 j[i] = d;
323 tl = YAP_TailOfTerm(tl);
324 }
325 if (tl != YAP_TermNil()) {
326 /* ERROR */
327 return FALSE;
328 }
329 return TRUE;
330}
331
332static YAP_Bool cp_float_matrix(YAP_Term tl, YAP_Term matrix) {
333 int *mat = (int *)YAP_BlobOfTerm(matrix);
334 int i, nelems = mat[MAT_SIZE];
335 double *j = matrix_double_data(mat, mat[MAT_NDIMS]);
336
337 for (i = 0; i < nelems; i++) {
338 YAP_Term th;
339 double d;
340
341 if (!YAP_IsPairTerm(tl)) {
342 return FALSE;
343 }
344 th = YAP_HeadOfTerm(tl);
345 if (YAP_IsIntTerm(th)) {
346 d = YAP_IntOfTerm(th);
347 } else if (!YAP_IsFloatTerm(th)) {
348 /* ERROR */
349 return FALSE;
350 } else {
351 d = YAP_FloatOfTerm(th);
352 }
353 j[i] = d;
354 tl = YAP_TailOfTerm(tl);
355 }
356 if (tl != YAP_TermNil()) {
357 /* ERROR */
358 return FALSE;
359 }
360 return TRUE;
361}
362
363static YAP_Bool set_int_matrix(YAP_Term matrix, YAP_Int set) {
364 int *mat = (int *)YAP_BlobOfTerm(matrix);
365 int i, nelems = mat[MAT_SIZE];
366 YAP_Int *j = matrix_long_data(mat, mat[MAT_NDIMS]);
367
368 for (i = 0; i < nelems; i++) {
369 j[i] = set;
370 }
371 return TRUE;
372}
373
374static YAP_Bool set_float_matrix(YAP_Term matrix, double set) {
375 int *mat = (int *)YAP_BlobOfTerm(matrix);
376 int i, nelems = mat[MAT_SIZE];
377 double *j = matrix_double_data(mat, mat[MAT_NDIMS]);
378
379 for (i = 0; i < nelems; i++) {
380 j[i] = set;
381 }
382 return TRUE;
383}
384
385static YAP_Bool new_ints_matrix(void) {
386 int ndims = YAP_IntOfTerm(YAP_ARG1);
387 YAP_Term tl = YAP_ARG2, out;
388 int dims[MAX_DIMS];
389 YAP_Term data;
390
391 if (!scan_dims(ndims, tl, dims))
392 return FALSE;
393 out = new_int_matrix(ndims, dims, NULL);
394 if (out == YAP_TermNil())
395 return FALSE;
396 data = YAP_ARG3;
397 if (!YAP_IsVarTerm(data) && !cp_int_matrix(data, out))
398 return FALSE;
399 return YAP_Unify(YAP_ARG4, out);
400}
401
402static YAP_Bool new_ints_matrix_set(void) {
403 int ndims = YAP_IntOfTerm(YAP_ARG1);
404 YAP_Term tl = YAP_ARG2, out, tset = YAP_ARG3;
405 int dims[MAX_DIMS];
406 YAP_Int set;
407
408 if (!YAP_IsIntTerm(tset)) {
409 return FALSE;
410 }
411 set = YAP_IntOfTerm(tset);
412 if (!scan_dims(ndims, tl, dims))
413 return FALSE;
414 size_t sz = 1, i;
415 for (i=0;i<ndims;i++)
416 sz*=dims[i];
417 YAP_RequiresExtraStack(sz*sizeof(YAP_CELL));
418 out = new_int_matrix(ndims, dims, NULL);
419 return set_int_matrix(out, set) &&
420 YAP_Unify(YAP_ARG4, out);
421}
422
423static YAP_Bool new_floats_matrix(void) {
424 int ndims = YAP_IntOfTerm(YAP_ARG1);
425 YAP_Term tl = YAP_ARG2, out, data = YAP_ARG3;
426 int dims[MAX_DIMS];
427 if (!scan_dims(ndims, tl, dims))
428 return FALSE;
429 out = new_float_matrix(ndims, dims, NULL);
430 if (out == YAP_TermNil())
431 return FALSE;
432 if (!YAP_IsVarTerm(data) && !cp_float_matrix(data, out))
433 return FALSE;
434 return YAP_Unify(YAP_ARG4, out);
435}
436
437static YAP_Bool new_floats_matrix_set(void) {
438 int ndims = YAP_IntOfTerm(YAP_ARG1);
439 YAP_Term tl = YAP_ARG2, out, tset = YAP_ARG3;
440 int dims[MAX_DIMS];
441 double set;
442
443 if (!YAP_IsFloatTerm(tset)) {
444 return FALSE;
445 }
446 set = YAP_FloatOfTerm(tset);
447 if (!scan_dims(ndims, tl, dims))
448 return FALSE;
449 out = new_float_matrix(ndims, dims, NULL);
450 if (!set_float_matrix(out, set))
451 return FALSE;
452 return YAP_Unify(YAP_ARG4, out);
453}
454
455static YAP_Term float_matrix_to_list(int *mat) {
456 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
457
458 /* prepare for worst case with double taking two cells */
459 if (YAP_RequiresExtraStack(6 * mat[MAT_SIZE])) {
460 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
461 data = matrix_double_data(mat, mat[MAT_NDIMS]);
462 }
463 return YAP_FloatsToList(data, mat[MAT_SIZE]);
464}
465
466static YAP_Term mk_int_list(int nelems, int *data) {
467 YAP_Term tn = YAP_TermNil();
468 YAP_Term tf = tn;
469 int i = 0;
470
471 for (i = nelems - 1; i >= 0; i--) {
472 tf = YAP_MkPairTerm(YAP_MkIntTerm(data[i]), tf);
473 if (tf == tn) {
474 /* error */
475 return tn;
476 }
477 }
478 return tf;
479}
480
481static YAP_Term mk_int_list2(int nelems, int base, int *data) {
482 YAP_Term tn = YAP_TermNil();
483 YAP_Term tf = tn;
484 int i = 0;
485
486 for (i = nelems - 1; i >= 0; i--) {
487 tf = YAP_MkPairTerm(YAP_MkIntTerm(data[i] + base), tf);
488 if (tf == tn) {
489 /* error */
490 return tn;
491 }
492 }
493 return tf;
494}
495
496static YAP_Term mk_rep_int_list(int nelems, int data) {
497 YAP_Term tn = YAP_TermNil();
498 YAP_Term tf = tn;
499 int i = 0;
500
501 for (i = nelems - 1; i >= 0; i--) {
502 tf = YAP_MkPairTerm(YAP_MkIntTerm(data), tf);
503 if (tf == tn) {
504 /* error */
505 return tn;
506 }
507 }
508 return tf;
509}
510
511static YAP_Term mk_long_list(int nelems, YAP_Int *data) {
512 YAP_Term tn = YAP_TermNil();
513 YAP_Term tf = tn;
514 int i = 0;
515
516 for (i = nelems - 1; i >= 0; i--) {
517 tf = YAP_MkPairTerm(YAP_MkIntTerm(data[i]), tf);
518 if (tf == tn) {
519 /* error */
520 return tn;
521 }
522 }
523 return tf;
524}
525
526static YAP_Term long_matrix_to_list(int *mat) {
527 YAP_Int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
528
529 /* prepare for worst case with longs evrywhere (3cells + 1) */
530 if (YAP_RequiresExtraStack(5 * mat[MAT_SIZE])) {
531 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
532 data = matrix_long_data(mat, mat[MAT_NDIMS]);
533 }
534 return mk_long_list(mat[MAT_SIZE], data);
535}
536
537static YAP_Term matrix_access(int *mat, int *indx) {
538 unsigned int off = matrix_get_offset(mat, indx);
539 if (mat[MAT_TYPE] == FLOAT_MATRIX)
540 return YAP_MkFloatTerm((matrix_double_data(mat, mat[MAT_NDIMS]))[off]);
541 else
542 return YAP_MkIntTerm((matrix_long_data(mat, mat[MAT_NDIMS]))[off]);
543}
544
545static void matrix_float_set(int *mat, int *indx, double nval) {
546 unsigned int off = 0;
547
548 off = matrix_get_offset(mat, indx);
549 (matrix_double_data(mat, mat[MAT_NDIMS]))[off] = nval;
550}
551
552static void matrix_long_set(int *mat, int *indx, YAP_Int nval) {
553 unsigned int off = matrix_get_offset(mat, indx);
554 (matrix_long_data(mat, mat[MAT_NDIMS]))[off] = nval;
555}
556
557static void matrix_float_set_all(int *mat, double nval) {
558 int i;
559 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
560
561 for (i = 0; i < mat[MAT_SIZE]; i++)
562 data[i] = nval;
563}
564
565static void matrix_long_set_all(int *mat, YAP_Int nval) {
566 int i;
567 YAP_Int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
568
569 if (nval == 0) {
570 memset((void *)data, 0, sizeof(YAP_Int) * mat[MAT_SIZE]);
571 } else {
572 for (i = 0; i < mat[MAT_SIZE]; i++)
573 data[i] = nval;
574 }
575}
576
577static void matrix_float_add(int *mat, int *indx, double nval) {
578 unsigned int off;
579 double *dat = matrix_double_data(mat, mat[MAT_NDIMS]);
580
581 off = matrix_get_offset(mat, indx);
582 dat[off] += nval;
583}
584
585static void matrix_long_add(int *mat, int *indx, YAP_Int nval) {
586 YAP_Int *dat = matrix_long_data(mat, mat[MAT_NDIMS]);
587 unsigned int off = matrix_get_offset(mat, indx);
588 dat[off] += nval;
589}
590
591static void matrix_inc(int *mat, int *indx) {
592 unsigned int off = matrix_get_offset(mat, indx);
593 if (mat[MAT_TYPE] == FLOAT_MATRIX)
594 (matrix_double_data(mat, mat[MAT_NDIMS])[off])++;
595 else
596 ((matrix_long_data(mat, mat[MAT_NDIMS]))[off])++;
597}
598
599static void matrix_dec(int *mat, int *indx) {
600 unsigned int off = matrix_get_offset(mat, indx);
601 if (mat[MAT_TYPE] == FLOAT_MATRIX)
602 (matrix_double_data(mat, mat[MAT_NDIMS])[off])--;
603 else
604 ((matrix_long_data(mat, mat[MAT_NDIMS]))[off])--;
605}
606
607static YAP_Term matrix_inc2(int *mat, int *indx) {
608 unsigned int off = matrix_get_offset(mat, indx);
609 if (mat[MAT_TYPE] == FLOAT_MATRIX) {
610 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
611 double d = data[off];
612 d++;
613 data[off] = d;
614 return YAP_MkFloatTerm(d);
615 } else {
616 YAP_Int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
617 YAP_Int d = data[off];
618 d++;
619 data[off] = d;
620 return YAP_MkIntTerm(d);
621 }
622}
623
624static YAP_Term matrix_dec2(int *mat, int *indx) {
625 unsigned int off = matrix_get_offset(mat, indx);
626 if (mat[MAT_TYPE] == FLOAT_MATRIX) {
627 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
628 double d = data[off];
629 d--;
630 data[off] = d;
631 return YAP_MkFloatTerm(d);
632 } else {
633 YAP_Int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
634 YAP_Int d = data[off];
635 d--;
636 data[off] = d;
637 return YAP_MkIntTerm(d);
638 }
639}
640
641
642static YAP_Bool matrix_set2(void) {
643 int dims[MAX_DIMS], *mat;
644 YAP_Term tf, t = YAP_ARG1;
645
646 mat = (int *)YAP_BlobOfTerm(YAP_ArgOfTerm(1, t));
647 if (!mat) {
648 /* Error */
649 return FALSE;
650 }
651 if (!scan_dims_args(mat[MAT_NDIMS], t, dims)) {
652 /* Error */
653 return FALSE;
654 }
655 tf = YAP_ARG2;
656 if (mat[MAT_TYPE] == INT_MATRIX) {
657 if (YAP_IsIntTerm(tf)) {
658 matrix_long_set(mat, dims, YAP_IntOfTerm(tf));
659 } else if (YAP_IsFloatTerm(tf)) {
660 matrix_long_set(mat, dims, YAP_FloatOfTerm(tf));
661 } else {
662 /* Error */
663 return FALSE;
664 }
665 } else {
666 if (YAP_IsIntTerm(tf)) {
667 matrix_float_set(mat, dims, YAP_IntOfTerm(tf));
668 } else if (YAP_IsFloatTerm(tf)) {
669 matrix_float_set(mat, dims, YAP_FloatOfTerm(tf));
670 } else {
671 /* Error */
672 return FALSE;
673 }
674 }
675 return TRUE;
676}
677
678static YAP_Bool matrix_set_all(void) {
679 int *mat;
680 YAP_Term tf;
681
682 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
683 if (!mat) {
684 if (YAP_IsApplTerm(YAP_ARG1) && YAP_FunctorOfTerm(YAP_ARG1) == FunctorFloats) {
685 double * v = YAP_PointerOfTerm(YAP_ArgOfTerm(1,YAP_ARG1));
686 size_t sz = YAP_IntOfTerm(YAP_ArgOfTerm(2,YAP_ARG1));
687 double val = YAP_FloatOfTerm(YAP_ARG2);
688 if (val == 0.0) memset(v, 0,sz*sizeof(double));
689 else {
690 int j;
691 for (j=0; j< sz; j++)
692 v[j] = val;
693
694 }
695
696 }
697 /* Error */
698 return FALSE;
699 }
700 tf = YAP_ARG2;
701 if (mat[MAT_TYPE] == INT_MATRIX) {
702 if (YAP_IsIntTerm(tf)) {
703 matrix_long_set_all(mat, YAP_IntOfTerm(tf));
704 } else if (YAP_IsFloatTerm(tf)) {
705 matrix_long_set_all(mat, YAP_FloatOfTerm(tf));
706 } else {
707 /* Error */
708 return FALSE;
709 }
710 } else {
711 if (YAP_IsIntTerm(tf)) {
712 matrix_float_set_all(mat, YAP_IntOfTerm(tf));
713 } else if (YAP_IsFloatTerm(tf)) {
714 matrix_float_set_all(mat, YAP_FloatOfTerm(tf));
715 } else {
716 /* Error */
717 return FALSE;
718 }
719 }
720 return TRUE;
721}
722
723static YAP_Bool matrix_add(void) {
724 int dims[MAX_DIMS], *mat;
725 YAP_Term tf;
726
727 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
728 if (!mat) {
729 /* Error */
730 return FALSE;
731 }
732 if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, dims)) {
733 /* Error */
734 return FALSE;
735 }
736 tf = YAP_ARG3;
737 if (mat[MAT_TYPE] == INT_MATRIX) {
738 if (YAP_IsIntTerm(tf)) {
739 matrix_long_add(mat, dims, YAP_IntOfTerm(tf));
740 } else if (YAP_IsFloatTerm(tf)) {
741 matrix_long_add(mat, dims, YAP_FloatOfTerm(tf));
742 } else {
743 /* Error */
744 return FALSE;
745 }
746 } else {
747 if (YAP_IsIntTerm(tf)) {
748 matrix_float_add(mat, dims, YAP_IntOfTerm(tf));
749 } else if (YAP_IsFloatTerm(tf)) {
750 matrix_float_add(mat, dims, YAP_FloatOfTerm(tf));
751 } else {
752 /* Error */
753 return FALSE;
754 }
755 }
756 return TRUE;
757}
758
759static YAP_Bool matrix_get(void) {
760 M mat;
761 YAP_Term tf;
762 int offset;
763
764 if (!GET_MATRIX(YAP_ARG1 ,&mat) ||
765 !(offset = GET_OFFSET(YAP_ARG2,&mat))) {
766 /* Error */
767 return false;
768 }
769 if (mat.floats)
770 tf = YAP_MkFloatTerm(mat.data[offset]);
771 else
772 tf = YAP_MkIntTerm(mat.ls[offset]);
773 return YAP_Unify(tf, YAP_ARG3);
774}
775
776static YAP_Bool matrix_set(void) {
777 M mat;
778 YAP_Term tf;
779 int offset;
780
781 if (!GET_MATRIX(YAP_ARG1 ,&mat) ||
782 !(offset = GET_OFFSET(YAP_ARG2,&mat))) {
783 /* Error */
784 return false;
785 }
786 if (mat.floats) {
787 if (YAP_IsIntTerm(YAP_ARG3))
788 mat.data[offset] = YAP_IntOfTerm(YAP_ARG3);
789 else
790 mat.data[offset] = YAP_FloatOfTerm(YAP_ARG3);
791 } else {
792 if (YAP_IsIntTerm(YAP_ARG3))
793 mat.ls[offset] = YAP_IntOfTerm(YAP_ARG3);
794 else
795 mat.ls[offset] = YAP_FloatOfTerm(YAP_ARG3);
796 }
797 return true;
798}
799
800static YAP_Bool do_matrix_access2(void) {
801 int dims[MAX_DIMS], *mat;
802 YAP_Term tf, t = YAP_ARG1;
803
804 mat = (int *)YAP_BlobOfTerm(YAP_ArgOfTerm(1, t));
805 if (!mat) {
806 /* Error */
807 return FALSE;
808 }
809 if (!scan_dims_args(mat[MAT_NDIMS], t, dims)) {
810 /* Error */
811 return FALSE;
812 }
813 tf = matrix_access(mat, dims);
814 return YAP_Unify(tf, YAP_ARG2);
815
816}
817
818static YAP_Bool do_matrix_inc(void) {
819 M mat;
820 int offset;
821 if (!GET_MATRIX(YAP_ARG1 ,&mat) ||
822 !(offset = GET_OFFSET(YAP_ARG2,&mat))) {
823 /* Error */
824 return false;
825 }
826 if (mat.floats)
827 mat.data[offset] += 1.0;
828 else
829 mat.ls[offset] += 1;
830 return true;
831}
832
833static YAP_Bool do_matrix_dec(void) {
834 M mat;
835 int offset;
836 if (!GET_MATRIX(YAP_ARG1 ,&mat) ||
837 !(offset = GET_OFFSET(YAP_ARG2,&mat))) {
838 /* Error */
839 return false;
840 }
841 if (mat.floats)
842 mat.data[offset] -= 1.0;
843 else
844 mat.ls[offset] -= 1;
845 return true;
846}
847
848static YAP_Bool do_matrix_inc2(void) {
849 int dims[MAX_DIMS], *mat;
850
851 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
852 if (!mat) {
853 /* Error */
854 return FALSE;
855 }
856 if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, dims)) {
857 /* Error */
858 return FALSE;
859 }
860 return YAP_Unify(matrix_inc2(mat, dims), YAP_ARG3);
861}
862
863static YAP_Bool do_matrix_dec2(void) {
864 int dims[MAX_DIMS], *mat;
865
866 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
867 if (!mat) {
868 /* Error */
869 return FALSE;
870 }
871 if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, dims)) {
872 /* Error */
873 return FALSE;
874 }
875 return YAP_Unify(matrix_dec2(mat, dims), YAP_ARG3);
876}
877
878static YAP_Bool matrix_to_list(void) {
879 int *mat;
880 YAP_Term tf;
881
882 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
883 if (!mat) {
884 /* Error */
885 return FALSE;
886 }
887 if (mat[MAT_TYPE] == INT_MATRIX)
888 tf = long_matrix_to_list(mat);
889 else
890 tf = float_matrix_to_list(mat);
891 return YAP_Unify(YAP_ARG2, tf);
892}
893
894static YAP_Bool matrix_set_base(void) {
895 int *mat;
896
897 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
898 if (!mat) {
899 /* Error */
900 return FALSE;
901 }
902 mat[MAT_BASE] = YAP_IntOfTerm(YAP_ARG2);
903 return TRUE;
904}
905
906static YAP_Bool matrix_dims(void) {
907 int *mat;
908 YAP_Term tf;
909
910 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
911 if (!mat) {
912 /* Error */
913 return FALSE;
914 }
915 tf = mk_int_list(mat[MAT_NDIMS], mat + MAT_DIMS);
916 return YAP_Unify(YAP_ARG2, tf);
917}
918
919static YAP_Bool matrix_dims3(void) {
920 int *mat;
921 YAP_Term tf, tof;
922
923 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
924 if (!mat) {
925 /* Error */
926 return FALSE;
927 }
928 tf = mk_int_list(mat[MAT_NDIMS], mat + MAT_DIMS);
929 tof = mk_rep_int_list(mat[MAT_NDIMS], mat[MAT_BASE]);
930 return YAP_Unify(YAP_ARG2, tf) && YAP_Unify(YAP_ARG3, tof);
931}
932
933static YAP_Bool matrix_size(void) {
934 int *mat;
935
936 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
937 if (!mat) {
938 /* Error */
939 return FALSE;
940 }
941 return YAP_Unify(YAP_ARG2, YAP_MkIntTerm(mat[MAT_SIZE]));
942}
943
944static YAP_Bool matrix_ndims(void) {
945 int *mat;
946
947 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
948 if (!mat) {
949 /* Error */
950 return FALSE;
951 }
952 return YAP_Unify(YAP_ARG2, YAP_MkIntTerm(mat[MAT_NDIMS]));
953}
954
955static YAP_Bool matrix_type(void) {
956 int *mat;
957 YAP_Term tf;
958
959 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
960 if (!mat) {
961 /* not an error, it may be called on a term matrix */
962 return FALSE;
963 }
964 if (mat[MAT_TYPE] == INT_MATRIX) {
965 tf = YAP_MkIntTerm(0);
966 } else {
967 tf = YAP_MkIntTerm(1);
968 }
969 return YAP_Unify(YAP_ARG2, tf);
970}
971
972static YAP_Bool matrix_arg_to_offset(void) {
973 int indx[MAX_DIMS], *mat;
974 unsigned int off;
975
976 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
977 if (!mat) {
978 /* Error */
979 return FALSE;
980 }
981 if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, indx)) {
982 /* Error */
983 return FALSE;
984 }
985 off = matrix_get_offset(mat, indx);
986
987 return YAP_Unify(YAP_ARG3, YAP_MkIntTerm(off));
988}
989
990static YAP_Bool matrix_offset_to_arg(void) {
991 int indx[MAX_DIMS], *mat;
992 unsigned int off;
993 YAP_Term ti, tf;
994
995 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
996 if (!mat) {
997 /* Error */
998 return FALSE;
999 }
1000 if (!YAP_IsIntTerm(ti = YAP_ARG2)) {
1001 /* Error */
1002 return FALSE;
1003 }
1004 off = YAP_IntOfTerm(ti);
1005 matrix_get_index(mat, off, indx);
1006 tf = mk_int_list2(mat[MAT_NDIMS], mat[MAT_BASE], indx);
1007 return YAP_Unify(YAP_ARG3, tf);
1008}
1009
1010static unsigned int scan_max_long(int sz, YAP_Int *data) {
1011 int i, off = 0;
1012 YAP_Int max = data[0];
1013 for (i = 1; i < sz; i++) {
1014 if (data[i] > max) {
1015 off = i;
1016 max = data[i];
1017 }
1018 }
1019 return off;
1020}
1021
1022static unsigned int scan_max_float(int sz, double *data) {
1023 int i, off = 0;
1024 double max = data[0];
1025 for (i = 1; i < sz; i++) {
1026 if (data[i] > max) {
1027 max = data[i];
1028 off = i;
1029 }
1030 }
1031 return off;
1032}
1033
1034static unsigned int scan_min_long(int sz, YAP_Int *data) {
1035 int i, off = 0;
1036 YAP_Int max = data[0];
1037 for (i = 1; i < sz; i++) {
1038 if (data[i] < max) {
1039 max = data[i];
1040 off = i;
1041 }
1042 }
1043 return off;
1044}
1045
1046static unsigned int scan_min_float(int sz, double *data) {
1047 int i, off = 0;
1048 double max = data[0];
1049 for (i = 1; i < sz; i++) {
1050 if (data[i] < max) {
1051 max = data[i];
1052 off = i;
1053 }
1054 }
1055 return off;
1056}
1057
1058static YAP_Bool matrix_max(void) {
1059 int *mat;
1060 unsigned int off;
1061 YAP_Term tf;
1062
1063 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1064 if (!mat) {
1065 /* Error */
1066 return FALSE;
1067 }
1068 if (mat[MAT_TYPE] == INT_MATRIX) {
1069 YAP_Int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
1070 off = scan_max_long(mat[MAT_SIZE], data);
1071 tf = YAP_MkIntTerm(data[off]);
1072 } else {
1073 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
1074 off = scan_max_float(mat[MAT_SIZE], data);
1075 tf = YAP_MkFloatTerm(data[off]);
1076 }
1077 return YAP_Unify(YAP_ARG2, tf);
1078}
1079
1080static YAP_Bool matrix_maxarg(void) {
1081 int indx[MAX_DIMS], *mat;
1082 unsigned int off;
1083 YAP_Term tf;
1084
1085 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1086 if (!mat) {
1087 /* Error */
1088 return FALSE;
1089 }
1090 if (mat[MAT_TYPE] == INT_MATRIX) {
1091 YAP_Int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
1092 off = scan_max_long(mat[MAT_SIZE], data);
1093 } else {
1094 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
1095 off = scan_max_float(mat[MAT_SIZE], data);
1096 }
1097 matrix_get_index(mat, off, indx);
1098 tf = mk_int_list(mat[MAT_NDIMS], indx);
1099 return YAP_Unify(YAP_ARG2, tf);
1100}
1101
1102static YAP_Bool matrix_min(void) {
1103 int *mat;
1104 unsigned int off;
1105 YAP_Term tf;
1106
1107 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1108 if (!mat) {
1109 /* Error */
1110 return FALSE;
1111 }
1112 if (mat[MAT_TYPE] == INT_MATRIX) {
1113 YAP_Int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
1114 off = scan_min_long(mat[MAT_SIZE], data);
1115 tf = YAP_MkIntTerm(data[off]);
1116 } else {
1117 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
1118 off = scan_min_float(mat[MAT_SIZE], data);
1119 tf = YAP_MkFloatTerm(data[off]);
1120 }
1121 return YAP_Unify(YAP_ARG2, tf);
1122}
1123
1124static YAP_Bool matrix_log_all(void) {
1125 int *mat;
1126
1127 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1128 if (!mat) {
1129 /* Error */
1130 return FALSE;
1131 }
1132 if (mat[MAT_TYPE] == INT_MATRIX) {
1133 return FALSE;
1134 } else {
1135 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
1136 int i;
1137
1138 for (i = 0; i < mat[MAT_SIZE]; i++) {
1139 data[i] = log(data[i]);
1140 }
1141 }
1142 return TRUE;
1143}
1144
1145static YAP_Bool matrix_log_all2(void) {
1146 int *mat;
1147
1148 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1149 if (!mat) {
1150 /* Error */
1151 return FALSE;
1152 }
1153 if (mat[MAT_TYPE] == INT_MATRIX) {
1154 YAP_Term out;
1155 YAP_Int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
1156 double *ndata;
1157 int i;
1158 int *nmat;
1159
1160 if (!YAP_IsVarTerm(YAP_ARG2)) {
1161 out = YAP_ARG2;
1162 } else {
1163 out = new_float_matrix(mat[MAT_NDIMS], mat + MAT_DIMS, NULL);
1164 if (out == YAP_TermNil())
1165 return FALSE;
1166 }
1167 nmat = (int *)YAP_BlobOfTerm(out);
1168 ndata = matrix_double_data(nmat, mat[MAT_NDIMS]);
1169 for (i = 0; i < mat[MAT_SIZE]; i++) {
1170 ndata[i] = log((double)data[i]);
1171 }
1172 if (YAP_IsVarTerm(YAP_ARG2)) {
1173 return YAP_Unify(YAP_ARG2, out);
1174 }
1175 } else {
1176 YAP_Term out;
1177 double *data = matrix_double_data(mat, mat[MAT_NDIMS]), *ndata;
1178 int i;
1179 int *nmat;
1180
1181 if (!YAP_IsVarTerm(YAP_ARG2)) {
1182 out = YAP_ARG2;
1183 } else {
1184 out = new_float_matrix(mat[MAT_NDIMS], mat + MAT_DIMS, NULL);
1185 if (out == YAP_TermNil())
1186 return FALSE;
1187 }
1188 nmat = (int *)YAP_BlobOfTerm(out);
1189 ndata = matrix_double_data(nmat, mat[MAT_NDIMS]);
1190 for (i = 0; i < mat[MAT_SIZE]; i++) {
1191 ndata[i] = log(data[i]);
1192 }
1193 if (YAP_IsVarTerm(YAP_ARG2)) {
1194 return YAP_Unify(YAP_ARG2, out);
1195 }
1196 }
1197 return TRUE;
1198}
1199
1200static YAP_Bool matrix_exp_all(void) {
1201 int *mat;
1202
1203 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1204 if (!mat) {
1205 /* Error */
1206 return FALSE;
1207 }
1208 if (mat[MAT_TYPE] == INT_MATRIX) {
1209 return FALSE;
1210 } else {
1211 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
1212 int i;
1213
1214 for (i = 0; i < mat[MAT_SIZE]; i++) {
1215 data[i] = exp(data[i]);
1216 }
1217 }
1218 return TRUE;
1219}
1220
1221static YAP_Bool matrix_exp2_all(void) {
1222 int *mat;
1223
1224 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1225 if (!mat) {
1226 /* Error */
1227 return FALSE;
1228 }
1229 if (mat[MAT_TYPE] == INT_MATRIX) {
1230 return FALSE;
1231 } else {
1232 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
1233 int i;
1234 double max = data[0];
1235
1236 for (i = 1; i < mat[MAT_SIZE]; i++) {
1237 if (data[i] > max)
1238 max = data[i];
1239 }
1240
1241 for (i = 0; i < mat[MAT_SIZE]; i++) {
1242 data[i] = exp(data[i] - max);
1243 }
1244 }
1245 return TRUE;
1246}
1247
1248static YAP_Bool matrix_exp_all2(void) {
1249 int *mat;
1250
1251 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1252 if (!mat) {
1253 /* Error */
1254 return FALSE;
1255 }
1256 if (mat[MAT_TYPE] == INT_MATRIX) {
1257 YAP_Term out;
1258 YAP_Int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
1259 double *ndata;
1260 int i;
1261 int *nmat;
1262
1263 if (!YAP_IsVarTerm(YAP_ARG2)) {
1264 out = YAP_ARG2;
1265 } else {
1266 out = new_float_matrix(mat[MAT_NDIMS], mat + MAT_DIMS, NULL);
1267 if (out == YAP_TermNil())
1268 return FALSE;
1269 }
1270 nmat = (int *)YAP_BlobOfTerm(out);
1271 ndata = matrix_double_data(nmat, mat[MAT_NDIMS]);
1272 for (i = 0; i < mat[MAT_SIZE]; i++) {
1273 ndata[i] = exp((double)data[i]);
1274 }
1275 if (YAP_IsVarTerm(YAP_ARG2)) {
1276 return YAP_Unify(YAP_ARG2, out);
1277 }
1278 } else {
1279 YAP_Term out;
1280 double *data = matrix_double_data(mat, mat[MAT_NDIMS]), *ndata;
1281 int i;
1282 int *nmat;
1283
1284 if (!YAP_IsVarTerm(YAP_ARG2)) {
1285 out = YAP_ARG2;
1286 } else {
1287 out = new_float_matrix(mat[MAT_NDIMS], mat + MAT_DIMS, NULL);
1288 if (out == YAP_TermNil())
1289 return FALSE;
1290 }
1291 nmat = (int *)YAP_BlobOfTerm(out);
1292 ndata = matrix_double_data(nmat, mat[MAT_NDIMS]);
1293 for (i = 0; i < mat[MAT_SIZE]; i++) {
1294 ndata[i] = exp(data[i]);
1295 }
1296 if (YAP_IsVarTerm(YAP_ARG2)) {
1297 return YAP_Unify(YAP_ARG2, out);
1298 }
1299 }
1300 return TRUE;
1301}
1302
1303static YAP_Bool matrix_minarg(void) {
1304 int indx[MAX_DIMS], *mat;
1305 unsigned int off;
1306 YAP_Term tf;
1307
1308 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1309 if (!mat) {
1310 /* Error */
1311 return FALSE;
1312 }
1313 if (mat[MAT_TYPE] == INT_MATRIX) {
1314 YAP_Int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
1315 off = scan_min_long(mat[MAT_SIZE], data);
1316 } else {
1317 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
1318 off = scan_min_float(mat[MAT_SIZE], data);
1319 }
1320 matrix_get_index(mat, off, indx);
1321 tf = mk_int_list(mat[MAT_NDIMS], indx);
1322 return YAP_Unify(YAP_ARG2, tf);
1323}
1324
1325static YAP_Bool matrix_sum(void) {
1326 int *mat;
1327 YAP_Term tf;
1328
1329 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1330 if (!mat) {
1331 /* Error */
1332 return FALSE;
1333 }
1334 if (mat[MAT_TYPE] == INT_MATRIX) {
1335 YAP_Int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
1336 int i;
1337 YAP_Int sum = 0;
1338
1339 for (i = 0; i < mat[MAT_SIZE]; i++) {
1340 sum += data[i];
1341 }
1342 tf = YAP_MkIntTerm(sum);
1343 } else {
1344 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
1345 int i;
1346 double sum = 0.0;
1347 // function KahanSum(input)
1348 double c = 0.0; // A running compensation for lost low-order bits.
1349 for (i = 0; i < mat[MAT_SIZE]; i++) {
1350 double y = data[i] - c; // So far, so good: c is zero.
1351 double t =
1352 sum +
1353 y; // Alas, sum is big, y small, so low-order digits of y are lost.
1354 c = (t - sum) - y; // (t - sum) cancels the high-order part of y;
1355 // subtracting y recovers negative (low part of y)
1356 sum = t; // Algebraically, c should always be zero. Beware
1357 // overly-aggressive optimizing compilers!
1358 }
1359 tf = YAP_MkFloatTerm(sum);
1360 }
1361 return YAP_Unify(YAP_ARG2, tf);
1362}
1363
1364static void add_int_lines(int total, int nlines, YAP_Int *mat0,
1365 YAP_Int *matf) {
1366 int ncols = total / nlines, i;
1367 for (i = 0; i < ncols; i++) {
1368 YAP_Int sum = 0;
1369 int j;
1370
1371 for (j = i; j < total; j += ncols) {
1372 sum += mat0[j];
1373 }
1374 matf[i] = sum;
1375 }
1376}
1377
1378static void add_double_lines(int total, int nlines, double *mat0,
1379 double *matf) {
1380 int ncols = total / nlines, i;
1381 for (i = 0; i < ncols; i++) {
1382 double sum = 0;
1383 int j;
1384
1385 for (j = i; j < total; j += ncols) {
1386 sum += mat0[j];
1387 }
1388 matf[i] = sum;
1389 }
1390}
1391
1392static YAP_Bool matrix_agg_lines(void) {
1393 int *mat;
1394 YAP_Term tf;
1395 YAP_Term top = YAP_ARG2;
1396 op_type op;
1397
1398 if (!YAP_IsIntTerm(top)) {
1399 return FALSE;
1400 }
1401 op = YAP_IntOfTerm(top);
1402 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1403 if (!mat) {
1404 /* Error */
1405 return FALSE;
1406 }
1407 /* create a new array without first dimension */
1408 if (mat[MAT_TYPE] == INT_MATRIX) {
1409 YAP_Int *data, *ndata;
1410 int dims = mat[MAT_NDIMS];
1411 int *nmat;
1412
1413 tf = new_int_matrix(dims - 1, mat + (MAT_DIMS + 1), NULL);
1414 if (tf == YAP_TermNil())
1415 return FALSE;
1416 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1417 nmat = (int *)YAP_BlobOfTerm(tf);
1418 data = matrix_long_data(mat, dims);
1419 ndata = matrix_long_data(nmat, dims - 1);
1420 if (op == MAT_PLUS) {
1421 add_int_lines(mat[MAT_SIZE], mat[MAT_DIMS], data, ndata);
1422 } else
1423 return FALSE;
1424 } else {
1425 double *data, *ndata;
1426 int dims = mat[MAT_NDIMS];
1427 int *nmat;
1428
1429 tf = new_float_matrix(dims - 1, mat + (MAT_DIMS + 1), NULL);
1430 nmat = (int *)YAP_BlobOfTerm(tf);
1431 if (tf == YAP_TermNil())
1432 return FALSE;
1433 data = matrix_double_data(mat, dims);
1434 ndata = matrix_double_data(nmat, dims - 1);
1435 if (op == MAT_PLUS) {
1436 add_double_lines(mat[MAT_SIZE], mat[MAT_DIMS], data, ndata);
1437 } else
1438 return FALSE;
1439 }
1440 return YAP_Unify(YAP_ARG3, tf);
1441}
1442
1443static void add_int_cols(int total, int nlines, YAP_Int *mat0,
1444 YAP_Int *matf) {
1445 int ncols = total / nlines, i, j = 0;
1446 for (i = 0; i < nlines; i++) {
1447 YAP_Int sum = 0;
1448 int max = (i + 1) * ncols;
1449
1450 for (; j < max; j++) {
1451 sum += mat0[j];
1452 }
1453 matf[i] = sum;
1454 }
1455}
1456
1457static void add_double_cols(int total, int nlines, double *mat0, double *matf) {
1458 int ncols = total / nlines, i, j = 0;
1459 for (i = 0; i < nlines; i++) {
1460 double sum = 0;
1461 int max = (i + 1) * ncols;
1462
1463 for (; j < max; j++) {
1464 sum += mat0[j];
1465 }
1466 matf[i] = sum;
1467 }
1468}
1469
1470static YAP_Bool matrix_agg_cols(void) {
1471 int *mat;
1472 YAP_Term tf;
1473 YAP_Term top = YAP_ARG2;
1474 op_type op;
1475
1476 if (!YAP_IsIntTerm(top)) {
1477 return FALSE;
1478 }
1479 op = YAP_IntOfTerm(top);
1480 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1481 if (!mat) {
1482 /* Error */
1483 return FALSE;
1484 }
1485 /* create a new array without first dimension */
1486 if (mat[MAT_TYPE] == INT_MATRIX) {
1487 YAP_Int *data, *ndata;
1488 int dims = mat[MAT_NDIMS];
1489 int *nmat;
1490
1491 tf = new_int_matrix(1, mat + MAT_DIMS, NULL);
1492 if (tf == YAP_TermNil())
1493 return FALSE;
1494 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1495 nmat = (int *)YAP_BlobOfTerm(tf);
1496 data = matrix_long_data(mat, dims);
1497 ndata = matrix_long_data(nmat, 1);
1498 if (op == MAT_PLUS) {
1499 add_int_cols(mat[MAT_SIZE], mat[MAT_DIMS], data, ndata);
1500 } else
1501 return FALSE;
1502 } else {
1503 double *data, *ndata;
1504 int dims = mat[MAT_NDIMS];
1505 int *nmat;
1506
1507 tf = new_float_matrix(1, mat + MAT_DIMS, NULL);
1508 if (tf == YAP_TermNil())
1509 return FALSE;
1510 nmat = (int *)YAP_BlobOfTerm(tf);
1511 data = matrix_double_data(mat, dims);
1512 ndata = matrix_double_data(nmat, 1);
1513 if (op == MAT_PLUS) {
1514 add_double_cols(mat[MAT_SIZE], mat[MAT_DIMS], data, ndata);
1515 } else
1516 return FALSE;
1517 }
1518 return YAP_Unify(YAP_ARG3, tf);
1519}
1520
1521static void div_int_by_lines(int total, int nlines, YAP_Int *mat1,
1522 YAP_Int *mat2, double *ndata) {
1523 int ncols = total / nlines, i;
1524 for (i = 0; i < total; i++) {
1525 ndata[i] = ((double)mat1[i]) / mat2[i % ncols];
1526 }
1527}
1528
1529static void div_int_by_dlines(int total, int nlines, YAP_Int *mat1,
1530 double *mat2, double *ndata) {
1531 int ncols = total / nlines, i;
1532 for (i = 0; i < total; i++) {
1533 ndata[i] = mat1[i] / mat2[i % ncols];
1534 }
1535}
1536
1537static void div_float_long_by_lines(int total, int nlines, double *mat1,
1538 YAP_Int *mat2, double *ndata) {
1539 int ncols = total / nlines, i;
1540 for (i = 0; i < total; i++) {
1541 ndata[i] = mat1[i] / mat2[i % ncols];
1542 }
1543}
1544
1545static void div_float_by_lines(int total, int nlines, double *mat1,
1546 double *mat2, double *ndata) {
1547 int ncols = total / nlines, i;
1548 for (i = 0; i < total; i++) {
1549 ndata[i] = mat1[i] / mat2[i % ncols];
1550 }
1551}
1552
1553static YAP_Bool matrix_op_to_lines(void) {
1554 int *mat1, *mat2;
1555 YAP_Term top = YAP_ARG3;
1556 op_type op;
1557 YAP_Term tf;
1558
1559 if (!YAP_IsIntTerm(top)) {
1560 return FALSE;
1561 }
1562 op = YAP_IntOfTerm(top);
1563 mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
1564 if (!mat1) {
1565 /* Error */
1566 return FALSE;
1567 }
1568 mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
1569 if (!mat2) {
1570 /* Error */
1571 return FALSE;
1572 }
1573 /* create a new array without first dimension */
1574 if (mat1[MAT_TYPE] == INT_MATRIX) {
1575 YAP_Int *data1;
1576 int dims = mat1[MAT_NDIMS];
1577 int *nmat;
1578 data1 = matrix_long_data(mat1, dims);
1579
1580 if (mat2[MAT_TYPE] == INT_MATRIX) {
1581 YAP_Int *data2 = matrix_long_data(mat2, dims - 1);
1582 if (op == MAT_DIV) {
1583 double *ndata;
1584
1585 tf = new_float_matrix(dims, mat1 + MAT_DIMS, NULL);
1586 if (tf == YAP_TermNil())
1587 return FALSE;
1588 nmat = YAP_BlobOfTerm(tf);
1589 ndata = matrix_double_data(nmat, dims);
1590 div_int_by_lines(mat1[MAT_SIZE], mat1[MAT_DIMS], data1, data2, ndata);
1591 } else {
1592 return FALSE;
1593 }
1594 } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
1595 double *data2 = matrix_double_data(mat2, dims - 1);
1596 if (op == MAT_DIV) {
1597 double *ndata;
1598
1599 tf = new_float_matrix(dims, mat1 + MAT_DIMS, NULL);
1600 if (tf == YAP_TermNil())
1601 return FALSE;
1602 nmat = YAP_BlobOfTerm(tf);
1603 ndata = matrix_double_data(nmat, dims);
1604 div_int_by_dlines(mat1[MAT_SIZE], mat1[MAT_DIMS], data1, data2, ndata);
1605 } else {
1606 return FALSE;
1607 }
1608 } else {
1609 return FALSE;
1610 }
1611 } else {
1612 double *data1, *ndata;
1613 int dims = mat1[MAT_NDIMS];
1614 int *nmat;
1615
1616 data1 = matrix_double_data(mat1, dims);
1617 tf = new_float_matrix(dims, mat1 + MAT_DIMS, NULL);
1618 nmat = YAP_BlobOfTerm(tf);
1619 if (tf == YAP_TermNil())
1620 return FALSE;
1621 ndata = matrix_double_data(nmat, dims);
1622 if (mat2[MAT_TYPE] == INT_MATRIX) {
1623 YAP_Int *data2 = matrix_long_data(mat2, dims - 1);
1624 if (op == MAT_DIV) {
1625 div_float_long_by_lines(mat1[MAT_SIZE], mat1[MAT_DIMS], data1, data2,
1626 ndata);
1627 } else {
1628 return FALSE;
1629 }
1630 } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
1631 double *data2 = matrix_double_data(mat2, dims - 1);
1632 if (op == MAT_DIV) {
1633 div_float_by_lines(mat1[MAT_SIZE], mat1[MAT_DIMS], data1, data2, ndata);
1634 } else {
1635 return FALSE;
1636 }
1637 } else {
1638 return FALSE;
1639 }
1640 }
1641 return YAP_Unify(YAP_ARG4, tf);
1642}
1643
1644static void matrix_long_add_data(YAP_Int *nmat, int siz, YAP_Int mat1[],
1645 YAP_Int mat2[]) {
1646 int i;
1647
1648 for (i = 0; i < siz; i++) {
1649 nmat[i] = mat1[i] + mat2[i];
1650 }
1651}
1652
1653static void matrix_long_double_add_data(double *nmat, int siz, YAP_Int mat1[],
1654 double mat2[]) {
1655 int i;
1656
1657 for (i = 0; i < siz; i++) {
1658 nmat[i] = mat1[i] + mat2[i];
1659 }
1660}
1661
1662static void matrix_double_add_data(double *nmat, int siz, double mat1[],
1663 double mat2[]) {
1664 int i;
1665
1666 for (i = 0; i < siz; i++) {
1667 nmat[i] = mat1[i] + mat2[i];
1668 }
1669}
1670
1671static void matrix_long_sub_data(YAP_Int *nmat, int siz, YAP_Int mat1[],
1672 YAP_Int mat2[]) {
1673 int i;
1674
1675 for (i = 0; i < siz; i++) {
1676 nmat[i] = mat1[i] - mat2[i];
1677 }
1678}
1679
1680static void matrix_long_double_sub_data(double *nmat, int siz, YAP_Int mat1[],
1681 double mat2[]) {
1682 int i;
1683
1684 for (i = 0; i < siz; i++) {
1685 nmat[i] = mat1[i] - mat2[i];
1686 }
1687}
1688
1689static void matrix_long_double_rsub_data(double *nmat, int siz, double mat1[],
1690 YAP_Int mat2[]) {
1691 int i;
1692
1693 for (i = 0; i < siz; i++) {
1694 nmat[i] = mat2[i] - mat1[i];
1695 }
1696}
1697
1698static void matrix_double_sub_data(double *nmat, int siz, double mat1[],
1699 double mat2[]) {
1700 int i;
1701
1702 for (i = 0; i < siz; i++) {
1703 nmat[i] = mat1[i] - mat2[i];
1704 }
1705}
1706
1707static void matrix_long_mult_data(YAP_Int *nmat, int siz, YAP_Int mat1[],
1708 YAP_Int mat2[]) {
1709 int i;
1710
1711 for (i = 0; i < siz; i++) {
1712 nmat[i] = mat1[i] * mat2[i];
1713 }
1714}
1715
1716static void matrix_long_double_mult_data(double *nmat, int siz, YAP_Int mat1[],
1717 double mat2[]) {
1718 int i;
1719
1720 for (i = 0; i < siz; i++) {
1721 nmat[i] = mat1[i] * mat2[i];
1722 }
1723}
1724
1725static void matrix_double_mult_data(double *nmat, int siz, double mat1[],
1726 double mat2[]) {
1727 int i;
1728
1729 for (i = 0; i < siz; i++) {
1730 nmat[i] = mat1[i] * mat2[i];
1731 }
1732}
1733
1734static void matrix_long_div_data(YAP_Int *nmat, int siz, YAP_Int mat1[],
1735 YAP_Int mat2[]) {
1736 int i;
1737
1738 for (i = 0; i < siz; i++) {
1739 nmat[i] = mat1[i] / mat2[i];
1740 }
1741}
1742
1743static void matrix_long_double_div_data(double *nmat, int siz, YAP_Int mat1[],
1744 double mat2[]) {
1745 int i;
1746
1747 for (i = 0; i < siz; i++) {
1748 nmat[i] = mat1[i] / mat2[i];
1749 }
1750}
1751
1752static void matrix_long_double_div2_data(double *nmat, int siz, double mat1[],
1753 YAP_Int mat2[]) {
1754 int i;
1755
1756 for (i = 0; i < siz; i++) {
1757 nmat[i] = mat1[i] / mat2[i];
1758 }
1759}
1760
1761static void matrix_double_div_data(double *nmat, int siz, double mat1[],
1762 double mat2[]) {
1763 int i;
1764
1765 for (i = 0; i < siz; i++) {
1766 nmat[i] = mat1[i] / mat2[i];
1767 }
1768}
1769
1770static void matrix_long_zdiv_data(YAP_Int *nmat, int siz, YAP_Int mat1[],
1771 YAP_Int mat2[]) {
1772 int i;
1773
1774 for (i = 0; i < siz; i++) {
1775 if (mat1[i] == 0)
1776 nmat[i] = 0;
1777 else
1778 nmat[i] = mat1[i] / mat2[i];
1779 }
1780}
1781
1782static void matrix_long_double_zdiv_data(double *nmat, int siz, YAP_Int mat1[],
1783 double mat2[]) {
1784 int i;
1785
1786 for (i = 0; i < siz; i++) {
1787 if (mat1[i] == 0)
1788 nmat[i] = 0;
1789 else
1790 nmat[i] = mat1[i] / mat2[i];
1791 }
1792}
1793
1794static void matrix_long_double_zdiv2_data(double *nmat, int siz, double mat1[],
1795 YAP_Int mat2[]) {
1796 int i;
1797
1798 for (i = 0; i < siz; i++) {
1799 if (mat1[i] == 0.0)
1800 nmat[i] = 0;
1801 else
1802 nmat[i] = mat1[i] / mat2[i];
1803 }
1804}
1805
1806static void matrix_double_zdiv_data(double *nmat, int siz, double mat1[],
1807 double mat2[]) {
1808 int i;
1809
1810 for (i = 0; i < siz; i++) {
1811 if (mat1[i] == 0.0) {
1812 nmat[i] = 0.0;
1813 } else {
1814 nmat[i] = mat1[i] / mat2[i];
1815 }
1816 }
1817}
1818
1819static YAP_Bool matrix_op(void) {
1820 int *mat1, *mat2;
1821 YAP_Term top = YAP_ARG3;
1822 op_type op;
1823 YAP_Term tf = YAP_ARG4;
1824 int create = TRUE;
1825
1826 if (!YAP_IsIntTerm(top)) {
1827 return FALSE;
1828 }
1829 op = YAP_IntOfTerm(top);
1830 mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
1831 if (!mat1) {
1832 /* Error */
1833 return FALSE;
1834 }
1835 mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
1836 if (!mat2) {
1837 /* Error */
1838 return FALSE;
1839 }
1840 if (tf == YAP_ARG1 || tf == YAP_ARG2) {
1841 create = FALSE;
1842 }
1843 if (mat1[MAT_TYPE] == INT_MATRIX) {
1844 YAP_Int *data1;
1845 int dims = mat1[MAT_NDIMS];
1846 int *nmat;
1847 data1 = matrix_long_data(mat1, dims);
1848
1849 if (mat2[MAT_TYPE] == INT_MATRIX) {
1850 YAP_Int *data2;
1851 YAP_Int *ndata;
1852
1853 if (create)
1854 tf = new_int_matrix(dims, mat1 + MAT_DIMS, NULL);
1855 if (tf == YAP_TermNil()) {
1856 return FALSE;
1857 } else {
1858 /* there may have been an overflow */
1859 mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
1860 data1 = matrix_long_data(mat1, dims);
1861 mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
1862 data2 = matrix_long_data(mat2, dims);
1863 }
1864 nmat = YAP_BlobOfTerm(tf);
1865 ndata = matrix_long_data(nmat, dims);
1866 switch (op) {
1867 case MAT_PLUS:
1868 matrix_long_add_data(ndata, mat1[MAT_SIZE], data1, data2);
1869 break;
1870 case MAT_SUB:
1871 matrix_long_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
1872 break;
1873 case MAT_TIMES:
1874 matrix_long_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
1875 break;
1876 case MAT_DIV:
1877 matrix_long_div_data(ndata, mat1[MAT_SIZE], data1, data2);
1878 break;
1879 case MAT_ZDIV:
1880 matrix_long_zdiv_data(ndata, mat1[MAT_SIZE], data1, data2);
1881 break;
1882 default:
1883 return FALSE;
1884 }
1885 } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
1886 double *data2;
1887 double *ndata;
1888
1889 if (create)
1890 tf = new_float_matrix(dims, mat1 + MAT_DIMS, NULL);
1891 if (tf == YAP_TermNil()) {
1892 return FALSE;
1893 } else {
1894 /* there may have been an overflow */
1895 mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
1896 data1 = matrix_long_data(mat1, dims);
1897 mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
1898 data2 = matrix_double_data(mat2, dims);
1899 }
1900 nmat = YAP_BlobOfTerm(tf);
1901 ndata = matrix_double_data(nmat, dims);
1902 switch (op) {
1903 case MAT_PLUS:
1904 matrix_long_double_add_data(ndata, mat1[MAT_SIZE], data1, data2);
1905 break;
1906 case MAT_SUB:
1907 matrix_long_double_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
1908 break;
1909 case MAT_TIMES:
1910 matrix_long_double_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
1911 break;
1912 case MAT_DIV:
1913 matrix_long_double_div_data(ndata, mat1[MAT_SIZE], data1, data2);
1914 break;
1915 case MAT_ZDIV:
1916 matrix_long_double_zdiv_data(ndata, mat1[MAT_SIZE], data1, data2);
1917 break;
1918 default:
1919 return FALSE;
1920 }
1921 } else {
1922 return FALSE;
1923 }
1924 } else {
1925 double *data1;
1926 int dims = mat1[MAT_NDIMS];
1927 int *nmat;
1928 data1 = matrix_double_data(mat1, dims);
1929
1930 if (mat2[MAT_TYPE] == INT_MATRIX) {
1931 YAP_Int *data2;
1932 double *ndata;
1933
1934 if (create)
1935 tf = new_float_matrix(dims, mat1 + MAT_DIMS, NULL);
1936 if (tf == YAP_TermNil()) {
1937 return FALSE;
1938 } else {
1939 /* there may have been an overflow */
1940 mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
1941 data1 = matrix_double_data(mat1, dims);
1942 mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
1943 data2 = matrix_long_data(mat2, dims);
1944 }
1945 nmat = YAP_BlobOfTerm(tf);
1946 ndata = matrix_double_data(nmat, dims);
1947 switch (op) {
1948 case MAT_PLUS:
1949 matrix_long_double_add_data(ndata, mat1[MAT_SIZE], data2, data1);
1950 break;
1951 case MAT_SUB:
1952 matrix_long_double_rsub_data(ndata, mat1[MAT_SIZE], data1, data2);
1953 break;
1954 case MAT_TIMES:
1955 matrix_long_double_mult_data(ndata, mat1[MAT_SIZE], data2, data1);
1956 break;
1957 case MAT_DIV:
1958 matrix_long_double_div2_data(ndata, mat1[MAT_SIZE], data1, data2);
1959 break;
1960 case MAT_ZDIV:
1961 matrix_long_double_zdiv2_data(ndata, mat1[MAT_SIZE], data1, data2);
1962 break;
1963 default:
1964 return FALSE;
1965 }
1966 } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
1967 double *data2;
1968 double *ndata;
1969
1970 if (create)
1971 tf = new_float_matrix(dims, mat1 + MAT_DIMS, NULL);
1972 if (tf == YAP_TermNil()) {
1973 return FALSE;
1974 } else {
1975 /* there may have been an overflow */
1976 mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
1977 data1 = matrix_double_data(mat1, dims);
1978 mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
1979 data2 = matrix_double_data(mat2, dims);
1980 }
1981 nmat = YAP_BlobOfTerm(tf);
1982 ndata = matrix_double_data(nmat, dims);
1983 switch (op) {
1984 case MAT_PLUS:
1985 matrix_double_add_data(ndata, mat1[MAT_SIZE], data1, data2);
1986 break;
1987 case MAT_SUB:
1988 matrix_double_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
1989 break;
1990 case MAT_TIMES:
1991 matrix_double_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
1992 break;
1993 case MAT_DIV:
1994 matrix_double_div_data(ndata, mat1[MAT_SIZE], data1, data2);
1995 break;
1996 case MAT_ZDIV:
1997 matrix_double_zdiv_data(ndata, mat1[MAT_SIZE], data1, data2);
1998 break;
1999 default:
2000 return FALSE;
2001 }
2002 } else {
2003 return FALSE;
2004 }
2005 }
2006 return YAP_Unify(YAP_ARG4, tf);
2007}
2008
2009static void add_int_by_cols(int total, int nlines, YAP_Int *mat1,
2010 YAP_Int *mat2, YAP_Int *ndata) {
2011 int i, ncols = total / nlines;
2012 for (i = 0; i < total; i++) {
2013 ndata[i] = mat1[i] + mat2[i / ncols];
2014 }
2015}
2016
2017static void add_int_by_dcols(int total, int nlines, YAP_Int *mat1,
2018 double *mat2, double *ndata) {
2019 int i, ncols = total / nlines;
2020 for (i = 0; i < total; i++) {
2021 ndata[i] = mat1[i] + mat2[i / ncols];
2022 }
2023}
2024
2025static void add_double_by_cols(int total, int nlines, double *mat1,
2026 double *mat2, double *ndata) {
2027 int i;
2028 int ncols = total / nlines;
2029
2030 for (i = 0; i < total; i++) {
2031 ndata[i] = mat1[i] + mat2[i / ncols];
2032 }
2033}
2034
2035static YAP_Bool matrix_op_to_cols(void) {
2036 int *mat1, *mat2;
2037 YAP_Term top = YAP_ARG3;
2038 op_type op;
2039 YAP_Term tf;
2040
2041 if (!YAP_IsIntTerm(top)) {
2042 return FALSE;
2043 }
2044 op = YAP_IntOfTerm(top);
2045 mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
2046 if (!mat1) {
2047 /* Error */
2048 return FALSE;
2049 }
2050 mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
2051 if (!mat2) {
2052 /* Error */
2053 return FALSE;
2054 }
2055 if (mat1[MAT_TYPE] == INT_MATRIX) {
2056 YAP_Int *data1;
2057 int dims = mat1[MAT_NDIMS];
2058 int *nmat;
2059 data1 = matrix_long_data(mat1, dims);
2060
2061 if (mat2[MAT_TYPE] == INT_MATRIX) {
2062 YAP_Int *data2 = matrix_long_data(mat2, 1);
2063 if (op == MAT_PLUS) {
2064 YAP_Int *ndata;
2065
2066 tf = new_int_matrix(dims, mat1 + MAT_DIMS, NULL);
2067 if (tf == YAP_TermNil())
2068 return FALSE;
2069 nmat = YAP_BlobOfTerm(tf);
2070 ndata = matrix_long_data(nmat, dims);
2071 add_int_by_cols(mat1[MAT_SIZE], mat1[MAT_DIMS], data1, data2, ndata);
2072 } else {
2073 return FALSE;
2074 }
2075 } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
2076 double *data2 = matrix_double_data(mat2, 1);
2077 if (op == MAT_PLUS) {
2078 double *ndata;
2079
2080 tf = new_float_matrix(dims, mat1 + MAT_DIMS, NULL);
2081 if (tf == YAP_TermNil())
2082 return FALSE;
2083 nmat = YAP_BlobOfTerm(tf);
2084 ndata = matrix_double_data(nmat, dims);
2085 add_int_by_dcols(mat1[MAT_SIZE], mat1[MAT_DIMS], data1, data2, ndata);
2086 } else {
2087 return FALSE;
2088 }
2089 } else {
2090 return FALSE;
2091 }
2092 } else {
2093 double *data1, *data2, *ndata;
2094 int dims = mat1[MAT_NDIMS];
2095 int *nmat;
2096
2097 if (mat2[MAT_TYPE] != FLOAT_MATRIX)
2098 return FALSE;
2099 tf = new_float_matrix(dims, mat1 + MAT_DIMS, NULL);
2100 if (tf == YAP_TermNil())
2101 return FALSE;
2102 nmat = YAP_BlobOfTerm(tf);
2103 data1 = matrix_double_data(mat1, dims);
2104 data2 = matrix_double_data(mat2, 1);
2105 ndata = matrix_double_data(nmat, dims);
2106 if (op == MAT_PLUS) {
2107 add_double_by_cols(mat1[MAT_SIZE], mat1[MAT_DIMS], data1, data2, ndata);
2108 } else
2109 return FALSE;
2110 }
2111 return YAP_Unify(YAP_ARG4, tf);
2112}
2113
2114static YAP_Bool matrix_op_to_all(void) {
2115 int *mat;
2116 YAP_Term tf = 0;
2117 YAP_Term top = YAP_ARG2;
2118 op_type op;
2119 int create = FALSE;
2120
2121 if (!YAP_IsIntTerm(top)) {
2122 return FALSE;
2123 }
2124 op = YAP_IntOfTerm(top);
2125 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2126 if (!mat) {
2127 /* Error */
2128 return FALSE;
2129 }
2130 if (YAP_IsVarTerm(YAP_ARG4)) {
2131 create = TRUE;
2132 }
2133 /* create a new array with same dimensions */
2134 if (mat[MAT_TYPE] == INT_MATRIX) {
2135 YAP_Int *data;
2136 int dims = mat[MAT_NDIMS];
2137 int *nmat;
2138 YAP_Term tnum = YAP_ARG3;
2139
2140 if (YAP_IsIntTerm(tnum)) {
2141 YAP_Int num;
2142 YAP_Int *ndata;
2143
2144 num = YAP_IntOfTerm(tnum);
2145 data = matrix_long_data(mat, dims);
2146 if (create) {
2147 tf = new_int_matrix(dims, mat + (MAT_DIMS), NULL);
2148 if (tf == YAP_TermNil())
2149 return FALSE;
2150 nmat = (int *)YAP_BlobOfTerm(tf);
2151 ndata = matrix_long_data(nmat, dims);
2152 } else {
2153 nmat = mat;
2154 ndata = data;
2155 }
2156 if (op == MAT_PLUS) {
2157 int i;
2158
2159 for (i = 0; i < mat[MAT_SIZE]; i++) {
2160 ndata[i] = data[i] + num;
2161 }
2162 } else if (op == MAT_TIMES) {
2163 int i;
2164
2165 for (i = 0; i < mat[MAT_SIZE]; i++) {
2166 ndata[i] = data[i] * num;
2167 }
2168 } else {
2169 return FALSE;
2170 }
2171 } else if (YAP_IsFloatTerm(tnum)) {
2172 double num;
2173 double *ndata;
2174
2175 num = YAP_FloatOfTerm(tnum);
2176 if (create) {
2177 tf = new_float_matrix(dims, mat + (MAT_DIMS), NULL);
2178 if (tf == YAP_TermNil())
2179 return FALSE;
2180 nmat = (int *)YAP_BlobOfTerm(tf);
2181 ndata = matrix_double_data(nmat, dims);
2182 } else {
2183 return FALSE;
2184 }
2185 data = matrix_long_data(mat, dims);
2186 if (op == MAT_PLUS) {
2187 int i;
2188
2189 for (i = 0; i < mat[MAT_SIZE]; i++) {
2190 ndata[i] = data[i] + num;
2191 }
2192 } else if (op == MAT_SUB) {
2193 int i;
2194
2195 for (i = 0; i < mat[MAT_SIZE]; i++) {
2196 ndata[i] = num - data[i];
2197 }
2198 } else if (op == MAT_TIMES) {
2199 int i;
2200
2201 for (i = 0; i < mat[MAT_SIZE]; i++) {
2202 ndata[i] = data[i] * num;
2203 }
2204 } else if (op == MAT_DIV) {
2205 int i;
2206
2207 for (i = 0; i < mat[MAT_SIZE]; i++) {
2208 ndata[i] = data[i] / num;
2209 }
2210 }
2211 } else {
2212 return FALSE;
2213 }
2214 } else {
2215 double *data, *ndata;
2216 int dims = mat[MAT_NDIMS];
2217 int *nmat;
2218 YAP_Term tnum = YAP_ARG3;
2219 double num;
2220
2221 if (YAP_IsFloatTerm(tnum)) {
2222 num = YAP_FloatOfTerm(tnum);
2223 } else if (!YAP_IntOfTerm(tnum)) {
2224 return FALSE;
2225 } else {
2226 if (!create)
2227 return FALSE;
2228 num = (double)YAP_IntOfTerm(tnum);
2229 }
2230 data = matrix_double_data(mat, dims);
2231 if (create) {
2232 tf = new_float_matrix(dims, mat + (MAT_DIMS), NULL);
2233 if (tf == YAP_TermNil())
2234 return FALSE;
2235 nmat = (int *)YAP_BlobOfTerm(tf);
2236 ndata = matrix_double_data(nmat, dims);
2237 } else {
2238 nmat = mat;
2239 ndata = data;
2240 }
2241 switch (op) {
2242 case MAT_PLUS: {
2243 int i;
2244
2245 for (i = 0; i < mat[MAT_SIZE]; i++) {
2246 ndata[i] = data[i] + num;
2247 }
2248 } break;
2249 case MAT_SUB: {
2250 int i;
2251
2252 for (i = 0; i < mat[MAT_SIZE]; i++) {
2253 ndata[i] = num - data[i];
2254 }
2255 } break;
2256 case MAT_TIMES: {
2257 int i;
2258
2259 for (i = 0; i < mat[MAT_SIZE]; i++) {
2260 ndata[i] = data[i] * num;
2261 }
2262 } break;
2263 case MAT_DIV: {
2264 int i;
2265
2266 for (i = 0; i < mat[MAT_SIZE]; i++) {
2267 ndata[i] = data[i] / num;
2268 }
2269 } break;
2270 default:
2271 return FALSE;
2272 }
2273 }
2274 if (create)
2275 return YAP_Unify(YAP_ARG4, tf);
2276 return YAP_Unify(YAP_ARG4, YAP_ARG1);
2277}
2278
2279/* given a matrix M and a set of dims, build a new reordered matrix to follow
2280 the new order
2281*/
2282static YAP_Bool matrix_transpose(void) {
2283 int ndims, i, *dims, *dimsn;
2284 int conv[MAX_DIMS], indx[MAX_DIMS], nindx[MAX_DIMS];
2285 YAP_Term tconv, tf;
2286 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2287 if (!mat) {
2288 /* Error */
2289 return FALSE;
2290 }
2291 ndims = mat[MAT_NDIMS];
2292 if (mat[MAT_TYPE] == INT_MATRIX) {
2293 /* create a new matrix with the same size */
2294 tf = new_int_matrix(ndims, mat + MAT_DIMS, NULL);
2295 if (tf == YAP_TermNil())
2296 return FALSE;
2297 } else {
2298 /* create a new matrix with the same size */
2299 tf = new_float_matrix(ndims, mat + MAT_DIMS, NULL);
2300 if (tf == YAP_TermNil())
2301 return FALSE;
2302 }
2303 /* just in case there was an overflow */
2304 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2305 nmat = (int *)YAP_BlobOfTerm(tf);
2306 dims = mat + MAT_DIMS;
2307 dimsn = nmat + MAT_DIMS;
2308 /* we now have our target matrix, let us grab our conversion matrix */
2309 tconv = YAP_ARG2;
2310 for (i = 0; i < ndims; i++) {
2311 YAP_Term th;
2312 YAP_Int j;
2313
2314 if (!YAP_IsPairTerm(tconv))
2315 return FALSE;
2316 th = YAP_HeadOfTerm(tconv);
2317 if (!YAP_IsIntTerm(th))
2318 return FALSE;
2319 conv[i] = j = YAP_IntOfTerm(th);
2320 dimsn[i] = dims[j];
2321 tconv = YAP_TailOfTerm(tconv);
2322 }
2323 /*
2324 we now got all the dimensions set up, so what we need to do
2325 next is to copy the elements to the new matrix.
2326 */
2327 if (mat[MAT_TYPE] == INT_MATRIX) {
2328 YAP_Int *data = matrix_long_data(mat, ndims);
2329 /* create a new matrix with the same size */
2330 for (i = 0; i < mat[MAT_SIZE]; i++) {
2331 YAP_Int x = data[i];
2332 int j;
2333 /*
2334 not very efficient, we could try to take advantage of the fact
2335 that we usually only change an index at a time
2336 */
2337 matrix_get_index(mat, i, indx);
2338 for (j = 0; j < ndims; j++) {
2339 nindx[j] = indx[conv[j]];
2340 }
2341 matrix_long_set(nmat, nindx, x);
2342 }
2343 } else {
2344 double *data = matrix_double_data(mat, ndims);
2345 /* create a new matrix with the same size */
2346 for (i = 0; i < mat[MAT_SIZE]; i++) {
2347 double x = data[i];
2348 long j;
2349
2350 matrix_get_index(mat, i, indx);
2351 for (j = 0; j < ndims; j++)
2352 nindx[j] = indx[conv[j]];
2353 matrix_float_set(nmat, nindx, x);
2354 }
2355 }
2356 return YAP_Unify(YAP_ARG3, tf);
2357}
2358
2359/* given a matrix M and a set of dims, fold one of the dimensions of the
2360 matrix on one of the elements
2361*/
2362static YAP_Bool matrix_select(void) {
2363 int ndims, i, j, newdims, prdim, leftarg, *dims, indx[MAX_DIMS];
2364 int nindx[MAX_DIMS];
2365 YAP_Term tpdim, tdimarg, tf;
2366 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2367 if (!mat) {
2368 /* Error */
2369 return FALSE;
2370 }
2371 /* we now have our target matrix, let us grab our conversion arguments */
2372 tpdim = YAP_ARG2;
2373 ndims = mat[MAT_NDIMS];
2374 dims = mat + MAT_DIMS;
2375 if (!YAP_IsIntTerm(tpdim)) {
2376 return FALSE;
2377 }
2378 prdim = YAP_IntOfTerm(tpdim);
2379 tdimarg = YAP_ARG3;
2380 if (!YAP_IsIntTerm(tdimarg)) {
2381 return FALSE;
2382 }
2383 leftarg = YAP_IntOfTerm(tdimarg);
2384 for (i = 0, j = 0; i < ndims; i++) {
2385 if (i != prdim) {
2386 nindx[j] = dims[i];
2387 j++;
2388 }
2389 }
2390 newdims = ndims - 1;
2391 if (mat[MAT_TYPE] == INT_MATRIX) {
2392 YAP_Int *data, *ndata;
2393
2394 /* create a new matrix with the same size */
2395 tf = new_int_matrix(newdims, nindx, NULL);
2396 if (tf == YAP_TermNil())
2397 return FALSE;
2398 /* in case the matrix moved */
2399 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2400 nmat = (int *)YAP_BlobOfTerm(tf);
2401 data = matrix_long_data(mat, ndims);
2402 ndata = matrix_long_data(nmat, newdims);
2403 /* create a new matrix with smaller size */
2404 for (i = 0; i < nmat[MAT_SIZE]; i++) {
2405 int j, k;
2406 /*
2407 not very efficient, we could try to take advantage of the fact
2408 that we usually only change an index at a time
2409 */
2410 matrix_get_index(nmat, i, indx);
2411 for (j = 0, k = 0; j < newdims; j++, k++) {
2412 if (j == prdim) {
2413 nindx[k] = leftarg;
2414 k++;
2415 }
2416 nindx[k] = indx[j];
2417 }
2418 if (k == prdim) {
2419 nindx[k] = leftarg;
2420 }
2421 ndata[i] = data[matrix_get_offset(mat, nindx)];
2422 }
2423 } else {
2424 double *data, *ndata;
2425
2426 /* create a new matrix with the same size */
2427 tf = new_float_matrix(newdims, nindx, NULL);
2428 if (tf == YAP_TermNil())
2429 return FALSE;
2430 /* in case the matrix moved */
2431 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2432 nmat = (int *)YAP_BlobOfTerm(tf);
2433 data = matrix_double_data(mat, ndims);
2434 ndata = matrix_double_data(nmat, newdims);
2435 /* create a new matrix with the same size */
2436 for (i = 0; i < nmat[MAT_SIZE]; i++) {
2437 int j, k;
2438 /*
2439 not very efficient, we could try to take advantage of the fact
2440 that we usually only change an index at a time
2441 */
2442 matrix_get_index(nmat, i, indx);
2443 for (j = 0, k = 0; j < newdims; j++, k++) {
2444 if (j == prdim) {
2445 nindx[k] = leftarg;
2446 k++;
2447 }
2448 nindx[k] = indx[j];
2449 }
2450 if (k == prdim) {
2451 nindx[k] = leftarg;
2452 }
2453 ndata[i] = data[matrix_get_offset(mat, nindx)];
2454 }
2455 }
2456 return YAP_Unify(YAP_ARG4, tf);
2457}
2458
2459/* given a matrix M and a set of N-1 dims, get the first dimension
2460 */
2461static YAP_Bool matrix_column(void) {
2462 int size, i, ndims, newdims[1];
2463 int indx[MAX_DIMS];
2464 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2465 YAP_Term tconv, tf;
2466
2467 if (!mat) {
2468 /* Error */
2469 return FALSE;
2470 }
2471 ndims = mat[MAT_NDIMS];
2472 /* we now have our target matrix, let us grab our conversion arguments */
2473 tconv = YAP_ARG2;
2474 for (i = 1; i < ndims; i++) {
2475 YAP_Term th;
2476
2477 if (!YAP_IsPairTerm(tconv))
2478 return FALSE;
2479 th = YAP_HeadOfTerm(tconv);
2480 if (!YAP_IsIntTerm(th))
2481 return FALSE;
2482 indx[i] = YAP_IntOfTerm(th);
2483 tconv = YAP_TailOfTerm(tconv);
2484 }
2485 if (tconv != YAP_TermNil())
2486 return FALSE;
2487 newdims[0] = size = mat[MAT_DIMS];
2488 if (mat[MAT_TYPE] == INT_MATRIX) {
2489 YAP_Int *data, *ndata;
2490
2491 /* create a new matrix with the same size */
2492 tf = new_int_matrix(1, newdims, NULL);
2493 if (tf == YAP_TermNil())
2494 return FALSE;
2495 /* in case the matrix moved */
2496 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2497 nmat = (int *)YAP_BlobOfTerm(tf);
2498 data = matrix_long_data(mat, ndims);
2499 ndata = matrix_long_data(nmat, 1);
2500 /* create a new matrix with smaller size */
2501 for (i = 0; i < size; i++) {
2502 indx[0] = i;
2503 ndata[i] = data[matrix_get_offset(mat, indx)];
2504 }
2505 } else {
2506 double *data, *ndata;
2507
2508 /* create a new matrix with the same size */
2509 tf = new_float_matrix(1, newdims, NULL);
2510 if (tf == YAP_TermNil())
2511 return FALSE;
2512 /* in case the matrix moved */
2513 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2514 nmat = (int *)YAP_BlobOfTerm(tf);
2515 data = matrix_double_data(mat, ndims);
2516 ndata = matrix_double_data(nmat, 1);
2517 /* create a new matrix with smaller size */
2518 for (i = 0; i < size; i++) {
2519 indx[0] = i;
2520 ndata[i] = data[matrix_get_offset(mat, indx)];
2521 }
2522 }
2523 return YAP_Unify(YAP_ARG3, tf);
2524}
2525
2526/* given a matrix M and a set of dims, sum out one of the dimensions
2527 */
2528static YAP_Bool matrix_sum_out(void) {
2529 int ndims, i, j, newdims, prdim;
2530 int indx[MAX_DIMS], nindx[MAX_DIMS];
2531 YAP_Term tpdim, tf;
2532 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2533 if (!mat) {
2534 /* Error */
2535 return FALSE;
2536 }
2537 /* we now have our target matrix, let us grab our conversion arguments */
2538 tpdim = YAP_ARG2;
2539 ndims = mat[MAT_NDIMS];
2540 if (!YAP_IsIntTerm(tpdim)) {
2541 return FALSE;
2542 }
2543 prdim = YAP_IntOfTerm(tpdim);
2544 newdims = ndims - 1;
2545 for (i = 0, j = 0; i < ndims; i++) {
2546 if (i != prdim) {
2547 nindx[j] = (mat + MAT_DIMS)[i];
2548 j++;
2549 }
2550 }
2551 if (mat[MAT_TYPE] == INT_MATRIX) {
2552 YAP_Int *data, *ndata;
2553
2554 /* create a new matrix with the same size */
2555 tf = new_int_matrix(newdims, nindx, NULL);
2556 if (tf == YAP_TermNil())
2557 return FALSE;
2558 /* in case the matrix moved */
2559 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2560 nmat = (int *)YAP_BlobOfTerm(tf);
2561 data = matrix_long_data(mat, ndims);
2562 ndata = matrix_long_data(nmat, newdims);
2563 /* create a new matrix with smaller size */
2564 for (i = 0; i < nmat[MAT_SIZE]; i++)
2565 ndata[i] = 0;
2566 for (i = 0; i < mat[MAT_SIZE]; i++) {
2567 int j, k;
2568 /*
2569 not very efficient, we could try to take advantage of the fact
2570 that we usually only change an index at a time
2571 */
2572 matrix_get_index(mat, i, indx);
2573 for (j = 0, k = 0; j < ndims; j++) {
2574 if (j != prdim) {
2575 nindx[k++] = indx[j];
2576 }
2577 }
2578 ndata[matrix_get_offset(nmat, nindx)] += data[i];
2579 }
2580 } else {
2581 double *data, *ndata;
2582
2583 /* create a new matrix with the same size */
2584 tf = new_float_matrix(newdims, nindx, NULL);
2585 if (tf == YAP_TermNil())
2586 return FALSE;
2587 /* in case the matrix moved */
2588 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2589 nmat = (int *)YAP_BlobOfTerm(tf);
2590 data = matrix_double_data(mat, ndims);
2591 ndata = matrix_double_data(nmat, newdims);
2592 /* create a new matrix with smaller size */
2593 for (i = 0; i < nmat[MAT_SIZE]; i++)
2594 ndata[i] = 0.0;
2595 for (i = 0; i < mat[MAT_SIZE]; i++) {
2596 int j, k;
2597 /*
2598 not very efficient, we could try to take advantage of the fact
2599 that we usually only change an index at a time
2600 */
2601 matrix_get_index(mat, i, indx);
2602 for (j = 0, k = 0; j < ndims; j++) {
2603 if (j != prdim) {
2604 nindx[k++] = indx[j];
2605 }
2606 }
2607 ndata[matrix_get_offset(nmat, nindx)] += data[i];
2608 }
2609 }
2610 return YAP_Unify(YAP_ARG3, tf);
2611}
2612
2613/* given a matrix M and a set of dims, sum out one of the dimensions
2614 */
2615static YAP_Bool matrix_sum_out_several(void) {
2616 int ndims, i, *dims, newdims;
2617 int indx[MAX_DIMS], nindx[MAX_DIMS], conv[MAX_DIMS];
2618 YAP_Term tf, tconv;
2619 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2620 if (!mat) {
2621 /* Error */
2622 return FALSE;
2623 }
2624 ndims = mat[MAT_NDIMS];
2625 dims = mat + MAT_DIMS;
2626 /* we now have our target matrix, let us grab our conversion arguments */
2627 tconv = YAP_ARG2;
2628 for (i = 0, newdims = 0; i < ndims; i++) {
2629 YAP_Term th;
2630
2631 if (!YAP_IsPairTerm(tconv))
2632 return FALSE;
2633 th = YAP_HeadOfTerm(tconv);
2634 if (!YAP_IsIntTerm(th))
2635 return FALSE;
2636 conv[i] = YAP_IntOfTerm(th);
2637 if (!conv[i]) {
2638 nindx[newdims++] = dims[i];
2639 }
2640 tconv = YAP_TailOfTerm(tconv);
2641 }
2642 if (mat[MAT_TYPE] == INT_MATRIX) {
2643 YAP_Int *data, *ndata;
2644
2645 /* create a new matrix with the same size */
2646 tf = new_int_matrix(newdims, nindx, NULL);
2647 if (tf == YAP_TermNil())
2648 return FALSE;
2649 /* in case the matrix moved */
2650 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2651 nmat = (int *)YAP_BlobOfTerm(tf);
2652 data = matrix_long_data(mat, ndims);
2653 ndata = matrix_long_data(nmat, newdims);
2654 /* create a new matrix with smaller size */
2655 for (i = 0; i < nmat[MAT_SIZE]; i++)
2656 ndata[i] = 0;
2657 for (i = 0; i < mat[MAT_SIZE]; i++) {
2658 int j, k;
2659 /*
2660 not very efficient, we could try to take advantage of the fact
2661 that we usually only change an index at a time
2662 */
2663 matrix_get_index(mat, i, indx);
2664 for (j = 0, k = 0; j < ndims; j++) {
2665 if (!conv[j]) {
2666 nindx[k++] = indx[j];
2667 }
2668 }
2669 ndata[matrix_get_offset(nmat, nindx)] =
2670 log(exp(ndata[matrix_get_offset(nmat, nindx)]) + exp(data[i]));
2671 }
2672 } else {
2673 double *data, *ndata;
2674
2675 /* create a new matrix with the same size */
2676 tf = new_float_matrix(newdims, nindx, NULL);
2677 if (tf == YAP_TermNil())
2678 return FALSE;
2679 /* in case the matrix moved */
2680 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2681 nmat = (int *)YAP_BlobOfTerm(tf);
2682 data = matrix_double_data(mat, ndims);
2683 ndata = matrix_double_data(nmat, newdims);
2684 /* create a new matrix with smaller size */
2685 for (i = 0; i < nmat[MAT_SIZE]; i++)
2686 ndata[i] = 0.0;
2687 for (i = 0; i < mat[MAT_SIZE]; i++) {
2688 int j, k;
2689 /*
2690 not very efficient, we could try to take advantage of the fact
2691 that we usually only change an index at a time
2692 */
2693 matrix_get_index(mat, i, indx);
2694 for (j = 0, k = 0; j < ndims; j++) {
2695 if (!conv[j]) {
2696 nindx[k++] = indx[j];
2697 }
2698 }
2699 ndata[matrix_get_offset(nmat, nindx)] =
2700 log(exp(ndata[matrix_get_offset(nmat, nindx)]) + exp(data[i]));
2701 }
2702 }
2703 return YAP_Unify(YAP_ARG3, tf);
2704}
2705
2706/* given a matrix M and a set of dims, sum out one of the dimensions
2707 */
2708static YAP_Bool matrix_sum_out_logs(void) {
2709 int ndims, i, j, *dims, newdims, prdim;
2710 int nindx[MAX_DIMS];
2711 YAP_Term tpdim, tf;
2712 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2713 if (!mat) {
2714 /* Error */
2715 return FALSE;
2716 }
2717 /* we now have our target matrix, let us grab our conversion arguments */
2718 tpdim = YAP_ARG2;
2719 ndims = mat[MAT_NDIMS];
2720 dims = mat + MAT_DIMS;
2721 if (!YAP_IsIntTerm(tpdim)) {
2722 return FALSE;
2723 }
2724 prdim = YAP_IntOfTerm(tpdim);
2725 newdims = ndims - 1;
2726 for (i = 0, j = 0; i < ndims; i++) {
2727 if (i != prdim) {
2728 nindx[j] = dims[i];
2729 j++;
2730 }
2731 }
2732 if (mat[MAT_TYPE] == INT_MATRIX) {
2733 YAP_Int *data, *ndata;
2734 int d = 1, j = 0, dd = 1;
2735
2736 /* create a new matrix with the same size */
2737 tf = new_int_matrix(newdims, nindx, NULL);
2738 if (tf == YAP_TermNil())
2739 return FALSE;
2740 /* in case the matrix moved */
2741 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2742 nmat = (int *)YAP_BlobOfTerm(tf);
2743 data = matrix_long_data(mat, ndims);
2744 ndata = matrix_long_data(nmat, newdims);
2745 while (j < prdim) {
2746 d = d * dims[j];
2747 j++;
2748 }
2749 dd = d * dims[prdim];
2750 for (i = 0; i < nmat[MAT_SIZE]; i++) {
2751 int j = i % d + (i / dd) * d;
2752 ndata[j] = exp(data[i]);
2753 }
2754 for (; i < mat[MAT_SIZE]; i++) {
2755 int j = i % d + (i / dd) * d;
2756 ndata[j] += exp(data[i]);
2757 }
2758 for (i = 0; i < nmat[MAT_SIZE]; i++) {
2759 ndata[i] = log(ndata[i]);
2760 }
2761 } else {
2762 double *data, *ndata;
2763 int d = 1, j = 0, dd = 1;
2764
2765 /* create a new matrix with the same size */
2766 tf = new_float_matrix(newdims, nindx, NULL);
2767 if (tf == YAP_TermNil())
2768 return FALSE;
2769 /* in case the matrix moved */
2770 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2771 nmat = (int *)YAP_BlobOfTerm(tf);
2772 data = matrix_double_data(mat, ndims);
2773 ndata = matrix_double_data(nmat, newdims);
2774
2775 j = ndims - 1;
2776 while (j > prdim) {
2777 d = d * dims[j];
2778 j--;
2779 }
2780 dd = d * dims[prdim];
2781 memset(ndata, 0, sizeof(double) * nmat[MAT_SIZE]);
2782 for (i = 0; i < mat[MAT_SIZE]; i++) {
2783 YAP_Int k = i % d + (i / dd) * d;
2784 ndata[k] += exp(data[i]);
2785 }
2786 for (i = 0; i < nmat[MAT_SIZE]; i++) {
2787 ndata[i] = log(ndata[i]);
2788 }
2789 }
2790 return YAP_Unify(YAP_ARG3, tf);
2791}
2792
2793/* given a matrix M and a set of dims, sum out one of the dimensions
2794 */
2795static YAP_Bool matrix_sum_out_logs_several(void) {
2796 int ndims, i, *dims, newdims;
2797 int indx[MAX_DIMS], nindx[MAX_DIMS], conv[MAX_DIMS];
2798 YAP_Term tf, tconv;
2799 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2800
2801 if (!mat) {
2802 /* Error */
2803 return FALSE;
2804 }
2805 ndims = mat[MAT_NDIMS];
2806 dims = mat + MAT_DIMS;
2807 /* we now have our target matrix, let us grab our conversion arguments */
2808 tconv = YAP_ARG2;
2809 for (i = 0, newdims = 0; i < ndims; i++) {
2810 YAP_Term th;
2811
2812 if (!YAP_IsPairTerm(tconv))
2813 return FALSE;
2814 th = YAP_HeadOfTerm(tconv);
2815 if (!YAP_IsIntTerm(th))
2816 return FALSE;
2817 conv[i] = YAP_IntOfTerm(th);
2818 if (!conv[i]) {
2819 nindx[newdims++] = dims[i];
2820 }
2821 tconv = YAP_TailOfTerm(tconv);
2822 }
2823 if (mat[MAT_TYPE] == INT_MATRIX) {
2824 YAP_Int *data, *ndata;
2825
2826 /* create a new matrix with the same size */
2827 tf = new_int_matrix(newdims, nindx, NULL);
2828 if (tf == YAP_TermNil())
2829 return FALSE;
2830 /* in case the matrix moved */
2831 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2832 nmat = (int *)YAP_BlobOfTerm(tf);
2833 data = matrix_long_data(mat, ndims);
2834 ndata = matrix_long_data(nmat, newdims);
2835 /* create a new matrix with smaller size */
2836 for (i = 0; i < nmat[MAT_SIZE]; i++)
2837 ndata[i] = 0;
2838 for (i = 0; i < mat[MAT_SIZE]; i++) {
2839 int j, k;
2840 /*
2841 not very efficient, we could try to take advantage of the fact
2842 that we usually only change an index at a time
2843 */
2844 matrix_get_index(mat, i, indx);
2845 for (j = 0, k = 0; j < ndims; j++) {
2846 if (!conv[j]) {
2847 nindx[k++] = indx[j];
2848 }
2849 }
2850 ndata[matrix_get_offset(nmat, nindx)] += exp(data[i]);
2851 }
2852 } else {
2853 double *data, *ndata;
2854
2855 /* create a new matrix with the same size */
2856 tf = new_float_matrix(newdims, nindx, NULL);
2857 if (tf == YAP_TermNil())
2858 return FALSE;
2859 /* in case the matrix moved */
2860 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2861 nmat = (int *)YAP_BlobOfTerm(tf);
2862 data = matrix_double_data(mat, ndims);
2863 ndata = matrix_double_data(nmat, newdims);
2864 /* create a new matrix with smaller size */
2865 for (i = 0; i < nmat[MAT_SIZE]; i++)
2866 ndata[i] = 0.0;
2867 for (i = 0; i < mat[MAT_SIZE]; i++) {
2868 int j, k;
2869 /*
2870 not very efficient, we could try to take advantage of the fact
2871 that we usually only change an index at a time
2872 */
2873 matrix_get_index(mat, i, indx);
2874 for (j = 0, k = 0; j < ndims; j++) {
2875 if (!conv[j]) {
2876 nindx[k++] = indx[j];
2877 }
2878 }
2879 ndata[matrix_get_offset(nmat, nindx)] += exp(data[i]);
2880 }
2881 for (i = 0; i < nmat[MAT_SIZE]; i++) {
2882 ndata[i] = log(ndata[i]);
2883 }
2884 }
2885 return YAP_Unify(YAP_ARG3, tf);
2886}
2887
2888/* given a matrix M and a set of dims, build a matrix to follow
2889 the new order
2890*/
2891static YAP_Bool matrix_expand(void) {
2892 int ndims, i, *dims, newdims = 0, olddims = 0;
2893 int new[MAX_DIMS], indx[MAX_DIMS], nindx[MAX_DIMS];
2894 YAP_Term tconv, tf;
2895 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2896 if (!mat) {
2897 /* Error */
2898 return FALSE;
2899 }
2900 /* we now have our target matrix, let us grab our conversion matrix */
2901 tconv = YAP_ARG2;
2902 ndims = mat[MAT_NDIMS];
2903 dims = mat + MAT_DIMS;
2904 for (i = 0; i < MAX_DIMS; i++) {
2905 YAP_Term th;
2906 YAP_Int j;
2907
2908 if (!YAP_IsPairTerm(tconv)) {
2909 if (tconv != YAP_TermNil())
2910 return FALSE;
2911 break;
2912 }
2913 th = YAP_HeadOfTerm(tconv);
2914 if (!YAP_IsIntTerm(th))
2915 return FALSE;
2916 newdims++;
2917 j = YAP_IntOfTerm(th);
2918 if (j == 0) {
2919 new[i] = 0;
2920 nindx[i] = dims[olddims];
2921 olddims++;
2922 } else {
2923 new[i] = 1;
2924 nindx[i] = j;
2925 }
2926 tconv = YAP_TailOfTerm(tconv);
2927 }
2928 if (olddims != ndims)
2929 return FALSE;
2930 if (mat[MAT_TYPE] == INT_MATRIX) {
2931 YAP_Int *data, *ndata;
2932
2933 /* create a new matrix with the same size */
2934 tf = new_int_matrix(newdims, nindx, NULL);
2935 if (tf == YAP_TermNil())
2936 return FALSE;
2937 /* in case the matrix moved */
2938 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2939 nmat = (int *)YAP_BlobOfTerm(tf);
2940 data = matrix_long_data(mat, ndims);
2941 ndata = matrix_long_data(nmat, newdims);
2942 /* create a new matrix with the same size */
2943 for (i = 0; i < nmat[MAT_SIZE]; i++) {
2944 int j, k = 0;
2945 /*
2946 not very efficient, we could try to take advantage of the fact
2947 that we usually only change an index at a time
2948 */
2949 matrix_get_index(nmat, i, indx);
2950 for (j = 0; j < newdims; j++) {
2951 if (!new[j])
2952 nindx[k++] = indx[j];
2953 }
2954 ndata[i] = data[matrix_get_offset(mat, nindx)];
2955 }
2956 } else {
2957 double *data, *ndata;
2958
2959 /* create a new matrix with the same size */
2960 tf = new_float_matrix(newdims, nindx, NULL);
2961 if (tf == YAP_TermNil())
2962 return FALSE;
2963 /* in case the matrix moved */
2964 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2965 nmat = (int *)YAP_BlobOfTerm(tf);
2966 data = matrix_double_data(mat, ndims);
2967 ndata = matrix_double_data(nmat, newdims);
2968 /* create a new matrix with the same size */
2969 for (i = 0; i < newdims; i++)
2970 indx[i] = 0;
2971 for (i = 0; i < nmat[MAT_SIZE]; i++) {
2972 int j, k = 0;
2973 /*
2974 not very efficient, we could try to take advantage of the fact
2975 that we usually only change an index at a time
2976 */
2977 for (j = 0; j < newdims; j++) {
2978 if (!new[j])
2979 nindx[k++] = indx[j];
2980 }
2981 ndata[i] = data[matrix_get_offset(mat, nindx)];
2982 matrix_next_index(nmat + MAT_DIMS, newdims, indx);
2983 }
2984 }
2985 return YAP_Unify(YAP_ARG3, tf);
2986}
2987
2988/* given a matrix M and a set of dims, build contract a matrix to follow
2989 the new order
2990*/
2991static YAP_Bool matrix_set_all_that_disagree(void) {
2992 int ndims, i, *dims;
2993 int indx[MAX_DIMS];
2994 YAP_Term tf;
2995 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2996 int dim = YAP_IntOfTerm(YAP_ARG2);
2997 int pos = YAP_IntOfTerm(YAP_ARG3);
2998
2999 if (!mat) {
3000 /* Error */
3001 return FALSE;
3002 }
3003 ndims = mat[MAT_NDIMS];
3004 dims = mat + MAT_DIMS;
3005 if (mat[MAT_TYPE] == INT_MATRIX) {
3006 YAP_Int *data, *ndata, val;
3007
3008 /* create a new matrix with the same size */
3009 tf = new_int_matrix(ndims, dims, NULL);
3010 if (tf == YAP_TermNil())
3011 return FALSE;
3012 /* in case the matrix moved */
3013 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
3014 nmat = (int *)YAP_BlobOfTerm(tf);
3015 data = matrix_long_data(mat, ndims);
3016 ndata = matrix_long_data(nmat, ndims);
3017 if (!YAP_IsIntTerm(YAP_ARG4))
3018 return FALSE;
3019 val = YAP_IntOfTerm(YAP_ARG4);
3020 /* create a new matrix with the same size */
3021 for (i = 0; i < nmat[MAT_SIZE]; i++) {
3022
3023 /*
3024 not very efficient, we could try to take advantage of the fact
3025 that we usually only change an index at a time
3026 */
3027 matrix_get_index(mat, i, indx);
3028 if (indx[dim] != pos)
3029 ndata[i] = val;
3030 else
3031 ndata[i] = data[i];
3032 }
3033 } else {
3034 double *data, *ndata, val;
3035
3036 /* create a new matrix with the same size */
3037 tf = new_float_matrix(ndims, dims, NULL);
3038 if (tf == YAP_TermNil())
3039 return FALSE;
3040 /* in case the matrix moved */
3041 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
3042 nmat = (int *)YAP_BlobOfTerm(tf);
3043 data = matrix_double_data(mat, ndims);
3044 ndata = matrix_double_data(nmat, ndims);
3045 if (YAP_IsFloatTerm(YAP_ARG4))
3046 val = YAP_FloatOfTerm(YAP_ARG4);
3047 else if (YAP_IsIntTerm(YAP_ARG4))
3048 val = YAP_IntOfTerm(YAP_ARG4);
3049 else
3050 return FALSE;
3051 /* create a new matrix with the same size */
3052 for (i = 0; i < nmat[MAT_SIZE]; i++) {
3053
3054 /*
3055 not very efficient, we could try to take advantage of the fact
3056 that we usually only change an index at a time
3057 */
3058 matrix_get_index(mat, i, indx);
3059 if (indx[dim] != pos)
3060 ndata[i] = val;
3061 else
3062 ndata[i] = data[i];
3063 }
3064 }
3065 return YAP_Unify(YAP_ARG5, tf);
3066}
3067
3068static YAP_Bool matrix_m(void) {
3069 int ndims, i, size;
3070 YAP_Term tm, *tp;
3071 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
3072
3073 if (!mat) {
3074 return YAP_Unify(YAP_ARG1, YAP_ARG2);
3075 }
3076 ndims = mat[MAT_NDIMS];
3077 size = mat[MAT_SIZE];
3078 tm = YAP_MkNewApplTerm(FunctorM, 5);
3079 tp = YAP_ArgsOfTerm(tm);
3080 tp[0] = mk_int_list(ndims, mat + MAT_DIMS);
3081 tp[1] = YAP_MkIntTerm(ndims);
3082 tp[2] = YAP_MkIntTerm(size);
3083 tp[3] = mk_rep_int_list(ndims, mat[MAT_BASE]);
3084 tp[4] = YAP_MkNewApplTerm(YAP_MkFunctor(AtomC, size), size);
3085 tp = YAP_ArgsOfTerm(tp[3]);
3086 if (mat[MAT_TYPE] == INT_MATRIX) {
3087 YAP_Int *data;
3088
3089 /* in case the matrix moved */
3090 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
3091 data = matrix_long_data(mat, ndims);
3092 for (i = 0; i < mat[MAT_SIZE]; i++) {
3093 tp[i] = YAP_MkIntTerm(data[i]);
3094 }
3095 } else {
3096 double *data;
3097
3098 /* in case the matrix moved */
3099 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
3100 data = matrix_double_data(mat, ndims);
3101 for (i = 0; i < mat[MAT_SIZE]; i++) {
3102 tp[i] = YAP_MkFloatTerm(data[i]);
3103 }
3104 }
3105 return YAP_Unify(YAP_ARG2, tm);
3106}
3107
3108static YAP_Bool is_matrix(void) {
3109 YAP_Term t = YAP_ARG1;
3110 int *mat = (int *)YAP_BlobOfTerm(t);
3111
3112 if (!mat) {
3113 if (!YAP_IsApplTerm(t))
3114 return FALSE;
3115 return YAP_FunctorOfTerm(t) == FunctorM;
3116 }
3117 return TRUE;
3118}
3119
3120static YAP_Bool get_float_from_address(void) {
3121 YAP_Float *fp = (YAP_Float *)YAP_IntOfTerm(YAP_ARG1);
3122 YAP_Int off = YAP_IntOfTerm(YAP_ARG2);
3123
3124 return YAP_Unify(YAP_ARG3, YAP_MkFloatTerm(fp[off]));
3125}
3126
3127static YAP_Bool set_float_from_address(void) {
3128 YAP_Float f, *fp = (YAP_Float *)YAP_IntOfTerm(YAP_ARG1);
3129 YAP_Int off = YAP_IntOfTerm(YAP_ARG2);
3130 if (YAP_IsIntTerm(YAP_ARG3)) {
3131 f = YAP_IntOfTerm(YAP_ARG3);
3132 }else {
3133 f = YAP_FloatOfTerm(YAP_ARG3);
3134 }
3135 fp[off] = f;
3136
3137 return true;
3138}
3139
3140
3141static YAP_Bool address_op_to_all(void) {
3142 YAP_Term top = YAP_ARG3;
3143 op_type op;
3144
3145 YAP_Float *fp = (YAP_Float *)YAP_IntOfTerm(YAP_ARG1);
3146 YAP_Int sz = YAP_IntOfTerm(YAP_ARG2);
3147 if (!YAP_IsIntTerm(top)) {
3148 return FALSE;
3149 }
3150 op = YAP_IntOfTerm(top);
3151 YAP_Term tnum = YAP_ARG4;
3152 double num;
3153
3154 if (YAP_IsFloatTerm(tnum)) {
3155 num = YAP_FloatOfTerm(tnum);
3156 } else if (!YAP_IntOfTerm(tnum)) {
3157 return FALSE;
3158 } else {
3159 num = (double)YAP_IntOfTerm(tnum);
3160 }
3161 switch (op) {
3162 case MAT_PLUS: {
3163 int i;
3164
3165 for (i = 0; i < sz; i++) {
3166 fp[i] = fp[i] + num;
3167 }
3168 } break;
3169 case MAT_SUB: {
3170 int i;
3171
3172 for (i = 0; i < sz; i++) {
3173 fp[i] = num - fp[i];
3174 }
3175 } break;
3176 case MAT_TIMES: {
3177 int i;
3178
3179 for (i = 0; i < sz; i++) {
3180 fp[i] = fp[i] * num;
3181 }
3182 } break;
3183 case MAT_DIV: {
3184 int i;
3185
3186 for (i = 0; i < sz; i++) {
3187 fp[i] = fp[i] / num;
3188 }
3189 } break;
3190 default:
3191 return FALSE;
3192 }
3193 return true;
3194}
3195
3196static YAP_Bool address_to_list(void) {
3197 YAP_Term t = YAP_TermNil();
3198 int i;
3199 YAP_Float *fp = (YAP_Float *)YAP_IntOfTerm(YAP_ARG1);
3200 YAP_Int sz = YAP_IntOfTerm(YAP_ARG2);
3201 for (i = 0; i< sz; i++)
3202 t = YAP_MkPairTerm(YAP_MkFloatTerm(fp[sz-i-1]),t);
3203 return YAP_Unify(YAP_ARG3, t);
3204}
3205
3206static YAP_Bool address_to_max(void) {
3207 int i;
3208 YAP_Float *fp = (YAP_Float *)YAP_IntOfTerm(YAP_ARG1);
3209 YAP_Int sz = YAP_IntOfTerm(YAP_ARG2);
3210 YAP_Float max = fp[0];
3211 for (i = 1; i< sz; i++)
3212 if (fp[i] > max) max = fp[i];
3213 return YAP_Unify(YAP_ARG3, YAP_MkFloatTerm(max));
3214}
3215
3216static YAP_Bool address_to_min(void) {
3217 int i;
3218 YAP_Float *fp = (YAP_Float *)YAP_IntOfTerm(YAP_ARG1);
3219 YAP_Int sz = YAP_IntOfTerm(YAP_ARG2);
3220 YAP_Float min = fp[0];
3221 for (i = 1; i< sz; i++)
3222 if (fp[i] < min) min = fp[i];
3223 return YAP_Unify(YAP_ARG3, YAP_MkFloatTerm(min));
3224}
3225
3226static YAP_Bool address_set_all(void) {
3227 int i;
3228 YAP_Float f, *fp = (YAP_Float *)YAP_IntOfTerm(YAP_ARG1);
3229 YAP_Int sz = YAP_IntOfTerm(YAP_ARG2);
3230 if (YAP_IsFloatTerm(YAP_ARG3)) {
3231 f = YAP_FloatOfTerm(YAP_ARG3);
3232 } else {
3233 f = YAP_IntOfTerm(YAP_ARG3);
3234 }
3235 for (i = 0; i< sz; i++)
3236 fp[i] = f;
3237 return true;
3238}
3239
3240static YAP_Bool address_to_sum(void) {
3241 YAP_Term t = YAP_TermNil();
3242 int i;
3243 YAP_Float *data = (YAP_Float *)YAP_IntOfTerm(YAP_ARG1);
3244 YAP_Int sz = YAP_IntOfTerm(YAP_ARG2);
3245 double sum = 0.0;
3246 // function KahanSum(input)
3247 double c = 0.0; // A running compensation for lost low-order bits.
3248 for (i = 0; i < sz; i++) {
3249 double y = data[i] - c; // So far, so good: c is zero.
3250 double t =
3251 sum +
3252 y; // Alas, sum is big, y small, so low-order digits of y are lost.
3253 c = (t - sum) - y; // (t - sum) cancels the high-order part of y;
3254 // subtracting y recovers negative (low part of y)
3255 sum = t; // Algebraically, c should always be zero. Beware
3256 // overly-aggressive optimizing compilers!
3257 }
3258 t = YAP_MkFloatTerm(sum);
3259 return YAP_Unify(YAP_ARG3, t);
3260}
3261
3262
3263
3264X_API void init_matrix(void);
3265
3266X_API void init_matrix(void) {
3267 AtomC = YAP_LookupAtom("c");
3268 FunctorM = YAP_MkFunctor(YAP_LookupAtom("$matrix"), 5);
3269 FunctorFloats = YAP_MkFunctor(YAP_LookupAtom("floats"), 2);
3270
3271 YAP_UserCPredicate("new_ints_matrix", new_ints_matrix, 4);
3272 YAP_UserCPredicate("new_ints_matrix_set", new_ints_matrix_set, 4);
3273 YAP_UserCPredicate("new_floats_matrix", new_floats_matrix, 4);
3274 YAP_UserCPredicate("new_floats_matrix_set", new_floats_matrix_set, 4);
3275 YAP_UserCPredicate("matrix_set", matrix_set, 3);
3276 YAP_UserCPredicate("matrix_set", matrix_set2, 2);
3277 YAP_UserCPredicate("matrix_set_all", matrix_set_all, 2);
3278 YAP_UserCPredicate("matrix_add", matrix_add, 3);
3279 YAP_UserCPredicate("matrix_get", matrix_get, 3);
3280 YAP_UserCPredicate("matrix_get", do_matrix_access2, 2);
3281 YAP_UserCPredicate("matrix_inc", do_matrix_inc, 2);
3282 YAP_UserCPredicate("matrix_dec", do_matrix_dec, 2);
3283 YAP_UserCPredicate("matrix_inc", do_matrix_inc2, 3);
3284 YAP_UserCPredicate("matrix_dec", do_matrix_dec2, 3);
3285 YAP_UserCPredicate("matrixn_to_list", matrix_to_list, 2);
3286 YAP_UserCPredicate("matrixn_set_base", matrix_set_base, 2);
3287 YAP_UserCPredicate("matrixn_dims", matrix_dims, 2);
3288 YAP_UserCPredicate("matrixn_dims", matrix_dims3, 3);
3289 YAP_UserCPredicate("matrixn_ndims", matrix_ndims, 2);
3290 YAP_UserCPredicate("matrixn_size", matrix_size, 2);
3291 YAP_UserCPredicate("matrix_type_as_number", matrix_type, 2);
3292 YAP_UserCPredicate("matrixn_arg_to_offset", matrix_arg_to_offset, 3);
3293 YAP_UserCPredicate("matrixn_offset_to_arg", matrix_offset_to_arg, 3);
3294 YAP_UserCPredicate("matrixn_max", matrix_max, 2);
3295 YAP_UserCPredicate("matrixn_maxarg", matrix_maxarg, 2);
3296 YAP_UserCPredicate("matrixn_min", matrix_min, 2);
3297 YAP_UserCPredicate("matrixn_minarg", matrix_minarg, 2);
3298 YAP_UserCPredicate("matrix_sum", matrix_sum, 2);
3299 YAP_UserCPredicate("matrix_shuffle", matrix_transpose, 3);
3300 YAP_UserCPredicate("matrix_expand", matrix_expand, 3);
3301 YAP_UserCPredicate("matrix_select", matrix_select, 4);
3302 YAP_UserCPredicate("matrix_column", matrix_column, 3);
3303 YAP_UserCPredicate("matrix_to_logs", matrix_log_all, 1);
3304 YAP_UserCPredicate("matrix_to_exps", matrix_exp_all, 1);
3305 YAP_UserCPredicate("matrix_to_exps2", matrix_exp2_all, 1);
3306 YAP_UserCPredicate("matrixn_to_logs", matrix_log_all2, 2);
3307 YAP_UserCPredicate("matrixn_to_exps", matrix_exp_all2, 2);
3308 YAP_UserCPredicate("matrix_sum_out", matrix_sum_out, 3);
3309 YAP_UserCPredicate("matrix_sum_out_several", matrix_sum_out_several, 3);
3310 YAP_UserCPredicate("matrix_sum_logs_out", matrix_sum_out_logs, 3);
3311 YAP_UserCPredicate("matrix_sum_logs_out_several", matrix_sum_out_logs_several,
3312 3);
3313 YAP_UserCPredicate("matrix_set_all_that_disagree",
3314 matrix_set_all_that_disagree, 5);
3315 YAP_UserCPredicate("do_matrix_op", matrix_op, 4);
3316 YAP_UserCPredicate("do_matrix_agg_lines", matrix_agg_lines, 3);
3317 YAP_UserCPredicate("do_matrix_agg_cols", matrix_agg_cols, 3);
3318 YAP_UserCPredicate("do_matrix_op_to_all", matrix_op_to_all, 4);
3319 YAP_UserCPredicate("do_matrix_op_to_lines", matrix_op_to_lines, 4);
3320 YAP_UserCPredicate("do_matrix_op_to_cols", matrix_op_to_cols, 4);
3321 YAP_UserCPredicate("matrix_m", matrix_m, 2);
3322 YAP_UserCPredicate("matrix", is_matrix, 1);
3323 YAP_UserCPredicate("get_float_from_address", get_float_from_address, 3);
3324 YAP_UserCPredicate("set_float_from_address", set_float_from_address, 3);
3325 YAP_UserCPredicate("address_to_list", address_to_list, 3);
3326 YAP_UserCPredicate("address_to_sum", address_to_sum, 3);
3327 YAP_UserCPredicate("address_to_min", address_to_min, 3);
3328 YAP_UserCPredicate("address_to_max", address_to_max, 3);
3329 YAP_UserCPredicate("address_set_all", address_set_all, 3);
3330 YAP_UserCPredicate("address_op_to_all", address_op_to_all, 3);}
3331
3332#ifdef _WIN32
3333
3334int WINAPI win_matrixs(HANDLE, DWORD, LPVOID);
3335
3336int WINAPI win_matrixs(HANDLE hinst, DWORD reason, LPVOID reserved) {
3337 switch (reason) {
3338 case DLL_PROCESS_ATTACH:
3339 break;
3340 case DLL_PROCESS_DETACH:
3341 break;
3342 case DLL_THREAD_ATTACH:
3343 break;
3344 case DLL_THREAD_DETACH:
3345 break;
3346 }
3347 return 1;
3348}
3349#endif
m
A matrix.
Definition: matrix.c:68