source: git/factory/cfNewtonPolygon.cc @ 72ebdb

spielwiese
Last change on this file since 72ebdb was e4fe2b, checked in by Oleksandr Motsak <motsak@…>, 13 years ago
FIX: Fixed huge BUG in cf_gmp.h CHG: starting to cleanup factory
  • 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  if (F.inCoeffDomain())
177  {
178    int* result= new int [1];
179    result [0]= 0;
180    sizeOfOutput= 1;
181    return result;
182  }
183  sizeOfOutput= size (F);
184  int* result= new int [sizeOfOutput];
185  int j= 0;
186  for (CFIterator i= F; i.hasTerms(); i++, j++)
187    result [j]= i.exp();
188  return result;
189}
190
191//get points in Z^2 whose convex hull is the Newton polygon
192int ** getPoints (const CanonicalForm& F, int& n)
193{
194  n= size (F);
195  int ** points= new int* [n];
196  for (int i= 0; i < n; i++)
197    points [i]= new int [2];
198
199  int j= 0;
200  int * buf;
201  int bufSize;
202  if (F.isUnivariate() && F.level() == 1)
203  {
204    for (CFIterator i= F; i.hasTerms(); i++, j++)
205    {
206      points [j] [0]= i.exp();
207      points [j] [1]= 0;
208    }
209    return points;
210  }
211  for (CFIterator i= F; i.hasTerms(); i++)
212  {
213    buf= getDegrees (i.coeff(), bufSize);
214    for (int k= 0; k < bufSize; k++, j++)
215    {
216      points [j] [0]= i.exp();
217      points [j] [1]= buf [k];
218    }
219    delete [] buf;
220  }
221  return points;
222}
223
224// assumes a bivariate poly as input
225int ** newtonPolygon (const CanonicalForm& F, int& sizeOfNewtonPoly)
226{
227  int sizeF= size (F);
228  int ** points= new int* [sizeF];
229  for (int i= 0; i < sizeF; i++)
230    points [i]= new int [2];
231  int j= 0;
232  int * buf;
233  int bufSize;
234  for (CFIterator i= F; i.hasTerms(); i++)
235  {
236    buf= getDegrees (i.coeff(), bufSize);
237    for (int k= 0; k < bufSize; k++, j++)
238    {
239      points [j] [0]= i.exp();
240      points [j] [1]= buf [k];
241    }
242    delete [] buf;
243  }
244
245  int n= polygon (points, sizeF);
246
247  int ** result= new int* [n];
248  for (int i= 0; i < n; i++)
249  {
250    result [i]= new int [2];
251    result [i] [0]= points [i] [0];
252    result [i] [1]= points [i] [1];
253  }
254
255  sizeOfNewtonPoly= n;
256  for (int i= 0; i < sizeF; i++)
257    delete [] points[i];
258  delete [] points;
259
260  return result;
261}
262
263// assumes first sizePoints entries of points form a Newton polygon
264bool isInPolygon (int ** points, int sizePoints, int* point)
265{
266  int ** buf= new int* [sizePoints + 1];
267  for (int i= 0; i < sizePoints; i++)
268  {
269    buf [i]= new int [2];
270    buf [i] [0]= points [i] [0];
271    buf [i] [1]= points [i] [1];
272  }
273  buf [sizePoints]= new int [2];
274  buf [sizePoints] [0]= point [0];
275  buf [sizePoints] [1]= point [1];
276  int sizeBuf= sizePoints + 1;
277
278  swap (buf, 0, smallestPointIndex (buf, sizeBuf));
279  int * minusPoint= new int [2];
280  minusPoint [0]= buf[0] [0];
281  minusPoint [1]= buf[0] [1];
282  translate (buf, minusPoint, sizeBuf);
283  sort (buf, sizeBuf);
284  minusPoint[0]= - minusPoint[0];
285  minusPoint[1]= - minusPoint[1];
286  translate (buf, minusPoint, sizeBuf); //reverse translation
287  delete [] minusPoint;
288
289  if (buf [0] [0] == point [0] && buf [0] [1] == point [1])
290  {
291    for (int i= 0; i < sizeBuf; i++)
292      delete [] buf[i];
293    delete [] buf;
294    return false;
295  }
296
297  for (int i= 1; i < sizeBuf-1; i++)
298  {
299    if (buf [i] [0] == point [0] && buf [i] [1] == point [1])
300    {
301      bool result= !isConvex (buf, i);
302      for (int i= 0; i < sizeBuf; i++)
303        delete [] buf [i];
304      delete [] buf;
305      return result;
306    }
307  }
308
309  if (buf [sizeBuf - 1] [0] == point [0] && buf [sizeBuf-1] [1] == point [1])
310  {
311    buf [1] [0]= point [0];
312    buf [1] [1]= point [1];
313    buf [2] [0]= buf [0] [0];
314    buf [2] [1]= buf [0] [1];
315    buf [0] [0]= buf [sizeBuf-2] [0];
316    buf [0] [1]= buf [sizeBuf-2] [1];
317    bool result= !isConvex (buf, 1);
318    for (int i= 0; i < sizeBuf; i++)
319      delete [] buf [i];
320    delete [] buf;
321    return result;
322  }
323  for (int i= 0; i < sizeBuf; i++)
324    delete [] buf [i];
325  delete [] buf;
326
327  return false;
328}
329
330void lambda (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 lambdaInverse (int** points, int sizePoints)
337{
338  for (int i= 0; i < sizePoints; i++)
339    points [i] [1]= points [i] [1] + points [i] [0];
340}
341
342void tau (int** points, int sizePoints, int k)
343{
344  for (int i= 0; i < sizePoints; i++)
345    points [i] [1]= points [i] [1] + k;
346}
347
348void mu (int** points, int sizePoints)
349{
350  int tmp;
351  for (int i= 0; i < sizePoints; i++)
352  {
353    tmp= points [i] [0];
354    points [i] [0]= points [i] [1];
355    points [i] [1]= tmp;
356  }
357}
358
359void getMaxMin (int** points, int sizePoints, int& minDiff, int& minSum,
360                int& maxDiff, int& maxSum, int& maxX, int& maxY
361               )
362{
363  minDiff= points [0] [1] - points [0] [0];
364  minSum= points [0] [1] + points [0] [0];
365  maxDiff= points [0] [1] - points [0] [0];
366  maxSum= points [0] [1] + points [0] [0];
367  maxX= points [0] [1];
368  maxY= points [0] [0];
369  int diff, sum;
370  for (int i= 1; i < sizePoints; i++)
371  {
372    diff= points [i] [1] - points [i] [0];
373    sum= points [i] [1] + points [i] [0];
374    minDiff= tmin (minDiff, diff);
375    minSum= tmin (minSum, sum);
376    maxDiff= tmax (maxDiff, diff);
377    maxSum= tmax (maxSum, sum);
378    maxX= tmax (maxX, points [i] [1]);
379    maxY= tmax (maxY, points [i] [0]);
380  }
381}
382
383#ifdef HAVE_NTL
384void convexDense(int** points, int sizePoints, mat_ZZ& M, vec_ZZ& A)
385{
386  if (sizePoints < 3)
387  {
388    if (sizePoints == 2)
389    {
390      int maxX= (points [1] [1] < points [0] [1])?points [0] [1]:points [1] [1];
391      int maxY= (points [1] [0] < points [0] [0])?points [0] [0]:points [1] [0];
392      long u,v,g;
393      XGCD (g, u, v, maxX, maxY);
394      M.SetDims (2,2);
395      A.SetLength (2);
396      if (points [0] [1] != points [0] [0] && points [1] [0] != points [1] [1])
397      {
398        M (1,1)= -u;
399        M (1,2)= v;
400        M (2,1)= maxY/g;
401        M (2,2)= maxX/g;
402        A (1)= u*maxX;
403        A (2)= -(maxY/g)*maxX;
404      }
405      else
406      {
407        M (1,1)= u;
408        M (1,2)= v;
409        M (2,1)= -maxY/g;
410        M (2,2)= maxX/g;
411        A (1)= to_ZZ (0);
412        A (2)= to_ZZ (0);
413      }
414    }
415    else if (sizePoints == 1)
416    {
417      ident (M, 2);
418      A.SetLength (2);
419      A (1)= to_ZZ (0);
420      A (2)= to_ZZ (0);
421    }
422    return;
423  }
424  A.SetLength (2);
425  A (1)= to_ZZ (0);
426  A (2)= to_ZZ (0);
427  ident (M, 2);
428  mat_ZZ Mu;
429  Mu.SetDims (2, 2);
430  Mu (2,1)= to_ZZ (1);
431  Mu (1,2)= to_ZZ (1);
432  Mu (1,1)= to_ZZ (0);
433  Mu (2,2)= to_ZZ (0);
434  mat_ZZ Lambda;
435  Lambda.SetDims (2, 2);
436  Lambda (1,1)= to_ZZ (1);
437  Lambda (1,2)= to_ZZ (-1);
438  Lambda (2,2)= to_ZZ (1);
439  Lambda (2,1)= to_ZZ (0);
440  mat_ZZ InverseLambda;
441  InverseLambda.SetDims (2,2);
442  InverseLambda (1,1)= to_ZZ (1);
443  InverseLambda (1,2)= to_ZZ (1);
444  InverseLambda (2,2)= to_ZZ (1);
445  InverseLambda (2,1)= to_ZZ (0);
446  ZZ tmp;
447  int minDiff, minSum, maxDiff, maxSum, maxX, maxY, b, d, f, h;
448  getMaxMin (points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY);
449  do
450  {
451    if (maxX < maxY)
452    {
453      mu (points, sizePoints);
454      M= Mu*M;
455      tmp= A (1);
456      A (1)= A (2);
457      A (2)= tmp;
458    }
459    getMaxMin (points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY);
460    b= maxX - maxDiff;
461    d= maxX + maxY - maxSum;
462    f= maxY + minDiff;
463    h= minSum;
464    if (b + f > maxY)
465    {
466      lambda (points, sizePoints);
467      tau (points, sizePoints, maxY - f);
468      M= Lambda*M;
469      A [0] += (maxY-f);
470      maxX= maxX + maxY - b - f;
471    }
472    else if (d + h > maxY)
473    {
474      lambdaInverse (points, sizePoints);
475      tau (points, sizePoints, -h);
476      M= InverseLambda*M;
477      A [0] += (-h);
478      maxX= maxX + maxY - d - h;
479    }
480    else
481      return;
482  } while (1);
483}
484
485CanonicalForm
486compress (const CanonicalForm& F, mat_ZZ& M, vec_ZZ& A)
487{
488  int n;
489  int ** newtonPolyg= newtonPolygon (F, n);
490  convexDense (newtonPolyg, n, M, A);
491  CanonicalForm result= 0;
492  ZZ expX, expY;
493  Variable x= Variable (1);
494  Variable y= Variable (2);
495
496  ZZ minExpX, minExpY;
497
498  int k= 0;
499  Variable alpha;
500  for (CFIterator i= F; i.hasTerms(); i++)
501  {
502    if (i.coeff().inCoeffDomain() && hasFirstAlgVar (i.coeff(), alpha))
503    {
504      expX= i.exp()*M (1,2) + A (1);
505      expY= i.exp()*M (2,2) + A (2);
506      if (k == 0)
507      {
508        minExpY= expY;
509        minExpX= expX;
510        k= 1;
511      }
512      else
513      {
514        minExpY= (minExpY > expY) ? expY : minExpY;
515        minExpX= (minExpX > expX) ? expX : minExpX;
516      }
517      result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
518      continue;
519    }
520    CFIterator j= i.coeff();
521    if (k == 0)
522    {
523      expX= j.exp()*M (1,1) + i.exp()*M (1,2) + A (1);
524      expY= j.exp()*M (2,1) + i.exp()*M (2,2) + A (2);
525      minExpX= expX;
526      minExpY= expY;
527      result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
528      j++;
529      k= 1;
530    }
531
532    for (; j.hasTerms(); j++)
533    {
534      expX= j.exp()*M (1,1) + i.exp()*M (1,2) + A (1);
535      expY= j.exp()*M (2,1) + i.exp()*M (2,2) + A (2);
536      result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
537      minExpY= (minExpY > expY) ? expY : minExpY;
538      minExpX= (minExpX > expX) ? expX : minExpX;
539    }
540  }
541
542  if (to_long (minExpX) < 0)
543  {
544    result *= power (x,-to_long(minExpX));
545    result /= CanonicalForm (x, 0);
546  }
547  else
548    result /= power (x,to_long(minExpX));
549
550  if (to_long (minExpY) < 0)
551  {
552    result *= power (y,-to_long(minExpY));
553    result /= CanonicalForm (y, 0);
554  }
555  else
556    result /= power (y,to_long(minExpY));
557
558  CanonicalForm tmp= LC (result);
559  if (tmp.inPolyDomain() && degree (tmp) <= 0)
560  {
561    int d= degree (result);
562    Variable x= result.mvar();
563    result -= tmp*power (x, d);
564    result += Lc (tmp)*power (x, d);
565  }
566
567  for (int i= 0; i < n; i++)
568    delete [] newtonPolyg [i];
569  delete [] newtonPolyg;
570
571  M= inv (M);
572  return result;
573}
574
575CanonicalForm
576decompress (const CanonicalForm& F, const mat_ZZ& inverseM, const vec_ZZ& A)
577{
578  CanonicalForm result= 0;
579  ZZ expX, expY;
580  Variable x= Variable (1);
581  Variable y= Variable (2);
582  ZZ minExpX, minExpY;
583  if (F.isUnivariate() && F.level() == 1)
584  {
585    CFIterator i= F;
586    expX= (i.exp() - A (1))*inverseM (1,1) + (-A (2))*inverseM (1,2);
587    expY= (i.exp() - A (1))*inverseM (2,1) + (-A (2))*inverseM (2,2);
588    minExpX= expX;
589    minExpY= expY;
590    result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
591    i++;
592    for (; i.hasTerms(); i++)
593    {
594      expX= (i.exp() - A (1))*inverseM (1,1) + (-A (2))*inverseM (1,2);
595      expY= (i.exp() - A (1))*inverseM (2,1) + (-A (2))*inverseM (2,2);
596      result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
597      minExpY= (minExpY > expY) ? expY : minExpY;
598      minExpX= (minExpX > expX) ? expX : minExpX;
599    }
600
601    if (to_long (minExpX) < 0)
602    {
603      result *= power (x,-to_long(minExpX));
604      result /= CanonicalForm (x, 0);
605    }
606    else
607      result /= power (x,to_long(minExpX));
608
609    if (to_long (minExpY) < 0)
610    {
611      result *= power (y,-to_long(minExpY));
612      result /= CanonicalForm (y, 0);
613    }
614    else
615      result /= power (y,to_long(minExpY));
616
617    return result/ Lc (result); //normalize
618  }
619
620  int k= 0;
621  Variable alpha;
622  for (CFIterator i= F; i.hasTerms(); i++)
623  {
624    if (i.coeff().inCoeffDomain() && hasFirstAlgVar (i.coeff(), alpha))
625    {
626      expX= -A(1)*inverseM (1,1) + (i.exp() - A (2))*inverseM (1,2);
627      expY= -A(1)*inverseM (2,1) + (i.exp() - A (2))*inverseM (2,2);
628      if (k == 0)
629      {
630        minExpY= expY;
631        minExpX= expX;
632        k= 1;
633      }
634      else
635      {
636        minExpY= (minExpY > expY) ? expY : minExpY;
637        minExpX= (minExpX > expX) ? expX : minExpX;
638      }
639      result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
640      continue;
641    }
642    CFIterator j= i.coeff();
643    if (k == 0)
644    {
645      expX= (j.exp() - A (1))*inverseM (1,1) + (i.exp() - A (2))*inverseM (1,2);
646      expY= (j.exp() - A (1))*inverseM (2,1) + (i.exp() - A (2))*inverseM (2,2);
647      minExpX= expX;
648      minExpY= expY;
649      result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
650      j++;
651      k= 1;
652    }
653
654    for (; j.hasTerms(); j++)
655    {
656      expX= (j.exp() - A (1))*inverseM (1,1) + (i.exp() - A (2))*inverseM (1,2);
657      expY= (j.exp() - A (1))*inverseM (2,1) + (i.exp() - A (2))*inverseM (2,2);
658      result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
659      minExpY= (minExpY > expY) ? expY : minExpY;
660      minExpX= (minExpX > expX) ? expX : minExpX;
661    }
662  }
663
664  if (to_long (minExpX) < 0)
665  {
666    result *= power (x,-to_long(minExpX));
667    result /= CanonicalForm (x, 0);
668  }
669  else
670    result /= power (x,to_long(minExpX));
671
672  if (to_long (minExpY) < 0)
673  {
674    result *= power (y,-to_long(minExpY));
675    result /= CanonicalForm (y, 0);
676  }
677  else
678    result /= power (y,to_long(minExpY));
679
680  return result/Lc (result); //normalize
681}
682#endif
Note: See TracBrowser for help on using the repository browser.