source: git/factory/cfNewtonPolygon.cc @ 7c118d

jengelh-datetimespielwiese
Last change on this file since 7c118d was 7c118d, checked in by Martin Lee <martinlee84@…>, 10 years ago
fix: memory leak
  • Property mode set to 100644
File size: 21.1 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 **/
12/*****************************************************************************/
13
14#include "config.h"
15#include <stdlib.h>
16
17#include "canonicalform.h"
18#include "cf_iter.h"
19#include "cf_algorithm.h"
20#include "cfNewtonPolygon.h"
21#include "templates/ftmpl_functions.h"
22#include "algext.h"
23
24static
25void translate (int** points, int* point, int sizePoints) //make point to 0
26{
27  for (int i= 0; i < sizePoints; i++)
28  {
29    points[i] [0] -= point [0];
30    points[i] [1] -= point [1];
31  }
32}
33
34static
35int smallestPointIndex (int** points, int sizePoints)
36{
37  int min= 0;
38  for (int i= 1; i < sizePoints; i++)
39  {
40    if (points[i][0] < points[min][0] ||
41        (points[i] [0] == points[min] [0] && points[i] [1] < points[min] [1]))
42      min= i;
43  }
44  return min;
45}
46
47static
48void swap (int** points, int i, int j)
49{
50  int* tmp= points[i];
51  points[i]= points[j];
52  points[j]= tmp;
53}
54
55static
56bool isLess (int* point1, int* point2)
57{
58  int area= point1[0]*point2[1]- point1[1]*point2[0];
59  if (area > 0) return true;
60  if (area == 0)
61  {
62    return (abs (point1[0]) + abs (point1[1]) >
63            abs (point2[0]) + abs (point2[1]));
64  }
65  return false;
66}
67
68static
69void quickSort (int lo, int hi, int** points)
70{
71  int i= lo, j= hi;
72  int* point= new int [2];
73  point [0]= points [(lo+hi)/2] [0];
74  point [1]= points [(lo+hi)/2] [1];
75  while (i <= j)
76  {
77    while (isLess (points [i], point) && i < hi) i++;
78    while (isLess (point, points[j]) && j > lo) j--;
79    if (i <= j)
80    {
81      swap (points, i, j);
82      i++;
83      j--;
84    }
85  }
86  delete [] point;
87  if (lo < j) quickSort (lo, j, points);
88  if (i < hi) quickSort (i, hi, points);
89}
90
91static
92void sort (int** points, int sizePoints)
93{
94  quickSort (1, sizePoints - 1, points);
95}
96
97// check whether p2 is convex
98static
99bool isConvex (int* point1, int* point2, int* point3)
100{
101  int relArea= (point1[0] - point2[0])*(point3[1] - point2[1]) -
102               (point1[1] - point2[1])*(point3[0] - point2[0]);
103  if (relArea < 0)
104    return true;
105  if (relArea == 0)
106  {
107    return !(abs (point1[0] - point3[0]) + abs (point1[1] - point3[1]) >=
108             (abs (point2[0] - point1[0]) + abs (point2[1] - point1[1]) +
109             abs (point2[0] - point3[0]) + abs (point2[1] - point3[1])));
110  }
111  return false;
112}
113
114static
115bool isConvex (int** points, int i)
116{
117  return isConvex (points[i - 1], points [i], points [i + 1]);
118}
119
120int grahamScan (int** points, int sizePoints)
121{
122  swap (points, 0, smallestPointIndex (points, sizePoints));
123  int * minusPoint= new int [2];
124  minusPoint [0]= points[0] [0];
125  minusPoint [1]= points[0] [1];
126  translate (points, minusPoint, sizePoints);
127  sort (points, sizePoints);
128  minusPoint[0]= - minusPoint[0];
129  minusPoint[1]= - minusPoint[1];
130  translate (points, minusPoint, sizePoints); //reverse translation
131  delete [] minusPoint;
132  int i= 3, k= 3;
133  while (k < sizePoints)
134  {
135    swap (points, i, k);
136    while (!isConvex (points, i - 1))
137    {
138      swap (points, i - 1, i);
139      i--;
140    }
141    k++;
142    i++;
143  }
144  if (i + 1 <= sizePoints || i == sizePoints)
145  {
146    int relArea=
147    (points [i-2][0] - points [i-1][0])*(points [0][1] - points [i-1][1])-
148    (points [i-2][1] - points [i-1][1])*(points [0][0] - points [i-1][0]);
149    if (relArea == 0)
150    {
151      if (abs (points [i-2][0] - points [0][0]) +
152          abs (points [i-2][1] - points [0][1]) >=
153          abs (points [i-1][0] - points [i-2][0]) +
154          abs (points [i-1][1] - points [i-2][1]) +
155          abs (points [i-1][0] - points [0][0]) +
156          abs (points [i-1][1] - points [0][1]))
157          i--;
158    }
159  }
160  return i;
161}
162
163//points[i] [0] is x-coordinate, points [i] [1] is y-coordinate
164int polygon (int** points, int sizePoints)
165{
166  if (sizePoints < 3) return sizePoints;
167  return grahamScan (points, sizePoints);
168}
169
170static
171int* getDegrees (const CanonicalForm& F, int& sizeOfOutput)
172{
173  if (F.inCoeffDomain())
174  {
175    int* result= new int [1];
176    result [0]= 0;
177    sizeOfOutput= 1;
178    return result;
179  }
180  sizeOfOutput= size (F);
181  int* result= new int [sizeOfOutput];
182  int j= 0;
183  for (CFIterator i= F; i.hasTerms(); i++, j++)
184    result [j]= i.exp();
185  return result;
186}
187
188//get points in Z^2 whose convex hull is the Newton polygon
189int ** getPoints (const CanonicalForm& F, int& n)
190{
191  n= size (F);
192  int ** points= new int* [n];
193  for (int i= 0; i < n; i++)
194    points [i]= new int [2];
195
196  int j= 0;
197  int * buf;
198  int bufSize;
199  if (F.isUnivariate() && F.level() == 1)
200  {
201    for (CFIterator i= F; i.hasTerms(); i++, j++)
202    {
203      points [j] [0]= i.exp();
204      points [j] [1]= 0;
205    }
206    return points;
207  }
208  for (CFIterator i= F; i.hasTerms(); i++)
209  {
210    buf= getDegrees (i.coeff(), bufSize);
211    for (int k= 0; k < bufSize; k++, j++)
212    {
213      points [j] [0]= i.exp();
214      points [j] [1]= buf [k];
215    }
216    delete [] buf;
217  }
218  return points;
219}
220
221int **
222merge (int ** points1, int sizePoints1, int ** points2, int sizePoints2,
223       int& sizeResult)
224{
225  int i, j;
226  sizeResult= sizePoints1+sizePoints2;
227  for (i= 0; i < sizePoints1; i++)
228  {
229    for (j= 0; j < sizePoints2; j++)
230    {
231      if (points1[i][0] != points2[j][0])
232        continue;
233      else
234      {
235        if (points1[i][1] != points2[j][1])
236          continue;
237        else
238        {
239          points2[j][0]= -1;
240          points2[j][1]= -1;
241          sizeResult--;
242        }
243      }
244    }
245  }
246  if (sizeResult == 0)
247    return points1;
248
249  int ** result= new int *[sizeResult];
250  for (i= 0; i < sizeResult; i++)
251    result [i]= new int [2];
252
253  int k= 0;
254  for (i= 0; i < sizePoints1; i++, k++)
255  {
256    result[k][0]= points1[i][0];
257    result[k][1]= points1[i][1];
258  }
259  for (i= 0; i < sizePoints2; i++)
260  {
261    if (points2[i][0] < 0)
262      continue;
263    else
264    {
265      result[k][0]= points2[i][0];
266      result[k][1]= points2[i][1];
267      k++;
268    }
269  }
270  return result;
271}
272
273// assumes a bivariate poly as input
274int ** newtonPolygon (const CanonicalForm& F, int& sizeOfNewtonPoly)
275{
276  int sizeF= size (F);
277  int ** points= new int* [sizeF];
278  for (int i= 0; i < sizeF; i++)
279    points [i]= new int [2];
280  int j= 0;
281  int * buf;
282  int bufSize;
283  for (CFIterator i= F; i.hasTerms(); i++)
284  {
285    buf= getDegrees (i.coeff(), bufSize);
286    for (int k= 0; k < bufSize; k++, j++)
287    {
288      points [j] [0]= i.exp();
289      points [j] [1]= buf [k];
290    }
291    delete [] buf;
292  }
293
294  int n= polygon (points, sizeF);
295
296  int ** result= new int* [n];
297  for (int i= 0; i < n; i++)
298  {
299    result [i]= new int [2];
300    result [i] [0]= points [i] [0];
301    result [i] [1]= points [i] [1];
302  }
303
304  sizeOfNewtonPoly= n;
305  for (int i= 0; i < sizeF; i++)
306    delete [] points[i];
307  delete [] points;
308
309  return result;
310}
311
312// assumes a bivariate polys as input
313int ** newtonPolygon (const CanonicalForm& F, const CanonicalForm& G,
314                      int& sizeOfNewtonPoly)
315{
316  int sizeF= size (F);
317  int ** pointsF= new int* [sizeF];
318  for (int i= 0; i < sizeF; i++)
319    pointsF [i]= new int [2];
320  int j= 0;
321  int * buf;
322  int bufSize;
323  for (CFIterator i= F; i.hasTerms(); i++)
324  {
325    buf= getDegrees (i.coeff(), bufSize);
326    for (int k= 0; k < bufSize; k++, j++)
327    {
328      pointsF [j] [0]= i.exp();
329      pointsF [j] [1]= buf [k];
330    }
331    delete [] buf;
332  }
333
334  int sizeG= size (G);
335  int ** pointsG= new int* [sizeG];
336  for (int i= 0; i < sizeG; i++)
337    pointsG [i]= new int [2];
338  j= 0;
339  for (CFIterator i= G; i.hasTerms(); i++)
340  {
341    buf= getDegrees (i.coeff(), bufSize);
342    for (int k= 0; k < bufSize; k++, j++)
343    {
344      pointsG [j] [0]= i.exp();
345      pointsG [j] [1]= buf [k];
346    }
347    delete [] buf;
348  }
349
350  int sizePoints;
351  int ** points= merge (pointsF, sizeF, pointsG, sizeG, sizePoints);
352
353  int n= polygon (points, sizePoints);
354
355  int ** result= new int* [n];
356  for (int i= 0; i < n; i++)
357  {
358    result [i]= new int [2];
359    result [i] [0]= points [i] [0];
360    result [i] [1]= points [i] [1];
361  }
362
363  sizeOfNewtonPoly= n;
364  for (int i= 0; i < sizeF; i++)
365    delete [] pointsF[i];
366  delete [] pointsF;
367  for (int i= 0; i < sizeG; i++)
368    delete [] pointsG[i];
369  delete [] pointsG;
370
371  return result;
372}
373
374// assumes first sizePoints entries of points form a Newton polygon
375bool isInPolygon (int ** points, int sizePoints, int* point)
376{
377  int ** buf= new int* [sizePoints + 1];
378  for (int i= 0; i < sizePoints; i++)
379  {
380    buf [i]= new int [2];
381    buf [i] [0]= points [i] [0];
382    buf [i] [1]= points [i] [1];
383  }
384  buf [sizePoints]= new int [2];
385  buf [sizePoints] [0]= point [0];
386  buf [sizePoints] [1]= point [1];
387  int sizeBuf= sizePoints + 1;
388
389  swap (buf, 0, smallestPointIndex (buf, sizeBuf));
390  int * minusPoint= new int [2];
391  minusPoint [0]= buf[0] [0];
392  minusPoint [1]= buf[0] [1];
393  translate (buf, minusPoint, sizeBuf);
394  sort (buf, sizeBuf);
395  minusPoint[0]= - minusPoint[0];
396  minusPoint[1]= - minusPoint[1];
397  translate (buf, minusPoint, sizeBuf); //reverse translation
398  delete [] minusPoint;
399
400  if (buf [0] [0] == point [0] && buf [0] [1] == point [1])
401  {
402    for (int i= 0; i < sizeBuf; i++)
403      delete [] buf[i];
404    delete [] buf;
405    return false;
406  }
407
408  for (int i= 1; i < sizeBuf-1; i++)
409  {
410    if (buf [i] [0] == point [0] && buf [i] [1] == point [1])
411    {
412      bool result= !isConvex (buf, i);
413      for (int i= 0; i < sizeBuf; i++)
414        delete [] buf [i];
415      delete [] buf;
416      return result;
417    }
418  }
419
420  if (buf [sizeBuf - 1] [0] == point [0] && buf [sizeBuf-1] [1] == point [1])
421  {
422    buf [1] [0]= point [0];
423    buf [1] [1]= point [1];
424    buf [2] [0]= buf [0] [0];
425    buf [2] [1]= buf [0] [1];
426    buf [0] [0]= buf [sizeBuf-2] [0];
427    buf [0] [1]= buf [sizeBuf-2] [1];
428    bool result= !isConvex (buf, 1);
429    for (int i= 0; i < sizeBuf; i++)
430      delete [] buf [i];
431    delete [] buf;
432    return result;
433  }
434  for (int i= 0; i < sizeBuf; i++)
435    delete [] buf [i];
436  delete [] buf;
437
438  return false;
439}
440
441void lambda (int** points, int sizePoints)
442{
443  for (int i= 0; i < sizePoints; i++)
444    points [i] [1]= points [i] [1] - points [i] [0];
445}
446
447void lambdaInverse (int** points, int sizePoints)
448{
449  for (int i= 0; i < sizePoints; i++)
450    points [i] [1]= points [i] [1] + points [i] [0];
451}
452
453void tau (int** points, int sizePoints, int k)
454{
455  for (int i= 0; i < sizePoints; i++)
456    points [i] [1]= points [i] [1] + k;
457}
458
459void mu (int** points, int sizePoints)
460{
461  int tmp;
462  for (int i= 0; i < sizePoints; i++)
463  {
464    tmp= points [i] [0];
465    points [i] [0]= points [i] [1];
466    points [i] [1]= tmp;
467  }
468}
469
470void getMaxMin (int** points, int sizePoints, int& minDiff, int& minSum,
471                int& maxDiff, int& maxSum, int& maxX, int& maxY
472               )
473{
474  minDiff= points [0] [1] - points [0] [0];
475  minSum= points [0] [1] + points [0] [0];
476  maxDiff= points [0] [1] - points [0] [0];
477  maxSum= points [0] [1] + points [0] [0];
478  maxX= points [0] [1];
479  maxY= points [0] [0];
480  int diff, sum;
481  for (int i= 1; i < sizePoints; i++)
482  {
483    diff= points [i] [1] - points [i] [0];
484    sum= points [i] [1] + points [i] [0];
485    minDiff= tmin (minDiff, diff);
486    minSum= tmin (minSum, sum);
487    maxDiff= tmax (maxDiff, diff);
488    maxSum= tmax (maxSum, sum);
489    maxX= tmax (maxX, points [i] [1]);
490    maxY= tmax (maxY, points [i] [0]);
491  }
492}
493
494#ifdef HAVE_NTL
495void convexDense(int** points, int sizePoints, mat_ZZ& M, vec_ZZ& A)
496{
497  if (sizePoints < 3)
498  {
499    if (sizePoints == 2)
500    {
501      int maxX= (points [1] [1] < points [0] [1])?points [0] [1]:points [1] [1];
502      int maxY= (points [1] [0] < points [0] [0])?points [0] [0]:points [1] [0];
503      long u,v,g;
504      XGCD (g, u, v, maxX, maxY);
505      M.SetDims (2,2);
506      A.SetLength (2);
507      if (points [0] [1] != points [0] [0] && points [1] [0] != points [1] [1])
508      {
509        M (1,1)= -u;
510        M (1,2)= v;
511        M (2,1)= maxY/g;
512        M (2,2)= maxX/g;
513        A (1)= u*maxX;
514        A (2)= -(maxY/g)*maxX;
515      }
516      else
517      {
518        M (1,1)= u;
519        M (1,2)= v;
520        M (2,1)= -maxY/g;
521        M (2,2)= maxX/g;
522        A (1)= to_ZZ (0);
523        A (2)= to_ZZ (0);
524      }
525    }
526    else if (sizePoints == 1)
527    {
528      ident (M, 2);
529      A.SetLength (2);
530      A (1)= to_ZZ (0);
531      A (2)= to_ZZ (0);
532    }
533    return;
534  }
535  A.SetLength (2);
536  A (1)= to_ZZ (0);
537  A (2)= to_ZZ (0);
538  ident (M, 2);
539  mat_ZZ Mu;
540  Mu.SetDims (2, 2);
541  Mu (2,1)= to_ZZ (1);
542  Mu (1,2)= to_ZZ (1);
543  Mu (1,1)= to_ZZ (0);
544  Mu (2,2)= to_ZZ (0);
545  mat_ZZ Lambda;
546  Lambda.SetDims (2, 2);
547  Lambda (1,1)= to_ZZ (1);
548  Lambda (1,2)= to_ZZ (-1);
549  Lambda (2,2)= to_ZZ (1);
550  Lambda (2,1)= to_ZZ (0);
551  mat_ZZ InverseLambda;
552  InverseLambda.SetDims (2,2);
553  InverseLambda (1,1)= to_ZZ (1);
554  InverseLambda (1,2)= to_ZZ (1);
555  InverseLambda (2,2)= to_ZZ (1);
556  InverseLambda (2,1)= to_ZZ (0);
557  ZZ tmp;
558  int minDiff, minSum, maxDiff, maxSum, maxX, maxY, b, d, f, h;
559  getMaxMin (points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY);
560  do
561  {
562    if (maxX < maxY)
563    {
564      mu (points, sizePoints);
565      M= Mu*M;
566      tmp= A (1);
567      A (1)= A (2);
568      A (2)= tmp;
569    }
570    getMaxMin (points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY);
571    b= maxX - maxDiff;
572    d= maxX + maxY - maxSum;
573    f= maxY + minDiff;
574    h= minSum;
575    if (b + f > maxY)
576    {
577      lambda (points, sizePoints);
578      tau (points, sizePoints, maxY - f);
579      M= Lambda*M;
580      A [0] += (maxY-f);
581      maxX= maxX + maxY - b - f;
582    }
583    else if (d + h > maxY)
584    {
585      lambdaInverse (points, sizePoints);
586      tau (points, sizePoints, -h);
587      M= InverseLambda*M;
588      A [0] += (-h);
589      maxX= maxX + maxY - d - h;
590    }
591    else
592      return;
593  } while (1);
594}
595
596CanonicalForm
597compress (const CanonicalForm& F, mat_ZZ& M, vec_ZZ& A, bool computeMA)
598{
599  int n;
600  int ** newtonPolyg= NULL;
601  if (computeMA)
602  {
603    newtonPolyg= newtonPolygon (F, n);
604    convexDense (newtonPolyg, n, M, A);
605  }
606  CanonicalForm result= 0;
607  ZZ expX, expY;
608  Variable x= Variable (1);
609  Variable y= Variable (2);
610
611  ZZ minExpX, minExpY;
612
613  int k= 0;
614  Variable alpha;
615  for (CFIterator i= F; i.hasTerms(); i++)
616  {
617    if (i.coeff().inCoeffDomain() && hasFirstAlgVar (i.coeff(), alpha))
618    {
619      expX= i.exp()*M (1,2) + A (1);
620      expY= i.exp()*M (2,2) + A (2);
621      if (k == 0)
622      {
623        minExpY= expY;
624        minExpX= expX;
625        k= 1;
626      }
627      else
628      {
629        minExpY= (minExpY > expY) ? expY : minExpY;
630        minExpX= (minExpX > expX) ? expX : minExpX;
631      }
632      result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
633      continue;
634    }
635    CFIterator j= i.coeff();
636    if (k == 0)
637    {
638      expX= j.exp()*M (1,1) + i.exp()*M (1,2) + A (1);
639      expY= j.exp()*M (2,1) + i.exp()*M (2,2) + A (2);
640      minExpX= expX;
641      minExpY= expY;
642      result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
643      j++;
644      k= 1;
645    }
646
647    for (; j.hasTerms(); j++)
648    {
649      expX= j.exp()*M (1,1) + i.exp()*M (1,2) + A (1);
650      expY= j.exp()*M (2,1) + i.exp()*M (2,2) + A (2);
651      result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
652      minExpY= (minExpY > expY) ? expY : minExpY;
653      minExpX= (minExpX > expX) ? expX : minExpX;
654    }
655  }
656
657  if (to_long (minExpX) < 0)
658  {
659    result *= power (x,-to_long(minExpX));
660    result /= CanonicalForm (x, 0);
661  }
662  else
663    result /= power (x,to_long(minExpX));
664
665  if (to_long (minExpY) < 0)
666  {
667    result *= power (y,-to_long(minExpY));
668    result /= CanonicalForm (y, 0);
669  }
670  else
671    result /= power (y,to_long(minExpY));
672
673  CanonicalForm tmp= LC (result);
674  if (tmp.inPolyDomain() && degree (tmp) <= 0)
675  {
676    int d= degree (result);
677    Variable x= result.mvar();
678    result -= tmp*power (x, d);
679    result += Lc (tmp)*power (x, d);
680  }
681
682  if (computeMA)
683  {
684    for (int i= 0; i < n; i++)
685      delete [] newtonPolyg [i];
686    delete [] newtonPolyg;
687    M= inv (M);
688  }
689
690  return result;
691}
692
693CanonicalForm
694decompress (const CanonicalForm& F, const mat_ZZ& inverseM, const vec_ZZ& A)
695{
696  CanonicalForm result= 0;
697  ZZ expX, expY;
698  Variable x= Variable (1);
699  Variable y= Variable (2);
700  ZZ minExpX, minExpY;
701  if (F.isUnivariate() && F.level() == 1)
702  {
703    CFIterator i= F;
704    expX= (i.exp() - A (1))*inverseM (1,1) + (-A (2))*inverseM (1,2);
705    expY= (i.exp() - A (1))*inverseM (2,1) + (-A (2))*inverseM (2,2);
706    minExpX= expX;
707    minExpY= expY;
708    result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
709    i++;
710    for (; i.hasTerms(); i++)
711    {
712      expX= (i.exp() - A (1))*inverseM (1,1) + (-A (2))*inverseM (1,2);
713      expY= (i.exp() - A (1))*inverseM (2,1) + (-A (2))*inverseM (2,2);
714      result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
715      minExpY= (minExpY > expY) ? expY : minExpY;
716      minExpX= (minExpX > expX) ? expX : minExpX;
717    }
718
719    if (to_long (minExpX) < 0)
720    {
721      result *= power (x,-to_long(minExpX));
722      result /= CanonicalForm (x, 0);
723    }
724    else
725      result /= power (x,to_long(minExpX));
726
727    if (to_long (minExpY) < 0)
728    {
729      result *= power (y,-to_long(minExpY));
730      result /= CanonicalForm (y, 0);
731    }
732    else
733      result /= power (y,to_long(minExpY));
734
735    return result/ Lc (result); //normalize
736  }
737
738  int k= 0;
739  Variable alpha;
740  for (CFIterator i= F; i.hasTerms(); i++)
741  {
742    if (i.coeff().inCoeffDomain() && hasFirstAlgVar (i.coeff(), alpha))
743    {
744      expX= -A(1)*inverseM (1,1) + (i.exp() - A (2))*inverseM (1,2);
745      expY= -A(1)*inverseM (2,1) + (i.exp() - A (2))*inverseM (2,2);
746      if (k == 0)
747      {
748        minExpY= expY;
749        minExpX= expX;
750        k= 1;
751      }
752      else
753      {
754        minExpY= (minExpY > expY) ? expY : minExpY;
755        minExpX= (minExpX > expX) ? expX : minExpX;
756      }
757      result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
758      continue;
759    }
760    CFIterator j= i.coeff();
761    if (k == 0)
762    {
763      expX= (j.exp() - A (1))*inverseM (1,1) + (i.exp() - A (2))*inverseM (1,2);
764      expY= (j.exp() - A (1))*inverseM (2,1) + (i.exp() - A (2))*inverseM (2,2);
765      minExpX= expX;
766      minExpY= expY;
767      result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
768      j++;
769      k= 1;
770    }
771
772    for (; j.hasTerms(); j++)
773    {
774      expX= (j.exp() - A (1))*inverseM (1,1) + (i.exp() - A (2))*inverseM (1,2);
775      expY= (j.exp() - A (1))*inverseM (2,1) + (i.exp() - A (2))*inverseM (2,2);
776      result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY));
777      minExpY= (minExpY > expY) ? expY : minExpY;
778      minExpX= (minExpX > expX) ? expX : minExpX;
779    }
780  }
781
782  if (to_long (minExpX) < 0)
783  {
784    result *= power (x,-to_long(minExpX));
785    result /= CanonicalForm (x, 0);
786  }
787  else
788    result /= power (x,to_long(minExpX));
789
790  if (to_long (minExpY) < 0)
791  {
792    result *= power (y,-to_long(minExpY));
793    result /= CanonicalForm (y, 0);
794  }
795  else
796    result /= power (y,to_long(minExpY));
797
798  return result/Lc (result); //normalize
799}
800#endif
801
802//assumes the input is a Newton polygon of a bivariate polynomial which is
803//primitive wrt. x and y, i.e. there is at least one point of the polygon lying
804//on the x-axis and one lying on the y-axis
805int* getRightSide (int** polygon, int sizeOfPolygon, int& sizeOfOutput)
806{
807  int maxY= polygon [0][0];
808  int indexY= 0;
809  for (int i= 1; i < sizeOfPolygon; i++)
810  {
811    if (maxY < polygon [i][0])
812    {
813      maxY= polygon [i][0];
814      indexY= i;
815    }
816    else if (maxY == polygon [i][0])
817    {
818      if (polygon [indexY][1] < polygon[i][1])
819        indexY= i;
820    }
821    if (maxY > polygon [i][0])
822      break;
823  }
824
825  int count= -1;
826  for (int i= indexY; i < sizeOfPolygon; i++)
827  {
828    if (polygon[i][0] == 0)
829    {
830      count= i - indexY;
831      break;
832    }
833  }
834
835  int * result;
836  int index= 0;
837  if (count < 0)
838  {
839    result= new int [sizeOfPolygon - indexY];
840    sizeOfOutput= sizeOfPolygon - indexY;
841    count= sizeOfPolygon - indexY - 1;
842    result [0]= polygon[sizeOfPolygon - 1][0] - polygon [0] [0];
843    index= 1;
844  }
845  else
846  {
847    sizeOfOutput= count;
848    result= new int [count];
849  }
850
851  for (int i= indexY + count; i > indexY; i--, index++)
852    result [index]= polygon [i - 1] [0] - polygon [i] [0];
853
854  return result;
855}
856
857bool irreducibilityTest (const CanonicalForm& F)
858{
859  int sizeOfNewtonPolygon;
860  int ** newtonPolyg= newtonPolygon (F, sizeOfNewtonPolygon);
861  if (sizeOfNewtonPolygon == 3)
862  {
863    bool check1=
864        (newtonPolyg[0][0]==0 || newtonPolyg[1][0]==0 || newtonPolyg[2][0]==0);
865    if (check1)
866    {
867      bool check2=
868        (newtonPolyg[0][1]==0 || newtonPolyg[1][1]==0 || newtonPolyg[2][0]==0);
869      if (check2)
870      {
871        bool isRat= isOn (SW_RATIONAL);
872        if (isRat)
873          Off (SW_RATIONAL);
874        CanonicalForm tmp= gcd (newtonPolyg[0][0],newtonPolyg[0][1]);
875        tmp= gcd (tmp, newtonPolyg[1][0]);
876        tmp= gcd (tmp, newtonPolyg[1][1]);
877        tmp= gcd (tmp, newtonPolyg[2][0]);
878        tmp= gcd (tmp, newtonPolyg[2][1]);
879        if (isRat)
880          On (SW_RATIONAL);
881        for (int i= 0; i < sizeOfNewtonPolygon; i++)
882          delete [] newtonPolyg [i];
883        delete [] newtonPolyg;
884        return (tmp==1);
885      }
886    }
887  }
888  for (int i= 0; i < sizeOfNewtonPolygon; i++)
889    delete [] newtonPolyg [i];
890  delete [] newtonPolyg;
891  return false;
892}
Note: See TracBrowser for help on using the repository browser.