]> www.ginac.de Git - ginac.git/blob - ginac/polynomial/pgcd.cpp
1 #include "pgcd.h"
2 #include "collect_vargs.h"
3 #include "smod_helpers.h"
4 #include "euclid_gcd_wrap.h"
5 #include "eval_point_finder.h"
6 #include "newton_interpolate.h"
7 #include "divide_in_z_p.h"
9 namespace GiNaC
10 {
12 extern void
13 primpart_content(ex& pp, ex& c, ex e, const exvector& vars, const long p);
15 // Computes the GCD of two polynomials over a prime field.
16 // Based on Algorithm 7.2 from "Algorithms for Computer Algebra"
17 // A and B are considered as Z_p[x_n][x_0, \ldots, x_{n-1}], that is,
18 // as a polynomials in variables x_0, \ldots x_{n-1} having coefficients
19 // from the ring Z_p[x_n]
20 ex pgcd(const ex& A, const ex& B, const exvector& vars, const long p)
21 {
22         static const ex ex1(1);
23         if (A.is_zero())
24                 return B;
26         if (B.is_zero())
27                 return A;
29         if (is_a<numeric>(A))
30                 return ex1;
32         if (is_a<numeric>(B))
33                 return ex1;
35         // Checks for univariate polynomial
36         if (vars.size() == 1) {
37                 ex ret = euclid_gcd(A, B, vars, p); // Univariate GCD
38                 return ret;
39         }
40         const ex& mainvar(vars.back());
42         // gcd of the contents
43         ex H = 0, Hprev = 0; // GCD candidate
44         ex newton_poly = 1;  // for Newton Interpolation
46         // Contents and primparts of A and B
47         ex Aprim, contA;
48         primpart_content(Aprim, contA, A, vars, p);
49         ex Bprim, contB;
50         primpart_content(Bprim, contB, B, vars, p);
51         // gcd of univariate polynomials
52         const ex cont_gcd = euclid_gcd(contA, contB, mainvar, p);
54         exvector restvars = vars;
55         restvars.pop_back();
56         const ex AL = lcoeff_wrt(Aprim, restvars);
57         const ex BL = lcoeff_wrt(Bprim, restvars);
58         // gcd of univariate polynomials
59         const ex lc_gcd = euclid_gcd(AL, BL, mainvar, p);
61         // The estimate of degree of the gcd of Ab and Bb
62         int gcd_deg = std::min(degree(Aprim, mainvar), degree(Bprim, mainvar));
63         eval_point_finder::value_type b;
65         eval_point_finder find_eval_point(p);
66         const numeric pn(p);
67         do {
68                 // Find a `good' evaluation point b.
69                 bool has_more_pts = find_eval_point(b, lc_gcd, mainvar);
70                 // If there are no more possible evaluation points, bail out
71                 if (!has_more_pts)
72                         throw pgcd_failed();
74                 const numeric bn(b);
75                 // Evaluate the polynomials in b
76                 ex Ab = Aprim.subs(mainvar == bn).smod(pn);
77                 ex Bb = Bprim.subs(mainvar == bn).smod(pn);
78                 ex Cb = pgcd(Ab, Bb, restvars, p);
80                 // Set the correct the leading coefficient
81                 const cln::cl_I lcb_gcd =
82                         smod(to_cl_I(lc_gcd.subs(mainvar == bn)), p);
83                 const cln::cl_I Cblc = integer_lcoeff(Cb, restvars);
84                 const cln::cl_I correct_lc = smod(lcb_gcd*recip(Cblc, p), p);
85                 Cb = (Cb*numeric(correct_lc)).smod(pn);
87                 // Test for unlucky homomorphisms
88                 const int img_gcd_deg = Cb.degree(restvars.back());
89                 if (img_gcd_deg < gcd_deg) {
90                         // The degree decreased, previous homomorphisms were
91                         // bad, so we have to start it all over.
92                         H = Cb;
93                         newton_poly = mainvar - numeric(b);
94                         Hprev = 0;
95                         gcd_deg  = img_gcd_deg;
96                         continue;
97                 }
98                 if (img_gcd_deg > gcd_deg) {
99                         // The degree of images GCD is too high, this
100                         // evaluation point is bad. Skip it.
101                         continue;
102                 }
104                 // Image has the same degree as the previous one
105                 // (or at least not higher than the limit)
106                 Hprev = H;
107                 H = newton_interp(Cb, b, H, newton_poly, mainvar, p);
108                 newton_poly = newton_poly*(mainvar - b);
110                 // try to reduce the number of division tests.
111                 const ex H_lcoeff = lcoeff_wrt(H, restvars);
113                 if (H_lcoeff.is_equal(lc_gcd)) {
114                         if ((Hprev-H).expand().smod(pn).is_zero())
115                                 continue;
116                         ex C /* primitive part of H */, contH /* dummy */;
117                         primpart_content(C, contH, H, vars, p);
118                         // Normalize GCD so that leading coefficient is 1
119                         const cln::cl_I Clc = recip(integer_lcoeff(C, vars), p);
120                         C = (C*numeric(Clc)).expand().smod(pn);
122                         ex dummy1, dummy2;
124                         if (divide_in_z_p(Aprim, C, dummy1, vars, p) &&
125                                         divide_in_z_p(Bprim, C, dummy2, vars, p))
126                                 return (cont_gcd*C).expand().smod(pn);
127                         else if (img_gcd_deg == 0)
128                                 return cont_gcd;
129                         // else continue building the candidate
130                 }
131         } while(true);
132         throw pgcd_failed();
133 }
135 } // namespace GiNaC
{"ft_lang_label":"__label__en","ft_lang_prob":0.5261226,"math_prob":0.90655655,"size":3486,"snap":"2022-40-2023-06","text_gpt3_token_len":1167,"char_repetition_ratio":0.11631247,"word_repetition_ratio":0.0073260074,"special_character_ratio":0.3591509,"punctuation_ratio":0.20617111,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.9922055,"pos_list":[0],"im_url_duplicate_count":[null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2022-09-26T13:41:41Z\",\"WARC-Record-ID\":\"<urn:uuid:15f28f4b-bfcd-4a81-bcdd-a20b49991f75>\",\"Content-Length\":\"43881\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:188b3900-13f2-4106-bbaa-8017ac1d0ac1>\",\"WARC-Concurrent-To\":\"<urn:uuid:38dc9042-950c-4669-beb3-e049c351063c>\",\"WARC-IP-Address\":\"188.68.41.94\",\"WARC-Target-URI\":\"https://www.ginac.de/ginac.git/?p=ginac.git;a=blob;f=ginac/polynomial/pgcd.cpp;hb=45b1e47372090352ac5af655b32473df2abab23b\",\"WARC-Payload-Digest\":\"sha1:DCWIV34RTGLX2BO4WSGZAJNKSORWFHUF\",\"WARC-Block-Digest\":\"sha1:P2LOJC4ILU4BA5IPY5DBAIVGXLSHM5IG\",\"WARC-Identified-Payload-Type\":\"application/xhtml+xml\",\"warc_filename\":\"/cc_download/warc_2022/CC-MAIN-2022-40/CC-MAIN-2022-40_segments_1664030334871.54_warc_CC-MAIN-20220926113251-20220926143251-00610.warc.gz\"}"}