source: git/factory/cfNewtonPolygon.cc @ 1e4b53

spielwiese
Last change on this file since 1e4b53 was e243418, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: use convex dense compression in gcd
  • Property mode set to 100644
File size: 20.0 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 || i == 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
224int **
225merge (int ** points1, int sizePoints1, int ** points2, int sizePoints2,
226       int& sizeResult)
227{
228  int i, j;
229  sizeResult= sizePoints1+sizePoints2;
230  for (i= 0; i < sizePoints1; i++)
231  {
232    for (j= 0; j < sizePoints2; j++)
233    {
234      if (points1[i][0] != points2[j][0])
235        continue;
236      else
237      {
238        if (points1[i][1] != points2[j][1])
239          continue;
240        else
241        {
242          points2[j][0]= -1;
243          points2[j][1]= -1;
244          sizeResult--;
245        }
246      }
247    }
248  }
249  if (sizeResult == 0)
250    return points1;
251
252  int ** result= new int *[sizeResult];
253  for (i= 0; i < sizeResult; i++)
254    result [i]= new int [2];
255
256  int k= 0;
257  for (i= 0; i < sizePoints1; i++, k++)
258  {
259    result[k][0]= points1[i][0];
260    result[k][1]= points1[i][1];
261  }
262  for (i= 0; i < sizePoints2; i++)
263  {
264    if (points2[i][0] < 0)
265      continue;
266    else
267    {
268      result[k][0]= points2[i][0];
269      result[k][1]= points2[i][1];
270      k++;
271    }
272  }
273  return result;
274}
275
276// assumes a bivariate poly as input
277int ** newtonPolygon (const CanonicalForm& F, int& sizeOfNewtonPoly)
278{
279  int sizeF= size (F);
280  int ** points= new int* [sizeF];
281  for (int i= 0; i < sizeF; i++)
282    points [i]= new int [2];
283  int j= 0;
284  int * buf;
285  int bufSize;
286  for (CFIterator i= F; i.hasTerms(); i++)
287  {
288    buf= getDegrees (i.coeff(), bufSize);
289    for (int k= 0; k < bufSize; k++, j++)
290    {
291      points [j] [0]= i.exp();
292      points [j] [1]= buf [k];
293    }
294    delete [] buf;
295  }
296
297  int n= polygon (points, sizeF);
298
299  int ** result= new int* [n];
300  for (int i= 0; i < n; i++)
301  {
302    result [i]= new int [2];
303    result [i] [0]= points [i] [0];
304    result [i] [1]= points [i] [1];
305  }
306
307  sizeOfNewtonPoly= n;
308  for (int i= 0; i < sizeF; i++)
309    delete [] points[i];
310  delete [] points;
311
312  return result;
313}
314
315// assumes a bivariate polys as input
316int ** newtonPolygon (const CanonicalForm& F, const CanonicalForm& G,
317                      int& sizeOfNewtonPoly)
318{
319  int sizeF= size (F);
320  int ** pointsF= new int* [sizeF];
321  for (int i= 0; i < sizeF; i++)
322    pointsF [i]= new int [2];
323  int j= 0;
324  int * buf;
325  int bufSize;
326  for (CFIterator i= F; i.hasTerms(); i++)
327  {
328    buf= getDegrees (i.coeff(), bufSize);
329    for (int k= 0; k < bufSize; k++, j++)
330    {
331      pointsF [j] [0]= i.exp();
332      pointsF [j] [1]= buf [k];
333    }
334    delete [] buf;
335  }
336
337  int sizeG= size (G);
338  int ** pointsG= new int* [sizeG];
339  for (int i= 0; i < sizeG; i++)
340    pointsG [i]= new int [2];
341  j= 0;
342  for (CFIterator i= G; i.hasTerms(); i++)
343  {
344    buf= getDegrees (i.coeff(), bufSize);
345    for (int k= 0; k < bufSize; k++, j++)
346    {
347      pointsG [j] [0]= i.exp();
348      pointsG [j] [1]= buf [k];
349    }
350    delete [] buf;
351  }
352
353  int sizePoints;
354  int ** points= merge (pointsF, sizeF, pointsG, sizeG, sizePoints);
355
356  int n= polygon (points, sizePoints);
357
358  int ** result= new int* [n];
359  for (int i= 0; i < n; i++)
360  {
361    result [i]= new int [2];
362    result [i] [0]= points [i] [0];
363    result [i] [1]= points [i] [1];
364  }
365
366  sizeOfNewtonPoly= n;
367  for (int i= 0; i < sizeF; i++)
368    delete [] pointsF[i];
369  delete [] pointsF;
370  for (int i= 0; i < sizeG; i++)
371    delete [] pointsG[i];
372  delete [] pointsG;
373
374  return result;
375}
376
377// assumes first sizePoints entries of points form a Newton polygon
378bool isInPolygon (int ** points, int sizePoints, int* point)
379{
380  int ** buf= new int* [sizePoints + 1];
381  for (int i= 0; i < sizePoints; i++)
382  {
383    buf [i]= new int [2];
384    buf [i] [0]= points [i] [0];
385    buf [i] [1]= points [i] [1];
386  }
387  buf [sizePoints]= new int [2];
388  buf [sizePoints] [0]= point [0];
389  buf [sizePoints] [1]= point [1];
390  int sizeBuf= sizePoints + 1;
391
392  swap (buf, 0, smallestPointIndex (buf, sizeBuf));
393  int * minusPoint= new int [2];
394  minusPoint [0]= buf[0] [0];
395  minusPoint [1]= buf[0] [1];
396  translate (buf, minusPoint, sizeBuf);
397  sort (buf, sizeBuf);
398  minusPoint[0]= - minusPoint[0];
399  minusPoint[1]= - minusPoint[1];
400  translate (buf, minusPoint, sizeBuf); //reverse translation
401  delete [] minusPoint;
402
403  if (buf [0] [0] == point [0] && buf [0] [1] == point [1])
404  {
405    for (int i= 0; i < sizeBuf; i++)
406      delete [] buf[i];
407    delete [] buf;
408    return false;
409  }
410
411  for (int i= 1; i < sizeBuf-1; i++)
412  {
413    if (buf [i] [0] == point [0] && buf [i] [1] == point [1])
414    {
415      bool result= !isConvex (buf, i);
416      for (int i= 0; i < sizeBuf; i++)
417        delete [] buf [i];
418      delete [] buf;
419      return result;
420    }
421  }
422
423  if (buf [sizeBuf - 1] [0] == point [0] && buf [sizeBuf-1] [1] == point [1])
424  {
425    buf [1] [0]= point [0];
426    buf [1] [1]= point [1];
427    buf [2] [0]= buf [0] [0];
428    buf [2] [1]= buf [0] [1];
429    buf [0] [0]= buf [sizeBuf-2] [0];
430    buf [0] [1]= buf [sizeBuf-2] [1];
431    bool result= !isConvex (buf, 1);
432    for (int i= 0; i < sizeBuf; i++)
433      delete [] buf [i];
434    delete [] buf;
435    return result;
436  }
437  for (int i= 0; i < sizeBuf; i++)
438    delete [] buf [i];
439  delete [] buf;
440
441  return false;
442}
443
444void lambda (int** points, int sizePoints)
445{
446  for (int i= 0; i < sizePoints; i++)
447    points [i] [1]= points [i] [1] - points [i] [0];
448}
449
450void lambdaInverse (int** points, int sizePoints)
451{
452  for (int i= 0; i < sizePoints; i++)
453    points [i] [1]= points [i] [1] + points [i] [0];
454}
455
456void tau (int** points, int sizePoints, int k)
457{
458  for (int i= 0; i < sizePoints; i++)
459    points [i] [1]= points [i] [1] + k;
460}
461
462void mu (int** points, int sizePoints)
463{
464  int tmp;
465  for (int i= 0; i < sizePoints; i++)
466  {
467    tmp= points [i] [0];
468    points [i] [0]= points [i] [1];
469    points [i] [1]= tmp;
470  }
471}
472
473void getMaxMin (int** points, int sizePoints, int& minDiff, int& minSum,
474                int& maxDiff, int& maxSum, int& maxX, int& maxY
475               )
476{
477  minDiff= points [0] [1] - points [0] [0];
478  minSum= points [0] [1] + points [0] [0];
479  maxDiff= points [0] [1] - points [0] [0];
480  maxSum= points [0] [1] + points [0] [0];
481  maxX= points [0] [1];
482  maxY= points [0] [0];
483  int diff, sum;
484  for (int i= 1; i < sizePoints; i++)
485  {
486    diff= points [i] [1] - points [i] [0];
487    sum= points [i] [1] + points [i] [0];
488    minDiff= tmin (minDiff, diff);
489    minSum= tmin (minSum, sum);
490    maxDiff= tmax (maxDiff, diff);
491    maxSum= tmax (maxSum, sum);
492    maxX= tmax (maxX, points [i] [1]);
493    maxY= tmax (maxY, points [i] [0]);
494  }
495}
496
497#ifdef HAVE_NTL
498void convexDense(int** points, int sizePoints, mat_ZZ& M, vec_ZZ& A)
499{
500  if (sizePoints < 3)
501  {
502    if (sizePoints == 2)
503    {
504      int maxX= (points [1] [1] < points [0] [1])?points [0] [1]:points [1] [1];
505      int maxY= (points [1] [0] < points [0] [0])?points [0] [0]:points [1] [0];
506      long u,v,g;
507      XGCD (g, u, v, maxX, maxY);
508      M.SetDims (2,2);
509      A.SetLength (2);
510      if (points [0] [1] != points [0] [0] && points [1] [0] != points [1] [1])
511      {
512        M (1,1)= -u;
513        M (1,2)= v;
514        M (2,1)= maxY/g;
515        M (2,2)= maxX/g;
516        A (1)= u*maxX;
517        A (2)= -(maxY/g)*maxX;
518      }
519      else
520      {
521        M (1,1)= u;
522        M (1,2)= v;
523        M (2,1)= -maxY/g;
524        M (2,2)= maxX/g;
525        A (1)= to_ZZ (0);
526        A (2)= to_ZZ (0);
527      }
528    }
529    else if (sizePoints == 1)
530    {
531      ident (M, 2);
532      A.SetLength (2);
533      A (1)= to_ZZ (0);
534      A (2)= to_ZZ (0);
535    }
536    return;
537  }
538  A.SetLength (2);
539  A (1)= to_ZZ (0);
540  A (2)= to_ZZ (0);
541  ident (M, 2);
542  mat_ZZ Mu;
543  Mu.SetDims (2, 2);
544  Mu (2,1)= to_ZZ (1);
545  Mu (1,2)= to_ZZ (1);
546  Mu (1,1)= to_ZZ (0);
547  Mu (2,2)= to_ZZ (0);
548  mat_ZZ Lambda;
549  Lambda.SetDims (2, 2);
550  Lambda (1,1)= to_ZZ (1);
551  Lambda (1,2)= to_ZZ (-1);
552  Lambda (2,2)= to_ZZ (1);
553  Lambda (2,1)= to_ZZ (0);
554  mat_ZZ InverseLambda;
555  InverseLambda.SetDims (2,2);
556  InverseLambda (1,1)= to_ZZ (1);
557  InverseLambda (1,2)= to_ZZ (1);
558  InverseLambda (2,2)= to_ZZ (1);
559  InverseLambda (2,1)= to_ZZ (0);
560  ZZ tmp;
561  int minDiff, minSum, maxDiff, maxSum, maxX, maxY, b, d, f, h;
562  getMaxMin (points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY);
563  do
564  {
565    if (maxX < maxY)
566    {
567      mu (points, sizePoints);
568      M= Mu*M;
569      tmp= A (1);
570      A (1)= A (2);
571      A (2)= tmp;
572    }
573    getMaxMin (points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY);
574    b= maxX - maxDiff;
575    d= maxX + maxY - maxSum;
576    f= maxY + minDiff;
577    h= minSum;
578    if (b + f > maxY)
579    {
580      lambda (points, sizePoints);
581      tau (points, sizePoints, maxY - f);
582      M= Lambda*M;
583      A [0] += (maxY-f);
584      maxX= maxX + maxY - b - f;
585    }
586    else if (d + h > maxY)
587    {
588      lambdaInverse (points, sizePoints);
589      tau (points, sizePoints, -h);
590      M= InverseLambda*M;
591      A [0] += (-h);
592      maxX= maxX + maxY - d - h;
593    }
594    else
595      return;
596  } while (1);
597}
598
599CanonicalForm
600compress (const CanonicalForm& F, mat_ZZ& M, vec_ZZ& A, bool computeMA)
601{
602  int n;
603  int ** newtonPolyg;
604  if (computeMA)
605  {
606    newtonPolyg= newtonPolygon (F, n);
607    convexDense (newtonPolyg, n, M, A);
608  }
609  CanonicalForm result= 0;
610  ZZ expX, expY;
611  Variable x= Variable (1);
612  Variable y= Variable (2);
613
614  ZZ minExpX, minExpY;
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= i.exp()*M (1,2) + A (1);
623      expY= i.exp()*M (2,2) + A (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()*M (1,1) + i.exp()*M (1,2) + A (1);
642      expY= j.exp()*M (2,1) + i.exp()*M (2,2) + A (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()*M (1,1) + i.exp()*M (1,2) + A (1);
653      expY= j.exp()*M (2,1) + i.exp()*M (2,2) + A (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  CanonicalForm tmp= LC (result);
677  if (tmp.inPolyDomain() && degree (tmp) <= 0)
678  {
679    int d= degree (result);
680    Variable x= result.mvar();
681    result -= tmp*power (x, d);
682    result += Lc (tmp)*power (x, d);
683  }
684
685  if (computeMA)
686  {
687    for (int i= 0; i < n; i++)
688      delete [] newtonPolyg [i];
689    delete [] newtonPolyg;
690    M= inv (M);
691  }
692
693  return result;
694}
695
696CanonicalForm
697decompress (const CanonicalForm& F, const mat_ZZ& inverseM, const vec_ZZ& A)
698{
699  CanonicalForm result= 0;
700  ZZ expX, expY;
701  Variable x= Variable (1);
702  Variable y= Variable (2);
703  ZZ minExpX, minExpY;
704  if (F.isUnivariate() && F.level() == 1)
705  {
706    CFIterator i= F;
707    expX= (i.exp() - A (1))*inverseM (1,1) + (-A (2))*inverseM (1,2);
708    expY= (i.exp() - A (1))*inverseM (2,1) + (-A (2))*inverseM (2,2);
709    minExpX= expX;
710    minExpY= expY;
711    result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
712    i++;
713    for (; i.hasTerms(); i++)
714    {
715      expX= (i.exp() - A (1))*inverseM (1,1) + (-A (2))*inverseM (1,2);
716      expY= (i.exp() - A (1))*inverseM (2,1) + (-A (2))*inverseM (2,2);
717      result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
718      minExpY= (minExpY > expY) ? expY : minExpY;
719      minExpX= (minExpX > expX) ? expX : minExpX;
720    }
721
722    if (to_long (minExpX) < 0)
723    {
724      result *= power (x,-to_long(minExpX));
725      result /= CanonicalForm (x, 0);
726    }
727    else
728      result /= power (x,to_long(minExpX));
729
730    if (to_long (minExpY) < 0)
731    {
732      result *= power (y,-to_long(minExpY));
733      result /= CanonicalForm (y, 0);
734    }
735    else
736      result /= power (y,to_long(minExpY));
737
738    return result/ Lc (result); //normalize
739  }
740
741  int k= 0;
742  Variable alpha;
743  for (CFIterator i= F; i.hasTerms(); i++)
744  {
745    if (i.coeff().inCoeffDomain() && hasFirstAlgVar (i.coeff(), alpha))
746    {
747      expX= -A(1)*inverseM (1,1) + (i.exp() - A (2))*inverseM (1,2);
748      expY= -A(1)*inverseM (2,1) + (i.exp() - A (2))*inverseM (2,2);
749      if (k == 0)
750      {
751        minExpY= expY;
752        minExpX= expX;
753        k= 1;
754      }
755      else
756      {
757        minExpY= (minExpY > expY) ? expY : minExpY;
758        minExpX= (minExpX > expX) ? expX : minExpX;
759      }
760      result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
761      continue;
762    }
763    CFIterator j= i.coeff();
764    if (k == 0)
765    {
766      expX= (j.exp() - A (1))*inverseM (1,1) + (i.exp() - A (2))*inverseM (1,2);
767      expY= (j.exp() - A (1))*inverseM (2,1) + (i.exp() - A (2))*inverseM (2,2);
768      minExpX= expX;
769      minExpY= expY;
770      result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
771      j++;
772      k= 1;
773    }
774
775    for (; j.hasTerms(); j++)
776    {
777      expX= (j.exp() - A (1))*inverseM (1,1) + (i.exp() - A (2))*inverseM (1,2);
778      expY= (j.exp() - A (1))*inverseM (2,1) + (i.exp() - A (2))*inverseM (2,2);
779      result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
780      minExpY= (minExpY > expY) ? expY : minExpY;
781      minExpX= (minExpX > expX) ? expX : minExpX;
782    }
783  }
784
785  if (to_long (minExpX) < 0)
786  {
787    result *= power (x,-to_long(minExpX));
788    result /= CanonicalForm (x, 0);
789  }
790  else
791    result /= power (x,to_long(minExpX));
792
793  if (to_long (minExpY) < 0)
794  {
795    result *= power (y,-to_long(minExpY));
796    result /= CanonicalForm (y, 0);
797  }
798  else
799    result /= power (y,to_long(minExpY));
800
801  return result/Lc (result); //normalize
802}
803#endif
804
805//assumes the input is a Newton polygon of a bivariate polynomial which is
806//primitive wrt. x and y, i.e. there is at least one point of the polygon lying
807//on the x-axis and one lying on the y-axis
808int* getRightSide (int** polygon, int sizeOfPolygon, int& sizeOfOutput)
809{
810  int maxY= polygon [0][0];
811  int indexY= 0;
812  for (int i= 1; i < sizeOfPolygon; i++)
813  {
814    if (maxY < polygon [i][0])
815    {
816      maxY= polygon [i][0];
817      indexY= i;
818    }
819    else if (maxY == polygon [i][0])
820    {
821      if (polygon [indexY][1] < polygon[i][1])
822        indexY= i;
823    }
824    if (maxY > polygon [i][0])
825      break;
826  }
827
828  int count= -1;
829  for (int i= indexY; i < sizeOfPolygon; i++)
830  {
831    if (polygon[i][0] == 0)
832    {
833      count= i - indexY;
834      break;
835    }
836  }
837
838  int * result;
839  int index= 0;
840  if (count < 0)
841  {
842    result= new int [sizeOfPolygon - indexY];
843    sizeOfOutput= sizeOfPolygon - indexY;
844    count= sizeOfPolygon - indexY - 1;
845    result [0]= polygon[sizeOfPolygon - 1][0] - polygon [0] [0];
846    index= 1;
847  }
848  else
849  {
850    sizeOfOutput= count;
851    result= new int [count];
852  }
853
854  for (int i= indexY + count; i > indexY; i--, index++)
855    result [index]= polygon [i - 1] [0] - polygon [i] [0];
856
857  return result;
858}
Note: See TracBrowser for help on using the repository browser.