35 for (
int i= 0;
i < sizePoints;
i++)
46 for (
int i= 1;
i < sizePoints;
i++)
64bool isLess (
int* point1,
int* point2)
66 long area= point1[0]*point2[1]- point1[1]*point2[0];
67 if (area > 0)
return true;
70 return (
abs (point1[0]) +
abs (point1[1]) >
71 abs (point2[0]) +
abs (point2[1]));
80 int* point=
new int [2];
81 point [0]=
points [(lo+hi)/2] [0];
82 point [1]=
points [(lo+hi)/2] [1];
107bool isConvex (
int* point1,
int* point2,
int* point3)
109 long relArea= (point1[0] - point2[0])*(point3[1] - point2[1]) -
110 (point1[1] - point2[1])*(point3[0] - point2[0]);
115 return !(
abs (point1[0] - point3[0]) +
abs (point1[1] - point3[1]) >=
116 (
abs (point2[0] - point1[0]) +
abs (point2[1] - point1[1]) +
117 abs (point2[0] - point3[0]) +
abs (point2[1] - point3[1])));
131 int * minusPoint=
new int [2];
132 minusPoint [0]=
points[0] [0];
133 minusPoint [1]=
points[0] [1];
136 minusPoint[0]= - minusPoint[0];
137 minusPoint[1]= - minusPoint[1];
139 delete [] minusPoint;
141 while (
k < sizePoints)
152 if (
i + 1 <= sizePoints ||
i == sizePoints)
174 if (sizePoints < 3)
return sizePoints;
188 sizeOfOutput=
size (F);
189 int*
result=
new int [sizeOfOutput];
200 int **
points=
new int* [n];
201 for (
int i= 0;
i < n;
i++)
219 for (
int k= 0;
k < bufSize;
k++,
j++)
230merge (
int ** points1,
int sizePoints1,
int ** points2,
int sizePoints2,
234 sizeResult= sizePoints1+sizePoints2;
235 for (
i= 0;
i < sizePoints1;
i++)
237 for (
j= 0;
j < sizePoints2;
j++)
239 if (points1[
i][0] != points2[
j][0])
243 if (points1[
i][1] != points2[
j][1])
257 int **
result=
new int *[sizeResult];
258 for (
i= 0;
i < sizeResult;
i++)
262 for (
i= 0;
i < sizePoints1;
i++,
k++)
267 for (
i= 0;
i < sizePoints2;
i++)
269 if (points2[
i][0] < 0)
285 int **
points=
new int* [sizeF];
286 for (
int i= 0;
i < sizeF;
i++)
294 for (
int k= 0;
k < bufSize;
k++,
j++)
304 int **
result=
new int* [n];
305 for (
int i= 0;
i < n;
i++)
313 for (
int i= 0;
i < sizeF;
i++)
322 int& sizeOfNewtonPoly)
325 int ** pointsF=
new int* [sizeF];
326 for (
int i= 0;
i < sizeF;
i++)
327 pointsF [
i]=
new int [2];
334 for (
int k= 0;
k < bufSize;
k++,
j++)
336 pointsF [
j] [0]=
i.exp();
337 pointsF [
j] [1]=
buf [
k];
343 int ** pointsG=
new int* [sizeG];
344 for (
int i= 0;
i < sizeG;
i++)
345 pointsG [
i]=
new int [2];
350 for (
int k= 0;
k < bufSize;
k++,
j++)
352 pointsG [
j] [0]=
i.exp();
353 pointsG [
j] [1]=
buf [
k];
359 int **
points=
merge (pointsF, sizeF, pointsG, sizeG, sizePoints);
363 int **
result=
new int* [n];
364 for (
int i= 0;
i < n;
i++)
372 for (
int i= 0;
i < sizeF;
i++)
373 delete [] pointsF[
i];
375 for (
int i= 0;
i < sizeG;
i++)
376 delete [] pointsG[
i];
385 int **
buf=
new int* [sizePoints + 1];
386 for (
int i= 0;
i < sizePoints;
i++)
388 buf [
i]=
new int [2];
392 buf [sizePoints]=
new int [2];
393 buf [sizePoints] [0]= point [0];
394 buf [sizePoints] [1]= point [1];
395 int sizeBuf= sizePoints + 1;
398 int * minusPoint=
new int [2];
399 minusPoint [0]=
buf[0] [0];
400 minusPoint [1]=
buf[0] [1];
403 minusPoint[0]= - minusPoint[0];
404 minusPoint[1]= - minusPoint[1];
406 delete [] minusPoint;
408 if (
buf [0] [0] == point [0] &&
buf [0] [1] == point [1])
410 for (
int i= 0;
i < sizeBuf;
i++)
416 for (
int i= 1;
i < sizeBuf-1;
i++)
418 if (
buf [
i] [0] == point [0] &&
buf [
i] [1] == point [1])
421 for (
int i= 0;
i < sizeBuf;
i++)
428 if (
buf [sizeBuf - 1] [0] == point [0] &&
buf [sizeBuf-1] [1] == point [1])
430 buf [1] [0]= point [0];
431 buf [1] [1]= point [1];
434 buf [0] [0]=
buf [sizeBuf-2] [0];
435 buf [0] [1]=
buf [sizeBuf-2] [1];
437 for (
int i= 0;
i < sizeBuf;
i++)
442 for (
int i= 0;
i < sizeBuf;
i++)
451 for (
int i= 0;
i < sizePoints;
i++)
457 for (
int i= 0;
i < sizePoints;
i++)
463 for (
int i= 0;
i < sizePoints;
i++)
470 for (
int i= 0;
i < sizePoints;
i++)
479 int& maxDiff,
int& maxSum,
int& maxX,
int& maxY
489 for (
int i= 1;
i < sizePoints;
i++)
494 minSum=
tmin (minSum, sum);
496 maxSum=
tmax (maxSum, sum);
504 mpz_t * tmp=
new mpz_t[4];
506 mpz_init_set (tmp[0],
N[0]);
507 mpz_mul (tmp[0], tmp[0],
M[0]);
508 mpz_addmul (tmp[0],
N[1],
M[2]);
510 mpz_init_set (tmp[1],
N[0]);
511 mpz_mul (tmp[1], tmp[1],
M[1]);
512 mpz_addmul (tmp[1],
N[1],
M[3]);
514 mpz_init_set (tmp[2],
N[2]);
515 mpz_mul (tmp[2], tmp[2],
M[0]);
516 mpz_addmul (tmp[2],
N[3],
M[2]);
518 mpz_init_set (tmp[3],
N[2]);
519 mpz_mul (tmp[3], tmp[3],
M[1]);
520 mpz_addmul (tmp[3],
N[3],
M[3]);
522 mpz_set (
M[0], tmp[0]);
523 mpz_set (
M[1], tmp[1]);
524 mpz_set (
M[2], tmp[2]);
525 mpz_set (
M[3], tmp[3]);
538 mpz_init_set (det,
M[0]);
539 mpz_mul (det, det,
M[3]);
540 mpz_submul (det,
M[1],
M[2]);
543 mpz_init_set (tmp,
M[0]);
544 mpz_divexact (tmp, tmp, det);
545 mpz_set (
M[0],
M[3]);
546 mpz_divexact (
M[0],
M[0], det);
549 mpz_neg (
M[1],
M[1]);
550 mpz_divexact (
M[1],
M[1], det);
551 mpz_neg (
M[2],
M[2]);
552 mpz_divexact (
M[2],
M[2], det);
564 mpz_t u,
v,
g,maxX,maxY;
568 mpz_init_set_si (maxX,
570 mpz_init_set_si (maxY,
572 mpz_gcdext (
g, u,
v, maxX, maxY);
576 mpz_mul (
A[0],
A[0], maxX);
577 mpz_set (
M[2], maxY);
578 mpz_divexact (
M[2],
M[2],
g);
579 mpz_set (
A[1],
M[2]);
580 mpz_neg (
A[1],
A[1]);
581 mpz_mul (
A[1],
A[1], maxX);
585 mpz_set (
M[3], maxX);
586 mpz_divexact (
M[3],
M[3],
g);
592 mpz_set (
M[2], maxY);
593 mpz_divexact (
M[2],
M[2],
g);
594 mpz_neg (
M[2],
M[2]);
595 mpz_set (
M[3], maxX);
596 mpz_divexact (
M[3],
M[3],
g);
604 else if (sizePoints == 1)
606 mpz_set_si (
M[0], 1);
607 mpz_set_si (
M[3], 1);
611 mpz_set_si (
M[0], 1);
612 mpz_set_si (
M[3], 1);
614 mpz_t * Mu=
new mpz_t[4];
615 mpz_init_set_si (Mu[1], 1);
616 mpz_init_set_si (Mu[2], 1);
620 mpz_t * Lambda=
new mpz_t[4];
621 mpz_init_set_si (Lambda[0], 1);
622 mpz_init_set_si (Lambda[1], -1);
623 mpz_init_set_si (Lambda[3], 1);
624 mpz_init (Lambda[2]);
626 mpz_t * InverseLambda=
new mpz_t[4];
627 mpz_init_set_si (InverseLambda[0], 1);
628 mpz_init_set_si (InverseLambda[1], 1);
629 mpz_init_set_si (InverseLambda[3], 1);
630 mpz_init (InverseLambda[2]);
634 int minDiff, minSum, maxDiff, maxSum, maxX, maxY,
b, d,
f,
h;
635 getMaxMin(
points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY);
645 mpz_set (
A[0],
A[1]);
648 getMaxMin(
points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY);
650 d= maxX + maxY - maxSum;
661 mpz_add_ui (
A[0],
A[0], maxY-
f);
663 mpz_add_ui (
A[0],
A[0],
f-maxY);
664 maxX= maxX + maxY -
b -
f;
666 else if (d +
h > maxY)
674 mpz_add_ui (
A[0],
A[0], -
h);
676 mpz_sub_ui (
A[0],
A[0],
h);
677 maxX= maxX + maxY - d -
h;
688 mpz_clear (Lambda[0]);
689 mpz_clear (Lambda[1]);
690 mpz_clear (Lambda[2]);
691 mpz_clear (Lambda[3]);
694 mpz_clear (InverseLambda[0]);
695 mpz_clear (InverseLambda[1]);
696 mpz_clear (InverseLambda[2]);
697 mpz_clear (InverseLambda[3]);
698 delete [] InverseLambda;
709 int ** newtonPolyg=
NULL;
720 mpz_t expX, expY, minExpX, minExpY;
728 mpz_t * exps=
new mpz_t [2*
size (F)];
735 mpz_set (expX,
A[0]);
736 mpz_set (expY,
A[1]);
737 mpz_addmul_ui (expX,
M[1],
i.exp());
738 mpz_addmul_ui (expY,
M[3],
i.exp());
742 mpz_set (minExpX, expX);
743 mpz_set (minExpY, expY);
748 if (mpz_cmp (minExpY, expY) > 0)
749 mpz_set (minExpY, expY);
750 if (mpz_cmp (minExpX, expX) > 0)
751 mpz_set (minExpX, expX);
753 mpz_init_set (exps[
count], expX);
755 mpz_init_set (exps[
count], expY);
763 mpz_set (expX,
A[0]);
764 mpz_addmul_ui (expX,
M[1],
i.exp());
765 mpz_addmul_ui (expX,
M[0],
j.exp());
767 mpz_set (expY,
A[1]);
768 mpz_addmul_ui (expY,
M[3],
i.exp());
769 mpz_addmul_ui (expY,
M[2],
j.exp());
771 mpz_set (minExpX, expX);
772 mpz_set (minExpY, expY);
773 mpz_init_set (exps[
count], expX);
775 mpz_init_set (exps[
count], expY);
782 for (;
j.hasTerms();
j++)
784 mpz_set (expX,
A[0]);
785 mpz_addmul_ui (expX,
M[1],
i.exp());
786 mpz_addmul_ui (expX,
M[0],
j.exp());
788 mpz_set (expY,
A[1]);
789 mpz_addmul_ui (expY,
M[3],
i.exp());
790 mpz_addmul_ui (expY,
M[2],
j.exp());
792 mpz_init_set (exps[
count], expX);
794 mpz_init_set (exps[
count], expY);
797 if (mpz_cmp (minExpY, expY) > 0)
798 mpz_set (minExpY, expY);
799 if (mpz_cmp (minExpX, expX) > 0)
800 mpz_set (minExpX, expX);
805 int mExpX= mpz_get_si (minExpX);
806 int mExpY= mpz_get_si (minExpY);
817 for (;
j.hasTerms();
j++)
836 for (
int i= 0;
i < n;
i++)
837 delete [] newtonPolyg [
i];
838 delete [] newtonPolyg;
842 for(
int tt=exps_maxcount;tt>=0;tt--) mpz_clear(exps[tt]);
859 mpz_t expX, expY, minExpX, minExpY;
865 mpz_t * exps=
new mpz_t [2*
size(F)];
872 mpz_set_si (expX,
i.exp());
873 mpz_sub (expX, expX,
A[0]);
874 mpz_mul (expX, expX, inverseM[0]);
875 mpz_submul (expX, inverseM[1],
A[1]);
877 mpz_set_si (expY,
i.exp());
878 mpz_sub (expY, expY,
A[0]);
879 mpz_mul (expY, expY, inverseM[2]);
880 mpz_submul (expY, inverseM[3],
A[1]);
882 mpz_set (minExpX, expX);
883 mpz_set (minExpY, expY);
885 mpz_init_set (exps[
count], expX);
887 mpz_init_set (exps[
count], expY);
892 for (;
i.hasTerms();
i++)
894 mpz_set_si (expX,
i.exp());
895 mpz_sub (expX, expX,
A[0]);
896 mpz_mul (expX, expX, inverseM[0]);
897 mpz_submul (expX, inverseM[1],
A[1]);
899 mpz_set_si (expY,
i.exp());
900 mpz_sub (expY, expY,
A[0]);
901 mpz_mul (expY, expY, inverseM[2]);
902 mpz_submul (expY, inverseM[3],
A[1]);
904 mpz_init_set (exps[
count], expX);
906 mpz_init_set (exps[
count], expY);
910 if (mpz_cmp (minExpY, expY) > 0)
911 mpz_set (minExpY, expY);
912 if (mpz_cmp (minExpX, expX) > 0)
913 mpz_set (minExpX, expX);
916 int mExpX= mpz_get_si (minExpX);
917 int mExpY= mpz_get_si (minExpY);
919 for (
i= F;
i.hasTerms();
i++)
932 for(
int tt=max_exps;tt>=0;tt--) mpz_clear(exps[tt]);
945 mpz_set_si (expX,
i.exp());
946 mpz_sub (expX, expX,
A[1]);
947 mpz_mul (expX, expX, inverseM[1]);
948 mpz_submul (expX,
A[0], inverseM[0]);
950 mpz_set_si (expY,
i.exp());
951 mpz_sub (expY, expY,
A[1]);
952 mpz_mul (expY, expY, inverseM[3]);
953 mpz_submul (expY,
A[0], inverseM[2]);
957 mpz_set (minExpX, expX);
958 mpz_set (minExpY, expY);
963 if (mpz_cmp (minExpY, expY) > 0)
964 mpz_set (minExpY, expY);
965 if (mpz_cmp (minExpX, expX) > 0)
966 mpz_set (minExpX, expX);
968 mpz_init_set (exps[
count], expX);
970 mpz_init_set (exps[
count], expY);
978 mpz_set_si (expX,
j.exp());
979 mpz_sub (expX, expX,
A[0]);
980 mpz_mul (expX, expX, inverseM[0]);
981 mpz_set_si (tmp,
i.exp());
982 mpz_sub (tmp, tmp,
A[1]);
983 mpz_addmul (expX, tmp, inverseM[1]);
985 mpz_set_si (expY,
j.exp());
986 mpz_sub (expY, expY,
A[0]);
987 mpz_mul (expY, expY, inverseM[2]);
988 mpz_set_si (tmp,
i.exp());
989 mpz_sub (tmp, tmp,
A[1]);
990 mpz_addmul (expY, tmp, inverseM[3]);
992 mpz_set (minExpX, expX);
993 mpz_set (minExpY, expY);
995 mpz_init_set (exps[
count], expX);
997 mpz_init_set (exps[
count], expY);
1005 for (;
j.hasTerms();
j++)
1007 mpz_set_si (expX,
j.exp());
1008 mpz_sub (expX, expX,
A[0]);
1009 mpz_mul (expX, expX, inverseM[0]);
1010 mpz_set_si (tmp,
i.exp());
1011 mpz_sub (tmp, tmp,
A[1]);
1012 mpz_addmul (expX, tmp, inverseM[1]);
1014 mpz_set_si (expY,
j.exp());
1015 mpz_sub (expY, expY,
A[0]);
1016 mpz_mul (expY, expY, inverseM[2]);
1017 mpz_set_si (tmp,
i.exp());
1018 mpz_sub (tmp, tmp,
A[1]);
1019 mpz_addmul (expY, tmp, inverseM[3]);
1021 mpz_init_set (exps[
count], expX);
1023 mpz_init_set (exps[
count], expY);
1027 if (mpz_cmp (minExpY, expY) > 0)
1028 mpz_set (minExpY, expY);
1029 if (mpz_cmp (minExpX, expX) > 0)
1030 mpz_set (minExpX, expX);
1034 int mExpX= mpz_get_si (minExpX);
1035 int mExpY= mpz_get_si (minExpY);
1058 mpz_clear (minExpX);
1059 mpz_clear (minExpY);
1062 for(
int tt=max_exps;tt>=0;tt--) mpz_clear(exps[tt]);
1075 for (
int i= 1;
i < sizeOfPolygon;
i++)
1092 for (
int i= indexY;
i < sizeOfPolygon;
i++)
1105 result=
new int [sizeOfPolygon - indexY];
1106 sizeOfOutput= sizeOfPolygon - indexY;
1107 count= sizeOfPolygon - indexY - 1;
1113 sizeOfOutput=
count;
1128 int sizeOfNewtonPolygon;
1130 if (sizeOfNewtonPolygon == 3)
1133 (newtonPolyg[0][0]==0 || newtonPolyg[1][0]==0 || newtonPolyg[2][0]==0);
1137 (newtonPolyg[0][1]==0 || newtonPolyg[1][1]==0 || newtonPolyg[2][0]==0);
1144 tmp=
gcd (tmp, newtonPolyg[1][0]);
1145 tmp=
gcd (tmp, newtonPolyg[1][1]);
1146 tmp=
gcd (tmp, newtonPolyg[2][0]);
1147 tmp=
gcd (tmp, newtonPolyg[2][1]);
1150 for (
int i= 0;
i < sizeOfNewtonPolygon;
i++)
1151 delete [] newtonPolyg [
i];
1152 delete [] newtonPolyg;
1157 for (
int i= 0;
i < sizeOfNewtonPolygon;
i++)
1158 delete [] newtonPolyg [
i];
1159 delete [] newtonPolyg;
1168 int sizeOfNewtonPolygon;
1188 while (!
g.isOne() &&
i < sizeOfNewtonPolygon)
1190 g=
gcd (
g, newtonPolyg[
i][0]);
1191 g=
gcd (
g, newtonPolyg[
i][1]);
1205 for (
int i= 0;
i < sizeOfNewtonPolygon;
i++)
1206 delete [] newtonPolyg[
i];
1208 delete [] newtonPolyg;
1365 for (
int j= 0;
j < 3;
j++)
Rational abs(const Rational &a)
const CanonicalForm CFMap CFMap & N
static void quickSort(int lo, int hi, int **points)
void lambda(int **points, int sizePoints)
static bool isLess(int *point1, int *point2)
bool modularIrredTest(const CanonicalForm &F)
modular absolute irreducibility test as described in "Modular Las Vegas Algorithms for Polynomial Abs...
void lambdaInverse(int **points, int sizePoints)
void mpz_mat_inv(mpz_t *&M)
static bool isConvex(int *point1, int *point2, int *point3)
int * getRightSide(int **polygon, int sizeOfPolygon, int &sizeOfOutput)
get the y-direction slopes of all edges with positive slope in y-direction of a convex polygon with a...
void convexDense(int **points, int sizePoints, mpz_t *&M, mpz_t *&A)
Algorithm 5 as described in Convex-Dense Bivariate Polynomial Factorization by Berthomieu,...
static void sort(int **points, int sizePoints)
bool modularIrredTestWithShift(const CanonicalForm &F)
modular absolute irreducibility test with shift as described in "Modular Las Vegas Algorithms for Pol...
bool irreducibilityTest(const CanonicalForm &F)
computes the Newton polygon of F and checks if it satisfies the irreducibility criterion from S....
void getMaxMin(int **points, int sizePoints, int &minDiff, int &minSum, int &maxDiff, int &maxSum, int &maxX, int &maxY)
bool absIrredTest(const CanonicalForm &F)
absolute irreducibility test as described in "Modular Las Vegas Algorithms for Polynomial Absolute Fa...
bool isInPolygon(int **points, int sizePoints, int *point)
check if point is inside a polygon described by points
CanonicalForm decompress(const CanonicalForm &F, const mpz_t *inverseM, const mpz_t *A)
decompress a bivariate poly
void tau(int **points, int sizePoints, int k)
int ** merge(int **points1, int sizePoints1, int **points2, int sizePoints2, int &sizeResult)
void mu(int **points, int sizePoints)
void mpz_mat_mul(const mpz_t *N, mpz_t *&M)
int ** getPoints(const CanonicalForm &F, int &n)
static int * getDegrees(const CanonicalForm &F, int &sizeOfOutput)
int grahamScan(int **points, int sizePoints)
CanonicalForm compress(const CanonicalForm &F, mpz_t *&M, mpz_t *&A, bool computeMA)
compress a bivariate poly
int polygon(int **points, int sizePoints)
compute a polygon
static void translate(int **points, int *point, int sizePoints)
static int smallestPointIndex(int **points, int sizePoints)
This file provides functions to compute the Newton polygon of a bivariate polynomial.
CanonicalForm bCommonDen(const CanonicalForm &f)
CanonicalForm bCommonDen ( const CanonicalForm & f )
CanonicalForm maxNorm(const CanonicalForm &f)
CanonicalForm maxNorm ( const CanonicalForm & f )
declarations of higher level algorithms.
CFFList FACTORY_PUBLIC factorize(const CanonicalForm &f, bool issqrfree=false)
factorization over or
#define ASSERT(expression, message)
static const int SW_RATIONAL
set to 1 for computations over Q
#define GaloisFieldDomain
Interface to generate InternalCF's over various domains from intrinsic types or mpz_t's.
Iterators for CanonicalForm's.
int cf_getNumSmallPrimes()
int cf_getSmallPrime(int i)
generate random evaluation points
class to iterate through CanonicalForm's
generate random elements in F_p
class to generate random evaluation points
factory's class for variables
const CanonicalForm int const CFList const Variable & y
REvaluation E(1, terms.length(), IntRandom(25))
const Variable & v
< [in] a sqrfree bivariate poly
static int min(int a, int b)
some useful template functions.
template CanonicalForm tmax(const CanonicalForm &, const CanonicalForm &)
template CanonicalForm tmin(const CanonicalForm &, const CanonicalForm &)
Operations in GF, where GF is a finite field of size less than 2^16 represented by a root of Conway p...
STATIC_VAR coordinates * points
static BOOLEAN length(leftv result, leftv arg)
STATIC_VAR gmp_float * diff
static int index(p_Length length, p_Ord ord)
int status int void size_t count
int status int void * buf