source: git/factory/cfNewtonPolygon.cc @ 2072126

fieker-DuValspielwiese
Last change on this file since 2072126 was 2072126, checked in by Hans Schoenemann <hannes@…>, 13 years ago
add missing HAVE_NTL git-svn-id: file:///usr/local/Singular/svn/trunk@14233 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 16.4 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file cfNewtonPolygon.cc
5 *
6 * This file provides functions to compute the Newton polygon of a bivariate
7 * polynomial
8 *
9 * @author Martin Lee
10 *
11 * @internal
12 * @version \$Id$
13 *
14 **/
15/*****************************************************************************/
16
17#include <config.h>
18#include <stdlib.h>
19
20#include "canonicalform.h"
21#include "cf_iter.h"
22#include "cf_algorithm.h"
23#include "cfNewtonPolygon.h"
24#include "templates/ftmpl_functions.h"
25#include "algext.h"
26
27static
28void translate (int** points, int* point, int sizePoints) //make point to 0
29{
30  for (int i= 0; i < sizePoints; i++)
31  {
32    points[i] [0] -= point [0];
33    points[i] [1] -= point [1];
34  }
35}
36
37static
38int smallestPointIndex (int** points, int sizePoints)
39{
40  int min= 0;
41  for (int i= 1; i < sizePoints; i++)
42  {
43    if (points[i][0] < points[min][0] || 
44        (points[i] [0] == points[min] [0] && points[i] [1] < points[min] [1]))
45      min= i;
46  }
47  return min;
48}
49
50static
51void swap (int** points, int i, int j)
52{
53  int* tmp= points[i];
54  points[i]= points[j];
55  points[j]= tmp;
56}
57
58static
59bool isLess (int* point1, int* point2)
60{
61  int area= point1[0]*point2[1]- point1[1]*point2[0];
62  if (area > 0) return true;
63  if (area == 0)
64  {
65    return (abs (point1[0]) + abs (point1[1]) > 
66            abs (point2[0]) + abs (point2[1]));
67  }
68  return false;
69}
70
71static
72void quickSort (int lo, int hi, int** points)
73{
74  int i= lo, j= hi;
75  int* point= new int [2];
76  point [0]= points [(lo+hi)/2] [0];
77  point [1]= points [(lo+hi)/2] [1];
78  while (i <= j)
79  {
80    while (isLess (points [i], point) && i < hi) i++;
81    while (isLess (point, points[j]) && j > lo) j--;
82    if (i <= j)
83    {
84      swap (points, i, j);
85      i++;
86      j--;
87    }
88  }
89  delete [] point;
90  if (lo < j) quickSort (lo, j, points);
91  if (i < hi) quickSort (i, hi, points);
92}
93
94static
95void sort (int** points, int sizePoints)
96{
97  quickSort (1, sizePoints - 1, points);
98}
99
100// check whether p2 is convex
101static
102bool isConvex (int* point1, int* point2, int* point3)
103{
104  int relArea= (point1[0] - point2[0])*(point3[1] - point2[1]) -
105               (point1[1] - point2[1])*(point3[0] - point2[0]);
106  if (relArea < 0)
107    return true;
108  if (relArea == 0)
109  {
110    return !(abs (point1[0] - point3[0]) + abs (point1[1] - point3[1]) >=
111             (abs (point2[0] - point1[0]) + abs (point2[1] - point1[1]) +
112             abs (point2[0] - point3[0]) + abs (point2[1] - point3[1])));
113  }
114  return false;
115}
116
117static
118bool isConvex (int** points, int i)
119{
120  return isConvex (points[i - 1], points [i], points [i + 1]);
121}
122
123int grahamScan (int** points, int sizePoints)
124{
125  swap (points, 0, smallestPointIndex (points, sizePoints));
126  int * minusPoint= new int [2];
127  minusPoint [0]= points[0] [0];
128  minusPoint [1]= points[0] [1];
129  translate (points, minusPoint, sizePoints);
130  sort (points, sizePoints);
131  minusPoint[0]= - minusPoint[0];
132  minusPoint[1]= - minusPoint[1];
133  translate (points, minusPoint, sizePoints); //reverse translation
134  delete [] minusPoint;
135  int i= 3, k= 3;
136  while (k < sizePoints)
137  {
138    swap (points, i, k);
139    while (!isConvex (points, i - 1))
140    {
141      swap (points, i - 1, i);
142      i--;
143    }
144    k++;
145    i++;
146  }
147  if (i + 1 < sizePoints)
148  {
149    int relArea=
150    (points [i-2][0] - points [i-1][0])*(points [0][1] - points [i-1][1])-
151    (points [i-2][1] - points [i-1][1])*(points [0][0] - points [i-1][0]);
152    if (relArea == 0)
153    {
154      if (abs (points [i-2][0] - points [0][0]) +
155          abs (points [i-2][1] - points [0][1]) >=
156          abs (points [i-1][0] - points [i-2][0]) +
157          abs (points [i-1][1] - points [i-2][1]) +
158          abs (points [i-1][0] - points [0][0]) +
159          abs (points [i-1][1] - points [0][1]))
160          i--;
161    }
162  }
163  return i;
164}
165
166//points[i] [0] is x-coordinate, points [i] [1] is y-coordinate
167int polygon (int** points, int sizePoints)
168{
169  if (sizePoints < 3) return sizePoints;
170  return grahamScan (points, sizePoints);
171}
172
173static
174int* getDegrees (const CanonicalForm& F, int& sizeOfOutput)
175{
176  ASSERT (F.isUnivariate(), "univariate poly expected");
177  sizeOfOutput= size (F);
178  int* result= new int [sizeOfOutput];
179  int j= 0;
180  for (CFIterator i= F; i.hasTerms(); i++, j++)
181    result [j]= i.exp();
182  return result;
183}
184
185//get points in Z^2 whose convex hull is the Newton polygon
186int ** getPoints (const CanonicalForm& F, int& n)
187{
188  n= size (F);
189  int ** points= new int* [n];
190  for (int i= 0; i < n; i++)
191    points [i]= new int [2];
192
193  int j= 0;
194  int * buf;
195  int bufSize;
196  if (F.isUnivariate() && F.level() == 1)
197  {
198    for (CFIterator i= F; i.hasTerms(); i++, j++)
199    {
200      points [j] [0]= i.exp();
201      points [j] [1]= 0;
202    }
203    return points;
204  }
205  for (CFIterator i= F; i.hasTerms(); i++)
206  {
207    buf= getDegrees (i.coeff(), bufSize);
208    for (int k= 0; k < bufSize; k++, j++)
209    {
210      points [j] [0]= i.exp();
211      points [j] [1]= buf [k];
212    }
213    delete [] buf;
214  }
215  return points;
216}
217
218// assumes a bivariate poly as input
219int ** newtonPolygon (const CanonicalForm& F, int& sizeOfNewtonPoly)
220{
221  int sizeF= size (F);
222  int ** points= new int* [sizeF];
223  for (int i= 0; i < sizeF; i++)
224    points [i]= new int [2];
225  int j= 0;
226  int * buf;
227  int bufSize;
228  for (CFIterator i= F; i.hasTerms(); i++)
229  {
230    buf= getDegrees (i.coeff(), bufSize);
231    for (int k= 0; k < bufSize; k++, j++)
232    {
233      points [j] [0]= i.exp();
234      points [j] [1]= buf [k];
235    }
236    delete [] buf;
237  }
238
239  int n= polygon (points, sizeF);
240
241  int ** result= new int* [n];
242  for (int i= 0; i < n; i++)
243  {
244    result [i]= new int [2];
245    result [i] [0]= points [i] [0];
246    result [i] [1]= points [i] [1];
247  }
248
249  sizeOfNewtonPoly= n;
250  for (int i= n; i < sizeF; i++)
251    delete [] points[i];
252  delete [] points;
253
254  return result;
255}
256
257// assumes first sizePoints entries of points form a Newton polygon
258bool isInPolygon (int ** points, int sizePoints, int* point)
259{
260  int ** buf= new int* [sizePoints + 1];
261  for (int i= 0; i < sizePoints; i++)
262  {
263    buf [i]= new int [2];
264    buf [i] [0]= points [i] [0];
265    buf [i] [1]= points [i] [1];
266  }
267  buf [sizePoints]= new int [2];
268  buf [sizePoints] [0]= point [0];
269  buf [sizePoints] [1]= point [1];
270  int sizeBuf= sizePoints + 1;
271
272  swap (buf, 0, smallestPointIndex (buf, sizeBuf));
273  int * minusPoint= new int [2];
274  minusPoint [0]= buf[0] [0];
275  minusPoint [1]= buf[0] [1];
276  translate (buf, minusPoint, sizeBuf);
277  sort (buf, sizeBuf);
278  minusPoint[0]= - minusPoint[0];
279  minusPoint[1]= - minusPoint[1];
280  translate (buf, minusPoint, sizeBuf); //reverse translation
281  delete [] minusPoint;
282
283  if (buf [0] [0] == point [0] && buf [0] [1] == point [1])
284  {
285    for (int i= 0; i < sizeBuf; i++)
286      delete [] buf[i];
287    delete [] buf;
288    return false;
289  }
290
291  for (int i= 1; i < sizeBuf-1; i++)
292  {
293    if (buf [i] [0] == point [0] && buf [i] [1] == point [1])
294    {
295      bool result= !isConvex (buf, i);
296      for (int i= 0; i < sizeBuf; i++)
297        delete [] buf [i];
298      delete [] buf;
299      return result;
300    }
301  }
302 
303  if (buf [sizeBuf - 1] [0] == point [0] && buf [sizeBuf-1] [1] == point [1])
304  {
305    buf [1] [0]= point [0];
306    buf [1] [1]= point [1];
307    buf [2] [0]= buf [0] [0];
308    buf [2] [1]= buf [0] [1];
309    buf [0] [0]= buf [sizeBuf-2] [0];
310    buf [0] [1]= buf [sizeBuf-2] [1];
311    bool result= !isConvex (buf, 1);
312    for (int i= 0; i < sizeBuf; i++)
313      delete [] buf [i];
314    delete [] buf;
315    return result;
316  }
317  for (int i= 0; i < sizeBuf; i++)
318    delete [] buf [i];
319  delete [] buf;
320
321  return false;
322}
323
324void lambda (int** points, int sizePoints)
325{
326  for (int i= 0; i < sizePoints; i++)
327    points [i] [1]= points [i] [1] - points [i] [0];
328}
329
330void lambdaInverse (int** points, int sizePoints)
331{
332  for (int i= 0; i < sizePoints; i++)
333    points [i] [1]= points [i] [1] + points [i] [0];
334}
335
336void tau (int** points, int sizePoints, int k)
337{
338  for (int i= 0; i < sizePoints; i++)
339    points [i] [1]= points [i] [1] + k;
340}
341
342void mu (int** points, int sizePoints)
343{
344  int tmp;
345  for (int i= 0; i < sizePoints; i++)
346  {
347    tmp= points [i] [0];
348    points [i] [0]= points [i] [1];
349    points [i] [1]= tmp;
350  }
351}
352
353void getMaxMin (int** points, int sizePoints, int& minDiff, int& minSum,
354                int& maxDiff, int& maxSum, int& maxX, int& maxY
355               )
356{
357  minDiff= points [0] [1] - points [0] [0];
358  minSum= points [0] [1] + points [0] [0];
359  maxDiff= points [0] [1] - points [0] [0];
360  maxSum= points [0] [1] + points [0] [0];
361  maxX= points [0] [1];
362  maxY= points [0] [0];
363  int diff, sum;
364  for (int i= 1; i < sizePoints; i++)
365  {
366    diff= points [i] [1] - points [i] [0];
367    sum= points [i] [1] + points [i] [0];
368    minDiff= tmin (minDiff, diff);
369    minSum= tmin (minSum, sum);
370    maxDiff= tmax (maxDiff, diff);
371    maxSum= tmax (maxSum, sum);
372    maxX= tmax (maxX, points [i] [1]);
373    maxY= tmax (maxY, points [i] [0]);
374  }
375}
376
377#ifdef HAVE_NTL
378void convexDense(int** points, int sizePoints, mat_ZZ& M, vec_ZZ& A)
379{
380  if (sizePoints < 3)
381  {
382    if (sizePoints == 2)
383    {
384      int maxX= (points [1] [1] < points [0] [1])?points [0] [1]:points [1] [1];
385      int maxY= (points [1] [0] < points [0] [0])?points [0] [0]:points [1] [0];
386      long u,v,g;
387      XGCD (g, u, v, maxX, maxY);
388      M.SetDims (2,2);
389      A.SetLength (2);
390      if (points [0] [1] != points [0] [0] && points [1] [0] != points [1] [1])
391      {
392        M (1,1)= -u;
393        M (1,2)= v;
394        M (2,1)= maxY/g;
395        M (2,2)= maxX/g;
396        A (1)= u*maxX;
397        A (2)= -(maxY/g)*maxX;
398      }
399      else
400      {
401        M (1,1)= u;
402        M (1,2)= v;
403        M (2,1)= -maxY/g;
404        M (2,2)= maxX/g;
405        A (1)= to_ZZ (0);
406        A (2)= to_ZZ (0);
407      }
408    }
409    else if (sizePoints == 1)
410    {
411      ident (M, 2);
412      A.SetLength (2);
413      A (1)= to_ZZ (0);
414      A (2)= to_ZZ (0);
415    }
416    return;
417  }
418  A.SetLength (2);
419  A (1)= to_ZZ (0);
420  A (2)= to_ZZ (0);
421  ident (M, 2);
422  mat_ZZ Mu;
423  Mu.SetDims (2, 2);
424  Mu (2,1)= to_ZZ (1);
425  Mu (1,2)= to_ZZ (1);
426  Mu (1,1)= to_ZZ (0);
427  Mu (2,2)= to_ZZ (0);
428  mat_ZZ Lambda;
429  Lambda.SetDims (2, 2);
430  Lambda (1,1)= to_ZZ (1);
431  Lambda (1,2)= to_ZZ (-1);
432  Lambda (2,2)= to_ZZ (1);
433  Lambda (2,1)= to_ZZ (0);
434  mat_ZZ InverseLambda;
435  InverseLambda.SetDims (2,2);
436  InverseLambda (1,1)= to_ZZ (1);
437  InverseLambda (1,2)= to_ZZ (1);
438  InverseLambda (2,2)= to_ZZ (1);
439  InverseLambda (2,1)= to_ZZ (0);
440  ZZ tmp;
441  int minDiff, minSum, maxDiff, maxSum, maxX, maxY, b, d, f, h;
442  getMaxMin (points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY);
443  do
444  {
445    if (maxX < maxY)
446    {
447      mu (points, sizePoints);
448      M= Mu*M;
449      tmp= A (1);
450      A (1)= A (2);
451      A (2)= tmp;
452    }
453    getMaxMin (points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY);
454    b= maxX - maxDiff;
455    d= maxX + maxY - maxSum;
456    f= maxY + minDiff;
457    h= minSum;
458    if (b + f > maxY)
459    {
460      lambda (points, sizePoints);
461      tau (points, sizePoints, maxY - f);
462      M= Lambda*M;
463      A [0] += (maxY-f);
464      maxX= maxX + maxY - b - f;
465    }
466    else if (d + h > maxY)
467    {
468      lambdaInverse (points, sizePoints);
469      tau (points, sizePoints, -h);
470      M= InverseLambda*M;
471      A [0] += (-h);
472      maxX= maxX + maxY - d - h;
473    }
474    else
475      return;
476  } while (1);
477}
478
479CanonicalForm
480compress (const CanonicalForm& F, mat_ZZ& M, vec_ZZ& A)
481{
482  int n;
483  int ** newtonPolyg= newtonPolygon (F, n);
484  convexDense (newtonPolyg, n, M, A);
485  CanonicalForm result= 0;
486  ZZ expX, expY;
487  Variable x= Variable (1);
488  Variable y= Variable (2);
489  int sizeF= size (F);
490  int ** points= getPoints (F, sizeF);
491
492  ZZ minExpX, minExpY;
493
494  int k= 0;
495  Variable alpha;
496  for (CFIterator i= F; i.hasTerms(); i++)
497  {
498    if (i.coeff().inCoeffDomain() && hasFirstAlgVar (i.coeff(), alpha))
499    {
500      expX= i.exp()*M (1,2) + A (1);
501      expY= i.exp()*M (2,2) + A (2);
502      if (k == 0)
503      {
504        minExpY= expY;
505        minExpX= expX;
506        k= 1;
507      }
508      else
509      {
510        minExpY= (minExpY > expY) ? expY : minExpY;
511        minExpX= (minExpX > expX) ? expX : minExpX;
512      }
513      result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
514      continue;
515    }
516    CFIterator j= i.coeff();
517    if (k == 0)
518    {
519      expX= j.exp()*M (1,1) + i.exp()*M (1,2) + A (1);
520      expY= j.exp()*M (2,1) + i.exp()*M (2,2) + A (2);
521      minExpX= expX;
522      minExpY= expY;
523      result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
524      j++;
525      k= 1;
526    }
527
528    for (; j.hasTerms(); j++)
529    {
530      expX= j.exp()*M (1,1) + i.exp()*M (1,2) + A (1);
531      expY= j.exp()*M (2,1) + i.exp()*M (2,2) + A (2);
532      result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
533      minExpY= (minExpY > expY) ? expY : minExpY;
534      minExpX= (minExpX > expX) ? expX : minExpX;
535    }
536  }
537
538  if (to_long (minExpX) < 0)
539  {
540    result *= power (x,-to_long(minExpX));
541    result /= CanonicalForm (x, 0);
542  }
543  else
544    result /= power (x,to_long(minExpX));
545
546  if (to_long (minExpY) < 0)
547  {
548    result *= power (y,-to_long(minExpY));
549    result /= CanonicalForm (y, 0);
550  }
551  else
552    result /= power (y,to_long(minExpY));
553
554  CanonicalForm tmp= LC (result);
555  if (tmp.inPolyDomain() && degree (tmp) <= 0)
556  {
557    int d= degree (result);
558    Variable x= result.mvar();
559    result -= tmp*power (x, d);
560    result += Lc (tmp)*power (x, d);
561  }
562
563  for (int i= 0; i < n; i++)
564    delete [] newtonPolyg [i];
565  delete [] newtonPolyg;
566
567  M= inv (M);
568  return result;
569}
570
571CanonicalForm
572decompress (const CanonicalForm& F, const mat_ZZ& inverseM, const vec_ZZ& A)
573{
574  CanonicalForm result= 0;
575  ZZ expX, expY;
576  Variable x= Variable (1);
577  Variable y= Variable (2);
578  ZZ minExpX, minExpY;
579  if (F.isUnivariate() && F.level() == 1)
580  {
581    CFIterator i= F;
582    expX= (i.exp() - A (1))*inverseM (1,1) + (-A (2))*inverseM (1,2);
583    expY= (i.exp() - A (1))*inverseM (2,1) + (-A (2))*inverseM (2,2);
584    minExpX= expX;
585    minExpY= expY;
586    result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
587    i++;
588    for (; i.hasTerms(); i++)
589    {
590      expX= (i.exp() - A (1))*inverseM (1,1) + (-A (2))*inverseM (1,2);
591      expY= (i.exp() - A (1))*inverseM (2,1) + (-A (2))*inverseM (2,2);
592      result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
593      minExpY= (minExpY > expY) ? expY : minExpY;
594      minExpX= (minExpX > expX) ? expX : minExpX;
595    }
596
597    if (to_long (minExpX) < 0)
598    {
599      result *= power (x,-to_long(minExpX));
600      result /= CanonicalForm (x, 0);
601    }
602    else
603      result /= power (x,to_long(minExpX));
604
605    if (to_long (minExpY) < 0)
606    {
607      result *= power (y,-to_long(minExpY));
608      result /= CanonicalForm (y, 0);
609    }
610    else
611      result /= power (y,to_long(minExpY));
612
613    return result/ Lc (result); //normalize
614  }
615
616  int k= 0;
617  Variable alpha;
618  for (CFIterator i= F; i.hasTerms(); i++)
619  {
620    if (i.coeff().inCoeffDomain() && hasFirstAlgVar (i.coeff(), alpha))
621    {
622      expX= -A(1)*inverseM (1,1) + (i.exp() - A (2))*inverseM (1,2);
623      expY= -A(1)*inverseM (2,1) + (i.exp() - A (2))*inverseM (2,2);
624      if (k == 0)
625      {
626        minExpY= expY;
627        minExpX= expX;
628        k= 1;
629      }
630      else
631      {
632        minExpY= (minExpY > expY) ? expY : minExpY;
633        minExpX= (minExpX > expX) ? expX : minExpX;
634      }
635      result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
636      continue;
637    }
638    CFIterator j= i.coeff();
639    if (k == 0)
640    {
641      expX= (j.exp() - A (1))*inverseM (1,1) + (i.exp() - A (2))*inverseM (1,2);
642      expY= (j.exp() - A (1))*inverseM (2,1) + (i.exp() - A (2))*inverseM (2,2);
643      minExpX= expX;
644      minExpY= expY;
645      result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
646      j++;
647      k= 1;
648    }
649
650    for (; j.hasTerms(); j++)
651    {
652      expX= (j.exp() - A (1))*inverseM (1,1) + (i.exp() - A (2))*inverseM (1,2);
653      expY= (j.exp() - A (1))*inverseM (2,1) + (i.exp() - A (2))*inverseM (2,2);
654      result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
655      minExpY= (minExpY > expY) ? expY : minExpY;
656      minExpX= (minExpX > expX) ? expX : minExpX;
657    }
658  }
659
660  if (to_long (minExpX) < 0)
661  {
662    result *= power (x,-to_long(minExpX));
663    result /= CanonicalForm (x, 0);
664  }
665  else
666    result /= power (x,to_long(minExpX));
667
668  if (to_long (minExpY) < 0)
669  {
670    result *= power (y,-to_long(minExpY));
671    result /= CanonicalForm (y, 0);
672  }
673  else
674    result /= power (y,to_long(minExpY));
675
676  return result/Lc (result); //normalize
677}
678#endif
Note: See TracBrowser for help on using the repository browser.