source:git/factory/cfNewtonPolygon.cc@ace936

spielwiese
Last change on this file since ace936 was ace936, checked in by Martin Lee <martinlee84@…>, 10 years ago
chg: keep exponents positive in convex dense (de)compression
• Property mode set to `100644`
File size: 31.5 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#ifdef HAVE_CONFIG_H
15#include "config.h"
16#endif /* HAVE_CONFIG_H */
17
18#include "cf_assert.h"
19
20#include <stdlib.h>
21
22#include "canonicalform.h"
23#include "cf_iter.h"
24#include "cf_algorithm.h"
25#include "cf_primes.h"
26#include "cf_reval.h"
27#include "cf_factory.h"
28#include "gfops.h"
29#include "cfNewtonPolygon.h"
30#include "templates/ftmpl_functions.h"
31#include "algext.h"
32
33static
34void translate (int** points, int* point, int sizePoints) //make point to 0
35{
36  for (int i= 0; i < sizePoints; i++)
37  {
38    points[i] [0] -= point [0];
39    points[i] [1] -= point [1];
40  }
41}
42
43static
44int smallestPointIndex (int** points, int sizePoints)
45{
46  int min= 0;
47  for (int i= 1; i < sizePoints; i++)
48  {
49    if (points[i][0] < points[min][0] ||
50        (points[i] [0] == points[min] [0] && points[i] [1] < points[min] [1]))
51      min= i;
52  }
53  return min;
54}
55
56static
57void swap (int** points, int i, int j)
58{
59  int* tmp= points[i];
60  points[i]= points[j];
61  points[j]= tmp;
62}
63
64static
65bool isLess (int* point1, int* point2)
66{
67  long area= point1[0]*point2[1]- point1[1]*point2[0];
68  if (area > 0) return true;
69  if (area == 0)
70  {
71    return (abs (point1[0]) + abs (point1[1]) >
72            abs (point2[0]) + abs (point2[1]));
73  }
74  return false;
75}
76
77static
78void quickSort (int lo, int hi, int** points)
79{
80  int i= lo, j= hi;
81  int* point= new int [2];
82  point [0]= points [(lo+hi)/2] [0];
83  point [1]= points [(lo+hi)/2] [1];
84  while (i <= j)
85  {
86    while (isLess (points [i], point) && i < hi) i++;
87    while (isLess (point, points[j]) && j > lo) j--;
88    if (i <= j)
89    {
90      swap (points, i, j);
91      i++;
92      j--;
93    }
94  }
95  delete [] point;
96  if (lo < j) quickSort (lo, j, points);
97  if (i < hi) quickSort (i, hi, points);
98}
99
100static
101void sort (int** points, int sizePoints)
102{
103  quickSort (1, sizePoints - 1, points);
104}
105
106// check whether p2 is convex
107static
108bool isConvex (int* point1, int* point2, int* point3)
109{
110  long relArea= (point1[0] - point2[0])*(point3[1] - point2[1]) -
111                (point1[1] - point2[1])*(point3[0] - point2[0]);
112  if (relArea < 0)
113    return true;
114  if (relArea == 0)
115  {
116    return !(abs (point1[0] - point3[0]) + abs (point1[1] - point3[1]) >=
117             (abs (point2[0] - point1[0]) + abs (point2[1] - point1[1]) +
118             abs (point2[0] - point3[0]) + abs (point2[1] - point3[1])));
119  }
120  return false;
121}
122
123static
124bool isConvex (int** points, int i)
125{
126  return isConvex (points[i - 1], points [i], points [i + 1]);
127}
128
129int grahamScan (int** points, int sizePoints)
130{
131  swap (points, 0, smallestPointIndex (points, sizePoints));
132  int * minusPoint= new int [2];
133  minusPoint [0]= points[0] [0];
134  minusPoint [1]= points[0] [1];
135  translate (points, minusPoint, sizePoints);
136  sort (points, sizePoints);
137  minusPoint[0]= - minusPoint[0];
138  minusPoint[1]= - minusPoint[1];
139  translate (points, minusPoint, sizePoints); //reverse translation
140  delete [] minusPoint;
141  int i= 3, k= 3;
142  while (k < sizePoints)
143  {
144    swap (points, i, k);
145    while (!isConvex (points, i - 1))
146    {
147      swap (points, i - 1, i);
148      i--;
149    }
150    k++;
151    i++;
152  }
153  if (i + 1 <= sizePoints || i == sizePoints)
154  {
155    long relArea=
156    (points [i-2][0] - points [i-1][0])*(points [0][1] - points [i-1][1])-
157    (points [i-2][1] - points [i-1][1])*(points [0][0] - points [i-1][0]);
158    if (relArea == 0)
159    {
160      if (abs (points [i-2][0] - points [0][0]) +
161          abs (points [i-2][1] - points [0][1]) >=
162          abs (points [i-1][0] - points [i-2][0]) +
163          abs (points [i-1][1] - points [i-2][1]) +
164          abs (points [i-1][0] - points [0][0]) +
165          abs (points [i-1][1] - points [0][1]))
166          i--;
167    }
168  }
169  return i;
170}
171
172//points[i] [0] is x-coordinate, points [i] [1] is y-coordinate
173int polygon (int** points, int sizePoints)
174{
175  if (sizePoints < 3) return sizePoints;
176  return grahamScan (points, sizePoints);
177}
178
179static
180int* getDegrees (const CanonicalForm& F, int& sizeOfOutput)
181{
182  if (F.inCoeffDomain())
183  {
184    int* result= new int [1];
185    result [0]= 0;
186    sizeOfOutput= 1;
187    return result;
188  }
189  sizeOfOutput= size (F);
190  int* result= new int [sizeOfOutput];
191  int j= 0;
192  for (CFIterator i= F; i.hasTerms(); i++, j++)
193    result [j]= i.exp();
194  return result;
195}
196
197//get points in Z^2 whose convex hull is the Newton polygon
198int ** getPoints (const CanonicalForm& F, int& n)
199{
200  n= size (F);
201  int ** points= new int* [n];
202  for (int i= 0; i < n; i++)
203    points [i]= new int [2];
204
205  int j= 0;
206  int * buf;
207  int bufSize;
208  if (F.isUnivariate() && F.level() == 1)
209  {
210    for (CFIterator i= F; i.hasTerms(); i++, j++)
211    {
212      points [j] [0]= i.exp();
213      points [j] [1]= 0;
214    }
215    return points;
216  }
217  for (CFIterator i= F; i.hasTerms(); i++)
218  {
219    buf= getDegrees (i.coeff(), bufSize);
220    for (int k= 0; k < bufSize; k++, j++)
221    {
222      points [j] [0]= i.exp();
223      points [j] [1]= buf [k];
224    }
225    delete [] buf;
226  }
227  return points;
228}
229
230int **
231merge (int ** points1, int sizePoints1, int ** points2, int sizePoints2,
232       int& sizeResult)
233{
234  int i, j;
235  sizeResult= sizePoints1+sizePoints2;
236  for (i= 0; i < sizePoints1; i++)
237  {
238    for (j= 0; j < sizePoints2; j++)
239    {
240      if (points1[i][0] != points2[j][0])
241        continue;
242      else
243      {
244        if (points1[i][1] != points2[j][1])
245          continue;
246        else
247        {
248          points2[j][0]= -1;
249          points2[j][1]= -1;
250          sizeResult--;
251        }
252      }
253    }
254  }
255  if (sizeResult == 0)
256    return points1;
257
258  int ** result= new int *[sizeResult];
259  for (i= 0; i < sizeResult; i++)
260    result [i]= new int [2];
261
262  int k= 0;
263  for (i= 0; i < sizePoints1; i++, k++)
264  {
265    result[k][0]= points1[i][0];
266    result[k][1]= points1[i][1];
267  }
268  for (i= 0; i < sizePoints2; i++)
269  {
270    if (points2[i][0] < 0)
271      continue;
272    else
273    {
274      result[k][0]= points2[i][0];
275      result[k][1]= points2[i][1];
276      k++;
277    }
278  }
279  return result;
280}
281
282// assumes a bivariate poly as input
283int ** newtonPolygon (const CanonicalForm& F, int& sizeOfNewtonPoly)
284{
285  int sizeF= size (F);
286  int ** points= new int* [sizeF];
287  for (int i= 0; i < sizeF; i++)
288    points [i]= new int [2];
289  int j= 0;
290  int * buf;
291  int bufSize;
292  for (CFIterator i= F; i.hasTerms(); i++)
293  {
294    buf= getDegrees (i.coeff(), bufSize);
295    for (int k= 0; k < bufSize; k++, j++)
296    {
297      points [j] [0]= i.exp();
298      points [j] [1]= buf [k];
299    }
300    delete [] buf;
301  }
302
303  int n= polygon (points, sizeF);
304
305  int ** result= new int* [n];
306  for (int i= 0; i < n; i++)
307  {
308    result [i]= new int [2];
309    result [i] [0]= points [i] [0];
310    result [i] [1]= points [i] [1];
311  }
312
313  sizeOfNewtonPoly= n;
314  for (int i= 0; i < sizeF; i++)
315    delete [] points[i];
316  delete [] points;
317
318  return result;
319}
320
321// assumes a bivariate polys as input
322int ** newtonPolygon (const CanonicalForm& F, const CanonicalForm& G,
323                      int& sizeOfNewtonPoly)
324{
325  int sizeF= size (F);
326  int ** pointsF= new int* [sizeF];
327  for (int i= 0; i < sizeF; i++)
328    pointsF [i]= new int [2];
329  int j= 0;
330  int * buf;
331  int bufSize;
332  for (CFIterator i= F; i.hasTerms(); i++)
333  {
334    buf= getDegrees (i.coeff(), bufSize);
335    for (int k= 0; k < bufSize; k++, j++)
336    {
337      pointsF [j] [0]= i.exp();
338      pointsF [j] [1]= buf [k];
339    }
340    delete [] buf;
341  }
342
343  int sizeG= size (G);
344  int ** pointsG= new int* [sizeG];
345  for (int i= 0; i < sizeG; i++)
346    pointsG [i]= new int [2];
347  j= 0;
348  for (CFIterator i= G; i.hasTerms(); i++)
349  {
350    buf= getDegrees (i.coeff(), bufSize);
351    for (int k= 0; k < bufSize; k++, j++)
352    {
353      pointsG [j] [0]= i.exp();
354      pointsG [j] [1]= buf [k];
355    }
356    delete [] buf;
357  }
358
359  int sizePoints;
360  int ** points= merge (pointsF, sizeF, pointsG, sizeG, sizePoints);
361
362  int n= polygon (points, sizePoints);
363
364  int ** result= new int* [n];
365  for (int i= 0; i < n; i++)
366  {
367    result [i]= new int [2];
368    result [i] [0]= points [i] [0];
369    result [i] [1]= points [i] [1];
370  }
371
372  sizeOfNewtonPoly= n;
373  for (int i= 0; i < sizeF; i++)
374    delete [] pointsF[i];
375  delete [] pointsF;
376  for (int i= 0; i < sizeG; i++)
377    delete [] pointsG[i];
378  delete [] pointsG;
379
380  return result;
381}
382
383// assumes first sizePoints entries of points form a Newton polygon
384bool isInPolygon (int ** points, int sizePoints, int* point)
385{
386  int ** buf= new int* [sizePoints + 1];
387  for (int i= 0; i < sizePoints; i++)
388  {
389    buf [i]= new int [2];
390    buf [i] [0]= points [i] [0];
391    buf [i] [1]= points [i] [1];
392  }
393  buf [sizePoints]= new int [2];
394  buf [sizePoints] [0]= point [0];
395  buf [sizePoints] [1]= point [1];
396  int sizeBuf= sizePoints + 1;
397
398  swap (buf, 0, smallestPointIndex (buf, sizeBuf));
399  int * minusPoint= new int [2];
400  minusPoint [0]= buf[0] [0];
401  minusPoint [1]= buf[0] [1];
402  translate (buf, minusPoint, sizeBuf);
403  sort (buf, sizeBuf);
404  minusPoint[0]= - minusPoint[0];
405  minusPoint[1]= - minusPoint[1];
406  translate (buf, minusPoint, sizeBuf); //reverse translation
407  delete [] minusPoint;
408
409  if (buf [0] [0] == point [0] && buf [0] [1] == point [1])
410  {
411    for (int i= 0; i < sizeBuf; i++)
412      delete [] buf[i];
413    delete [] buf;
414    return false;
415  }
416
417  for (int i= 1; i < sizeBuf-1; i++)
418  {
419    if (buf [i] [0] == point [0] && buf [i] [1] == point [1])
420    {
421      bool result= !isConvex (buf, i);
422      for (int i= 0; i < sizeBuf; i++)
423        delete [] buf [i];
424      delete [] buf;
425      return result;
426    }
427  }
428
429  if (buf [sizeBuf - 1] [0] == point [0] && buf [sizeBuf-1] [1] == point [1])
430  {
431    buf [1] [0]= point [0];
432    buf [1] [1]= point [1];
433    buf [2] [0]= buf [0] [0];
434    buf [2] [1]= buf [0] [1];
435    buf [0] [0]= buf [sizeBuf-2] [0];
436    buf [0] [1]= buf [sizeBuf-2] [1];
437    bool result= !isConvex (buf, 1);
438    for (int i= 0; i < sizeBuf; i++)
439      delete [] buf [i];
440    delete [] buf;
441    return result;
442  }
443  for (int i= 0; i < sizeBuf; i++)
444    delete [] buf [i];
445  delete [] buf;
446
447  return false;
448}
449
450void lambda (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 lambdaInverse (int** points, int sizePoints)
457{
458  for (int i= 0; i < sizePoints; i++)
459    points [i] [1]= points [i] [1] + points [i] [0];
460}
461
462void tau (int** points, int sizePoints, int k)
463{
464  for (int i= 0; i < sizePoints; i++)
465    points [i] [1]= points [i] [1] + k;
466}
467
468void mu (int** points, int sizePoints)
469{
470  int tmp;
471  for (int i= 0; i < sizePoints; i++)
472  {
473    tmp= points [i] [0];
474    points [i] [0]= points [i] [1];
475    points [i] [1]= tmp;
476  }
477}
478
479void getMaxMin (int** points, int sizePoints, int& minDiff, int& minSum,
480                int& maxDiff, int& maxSum, int& maxX, int& maxY
481               )
482{
483  minDiff= points [0] [1] - points [0] [0];
484  minSum= points [0] [1] + points [0] [0];
485  maxDiff= points [0] [1] - points [0] [0];
486  maxSum= points [0] [1] + points [0] [0];
487  maxX= points [0] [1];
488  maxY= points [0] [0];
489  int diff, sum;
490  for (int i= 1; i < sizePoints; i++)
491  {
492    diff= points [i] [1] - points [i] [0];
493    sum= points [i] [1] + points [i] [0];
494    minDiff= tmin (minDiff, diff);
495    minSum= tmin (minSum, sum);
496    maxDiff= tmax (maxDiff, diff);
497    maxSum= tmax (maxSum, sum);
498    maxX= tmax (maxX, points [i] [1]);
499    maxY= tmax (maxY, points [i] [0]);
500  }
501}
502
503void mpz_mat_mul (const mpz_t* N, mpz_t*& M)
504{
505  mpz_t * tmp= new mpz_t[4];
506
507  mpz_init_set (tmp[0], N[0]);
508  mpz_mul (tmp[0], tmp[0], M[0]);
510
511  mpz_init_set (tmp[1], N[0]);
512  mpz_mul (tmp[1], tmp[1], M[1]);
514
515  mpz_init_set (tmp[2], N[2]);
516  mpz_mul (tmp[2], tmp[2], M[0]);
518
519  mpz_init_set (tmp[3], N[2]);
520  mpz_mul (tmp[3], tmp[3], M[1]);
522
523  mpz_set (M[0], tmp[0]);
524  mpz_set (M[1], tmp[1]);
525  mpz_set (M[2], tmp[2]);
526  mpz_set (M[3], tmp[3]);
527
528  mpz_clear (tmp[0]);
529  mpz_clear (tmp[1]);
530  mpz_clear (tmp[2]);
531  mpz_clear (tmp[3]);
532
533  delete [] tmp;
534}
535
536void mpz_mat_inv (mpz_t*& M)
537{
538  mpz_t det;
539  mpz_init_set (det, M[0]);
540  mpz_mul (det, det, M[3]);
541  mpz_submul (det, M[1], M[2]);
542
543  mpz_t tmp;
544  mpz_init_set (tmp, M[0]);
545  mpz_divexact (tmp, tmp, det);
546  mpz_set (M[0], M[3]);
547  mpz_divexact (M[0], M[0], det);
548  mpz_set (M[3], tmp);
549
550  mpz_neg (M[1], M[1]);
551  mpz_divexact (M[1], M[1], det);
552  mpz_neg (M[2], M[2]);
553  mpz_divexact (M[2], M[2], det);
554
555  mpz_clear (det);
556  mpz_clear (tmp);
557}
558
559void convexDense(int** points, int sizePoints, mpz_t*& M, mpz_t*& A)
560{
561  if (sizePoints < 3)
562  {
563    if (sizePoints == 2)
564    {
565      mpz_t u,v,g,maxX,maxY;
566      mpz_init (u);
567      mpz_init (v);
568      mpz_init (g);
569      mpz_init_set_si (maxX,
570                       (points[1][1] < points[0][1])?points[0][1]:points[1][1]);
571      mpz_init_set_si (maxY,
572                       (points[1][0] < points[0][0])?points[0][0]:points[1][0]);
573      mpz_gcdext (g, u, v, maxX, maxY);
574      if (points [0] [1] != points [0] [0] && points [1] [0] != points [1] [1])
575      {
576        mpz_set (A[0], u);
577        mpz_mul (A[0], A[0], maxX);
578        mpz_set (M[2], maxY);
579        mpz_divexact (M[2], M[2], g);
580        mpz_set (A[1], M[2]);
581        mpz_neg (A[1], A[1]);
582        mpz_mul (A[1], A[1], maxX);
583        mpz_neg (u, u);
584        mpz_set (M[0], u);
585        mpz_set (M[1], v);
586        mpz_set (M[3], maxX);
587        mpz_divexact (M[3], M[3], g);
588      }
589      else
590      {
591        mpz_set (M[0], u);
592        mpz_set (M[1], v);
593        mpz_set (M[2], maxY);
594        mpz_divexact (M[2], M[2], g);
595        mpz_neg (M[2], M[2]);
596        mpz_set (M[3], maxX);
597        mpz_divexact (M[3], M[3], g);
598      }
599      mpz_clear (u);
600      mpz_clear (v);
601      mpz_clear (g);
602      mpz_clear (maxX);
603      mpz_clear (maxY);
604    }
605    else if (sizePoints == 1)
606    {
607      mpz_set_si (M[0], 1);
608      mpz_set_si (M[3], 1);
609    }
610    return;
611  }
612  mpz_set_si (M[0], 1);
613  mpz_set_si (M[3], 1);
614
615  mpz_t * Mu= new mpz_t[4];
616  mpz_init_set_si (Mu[1], 1);
617  mpz_init_set_si (Mu[2], 1);
618  mpz_init (Mu[0]);
619  mpz_init (Mu[3]);
620
621  mpz_t * Lambda= new mpz_t[4];
622  mpz_init_set_si (Lambda[0], 1);
623  mpz_init_set_si (Lambda[1], -1);
624  mpz_init_set_si (Lambda[3], 1);
625  mpz_init (Lambda[2]);
626
627  mpz_t * InverseLambda= new mpz_t[4];
628  mpz_init_set_si (InverseLambda[0], 1);
629  mpz_init_set_si (InverseLambda[1], 1);
630  mpz_init_set_si (InverseLambda[3], 1);
631  mpz_init (InverseLambda[2]);
632
633  mpz_t tmp;
634  mpz_init (tmp);
635  int minDiff, minSum, maxDiff, maxSum, maxX, maxY, b, d, f, h;
636  getMaxMin(points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY);
637  do
638  {
639    if (maxX < maxY)
640    {
641      mu (points, sizePoints);
642
643      mpz_mat_mul (Mu, M);
644
645      mpz_set (tmp, A[0]);
646      mpz_set (A[0], A[1]);
647      mpz_set (A[1], tmp);
648    }
649    getMaxMin(points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY);
650    b= maxX - maxDiff;
651    d= maxX + maxY - maxSum;
652    f= maxY + minDiff;
653    h= minSum;
654    if (b + f > maxY)
655    {
656      lambda (points, sizePoints);
657      tau (points, sizePoints, maxY - f);
658
659      mpz_mat_mul (Lambda, M);
660
661      if (maxY-f > 0)
663      else
665      maxX= maxX + maxY - b - f;
666    }
667    else if (d + h > maxY)
668    {
669      lambdaInverse (points, sizePoints);
670      tau (points, sizePoints, -h);
671
672      mpz_mat_mul (InverseLambda, M);
673
674      if (h < 0)
676      else
677        mpz_sub_ui (A[0], A[0], h);
678      maxX= maxX + maxY - d - h;
679    }
680    else
681    {
682      mpz_clear (tmp);
683      mpz_clear (Mu[0]);
684      mpz_clear (Mu[1]);
685      mpz_clear (Mu[2]);
686      mpz_clear (Mu[3]);
687      delete [] Mu;
688
689      mpz_clear (Lambda[0]);
690      mpz_clear (Lambda[1]);
691      mpz_clear (Lambda[2]);
692      mpz_clear (Lambda[3]);
693      delete [] Lambda;
694
695      mpz_clear (InverseLambda[0]);
696      mpz_clear (InverseLambda[1]);
697      mpz_clear (InverseLambda[2]);
698      mpz_clear (InverseLambda[3]);
699      delete [] InverseLambda;
700
701      return;
702    }
703  } while (1);
704}
705
706CanonicalForm
707compress (const CanonicalForm& F, mpz_t*& M, mpz_t*& A, bool computeMA)
708{
709  int n;
710  int ** newtonPolyg= NULL;
711  if (computeMA)
712  {
713    newtonPolyg= newtonPolygon (F, n);
714    convexDense (newtonPolyg, n, M, A);
715  }
716
717  CanonicalForm result= 0;
718  Variable x= Variable (1);
719  Variable y= Variable (2);
720
721  mpz_t expX, expY, minExpX, minExpY;
722  mpz_init (expX);
723  mpz_init (expY);
724  mpz_init (minExpX);
725  mpz_init (minExpY);
726
727  int k= 0;
728  Variable alpha;
729  mpz_t * exps= new mpz_t [2*size (F)];
730  int count= 0;
731  for (CFIterator i= F; i.hasTerms(); i++)
732  {
733    if (i.coeff().inCoeffDomain() && hasFirstAlgVar (i.coeff(), alpha))
734    {
735      mpz_set (expX, A[0]);
736      mpz_set (expY, A[1]);
739
740      if (k == 0)
741      {
742        mpz_set (minExpX, expX);
743        mpz_set (minExpY, expY);
744        k= 1;
745      }
746      else
747      {
748        if (mpz_cmp (minExpY, expY) > 0)
749          mpz_set (minExpY, expY);
750        if (mpz_cmp (minExpX, expX) > 0)
751          mpz_set (minExpX, expX);
752      }
753      mpz_init_set (exps[count], expX);
754      count++;
755      mpz_init_set (exps[count], expY);
756      count++;
757      continue;
758    }
759    CFIterator j= i.coeff();
760    if (k == 0)
761    {
762      mpz_set (expX, A[0]);
765
766      mpz_set (expY, A[1]);
769
770      mpz_set (minExpX, expX);
771      mpz_set (minExpY, expY);
772      mpz_init_set (exps[count], expX);
773      count++;
774      mpz_init_set (exps[count], expY);
775      count++;
776      j++;
777      k= 1;
778    }
779
780    for (; j.hasTerms(); j++)
781    {
782      mpz_set (expX, A[0]);
785
786      mpz_set (expY, A[1]);
789
790      mpz_init_set (exps[count], expX);
791      count++;
792      mpz_init_set (exps[count], expY);
793      count++;
794      if (mpz_cmp (minExpY, expY) > 0)
795        mpz_set (minExpY, expY);
796      if (mpz_cmp (minExpX, expX) > 0)
797        mpz_set (minExpX, expX);
798    }
799  }
800
801  count= 0;
802  int mExpX= mpz_get_si (minExpX);
803  int mExpY= mpz_get_si (minExpY);
804  for (CFIterator i= F; i.hasTerms(); i++)
805  {
806    if (i.coeff().inCoeffDomain() && hasFirstAlgVar (i.coeff(), alpha))
807    {
808      result += i.coeff()*power (x, mpz_get_si (exps[count])-mExpX)*
809                power (y, mpz_get_si (exps[count+1])-mExpY);
810      count += 2;
811      continue;
812    }
813    CFIterator j= i.coeff();
814    for (; j.hasTerms(); j++)
815    {
816      result += j.coeff()*power (x, mpz_get_si (exps[count])-mExpX)*
817                power (y, mpz_get_si (exps[count+1])-mExpY);
818      count += 2;
819    }
820  }
821
822  CanonicalForm tmp= LC (result);
823  if (tmp.inPolyDomain() && degree (tmp) <= 0)
824  {
825    int d= degree (result);
826    Variable x= result.mvar();
827    result -= tmp*power (x, d);
828    result += Lc (tmp)*power (x, d);
829  }
830
831  if (computeMA)
832  {
833    for (int i= 0; i < n; i++)
834      delete [] newtonPolyg [i];
835    delete [] newtonPolyg;
836    mpz_mat_inv (M);
837  }
838
839  delete [] exps;
840  mpz_clear (expX);
841  mpz_clear (expY);
842  mpz_clear (minExpX);
843  mpz_clear (minExpY);
844
845  return result;
846}
847
848CanonicalForm
849decompress (const CanonicalForm& F, const mpz_t* inverseM, const mpz_t * A)
850{
851  CanonicalForm result= 0;
852  Variable x= Variable (1);
853  Variable y= Variable (2);
854
855  mpz_t expX, expY, minExpX, minExpY;
856  mpz_init (expX);
857  mpz_init (expY);
858  mpz_init (minExpX);
859  mpz_init (minExpY);
860
861  mpz_t * exps= new mpz_t [2*size(F)];
862  int count= 0;
863  if (F.isUnivariate() && F.level() == 1)
864  {
865    CFIterator i= F;
866
867    mpz_set_si (expX, i.exp());
868    mpz_sub (expX, expX, A[0]);
869    mpz_mul (expX, expX, inverseM[0]);
870    mpz_submul (expX, inverseM[1], A[1]);
871
872    mpz_set_si (expY, i.exp());
873    mpz_sub (expY, expY, A[0]);
874    mpz_mul (expY, expY, inverseM[2]);
875    mpz_submul (expY, inverseM[3], A[1]);
876
877    mpz_set (minExpX, expX);
878    mpz_set (minExpY, expY);
879
880    mpz_init_set (exps[count], expX);
881    mpz_init_set (exps[count+1], expY);
882    count += 2;
883
884    i++;
885    for (; i.hasTerms(); i++)
886    {
887      mpz_set_si (expX, i.exp());
888      mpz_sub (expX, expX, A[0]);
889      mpz_mul (expX, expX, inverseM[0]);
890      mpz_submul (expX, inverseM[1], A[1]);
891
892      mpz_set_si (expY, i.exp());
893      mpz_sub (expY, expY, A[0]);
894      mpz_mul (expY, expY, inverseM[2]);
895      mpz_submul (expY, inverseM[3], A[1]);
896
897      mpz_init_set (exps[count], expX);
898      mpz_init_set (exps[count+1], expY);
899      count += 2;
900
901      if (mpz_cmp (minExpY, expY) > 0)
902        mpz_set (minExpY, expY);
903      if (mpz_cmp (minExpX, expX) > 0)
904        mpz_set (minExpX, expX);
905    }
906
907    int mExpX= mpz_get_si (minExpX);
908    int mExpY= mpz_get_si (minExpY);
909    count= 0;
910    for (i= F; i.hasTerms(); i++)
911    {
912      result += i.coeff()*power (x, mpz_get_si (exps[count])-mExpX)*
913                power (y, mpz_get_si (exps[count+1])-mExpY);
914      count += 2;
915    }
916
917    mpz_clear (expX);
918    mpz_clear (expY);
919    mpz_clear (minExpX);
920    mpz_clear (minExpY);
921
922    delete [] exps;
923    return result/ Lc (result); //normalize
924  }
925
926  mpz_t tmp;
927  mpz_init (tmp);
928  int k= 0;
929  Variable alpha;
930  for (CFIterator i= F; i.hasTerms(); i++)
931  {
932    if (i.coeff().inCoeffDomain() && hasFirstAlgVar (i.coeff(), alpha))
933    {
934      mpz_set_si (expX, i.exp());
935      mpz_sub (expX, expX, A[1]);
936      mpz_mul (expX, expX, inverseM[1]);
937      mpz_submul (expX, A[0], inverseM[0]);
938
939      mpz_set_si (expY, i.exp());
940      mpz_sub (expY, expY, A[1]);
941      mpz_mul (expY, expY, inverseM[3]);
942      mpz_submul (expY, A[0], inverseM[2]);
943
944      if (k == 0)
945      {
946        mpz_set (minExpX, expX);
947        mpz_set (minExpY, expY);
948        k= 1;
949      }
950      else
951      {
952        if (mpz_cmp (minExpY, expY) > 0)
953          mpz_set (minExpY, expY);
954        if (mpz_cmp (minExpX, expX) > 0)
955          mpz_set (minExpX, expX);
956      }
957      mpz_init_set (exps[count], expX);
958      mpz_init_set (exps[count+1], expY);
959      count += 2;
960      continue;
961    }
962    CFIterator j= i.coeff();
963    if (k == 0)
964    {
965      mpz_set_si (expX, j.exp());
966      mpz_sub (expX, expX, A[0]);
967      mpz_mul (expX, expX, inverseM[0]);
968      mpz_set_si (tmp, i.exp());
969      mpz_sub (tmp, tmp, A[1]);
971
972      mpz_set_si (expY, j.exp());
973      mpz_sub (expY, expY, A[0]);
974      mpz_mul (expY, expY, inverseM[2]);
975      mpz_set_si (tmp, i.exp());
976      mpz_sub (tmp, tmp, A[1]);
978
979      mpz_set (minExpX, expX);
980      mpz_set (minExpY, expY);
981
982      mpz_init_set (exps[count], expX);
983      mpz_init_set (exps[count+1], expY);
984      count += 2;
985
986      j++;
987      k= 1;
988    }
989
990    for (; j.hasTerms(); j++)
991    {
992      mpz_set_si (expX, j.exp());
993      mpz_sub (expX, expX, A[0]);
994      mpz_mul (expX, expX, inverseM[0]);
995      mpz_set_si (tmp, i.exp());
996      mpz_sub (tmp, tmp, A[1]);
998
999      mpz_set_si (expY, j.exp());
1000      mpz_sub (expY, expY, A[0]);
1001      mpz_mul (expY, expY, inverseM[2]);
1002      mpz_set_si (tmp, i.exp());
1003      mpz_sub (tmp, tmp, A[1]);
1005
1006      mpz_init_set (exps[count], expX);
1007      mpz_init_set (exps[count+1], expY);
1008      count += 2;
1009
1010      if (mpz_cmp (minExpY, expY) > 0)
1011        mpz_set (minExpY, expY);
1012      if (mpz_cmp (minExpX, expX) > 0)
1013        mpz_set (minExpX, expX);
1014    }
1015  }
1016
1017  int mExpX= mpz_get_si (minExpX);
1018  int mExpY= mpz_get_si (minExpY);
1019  count= 0;
1020  for (CFIterator i= F; i.hasTerms(); i++)
1021  {
1022    if (i.coeff().inCoeffDomain() && hasFirstAlgVar (i.coeff(), alpha))
1023    {
1024      result += i.coeff()*power (x, mpz_get_si (exps[count])-mExpX)*
1025                power (y, mpz_get_si (exps[count+1])-mExpY);
1026      count += 2;
1027      continue;
1028    }
1029    for (CFIterator j= i.coeff(); j.hasTerms(); j++)
1030    {
1031      result += j.coeff()*power (x, mpz_get_si (exps[count])-mExpX)*
1032                power (y, mpz_get_si (exps[count+1])-mExpY);
1033      count += 2;
1034    }
1035  }
1036
1037  mpz_clear (expX);
1038  mpz_clear (expY);
1039  mpz_clear (minExpX);
1040  mpz_clear (minExpY);
1041  mpz_clear (tmp);
1042
1043  delete [] exps;
1044
1045  return result/Lc (result); //normalize
1046}
1047
1048//assumes the input is a Newton polygon of a bivariate polynomial which is
1049//primitive wrt. x and y, i.e. there is at least one point of the polygon lying
1050//on the x-axis and one lying on the y-axis
1051int* getRightSide (int** polygon, int sizeOfPolygon, int& sizeOfOutput)
1052{
1053  int maxY= polygon [0][0];
1054  int indexY= 0;
1055  for (int i= 1; i < sizeOfPolygon; i++)
1056  {
1057    if (maxY < polygon [i][0])
1058    {
1059      maxY= polygon [i][0];
1060      indexY= i;
1061    }
1062    else if (maxY == polygon [i][0])
1063    {
1064      if (polygon [indexY][1] < polygon[i][1])
1065        indexY= i;
1066    }
1067    if (maxY > polygon [i][0])
1068      break;
1069  }
1070
1071  int count= -1;
1072  for (int i= indexY; i < sizeOfPolygon; i++)
1073  {
1074    if (polygon[i][0] == 0)
1075    {
1076      count= i - indexY;
1077      break;
1078    }
1079  }
1080
1081  int * result;
1082  int index= 0;
1083  if (count < 0)
1084  {
1085    result= new int [sizeOfPolygon - indexY];
1086    sizeOfOutput= sizeOfPolygon - indexY;
1087    count= sizeOfPolygon - indexY - 1;
1088    result [0]= polygon[sizeOfPolygon - 1][0] - polygon [0] [0];
1089    index= 1;
1090  }
1091  else
1092  {
1093    sizeOfOutput= count;
1094    result= new int [count];
1095  }
1096
1097  for (int i= indexY + count; i > indexY; i--, index++)
1098    result [index]= polygon [i - 1] [0] - polygon [i] [0];
1099
1100  return result;
1101}
1102
1103bool irreducibilityTest (const CanonicalForm& F)
1104{
1105  ASSERT (getNumVars (F) == 2, "expected bivariate polynomial");
1106  ASSERT (getCharacteristic() == 0, "expected polynomial over integers or rationals");
1107
1108  int sizeOfNewtonPolygon;
1109  int ** newtonPolyg= newtonPolygon (F, sizeOfNewtonPolygon);
1110  if (sizeOfNewtonPolygon == 3)
1111  {
1112    bool check1=
1113        (newtonPolyg[0][0]==0 || newtonPolyg[1][0]==0 || newtonPolyg[2][0]==0);
1114    if (check1)
1115    {
1116      bool check2=
1117        (newtonPolyg[0][1]==0 || newtonPolyg[1][1]==0 || newtonPolyg[2][0]==0);
1118      if (check2)
1119      {
1120        bool isRat= isOn (SW_RATIONAL);
1121        if (isRat)
1122          Off (SW_RATIONAL);
1123        CanonicalForm tmp= gcd (newtonPolyg[0][0],newtonPolyg[0][1]); // maybe it's better to use plain intgcd
1124        tmp= gcd (tmp, newtonPolyg[1][0]);
1125        tmp= gcd (tmp, newtonPolyg[1][1]);
1126        tmp= gcd (tmp, newtonPolyg[2][0]);
1127        tmp= gcd (tmp, newtonPolyg[2][1]);
1128        if (isRat)
1129          On (SW_RATIONAL);
1130        for (int i= 0; i < sizeOfNewtonPolygon; i++)
1131          delete [] newtonPolyg [i];
1132        delete [] newtonPolyg;
1133        return (tmp==1);
1134      }
1135    }
1136  }
1137  for (int i= 0; i < sizeOfNewtonPolygon; i++)
1138    delete [] newtonPolyg [i];
1139  delete [] newtonPolyg;
1140  return false;
1141}
1142
1143bool absIrredTest (const CanonicalForm& F)
1144{
1145  ASSERT (getNumVars (F) == 2, "expected bivariate polynomial");
1146  ASSERT (factorize (F).length() <= 2, " expected irreducible polynomial");
1147
1148  int sizeOfNewtonPolygon;
1149  int ** newtonPolyg= newtonPolygon (F, sizeOfNewtonPolygon);
1150  bool isRat= isOn (SW_RATIONAL);
1151  if (isRat)
1152    Off (SW_RATIONAL);
1153  int p=getCharacteristic();
1154  int d=1;
1155  char bufGFName='Z';
1156  bool GF= (CFFactory::gettype()==GaloisFieldDomain);
1157  if (GF)
1158  {
1159    d= getGFDegree();
1160    bufGFName=gf_name;
1161  }
1162
1163  setCharacteristic(0);
1164
1165  CanonicalForm g= gcd (newtonPolyg[0][0], newtonPolyg[0][1]); //maybe it's better to use plain intgcd
1166
1167  int i= 1;
1168  while (!g.isOne() && i < sizeOfNewtonPolygon)
1169  {
1170    g= gcd (g, newtonPolyg[i][0]);
1171    g= gcd (g, newtonPolyg[i][1]);
1172    i++;
1173  }
1174
1175  bool result= g.isOne();
1176
1177  if (GF)
1178    setCharacteristic (p, d, bufGFName);
1179  else
1180    setCharacteristic(p);
1181
1182  if (isRat)
1183    On (SW_RATIONAL);
1184
1185  for (int i= 0; i < sizeOfNewtonPolygon; i++)
1186    delete [] newtonPolyg[i];
1187
1188  delete [] newtonPolyg;
1189
1190  return result;
1191}
1192
1193bool modularIrredTest (const CanonicalForm& F)
1194{
1195  ASSERT (getNumVars (F) == 2, "expected bivariate polynomial");
1196  ASSERT (factorize (F).length() <= 2, " expected irreducible polynomial");
1197
1198  bool isRat= isOn (SW_RATIONAL);
1199  if (isRat)
1200    Off (SW_RATIONAL);
1201
1202  if (isRat)
1203  {
1204    ASSERT (bCommonDen (F).isOne(), "poly over Z expected");
1205  }
1206
1207  CanonicalForm Fp, N= maxNorm (F);
1208  int tdeg= totaldegree (F);
1209
1210  int i= 0;
1211  //TODO: maybe it's better to choose the characteristic as large as possible
1212  //      as factorization over large finite field will be faster
1213  //TODO: handle those cases where our factory primes are not enough
1214  //TODO: factorize coefficients corresponding to the vertices of the Newton
1215  //      polygon and only try the obtained factors
1216  if (N < cf_getSmallPrime (cf_getNumSmallPrimes()-1))
1217  {
1218    while (i < cf_getNumSmallPrimes() && N > cf_getSmallPrime(i))
1219    {
1220      setCharacteristic (cf_getSmallPrime (i));
1221      Fp= F.mapinto();
1222      i++;
1223      if (totaldegree (Fp) == tdeg)
1224      {
1225        if (absIrredTest (Fp))
1226        {
1227          CFFList factors= factorize (Fp);
1228          if (factors.length() == 2 && factors.getLast().exp() == 1)
1229          {
1230            if (isRat)
1231              On (SW_RATIONAL);
1232            setCharacteristic (0);
1233            return true;
1234          }
1235        }
1236      }
1237      setCharacteristic (0);
1238    }
1239  }
1240  else
1241  {
1242    while (i < cf_getNumPrimes() && N > cf_getPrime (i))
1243    {
1244      setCharacteristic (cf_getPrime (i));
1245      Fp= F.mapinto();
1246      i++;
1247      if (totaldegree (Fp) == tdeg)
1248      {
1249        if (absIrredTest (Fp))
1250        {
1251          CFFList factors= factorize (Fp);
1252          if (factors.length() == 2 && factors.getLast().exp() == 1)
1253          {
1254            if (isRat)
1255              On (SW_RATIONAL);
1256            setCharacteristic (0);
1257            return true;
1258          }
1259        }
1260      }
1261      setCharacteristic (0);
1262    }
1263  }
1264
1265  if (isRat)
1266    On (SW_RATIONAL);
1267
1268  return false;
1269}
1270
1271bool
1272modularIrredTestWithShift (const CanonicalForm& F)
1273{
1274  ASSERT (getNumVars (F) == 2, "expected bivariate polynomial");
1275  ASSERT (factorize (F).length() <= 2, " expected irreducible polynomial");
1276
1277  bool isRat= isOn (SW_RATIONAL);
1278  if (isRat)
1279    Off (SW_RATIONAL);
1280
1281  if (isRat)
1282  {
1283    ASSERT (bCommonDen (F).isOne(), "poly over Z expected");
1284  }
1285
1286  Variable x= Variable (1);
1287  Variable y= Variable (2);
1288  CanonicalForm Fp;
1289  int tdeg= totaldegree (F);
1290
1291  REvaluation E;
1292
1293  setCharacteristic (2);
1294  Fp= F.mapinto();
1295
1296  E= REvaluation (1,2, FFRandom());
1297
1298  E.nextpoint();
1299
1300  Fp= Fp (x+E[1], x);
1301  Fp= Fp (y+E[2], y);
1302
1303  if (tdeg == totaldegree (Fp))
1304  {
1305    if (absIrredTest (Fp))
1306    {
1307      CFFList factors= factorize (Fp);
1308      if (factors.length() == 2 && factors.getLast().exp() == 1)
1309      {
1310        if (isRat)
1311          On (SW_RATIONAL);
1312        setCharacteristic (0);
1313        return true;
1314      }
1315    }
1316  }
1317
1318  E.nextpoint();
1319
1320  Fp= Fp (x+E[1], x);
1321  Fp= Fp (y+E[2], y);
1322
1323  if (tdeg == totaldegree (Fp))
1324  {
1325    if (absIrredTest (Fp))
1326    {
1327      CFFList factors= factorize (Fp);
1328      if (factors.length() == 2 && factors.getLast().exp() == 1)
1329      {
1330        if (isRat)
1331          On (SW_RATIONAL);
1332        setCharacteristic (0);
1333        return true;
1334      }
1335    }
1336  }
1337
1338  int i= 0;
1339  while (cf_getSmallPrime (i) < 102)
1340  {
1341    setCharacteristic (cf_getSmallPrime (i));
1342    i++;
1343    E= REvaluation (1, 2, FFRandom());
1344
1345    for (int j= 0; j < 3; j++)
1346    {
1347      Fp= F.mapinto();
1348      E.nextpoint();
1349      Fp= Fp (x+E[1], x);
1350      Fp= Fp (y+E[2], y);
1351
1352      if (tdeg == totaldegree (Fp))
1353      {
1354        if (absIrredTest (Fp))
1355        {
1356          CFFList factors= factorize (Fp);
1357          if (factors.length() == 2 && factors.getLast().exp() == 1)
1358          {
1359            if (isRat)
1360              On (SW_RATIONAL);
1361            setCharacteristic (0);
1362            return true;
1363          }
1364        }
1365      }
1366    }
1367  }
1368
1369  setCharacteristic (0);
1370  if (isRat)
1371    On (SW_RATIONAL);
1372
1373  return false;
1374}
1375
1376
Note: See TracBrowser for help on using the repository browser.