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