interval.cc
Go to the documentation of this file.
1 #include "kernel/mod2.h"
2 #include "Singular/blackbox.h"
4 #include "Singular/ipshell.h" // for iiCheckTypes
6 #include "Singular/mod_lib.h"
7 
8 /*
9  * CONSTRUCTORS & DESTRUCTORS
10  */
11 
12 /* interval */
13 
14 interval::interval(const ring r)
15 {
16  lower = n_Init(0, r->cf);
17  upper = n_Init(0, r->cf);
18  R = r;
19  R->ref++;
20 }
21 
22 interval::interval(number a, const ring r)
23 {
24  // dangerous: check if a is in coefs r->cf
25  lower = a;
26  upper = n_Copy(a, r->cf);
27  R = r;
28  R->ref++;
29 }
30 
31 interval::interval(number a, number b, const ring r)
32 {
33  lower = a;
34  upper = b;
35  R = r;
36  R->ref++;
37 }
38 
40 {
41  lower = n_Copy(I->lower, I->R->cf);
42  upper = n_Copy(I->upper, I->R->cf);
43  R = I->R;
44  R->ref++;
45 }
46 
48 {
49  n_Delete(&lower, R->cf);
50  n_Delete(&upper, R->cf);
51  R->ref--;
52 }
53 
55 {
56  if (R != r)
57  {
58  // check if coefficient fields are equal
59  // if not pass numbers to new cf-field
60  if (R->cf != r->cf)
61  {
62  nMapFunc fun = n_SetMap(R->cf, r->cf);
63  number lo = fun(lower, R->cf, r->cf),
64  up = fun(upper, R->cf, r->cf);
65  n_Delete(&lower, R->cf);
66  n_Delete(&upper, R->cf);
67  lower = lo;
68  upper = up;
69  }
70  R->ref--;
71  r->ref++;
72  R = r;
73  }
74 
75  return *this;
76 }
77 
78 /* box */
79 
81 {
82  R = currRing;
83  int i, n = R->N;
84  intervals = (interval**) omAlloc0(n * sizeof(interval*));
85  if (intervals != NULL)
86  {
87  for (i = 0; i < n; i++)
88  {
89  intervals[i] = new interval();
90  }
91  }
92  R->ref++;
93 }
94 
96 {
97  R = B->R;
98  int i, n = R->N;
99  intervals = (interval**) omAlloc0(n * sizeof(interval*));
100  if (intervals != NULL)
101  {
102  for (i = 0; i < n; i++)
103  {
104  intervals[i] = new interval(B->intervals[i]);
105  }
106  }
107  R->ref++;
108 }
109 
111 {
112  int i, n = R->N;
113  for (i = 0; i < n; i++)
114  {
115  delete intervals[i];
116  }
117  omFree((void**) intervals);
118  R->ref--;
119 }
120 
121 // note: does not copy input
123 {
124  if (0 <= i && i < R->N)
125  {
126  delete intervals[i];
127  intervals[i] = I;
128  }
129  return *this;
130 }
131 
132 /*
133  * TYPE IDs
134  */
135 
136 static int intervalID;
137 static int boxID;
138 
139 /*
140  * INTERVAL FUNCTIONS
141  */
142 
143 static void* interval_Init(blackbox*)
144 {
145  return (void*) new interval();
146 }
147 
148 // convert interval to string
149 static char* interval_String(blackbox*, void *d)
150 {
151  if (d == NULL)
152  {
153  // invalid object
154  return omStrDup("[?]");
155  }
156  else
157  {
158  interval* i = (interval*) d;
159 
160  // use n_Write since nothing better (?) exists
161  StringSetS("[");
162  n_Write(i->lower, i->R->cf);
163  StringAppendS(", ");
164  n_Write(i->upper, i->R->cf);
165  StringAppendS("]");
166 
167  return StringEndS();
168  }
169 }
170 
171 static void* interval_Copy(blackbox*, void *d)
172 {
173  return (void*) new interval((interval*) d);
174 }
175 
176 // destroy interval
177 static void interval_Destroy(blackbox*, void *d)
178 {
179  if (d != NULL)
180  delete (interval*) d;
181 }
182 
183 // assigning values to intervals
185 {
186  interval *RES;
187 
188  /*
189  * Allow assignments of the form
190  * I = a,
191  * I = a, b,
192  * I = J
193  * where a, b numbers or ints and J interval
194  */
195 
196  if (args->Typ() == intervalID)
197  {
198  RES = new interval((interval*) args->CopyD());
199  }
200  else
201  {
202  number n1, n2;
203 
204  if (args->Typ() == INT_CMD)
205  {
206  n1 = nInit((int)(long) args->Data());
207  }
208  else if (args->Typ() == NUMBER_CMD)
209  {
210  n1 = (number) args->CopyD();
211  }
212  else
213  {
214  WerrorS("Input not supported: first argument not int or number");
215  return TRUE;
216  }
217 
218  // check if second argument exists
219  if (args->next == NULL)
220  {
221  RES = new interval(n1);
222  }
223  else
224  {
225  if (args->next->Typ() == INT_CMD)
226  {
227  n2 = nInit((int)(long) args->next->Data());
228  }
229  else if (args->next->Typ() == NUMBER_CMD)
230  {
231  n2 = (number) args->next->CopyD();
232  }
233  else
234  {
235  WerrorS("Input not supported: second argument not int or number");
236  return TRUE;
237  }
238 
239  RES = new interval(n1, n2);
240  }
241  }
242 
243  // destroy data of result if it exists
244  if (result->Data() != NULL)
245  {
246  delete (interval*) result->Data();
247  }
248 
249  if (result->rtyp == IDHDL)
250  {
251  IDDATA((idhdl)result->data) = (char*) RES;
252  }
253  else
254  {
255  result->rtyp = intervalID;
256  result->data = (void*) RES;
257  }
258 
259  args->CleanUp();
260  return FALSE;
261 }
262 
264 {
265  if (arg != NULL && arg->Typ() == intervalID)
266  {
267  interval *I = (interval*) arg->Data();
268  result->rtyp = NUMBER_CMD;
269  result->data = (void*) n_Sub(I->upper, I->lower, I->R->cf);
270  arg->CleanUp();
271  return FALSE;
272  }
273 
274  WerrorS("syntax: length(<interval>)");
275  return TRUE;
276 }
277 
278 // interval -> interval procedures
279 
281 {
282  // must assume a is in ring I->R
283  number lo, up;
284  if (nGreaterZero(a))
285  {
286  lo = n_Mult(a, I->lower, I->R->cf);
287  up = n_Mult(a, I->upper, I->R->cf);
288  }
289  else
290  {
291  lo = n_Mult(a, I->upper, I->R->cf);
292  up = n_Mult(a, I->lower, I->R->cf);
293  }
294 
295  n_Normalize(lo, I->R->cf);
296  n_Normalize(up, I->R->cf);
297 
298  return new interval(lo, up, I->R);
299 }
300 
302 {
303  number lo, up;
304  number nums[4];
305  nums[0] = n_Mult(I->lower, J->lower, I->R->cf);
306  nums[1] = n_Mult(I->lower, J->upper, I->R->cf);
307  nums[2] = n_Mult(I->upper, J->lower, I->R->cf);
308  nums[3] = n_Mult(I->upper, J->upper, I->R->cf);
309 
310  int i, imax = 0, imin = 0;
311  for (i = 1; i < 4; i++)
312  {
313  if (n_Greater(nums[i], nums[imax], I->R->cf))
314  {
315  imax = i;
316  }
317  if (n_Greater(nums[imin], nums[i], I->R->cf))
318  {
319  imin = i;
320  }
321  }
322 
323  lo = n_Copy(nums[imin], I->R->cf);
324  up = n_Copy(nums[imax], I->R->cf);
325 
326  // delete products
327  for (i = 0; i < 4; i++)
328  {
329  n_Delete(&nums[i], I->R->cf);
330  }
331 
332  n_Normalize(lo, I->R->cf);
333  n_Normalize(up, I->R->cf);
334 
335  return new interval(lo, up, I->R);
336 }
337 
339 {
340  number lo = n_Add(I->lower, J->lower, I->R->cf),
341  up = n_Add(I->upper, J->upper, I->R->cf);
342 
343  n_Normalize(lo, I->R->cf);
344  n_Normalize(up, I->R->cf);
345 
346  return new interval(lo, up);
347 }
348 
350 {
351  number lo = n_Sub(I->lower, J->upper, I->R->cf),
352  up = n_Sub(I->upper, J->lower, I->R->cf);
353 
354  n_Normalize(lo, I->R->cf);
355  n_Normalize(up, I->R->cf);
356 
357  return new interval(lo, up, I->R);
358 }
359 
360 static bool intervalEqual(interval *I, interval *J)
361 {
362  assume(I->R == J->R);
363  return n_Equal(I->lower, J->lower, I->R->cf)
364  && n_Equal(I->upper, J->upper, I->R->cf);
365 }
366 
367 // ckeck if zero is contained in an interval
369 {
370  number n = n_Mult(I->lower, I->upper, I->R->cf);
371  bool result = !n_GreaterZero(n, I->R->cf);
372  // delete helper number
373  n_Delete(&n, I->R->cf);
374 
375  return result;
376 }
377 
379 {
380  if (p == 0)
381  {
382  return new interval(n_Init(1,I->R->cf), I->R);
383  }
384 
385  // no initialisation required (?)
386  number lo, up;
387 
388  n_Power(I->lower, p, &lo, I->R->cf);
389  n_Power(I->upper, p, &up, I->R->cf);
390 
391  // should work now
392  if (p % 2 == 1)
393  {
394  return new interval(lo, up, I->R);
395  }
396  else
397  {
398  // perform pointer swap if necessary
399  number tmp;
400  if (n_Greater(lo, up, I->R->cf))
401  {
402  tmp = up;
403  up = lo;
404  lo = tmp;
405  }
406 
407  if (intervalContainsZero(I))
408  {
409  n_Delete(&lo, I->R->cf);
410  lo = n_Init(0, I->R->cf);
411  }
412  return new interval(lo, up, I->R);
413  }
414 }
415 
416 /*
417  * BINARY OPERATIONS:
418  * Cases handled:
419  * I + J,
420  *
421  * I - J,
422  *
423  * a * J,
424  * I * a,
425  * I * J,
426  *
427  * a / J,
428  * I / a,
429  * I / J
430  *
431  * I ^ n,
432  *
433  * I == J
434  *
435  * where I, J, interval, a, b int or number, n int
436  */
437 
438 static BOOLEAN interval_Op2(int op, leftv result, leftv i1, leftv i2)
439 {
440  interval *RES;
441 
442  switch(op)
443  {
444  case '+':
445  {
446  if (i1->Typ() != intervalID || i2->Typ() != intervalID)
447  {
448  WerrorS("syntax: <interval> + <interval>");
449  return TRUE;
450  }
451  interval *I1, *I2;
452  I1 = (interval*) i1->Data();
453  I2 = (interval*) i2->Data();
454  if (I1->R != I2->R)
455  {
456  WerrorS("adding intervals defined in different rings not supported");
457  return TRUE;
458  }
459 
460  RES = intervalAdd(I1, I2);
461  break;
462  }
463  case '-':
464  {
465  if (i1->Typ() != intervalID || i2->Typ() != intervalID)
466  {
467  WerrorS("syntax: <interval> - <interval>");
468  return TRUE;
469  }
470  interval *I1, *I2;
471  I1 = (interval*) i1->Data();
472  I2 = (interval*) i2->Data();
473  if (I1->R != I2->R)
474  {
475  WerrorS("subtracting intervals defined in different rings not supported");
476  return TRUE;
477  }
478 
479  RES = intervalSubtract(I1, I2);
480  break;
481  }
482  case '*':
483  {
484  if (i1->Typ() == i2->Typ())
485  {
486  // both must be intervals
487  interval *I1, *I2;
488  I1 = (interval*) i1->Data();
489  I2 = (interval*) i2->Data();
490  if (I1->R != I2->R)
491  {
492  WerrorS("multiplying intervals defined in different rings not supported");
493  return TRUE;
494  }
495 
496  RES = intervalMultiply(I1, I2);
497  }
498  else
499  {
500  // one arg is scalar, one is interval
501  // give new names to reduce to one case
502  leftv iscalar, iinterv;
503  if (i1->Typ() != intervalID)
504  {
505  // i1 is scalar
506  iscalar = i1;
507  iinterv = i2;
508  }
509  else
510  {
511  // i2 is scalar
512  iscalar = i2;
513  iinterv = i1;
514  }
515  number n;
516 
517  switch (iscalar->Typ())
518  {
519  case INT_CMD:
520  { n = nInit((int)(long) iscalar->Data()); break; }
521  case NUMBER_CMD:
522  { n = (number) iscalar->CopyD(); break; }
523  default:
524  { WerrorS("first argument not int/number/interval"); return TRUE; }
525  }
526 
527  interval *I = (interval*) iinterv->Data();
528 
529  RES = intervalScalarMultiply(n, I);
530  // n no longer needed, delete it
531  nDelete(&n);
532  }
533 
534  break;
535  }
536  case '/':
537  {
538  if(i2->Typ() == intervalID)
539  {
540  interval *I2;
541  I2 = (interval*) i2->Data();
542 
543  // make sure I2 is invertible
544  if (intervalContainsZero(I2))
545  {
546  WerrorS("second interval contains zero");
547  return TRUE;
548  }
549 
550  number invlo, invup;
551  invlo = n_Invers(I2->lower, I2->R->cf);
552  invup = n_Invers(I2->upper, I2->R->cf);
553 
554  // inverse interval
555  interval *I2inv = new interval(invup, invlo, I2->R);
556 
557  if (i1->Typ() == intervalID)
558  {
559  interval *I1 = (interval*) i1->Data();
560  if (I1->R != I2->R)
561  {
562  WerrorS("dividing intervals from different rings not supported");
563  delete I2inv;
564  return TRUE;
565  }
566  RES = intervalMultiply(I1, I2inv);
567  }
568  else
569  {
570  // i1 is not an interval
571  number n;
572  switch (i1->Typ())
573  {
574  case INT_CMD:
575  {
576  n = nInit((int)(long) i1->Data());
577  break;
578  }
579  case NUMBER_CMD:
580  {
581  n = nCopy((number) i1->Data());
582  break;
583  }
584  default:
585  {
586  WerrorS("first argument not int/number/interval");
587  delete I2inv;
588  return TRUE;
589  }
590  }
591  RES = intervalScalarMultiply(n, I2inv);
592  nDelete(&n);
593  }
594 
595  delete I2inv;
596  break;
597  }
598  else
599  {
600  // i2 is not an interval
601  interval *I1 = (interval*) i1->Data();
602  number n;
603  switch(i2->Typ())
604  {
605  case INT_CMD:
606  { n = nInit((int)(long) i2->Data()); break; }
607  case NUMBER_CMD:
608  { n = nCopy((number) i2->Data()); break; }
609  default:
610  {
611  WerrorS("second argument not int/number/interval");
612  return TRUE;
613  }
614  }
615  // handle zero to prevent memory leak (?)
616  if (nIsZero(n))
617  {
618  WerrorS("<interval>/0 not supported");
619  return TRUE;
620  }
621 
622  number nInv = nInvers(n);
623  nDelete(&n);
624  RES = intervalScalarMultiply(nInv, I1);
625  nDelete(&nInv);
626 
627  break;
628  }
629 
630  break;
631  }
632  case '^':
633  {
634  if (i1->Typ() != intervalID || i2->Typ() != INT_CMD)
635  {
636  WerrorS("syntax: <interval> ^ <int>");
637  return TRUE;
638  }
639  int p = (int)(long) i2->Data();
640  if (p < 0)
641  {
642  WerrorS("<interval> ^ n not implemented for n < 0");
643  return TRUE;
644  }
645  interval *I = (interval*) i1->Data();
646 
647  RES = intervalPower(I, p);
648  break;
649  }
650  case EQUAL_EQUAL:
651  {
652  if (i1->Typ() != intervalID || i2->Typ() != intervalID)
653  {
654  WerrorS("syntax: <interval> == <interval>");
655  return TRUE;
656  }
657  interval *I1, *I2;
658  I1 = (interval*) i1->Data();
659  I2 = (interval*) i2->Data();
660 
661  result->rtyp = INT_CMD;
662  result->data = (void*) intervalEqual(I1, I2);
663  i1->CleanUp();
664  i2->CleanUp();
665  return FALSE;
666  }
667  case '[':
668  {
669  if (i1->Typ() != intervalID || i2->Typ() != INT_CMD)
670  {
671  WerrorS("syntax: <interval>[<int>]");
672  return TRUE;
673  }
674 
675  interval *I = (interval*) i1->Data();
676  int n = (int)(long) i2->Data();
677 
678  number out;
679  if (n == 1)
680  {
681  out = nCopy(I->lower);
682  }
683  else if (n == 2)
684  {
685  out = nCopy(I->upper);
686  }
687  else
688  {
689  WerrorS("Allowed indices are 1 and 2");
690  return TRUE;
691  }
692 
693  // delete number in result
694  if (result != NULL && result->Data() != NULL)
695  {
696  number r = (number) result->Data();
697  nDelete(&r);
698  }
699 
700  result->rtyp = NUMBER_CMD;
701  result->data = (void*) out;
702  i1->CleanUp();
703  i2->CleanUp();
704  return FALSE;
705  }
706  default:
707  {
708  // use default error
709  return blackboxDefaultOp2(op, result, i1, i2);
710  }
711  }
712 
713  // destroy data of result if it exists
714  if (result->Data() != NULL)
715  {
716  delete (interval*) result->Data();
717  }
718 
719  result->rtyp = intervalID;
720  result->data = (void*) RES;
721  i1->CleanUp();
722  i2->CleanUp();
723  return FALSE;
724 }
725 
726 static BOOLEAN interval_serialize(blackbox*, void *d, si_link f)
727 {
728  /*
729  * Format: "interval" setring number number
730  */
731  interval *I = (interval*) d;
732 
733  sleftv l, lo, up;
734  memset(&l, 0, sizeof(l));
735  memset(&lo, 0, sizeof(lo));
736  memset(&up, 0, sizeof(up));
737 
738  l.rtyp = STRING_CMD;
739  l.data = (void*) "interval";
740 
741  lo.rtyp = NUMBER_CMD;
742  lo.data = (void*) I->lower;
743 
744  up.rtyp = NUMBER_CMD;
745  up.data = (void*) I->upper;
746 
747  f->m->Write(f, &l);
748 
749  ring R = I->R;
750 
751  f->m->SetRing(f, R, TRUE);
752  f->m->Write(f, &lo);
753  f->m->Write(f, &up);
754 
755  // no idea
756  if (R != currRing)
757  f->m->SetRing(f, currRing, FALSE);
758 
759  return FALSE;
760 }
761 
762 static BOOLEAN interval_deserialize(blackbox**, void **d, si_link f)
763 {
764  leftv l;
765 
766  l = f->m->Read(f);
767  l->next = f->m->Read(f);
768 
769  number lo = (number) l->CopyD(),
770  up = (number) l->next->CopyD();
771 
772  l->CleanUp();
773 
774  *d = (void*) new interval(lo, up);
775  return FALSE;
776 }
777 
778 /*
779  * BOX FUNCTIONS
780  */
781 
782 static void* box_Init(blackbox*)
783 {
784  return (void*) new box();
785 }
786 
787 static void* box_Copy(blackbox*, void *d)
788 {
789  return (void*) new box((box*) d);
790 }
791 
792 static void box_Destroy(blackbox*, void *d)
793 {
794  if (d != NULL)
795  delete (box*) d;
796 }
797 
798 static char* box_String(blackbox*, void *d)
799 {
800  blackbox *b_i = getBlackboxStuff(intervalID);
801  box *B = (box*) d;
802  int i, n = B->R->N;
803 
804  if (B == NULL || B->intervals == NULL)
805  {
806  return omStrDup("ooo");
807  }
808 
809  StringSetS(interval_String(b_i, (void*) B->intervals[0]));
810 
811  for (i = 1; i < n; i++)
812  {
813  // interpret box as Cartesian product, hence use " x "
814  StringAppendS(" x ");
815  StringAppendS(interval_String(b_i, (void*) B->intervals[i]));
816  }
817  return StringEndS();
818 }
819 
820 // assigning values to intervals
822 {
823  box *RES;
824 
825  /*
826  * Allow assignments of the form
827  *
828  * B = C,
829  * B = l,
830  *
831  * where B, C boxes, l list of intervals
832  */
833 
834  if (args->Typ() == boxID)
835  {
836  box *B = (box*) args->Data();
837  RES = new box(B);
838  }
839  else if (args->Typ() == LIST_CMD)
840  {
841  RES = new box();
842  lists l = (lists) args->Data();
843 
844  int i, m = lSize(l), n = currRing->N;
845  // minimum
846  int M = m > (n-1) ? (n-1) : m;
847 
848  for (i = 0; i <= M; i++)
849  {
850  if (l->m[i].Typ() != intervalID)
851  {
852  WerrorS("list contains non-intervals");
853  delete RES;
854  args->CleanUp();
855  return TRUE;
856  }
857  RES->setInterval(i, (interval*) l->m[i].CopyD());
858 
859  // make sure rings of boxes and their intervals are consistent
860  // this is important for serialization
861  RES->intervals[i]->setRing(RES->R);
862  }
863  }
864  else
865  {
866  WerrorS("Input not supported: first argument not box, list, or interval");
867  return TRUE;
868  }
869 
870  // destroy data of result if it exists
871  if (result != NULL && result->Data() != NULL)
872  {
873  delete (box*) result->Data();
874  }
875 
876  if (result->rtyp == IDHDL)
877  {
878  IDDATA((idhdl)result->data) = (char*) RES;
879  }
880  else
881  {
882  result->rtyp = boxID;
883  result->data = (void*) RES;
884  }
885  args->CleanUp();
886 
887  return FALSE;
888 }
889 
890 static BOOLEAN box_Op2(int op, leftv result, leftv b1, leftv b2)
891 {
892  if (b1 == NULL || b1->Typ() != boxID)
893  {
894  Werror("first argument is not box but type(%d), second is type(%d)",
895  b1->Typ(), b2->Typ());
896  return TRUE;
897  }
898 
899  box *B1 = (box*) b1->Data();
900  int n = B1->R->N;
901 
902  box *RES;
903  switch(op)
904  {
905  case '[':
906  {
907  if (b2 == NULL || b2->Typ() != INT_CMD)
908  {
909  WerrorS("second argument not int");
910  return TRUE;
911  }
912  if (result->Data() != NULL)
913  {
914  delete (interval*) result->Data();
915  }
916 
917  int i = (int)(long) b2->Data();
918 
919  if (i < 1 || i > n)
920  {
921  WerrorS("index out of bounds");
922  return TRUE;
923  }
924 
925  // delete data of result
926  if (result->Data() != NULL)
927  {
928  delete (interval*) result->Data();
929  }
930 
931  result->rtyp = intervalID;
932  result->data = (void*) new interval(B1->intervals[i-1]);
933  b1->CleanUp();
934  b2->CleanUp();
935  return FALSE;
936  }
937  case '-':
938  {
939  if (b2 == NULL || b2->Typ() != boxID)
940  {
941  WerrorS("second argument not box");
942  return TRUE;
943  }
944 
945  box *B2 = (box*) b2->Data();
946  // maybe try to skip this initialisation
947  // copying def of box() results in segfault?
948  if (B1->R != B2->R)
949  {
950  WerrorS("subtracting boxes from different rings not supported");
951  return TRUE;
952  }
953  RES = new box();
954  int i;
955  for (i = 0; i < n; i++)
956  {
957  RES->setInterval(i, intervalSubtract(B1->intervals[i], B2->intervals[i]));
958  }
959 
960  if (result->Data() != NULL)
961  {
962  delete (box*) result->Data();
963  }
964 
965  result->rtyp = boxID;
966  result->data = (void*) RES;
967  b1->CleanUp();
968  b2->CleanUp();
969  return FALSE;
970  }
971  case EQUAL_EQUAL:
972  {
973  if (b2 == NULL || b2->Typ() != boxID)
974  {
975  WerrorS("second argument not box");
976  }
977  box *B2 = (box*) b2->Data();
978  int i;
979  bool res = true;
980  for (i = 0; i < n; i++)
981  {
982  if (!intervalEqual(B1->intervals[i], B2->intervals[i]))
983  {
984  res = false;
985  break;
986  }
987  }
988 
989  result->rtyp = INT_CMD;
990  result->data = (void*) res;
991  b1->CleanUp();
992  b2->CleanUp();
993  return FALSE;
994  }
995  default:
996  return blackboxDefaultOp2(op, result, b1, b2);
997  }
998 }
999 
1000 static BOOLEAN box_OpM(int op, leftv result, leftv args)
1001 {
1002  leftv a = args;
1003  switch(op)
1004  {
1005  case INTERSECT_CMD:
1006  {
1007  if (args->Typ() != boxID)
1008  {
1009  WerrorS("can only intersect boxes");
1010  return TRUE;
1011  }
1012  box *B = (box*) args->Data();
1013  int i, n = B->R->N;
1014  number lowerb[n], upperb[n];
1015 
1016  // do not copy, use same pointers, copy at the end
1017  for (i = 0; i < n; i++)
1018  {
1019  lowerb[i] = B->intervals[i]->lower;
1020  upperb[i] = B->intervals[i]->upper;
1021  }
1022 
1023  args = args->next;
1024  while(args != NULL)
1025  {
1026  if (args->Typ() != boxID)
1027  {
1028  WerrorS("can only intersect boxes");
1029  return TRUE;
1030  }
1031 
1032  B = (box*) args->Data();
1033  for (i = 0; i < n; i++)
1034  {
1035  if (nGreater(B->intervals[i]->lower, lowerb[i]))
1036  {
1037  lowerb[i] = B->intervals[i]->lower;
1038  }
1039  if (nGreater(upperb[i], B->intervals[i]->upper))
1040  {
1041  upperb[i] = B->intervals[i]->upper;
1042  }
1043 
1044  if (nGreater(lowerb[i], upperb[i]))
1045  {
1046  result->rtyp = INT_CMD;
1047  result->data = (void*) (-1);
1048  a->CleanUp();
1049  return FALSE;
1050  }
1051  }
1052  args = args->next;
1053  }
1054 
1055  // now copy the numbers
1056  box *RES = new box();
1057  for (i = 0; i < n; i++)
1058  {
1059  RES->setInterval(i, new interval(nCopy(lowerb[i]), nCopy(upperb[i])));
1060  }
1061 
1062  result->rtyp = boxID;
1063  result->data = (void*) RES;
1064  a->CleanUp();
1065  return FALSE;
1066  }
1067  default:
1068  return blackboxDefaultOpM(op, result, args);
1069  }
1070 }
1071 
1072 static BOOLEAN box_serialize(blackbox*, void *d, si_link f)
1073 {
1074  /*
1075  * Format: "box" setring intervals[1] .. intervals[N]
1076  */
1077  box *B = (box*) d;
1078  int N = B->R->N, i;
1079  sleftv l, iv;
1080  memset(&l, 0, sizeof(l));
1081  memset(&iv, 0, sizeof(iv));
1082 
1083  l.rtyp = STRING_CMD;
1084  l.data = (void*) "box";
1085 
1086  f->m->Write(f, &l);
1087  f->m->SetRing(f, B->R, TRUE);
1088 
1089  iv.rtyp = intervalID;
1090  for (i = 0; i < N; i++)
1091  {
1092  iv.data = (void*) B->intervals[i];
1093  f->m->Write(f, &iv);
1094  }
1095 
1096  if (currRing != B->R)
1097  f->m->SetRing(f, currRing, FALSE);
1098 
1099  return FALSE;
1100 }
1101 
1102 static BOOLEAN box_deserialize(blackbox**, void **d, si_link f)
1103 {
1104  leftv l;
1105 
1106  // read once to set ring
1107  l = f->m->Read(f);
1108 
1109  ring R = currRing;
1110  int i, N = R->N;
1111  box *B = new box();
1112 
1113  B->setInterval(0, (interval*) l->CopyD());
1114  l->CleanUp();
1115 
1116  for (i = 1; i < N; i++)
1117  {
1118  l = f->m->Read(f);
1119  B->setInterval(i, (interval*) l->CopyD());
1120  l->CleanUp();
1121  }
1122 
1123  *d = (void*) B;
1124  return FALSE;
1125 }
1126 
1128 {
1129  // check for proper types
1130  const short t[] = {3, (short) boxID, INT_CMD, (short) intervalID};
1131  if (!iiCheckTypes(args, t, 1))
1132  {
1133  return TRUE;
1134  }
1135  box *B = (box*) args->Data();
1136  int n = B->R->N,
1137  i = (int)(long) args->next->Data();
1138  interval *I = (interval*) args->next->next->Data();
1139 
1140  if (i < 1 || i > n)
1141  {
1142  WerrorS("boxSet: index out of range");
1143  return TRUE;
1144  }
1145 
1146  box *RES = new box(B);
1147 
1148  RES->setInterval(i-1, new interval(I));
1149  // ensure consistency
1150  RES->intervals[i-1]->setRing(RES->R);
1151 
1152  result->rtyp = boxID;
1153  result->data = (void*) RES;
1154  args->CleanUp();
1155  return FALSE;
1156 }
1157 
1158 /*
1159  * POLY FUNCTIONS
1160  */
1161 
1163 {
1164  const short t[] = {2, POLY_CMD, (short) boxID};
1165  if (!iiCheckTypes(args, t, 1))
1166  {
1167  return TRUE;
1168  }
1169 
1170  poly p = (poly) args->Data();
1171  box *B = (box*) args->next->Data();
1172  int i, pot, n = B->R->N;
1173 
1174  interval *tmp, *tmpPot, *tmpMonom, *RES = new interval();
1175 
1176  while(p != NULL)
1177  {
1178  tmpMonom = new interval(nInit(1));
1179 
1180  for (i = 1; i <= n; i++)
1181  {
1182  pot = pGetExp(p, i);
1183 
1184  tmpPot = intervalPower(B->intervals[i-1], pot);
1185  tmp = intervalMultiply(tmpMonom, tmpPot);
1186 
1187  delete tmpMonom;
1188  delete tmpPot;
1189 
1190  tmpMonom = tmp;
1191  }
1192 
1193  tmp = intervalScalarMultiply(p->coef, tmpMonom);
1194  delete tmpMonom;
1195  tmpMonom = tmp;
1196 
1197  tmp = intervalAdd(RES, tmpMonom);
1198  delete RES;
1199  delete tmpMonom;
1200 
1201  RES = tmp;
1202 
1203  p = p->next;
1204  }
1205 
1206  if (result->Data() != NULL)
1207  {
1208  delete (box*) result->Data();
1209  }
1210 
1211  result->rtyp = intervalID;
1212  result->data = (void*) RES;
1213  args->CleanUp();
1214  return FALSE;
1215 }
1216 
1217 /*
1218  * INIT MODULE
1219  */
1220 
1221 extern "C" int SI_MOD_INIT(interval)(SModulFunctions* psModulFunctions)
1222 {
1223  blackbox *b_iv = (blackbox*) omAlloc0(sizeof(blackbox)),
1224  *b_bx = (blackbox*) omAlloc0(sizeof(blackbox));
1225 
1226  b_iv->blackbox_Init = interval_Init;
1227  b_iv->blackbox_Copy = interval_Copy;
1228  b_iv->blackbox_destroy = interval_Destroy;
1229  b_iv->blackbox_String = interval_String;
1230  b_iv->blackbox_Assign = interval_Assign;
1231  b_iv->blackbox_Op2 = interval_Op2;
1232  b_iv->blackbox_serialize = interval_serialize;
1233  b_iv->blackbox_deserialize = interval_deserialize;
1234 
1235  intervalID = setBlackboxStuff(b_iv, "interval");
1236 
1237  b_bx->blackbox_Init = box_Init;
1238  b_bx->blackbox_Copy = box_Copy;
1239  b_bx->blackbox_destroy = box_Destroy;
1240  b_bx->blackbox_String = box_String;
1241  b_bx->blackbox_Assign = box_Assign;
1242  b_bx->blackbox_Op2 = box_Op2;
1243  b_bx->blackbox_OpM = box_OpM;
1244  b_bx->blackbox_serialize = box_serialize;
1245  b_bx->blackbox_deserialize = box_deserialize;
1246 
1247  boxID = setBlackboxStuff(b_bx, "box");
1248 
1249  // add additional functions
1250  psModulFunctions->iiAddCproc("rootisolation.lib", "length", FALSE, length);
1251  psModulFunctions->iiAddCproc("rootisolation.lib", "boxSet", FALSE, boxSet);
1252  psModulFunctions->iiAddCproc("rootisolation.lib", "evalPolyAtBox", FALSE,
1253  evalPolyAtBox);
1254 
1255  return MAX_TOK;
1256 }
1257 // vim: spell spelllang=en
static FORCE_INLINE number n_Sub(number a, number b, const coeffs r)
return the difference of &#39;a&#39; and &#39;b&#39;, i.e., a-b
Definition: coeffs.h:669
box & setInterval(int, interval *)
Definition: interval.cc:122
static bool intervalContainsZero(interval *I)
Definition: interval.cc:368
static FORCE_INLINE BOOLEAN n_Greater(number a, number b, const coeffs r)
ordered fields: TRUE iff &#39;a&#39; is larger than &#39;b&#39;; in Z/pZ: TRUE iff la > lb, where la and lb are the l...
Definition: coeffs.h:511
static BOOLEAN interval_serialize(blackbox *, void *d, si_link f)
Definition: interval.cc:726
~interval()
Definition: interval.cc:47
sleftv * m
Definition: lists.h:46
Class used for (list of) interpreter objects.
Definition: subexpr.h:82
Definition: tok.h:96
number upper
Definition: interval.h:9
Definition: lists.h:23
static BOOLEAN evalPolyAtBox(leftv result, leftv args)
Definition: interval.cc:1162
#define FALSE
Definition: auxiliary.h:94
number lower
Definition: interval.h:8
interval(const ring r=currRing)
Definition: interval.cc:14
Definition: tok.h:216
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition: coeffs.h:538
ring R
Definition: interval.h:24
static void interval_Destroy(blackbox *, void *d)
Definition: interval.cc:177
#define TRUE
Definition: auxiliary.h:98
box()
Definition: interval.cc:80
static FORCE_INLINE void n_Normalize(number &n, const coeffs r)
inplace-normalization of n; produces some canonical representation of n;
Definition: coeffs.h:578
void WerrorS(const char *s)
Definition: feFopen.cc:24
char * StringEndS()
Definition: reporter.cc:151
BOOLEAN blackboxDefaultOp2(int, leftv, leftv, leftv)
default procedure blackboxDefaultOp2, to be called as "default:" branch
Definition: blackbox.cc:81
static interval * intervalScalarMultiply(number a, interval *I)
Definition: interval.cc:280
int Typ()
Definition: subexpr.cc:1033
Definition: idrec.h:34
#define IDHDL
Definition: tok.h:31
static BOOLEAN box_Op2(int op, leftv result, leftv b1, leftv b2)
Definition: interval.cc:890
void * data
Definition: subexpr.h:88
static char * interval_String(blackbox *, void *d)
Definition: interval.cc:149
static FORCE_INLINE number n_Mult(number a, number b, const coeffs r)
return the product of &#39;a&#39; and &#39;b&#39;, i.e., a*b
Definition: coeffs.h:636
#define M
Definition: sirandom.c:24
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
CanonicalForm b
Definition: cfModGcd.cc:4044
~box()
Definition: interval.cc:110
CanonicalForm res
Definition: facAbsFact.cc:64
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:48
static interval * intervalSubtract(interval *I, interval *J)
Definition: interval.cc:349
#define nGreaterZero(n)
Definition: numbers.h:27
#define omFree(addr)
Definition: omAllocDecl.h:261
static BOOLEAN box_deserialize(blackbox **, void **d, si_link f)
Definition: interval.cc:1102
#define assume(x)
Definition: mod2.h:390
static FORCE_INLINE number n_Add(number a, number b, const coeffs r)
return the sum of &#39;a&#39; and &#39;b&#39;, i.e., a+b
Definition: coeffs.h:656
ring R
Definition: interval.h:10
void StringSetS(const char *st)
Definition: reporter.cc:128
static BOOLEAN boxSet(leftv result, leftv args)
Definition: interval.cc:1127
interval ** intervals
Definition: interval.h:23
void StringAppendS(const char *st)
Definition: reporter.cc:107
static BOOLEAN box_Assign(leftv result, leftv args)
Definition: interval.cc:821
number(* nMapFunc)(number a, const coeffs src, const coeffs dst)
maps "a", which lives in src, into dst
Definition: coeffs.h:73
static char * box_String(blackbox *, void *d)
Definition: interval.cc:798
static FORCE_INLINE void n_Write(number n, const coeffs r, const BOOLEAN bShortOut=TRUE)
Definition: coeffs.h:591
static FORCE_INLINE number n_Invers(number a, const coeffs r)
return the multiplicative inverse of &#39;a&#39;; raise an error if &#39;a&#39; is not invertible ...
Definition: coeffs.h:564
int m
Definition: cfEzgcd.cc:121
FILE * f
Definition: checklibs.c:9
int i
Definition: cfEzgcd.cc:125
static void * interval_Copy(blackbox *, void *d)
Definition: interval.cc:171
static FORCE_INLINE nMapFunc n_SetMap(const coeffs src, const coeffs dst)
set the mapping function pointers for translating numbers from src to dst
Definition: coeffs.h:721
int lSize(lists L)
Definition: lists.cc:25
#define nDelete(n)
Definition: numbers.h:16
leftv next
Definition: subexpr.h:86
#define nInvers(a)
Definition: numbers.h:33
static void * interval_Init(blackbox *)
Definition: interval.cc:143
static void * box_Init(blackbox *)
Definition: interval.cc:782
static FORCE_INLINE void n_Power(number a, int b, number *res, const coeffs r)
fill res with the power a^b
Definition: coeffs.h:632
#define nIsZero(n)
Definition: numbers.h:19
#define NULL
Definition: omList.c:12
static BOOLEAN box_OpM(int op, leftv result, leftv args)
Definition: interval.cc:1000
slists * lists
Definition: mpr_numeric.h:146
static bool intervalEqual(interval *I, interval *J)
Definition: interval.cc:360
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of &#39;n&#39;
Definition: coeffs.h:451
static interval * intervalPower(interval *I, int p)
Definition: interval.cc:378
b *CanonicalForm B
Definition: facBivar.cc:52
BOOLEAN iiCheckTypes(leftv args, const short *type_list, int report)
check a list of arguemys against a given field of types return TRUE if the types match return FALSE (...
Definition: ipshell.cc:6551
static void box_Destroy(blackbox *, void *d)
Definition: interval.cc:792
static int intervalID
Definition: interval.cc:136
static interval * intervalMultiply(interval *I, interval *J)
Definition: interval.cc:301
static BOOLEAN length(leftv result, leftv arg)
Definition: interval.cc:263
BOOLEAN blackboxDefaultOpM(int op, leftv res, leftv args)
default procedure blackboxDefaultOpM, to be called as "default:" branch
Definition: blackbox.cc:91
int rtyp
Definition: subexpr.h:91
#define nCopy(n)
Definition: numbers.h:15
static BOOLEAN box_serialize(blackbox *, void *d, si_link f)
Definition: interval.cc:1072
void CleanUp(ring r=currRing)
Definition: subexpr.cc:348
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
interval & setRing(ring)
Definition: interval.cc:54
void * Data()
Definition: subexpr.cc:1176
static FORCE_INLINE BOOLEAN n_Equal(number a, number b, const coeffs r)
TRUE iff &#39;a&#39; and &#39;b&#39; represent the same number; they may have different representations.
Definition: coeffs.h:460
Definition: tok.h:118
static BOOLEAN interval_Assign(leftv result, leftv args)
Definition: interval.cc:184
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete &#39;p&#39;
Definition: coeffs.h:455
static FORCE_INLINE BOOLEAN n_GreaterZero(number n, const coeffs r)
ordered fields: TRUE iff &#39;n&#39; is positive; in Z/pZ: TRUE iff 0 < m <= roundedBelow(p/2), where m is the long representing n in C: TRUE iff (Im(n) != 0 and Im(n) >= 0) or (Im(n) == 0 and Re(n) >= 0) in K(a)/<p(a)>: TRUE iff (n != 0 and (LC(n) > 0 or deg(n) > 0)) in K(t_1, ..., t_n): TRUE iff (LC(numerator(n) is a constant and > 0) or (LC(numerator(n) is not a constant) in Z/2^kZ: TRUE iff 0 < n <= 2^(k-1) in Z/mZ: TRUE iff the internal mpz is greater than zero in Z: TRUE iff n > 0
Definition: coeffs.h:494
int p
Definition: cfModGcd.cc:4019
static BOOLEAN interval_Op2(int op, leftv result, leftv i1, leftv i2)
Definition: interval.cc:438
#define IDDATA(a)
Definition: ipid.h:121
#define nInit(i)
Definition: numbers.h:24
int setBlackboxStuff(blackbox *bb, const char *n)
define a new type
Definition: blackbox.cc:126
int BOOLEAN
Definition: auxiliary.h:85
static void * box_Copy(blackbox *, void *d)
Definition: interval.cc:787
static int boxID
Definition: interval.cc:137
void Werror(const char *fmt,...)
Definition: reporter.cc:189
void * CopyD(int t)
Definition: subexpr.cc:739
#define omAlloc0(size)
Definition: omAllocDecl.h:211
int l
Definition: cfEzgcd.cc:93
return result
Definition: facAbsBiFact.cc:76
Definition: interval.h:21
#define nGreater(a, b)
Definition: numbers.h:28
blackbox * getBlackboxStuff(const int t)
return the structure to the type given by t
Definition: blackbox.cc:16
static BOOLEAN interval_deserialize(blackbox **, void **d, si_link f)
Definition: interval.cc:762
#define omStrDup(s)
Definition: omAllocDecl.h:263
static interval * intervalAdd(interval *I, interval *J)
Definition: interval.cc:338