Context Navigation

source: git/factory/facFqBivar.cc @ dd8047

spielwiese

Last change on this file since dd8047 was dd8047, checked in by Martin Lee <martinlee84@…>, 13 years ago
fix in bivariate factorization git-svn-id: file:///usr/local/Singular/svn/trunk@14095 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 38.0 KB

Line
1	/*****************************************************************************\
2	* Computer Algebra System SINGULAR
3	\*****************************************************************************/
4	/** @file facFqBivar.cc
5	*
6	* This file provides functions for factorizing a bivariate polynomial over
7	* \f$ F_{p} \f$ , \f$ F_{p}(\alpha ) \f$ or GF, based on "Modern Computer
8	* Algebra, Chapter 15" by J. von zur Gathen & J. Gerhard and "Factoring
9	* multivariate polynomials over a finite field" by L. Bernardin.
10	*
11	* ABSTRACT: In contrast to biFactorizer() in facFqFactorice.cc we evaluate and
12	* factorize the polynomial in both variables. So far factor recombination is
13	* done naive!
14	*
15	* @author Martin Lee
16	*
17	* @internal @version \$Id$
18	*
19	**/
20	/*****************************************************************************/
21
22	#include <config.h>
23
24	#include "assert.h"
25	#include "debug.h"
26	#include "timing.h"
27
28	#include "canonicalform.h"
29	#include "cf_defs.h"
30	#include "cf_map_ext.h"
31	#include "cf_random.h"
32	#include "facHensel.h"
33	#include "cf_map.h"
34	#include "cf_gcd_smallp.h"
35	#include "facFqBivar.h"
36
37
38	#ifdef HAVE_NTL
39	#include "NTLconvert.h"
40
41	TIMING_DEFINE_PRINT(fac_uni_factorizer);
42	TIMING_DEFINE_PRINT(fac_hensel_lift);
43	TIMING_DEFINE_PRINT(fac_factor_recombination);
44
45	CanonicalForm prodMod0 (const CFList& L, const CanonicalForm& M)
46	{
47	if (L.isEmpty())
48	return 1;
49	else if (L.length() == 1)
50	return mod (L.getFirst()(0, 1) , M);
51	else if (L.length() == 2)
52	return mod (L.getFirst()(0, 1)*L.getLast()(0, 1), M);
53	else
54	{
55	int l= L.length()/2;
56	CFListIterator i= L;
57	CFList tmp1, tmp2;
58	CanonicalForm buf1, buf2;
59	for (int j= 1; j <= l; j++, i++)
60	tmp1.append (i.getItem());
61	tmp2= Difference (L, tmp1);
62	buf1= prodMod0 (tmp1, M);
63	buf2= prodMod0 (tmp2, M);
64	return mod (buf1*buf2, M);
65	}
66	}
67
68	CanonicalForm evalPoint (const CanonicalForm& F, CanonicalForm & eval,
69	const Variable& alpha, CFList& list, const bool& GF,
70	bool& fail)
71	{
72	fail= false;
73	Variable x= Variable(2);
74	FFRandom genFF;
75	GFRandom genGF;
76	CanonicalForm random, mipo;
77	double bound;
78	int p= getCharacteristic ();
79	if (alpha.level() != 1)
80	{
81	mipo= getMipo (alpha);
82	int d= degree (mipo);
83	bound= ipower (p, d);
84	}
85	else if (GF)
86	{
87	int d= getGFDegree();
88	bound= ipower (p, d);
89	}
90	else
91	bound= p;
92
93	do
94	{
95	if (list.length() >= bound)
96	{
97	fail= true;
98	break;
99	}
100	if (list.isEmpty())
101	random= 0;
102	else if (GF)
103	random= genGF.generate();
104	else if (list.length() < p \|\| alpha == x)
105	random= genFF.generate();
106	else if (alpha != x && list.length() >= p)
107	{
108	AlgExtRandomF genAlgExt (alpha);
109	random= genAlgExt.generate();
110	}
111	if (find (list, random)) continue;
112	eval= F (random, x);
113	if (degree (eval) != degree (F, Variable (1)))
114	{ //leading coeff vanishes
115	if (!find (list, random))
116	list.append (random);
117	continue;
118	}
119	if (degree (gcd (deriv (eval, eval.mvar()), eval), eval.mvar()) > 0)
120	{ //evaluated polynomial is not squarefree
121	if (!find (list, random))
122	list.append (random);
123	continue;
124	}
125	} while (find (list, random));
126
127	return random;
128	}
129
130	CFList
131	uniFactorizer (const CanonicalForm& A, const Variable& alpha, const bool& GF)
132	{
133	Variable x= A.mvar();
134	ASSERT (A.isUnivariate(), "univariate polynomial expected");
135	CFFList factorsA;
136	ZZ p= to_ZZ (getCharacteristic());
137	ZZ_p::init (p);
138	if (GF)
139	{
140	Variable beta= rootOf (gf_mipo);
141	int k= getGFDegree();
142	char cGFName= gf_name;
143	setCharacteristic (getCharacteristic());
144	CanonicalForm buf= GF2FalphaRep (A, beta);
145	if (getCharacteristic() > 2)
146	{
147	ZZ_pX NTLMipo= convertFacCF2NTLZZpX (gf_mipo);
148	ZZ_pE::init (NTLMipo);
149	ZZ_pEX NTLA= convertFacCF2NTLZZ_pEX (buf, NTLMipo);
150	MakeMonic (NTLA);
151	vec_pair_ZZ_pEX_long NTLFactorsA= CanZass (NTLA);
152	ZZ_pE multi= to_ZZ_pE (1);
153	factorsA= convertNTLvec_pair_ZZpEX_long2FacCFFList (NTLFactorsA, multi,
154	x, beta);
155	}
156	else
157	{
158	GF2X NTLMipo= convertFacCF2NTLGF2X (gf_mipo);
159	GF2E::init (NTLMipo);
160	GF2EX NTLA= convertFacCF2NTLGF2EX (buf, NTLMipo);
161	MakeMonic (NTLA);
162	vec_pair_GF2EX_long NTLFactorsA= CanZass (NTLA);
163	GF2E multi= to_GF2E (1);
164	factorsA= convertNTLvec_pair_GF2EX_long2FacCFFList (NTLFactorsA, multi,
165	x, beta);
166	}
167	setCharacteristic (getCharacteristic(), k, cGFName);
168	for (CFFListIterator i= factorsA; i.hasItem(); i++)
169	{
170	buf= i.getItem().factor();
171	buf= Falpha2GFRep (buf);
172	i.getItem()= CFFactor (buf, i.getItem().exp());
173	}
174	}
175	else if (alpha.level() != 1)
176	{
177	if (getCharacteristic() > 2)
178	{
179	ZZ_pX NTLMipo= convertFacCF2NTLZZpX (getMipo (alpha));
180	ZZ_pE::init (NTLMipo);
181	ZZ_pEX NTLA= convertFacCF2NTLZZ_pEX (A, NTLMipo);
182	MakeMonic (NTLA);
183	vec_pair_ZZ_pEX_long NTLFactorsA= CanZass (NTLA);
184	ZZ_pE multi= to_ZZ_pE (1);
185	factorsA= convertNTLvec_pair_ZZpEX_long2FacCFFList (NTLFactorsA, multi,
186	x, alpha);
187	}
188	else
189	{
190	GF2X NTLMipo= convertFacCF2NTLGF2X (getMipo (alpha));
191	GF2E::init (NTLMipo);
192	GF2EX NTLA= convertFacCF2NTLGF2EX (A, NTLMipo);
193	MakeMonic (NTLA);
194	vec_pair_GF2EX_long NTLFactorsA= CanZass (NTLA);
195	GF2E multi= to_GF2E (1);
196	factorsA= convertNTLvec_pair_GF2EX_long2FacCFFList (NTLFactorsA, multi,
197	x, alpha);
198	}
199	}
200	else
201	{
202	if (getCharacteristic() > 2)
203	{
204	ZZ_pX NTLA= convertFacCF2NTLZZpX (A);
205	MakeMonic (NTLA);
206	vec_pair_ZZ_pX_long NTLFactorsA= CanZass (NTLA);
207	ZZ_p multi= to_ZZ_p (1);
208	factorsA= convertNTLvec_pair_ZZpX_long2FacCFFList (NTLFactorsA, multi,
209	x);
210	}
211	else
212	{
213	GF2X NTLA= convertFacCF2NTLGF2X (A);
214	vec_pair_GF2X_long NTLFactorsA= CanZass (NTLA);
215	GF2 multi= to_GF2 (1);
216	factorsA= convertNTLvec_pair_GF2X_long2FacCFFList (NTLFactorsA, multi,
217	x);
218	}
219	}
220	CFList uniFactors;
221	for (CFFListIterator i= factorsA; i.hasItem(); i++)
222	uniFactors.append (i.getItem().factor());
223	return uniFactors;
224	}
225
226	/// naive factor recombination as decribed in "Factoring
227	/// multivariate polynomials over a finite field" by L Bernardin.
228	CFList
229	extFactorRecombination (const CFList& factors, const CanonicalForm& F,
230	const CanonicalForm& M, const ExtensionInfo& info,
231	const DegreePattern& degs, const CanonicalForm& eval)
232	{
233	if (factors.length() == 0)
234	return CFList();
235
236	Variable alpha= info.getAlpha();
237	Variable beta= info.getBeta();
238	CanonicalForm gamma= info.getGamma();
239	CanonicalForm delta= info.getDelta();
240	int k= info.getGFDegree();
241
242	Variable y= F.mvar();
243	CFList source, dest;
244	if (degs.getLength() <= 1 \|\| factors.length() == 1)
245	return CFList(mapDown (F(y-eval, y), info, source, dest));
246
247	DEBOUTLN (cerr, "LC (F, 1)*prodMod (factors, M) == F " <<
248	(LC (F, 1)*prodMod (factors, M) == F));
249	int degMipoBeta;
250	if (!k && beta.level() == 1)
251	degMipoBeta= 1;
252	else if (!k && beta.level() != 1)
253	degMipoBeta= degree (getMipo (beta));
254
255	CFList T, S, Diff;
256	T= factors;
257
258	int s= 1;
259	CFList result;
260	CanonicalForm buf, buf2;
261	if (beta != Variable (1))
262	buf= mapDown (F, gamma, delta, alpha, source, dest);
263	else
264	buf= F;
265
266	CanonicalForm g, LCBuf= LC (buf, Variable (1));
267	int v [T.length()];
268	for (int i= 0; i < T.length(); i++)
269	v[i]= 0;
270
271	CFArray TT;
272	DegreePattern bufDegs1, bufDegs2;
273	bufDegs1= degs;
274	int subsetDeg;
275	TT= copy (factors);
276	bool nosubset= false;
277	bool recombination= false;
278	bool trueFactor= false;
279	while (T.length() >= 2*s)
280	{
281	while (nosubset == false)
282	{
283	if (T.length() == s)
284	{
285	if (recombination)
286	{
287	T.insert (LCBuf);
288	g= prodMod (T, M);
289	T.removeFirst();
290	g /= content(g);
291	g= g (y - eval, y);
292	g /= Lc (g);
293	appendTestMapDown (result, g, info, source, dest);
294	return result;
295	}
296	else
297	{
298	appendMapDown (result, F (y - eval, y), info, source, dest);
299	return result;
300	}
301	}
302	S= subset (v, s, TT, nosubset);
303	if (nosubset) break;
304	subsetDeg= subsetDegree (S);
305	// skip those combinations that are not possible
306	if (!degs.find (subsetDeg))
307	continue;
308	else
309	{
310	g= prodMod0 (S, M);
311	g= mod (g*LCBuf, M);
312	g /= content (g);
313	if (fdivides (LC (g), LCBuf))
314	{
315	S.insert (LCBuf);
316	g= prodMod (S, M);
317	S.removeFirst();
318	g /= content (g, Variable (1));
319	if (fdivides (g, buf))
320	{
321	buf2= g (y - eval, y);
322	buf2 /= Lc (buf2);
323
324	if (!k && beta == Variable (1))
325	{
326	if (degree (buf2, alpha) < degMipoBeta)
327	{
328	buf /= g;
329	LCBuf= LC (buf, Variable (1));
330	recombination= true;
331	appendTestMapDown (result, buf2, info, source, dest);
332	trueFactor= true;
333	}
334	}
335	else
336	{
337	if (!isInExtension (buf2, delta, k))
338	{
339	buf /= g;
340	LCBuf= LC (buf, Variable (1));
341	recombination= true;
342	appendTestMapDown (result, buf2, info, source, dest);
343	trueFactor= true;
344	}
345	}
346	if (trueFactor)
347	{
348	T= Difference (T, S);
349	// compute new possible degree pattern
350	bufDegs2= DegreePattern (T);
351	bufDegs1.intersect (bufDegs2);
352	bufDegs1.refine ();
353	if (T.length() < 2*s \|\| T.length() == s \|\|
354	bufDegs1.getLength() == 1)
355	{
356	if (recombination)
357	{
358	appendTestMapDown (result, buf (y - eval, y), info, source,
359	dest);
360	return result;
361	}
362	else
363	{
364	appendMapDown (result, F (y - eval, y), info, source, dest);
365	return result;
366	}
367	}
368	trueFactor= false;
369	TT= copy (T);
370	indexUpdate (v, s, T.length(), nosubset);
371	if (nosubset) break;
372	}
373	}
374	}
375	}
376	}
377	s++;
378	if (T.length() < 2*s \|\| T.length() == s)
379	{
380	if (recombination)
381	{
382	appendTestMapDown (result, buf (y - eval, y), info, source, dest);
383	return result;
384	}
385	else
386	{
387	appendMapDown (result, F (y - eval, y), info, source, dest);
388	return result;
389	}
390	}
391	int v [T.length()];
392	for (int i= 0; i < T.length(); i++)
393	v[i]= 0;
394	nosubset= false;
395	}
396	if (T.length() < 2*s)
397	appendMapDown (result, F (y - eval, y), info, source, dest);
398
399	return result;
400	}
401
402	/// naive factor recombination as decribed in "Factoring
403	/// multivariate polynomials over a finite field" by L Bernardin.
404	CFList
405	factorRecombination (const CFList& factors, const CanonicalForm& F,
406	const CanonicalForm& M, const DegreePattern& degs)
407	{
408	if (factors.length() == 0)
409	return CFList ();
410	if (degs.getLength() <= 1 \|\| factors.length() == 1)
411	return CFList(F);
412	DEBOUTLN (cerr, "LC (F, 1)*prodMod (factors, M) == F " <<
413	(LC (F, 1)*prodMod (factors, M) == F));
414	CFList T, S;
415
416	T= factors;
417	int s= 1;
418	CFList result;
419	CanonicalForm LCBuf= LC (F, Variable (1));
420	CanonicalForm g, buf= F;
421	int v [T.length()];
422	for (int i= 0; i < T.length(); i++)
423	v[i]= 0;
424	bool nosubset= false;
425	CFArray TT;
426	DegreePattern bufDegs1, bufDegs2;
427	bufDegs1= degs;
428	int subsetDeg;
429	TT= copy (factors);
430	bool recombination= false;
431	while (T.length() >= 2*s)
432	{
433	while (nosubset == false)
434	{
435	if (T.length() == s)
436	{
437	if (recombination)
438	{
439	T.insert (LCBuf);
440	g= prodMod (T, M);
441	T.removeFirst();
442	result.append (g/content (g, Variable (1)));
443	return result;
444	}
445	else
446	return CFList (F);
447	}
448	S= subset (v, s, TT, nosubset);
449	if (nosubset) break;
450	subsetDeg= subsetDegree (S);
451	// skip those combinations that are not possible
452	if (!degs.find (subsetDeg))
453	continue;
454	else
455	{
456	g= prodMod0 (S, M);
457	g= mod (g*LCBuf, M);
458	g /= content (g);
459	if (fdivides (LC(g), LCBuf))
460	{
461	S.insert (LCBuf);
462	g= prodMod (S, M);
463	S.removeFirst();
464	g /= content (g, Variable (1));
465	if (fdivides (g, buf))
466	{
467	recombination= true;
468	result.append (g);
469	buf /= g;
470	LCBuf= LC (buf, Variable(1));
471	T= Difference (T, S);
472
473	// compute new possible degree pattern
474	bufDegs2= DegreePattern (T);
475	bufDegs1.intersect (bufDegs2);
476	bufDegs1.refine ();
477	if (T.length() < 2*s \|\| T.length() == s \|\|
478	bufDegs1.getLength() == 1)
479	{
480	if (recombination)
481	{
482	result.append (buf);
483	return result;
484	}
485	else
486	return CFList (F);
487	}
488	TT= copy (T);
489	indexUpdate (v, s, T.length(), nosubset);
490	if (nosubset) break;
491	}
492	}
493	}
494	}
495	s++;
496	if (T.length() < 2*s \|\| T.length() == s)
497	{
498	if (recombination)
499	{
500	result.append (buf);
501	return result;
502	}
503	else
504	return CFList (F);
505	}
506	int v [T.length()];
507	for (int i= 0; i < T.length(); i++)
508	v[i]= 0;
509	nosubset= false;
510	}
511	if (T.length() < 2*s)
512	result.append (F);
513
514	return result;
515	}
516
517	Variable chooseExtension (const CanonicalForm & A, const Variable & alpha)
518	{
519	int p= getCharacteristic();
520	ZZ NTLp= to_ZZ (p);
521	ZZ_p::init (NTLp);
522	ZZ_pX NTLirredpoly;
523	int d= degree (A);
524	int m= 1;
525	if (alpha != Variable (1))
526	m= degree (getMipo (alpha));
527	int i= (int) floor((double) d/(double) m) - 1;
528	if (i < 2)
529	i= 2;
530	if (i > 8)
531	i= i/4;
532	BuildIrred (NTLirredpoly, i*m);
533	Variable x= A.mvar();
534	CanonicalForm newMipo= convertNTLZZpX2CF (NTLirredpoly, x);
535	return rootOf (newMipo);
536	}
537
538	CFList
539	earlyFactorDetection (CanonicalForm& F, CFList& factors,int& adaptedLiftBound,
540	DegreePattern& degs, bool& success, int deg)
541	{
542	DegreePattern bufDegs1= degs;
543	DegreePattern bufDegs2;
544	CFList result;
545	CFList T= factors;
546	CanonicalForm buf= F;
547	CanonicalForm LCBuf= LC (buf, Variable (1));
548	CanonicalForm g;
549	CanonicalForm M= power (F.mvar(), deg);
550	adaptedLiftBound= 0;
551	int d;
552	if (degree (LCBuf) == degree (F))
553	d= degree (F);
554	else
555	d= degree (F) + degree (LCBuf);
556	for (CFListIterator i= factors; i.hasItem(); i++)
557	{
558	if (!bufDegs1.find (degree (i.getItem(), 1)))
559	continue;
560	else
561	{
562	g= i.getItem() (0, 1);
563	g *= LCBuf;
564	g= mod (g, M);
565	if (fdivides (LC (g), LCBuf))
566	{
567	g= mulMod2 (i.getItem(), LCBuf, M);
568	g /= content (g, Variable (1));
569	if (fdivides (g, buf))
570	{
571	result.append (g);
572	buf /= g;
573	d -= degree (g) + degree (LC (g, Variable (1)));
574	LCBuf= LC (buf, Variable (1));
575	T= Difference (T, CFList (i.getItem()));
576
577	// compute new possible degree pattern
578	bufDegs2= DegreePattern (T);
579	bufDegs1.intersect (bufDegs2);
580	bufDegs1.refine ();
581	if (bufDegs1.getLength() <= 1)
582	{
583	result.append (buf);
584	break;
585	}
586	}
587	}
588	}
589	}
590	adaptedLiftBound= d + 1;
591	if (d < deg)
592	{
593	factors= T;
594	degs= bufDegs1;
595	F= buf;
596	success= true;
597	}
598	if (bufDegs1.getLength() <= 1)
599	degs= bufDegs1;
600	return result;
601	}
602
603	CFList
604	extEarlyFactorDetection (CanonicalForm& F, CFList& factors,
605	int& adaptedLiftBound, DegreePattern& degs,
606	bool& success, const ExtensionInfo& info,
607	const CanonicalForm& eval, int deg)
608	{
609	Variable alpha= info.getAlpha();
610	Variable beta= info.getBeta();
611	CanonicalForm gamma= info.getGamma();
612	CanonicalForm delta= info.getDelta();
613	int k= info.getGFDegree();
614	DegreePattern bufDegs1= degs, bufDegs2;
615	CFList result;
616	CFList T= factors;
617	Variable y= F.mvar();
618	CanonicalForm buf= F, LCBuf= LC (buf, Variable (1)), g, buf2;
619	CanonicalForm M= power (y, deg);
620	adaptedLiftBound= 0;
621	bool trueFactor= false;
622	int d;
623	if (degree (F) == degree (LCBuf))
624	d= degree (F);
625	else
626	d= degree (F) + degree (LCBuf);
627	CFList source, dest;
628	int degMipoBeta;
629	if (!k && beta.level() == 1)
630	degMipoBeta= 1;
631	else if (!k && beta.level() != 1)
632	degMipoBeta= degree (getMipo (beta));
633	for (CFListIterator i= factors; i.hasItem(); i++)
634	{
635	if (!bufDegs1.find (degree (i.getItem(), 1)))
636	continue;
637	else
638	{
639	g= i.getItem() (0, 1);
640	g *= LCBuf;
641	g= mod (g, M);
642	if (fdivides (LC (g), LCBuf))
643	{
644	g= mulMod2 (i.getItem(), LCBuf, M);
645	g /= content (g, Variable (1));
646	if (fdivides (g, buf))
647	{
648	buf2= g (y - eval, y);
649	buf2 /= Lc (buf2);
650
651	if (!k && beta == Variable (1))
652	{
653	if (degree (buf2, alpha) < degMipoBeta)
654	{
655	appendTestMapDown (result, buf2, info, source, dest);
656	buf /= g;
657	d -= degree (g) + degree (LC (g, Variable (1)));
658	LCBuf= LC (buf, Variable (1));
659	trueFactor= true;
660	}
661	}
662	else
663	{
664	if (!isInExtension (buf2, delta, k))
665	{
666	appendTestMapDown (result, buf2, info, source, dest);
667	buf /= g;
668	d -= degree (g) + degree (LC (g, Variable (1)));
669	LCBuf= LC (buf, Variable (1));
670	trueFactor= true;
671	}
672	}
673	if (trueFactor)
674	{
675	T= Difference (T, CFList (i.getItem()));
676
677	// compute new possible degree pattern
678	bufDegs2= DegreePattern (T);
679	bufDegs1.intersect (bufDegs2);
680	bufDegs1.refine ();
681	trueFactor= false;
682	if (bufDegs1.getLength() <= 1)
683	{
684	buf= buf (y - eval, y);
685	buf /= Lc (buf);
686	appendMapDown (result, buf, info, source, dest);
687	break;
688	}
689	}
690	}
691	}
692	}
693	}
694	adaptedLiftBound= d + 1;
695	if (adaptedLiftBound < deg)
696	{
697	success= true;
698	factors= T;
699	degs= bufDegs1;
700	F= buf;
701	}
702	if (bufDegs1.getLength() <= 1)
703	degs= bufDegs1;
704
705	return result;
706	}
707
708	CFList
709	henselLiftAndEarly (CanonicalForm& A, bool& earlySuccess, CFList&
710	earlyFactors, DegreePattern& degs, int& liftBound,
711	const CFList& uniFactors, const ExtensionInfo& info,
712	const CanonicalForm& eval)
713	{
714	Variable alpha= info.getAlpha();
715	Variable beta= info.getBeta();
716	CanonicalForm gamma= info.getGamma();
717	CanonicalForm delta= info.getDelta();
718	int k= info.getGFDegree();
719	bool extension= info.isInExtension();
720
721	Variable x= Variable (1);
722	Variable y= Variable (2);
723	CFArray Pi;
724	CFList diophant;
725	CFList bufUniFactors= uniFactors;
726	bufUniFactors.insert (LC (A, x));
727	CFMatrix M= CFMatrix (liftBound, bufUniFactors.length() - 1);
728	earlySuccess= false;
729	int newLiftBound= 0;
730	int smallFactorDeg= 11; //this is a tunable parameter
731	if (smallFactorDeg >= liftBound)
732	henselLift12 (A, bufUniFactors, liftBound, Pi, diophant, M);
733	else if (smallFactorDeg >= degree (A, y) + 1)
734	{
735	henselLift12 (A, bufUniFactors, degree (A, y) + 1, Pi, diophant, M);
736	if (!extension)
737	earlyFactors= earlyFactorDetection (A, bufUniFactors, newLiftBound,
738	degs, earlySuccess, degree (A, y) + 1);
739	else
740	earlyFactors= extEarlyFactorDetection (A, bufUniFactors,
741	newLiftBound, degs, earlySuccess, info, eval,
742	degree (A, y) + 1);
743	if (degs.getLength() > 1 && !earlySuccess)
744	{
745	if (newLiftBound > degree (A, y) + 1)
746	{
747	liftBound= newLiftBound;
748	bufUniFactors.insert (LC(A, x));
749	henselLiftResume12 (A, bufUniFactors, degree (A, y) + 1, liftBound,
750	Pi, diophant, M);
751	}
752	}
753	else if (earlySuccess)
754	liftBound= newLiftBound;
755	}
756	else if (smallFactorDeg < degree (A, y) + 1)
757	{
758	henselLift12 (A, bufUniFactors, smallFactorDeg, Pi, diophant, M);
759	if (!extension)
760	earlyFactors= earlyFactorDetection (A, bufUniFactors, newLiftBound,
761	degs, earlySuccess,
762	smallFactorDeg);
763	else
764	earlyFactors= extEarlyFactorDetection (A, bufUniFactors,
765	newLiftBound, degs, earlySuccess,
766	info, eval, smallFactorDeg);
767	if (degs.getLength() > 1 && !earlySuccess)
768	{
769	bufUniFactors.insert (LC (A, x));
770	henselLiftResume12 (A, bufUniFactors, smallFactorDeg, degree (A, y)
771	+ 1, Pi, diophant, M);
772	if (!extension)
773	earlyFactors= earlyFactorDetection (A, bufUniFactors, newLiftBound,
774	degs, earlySuccess, degree (A, y) + 1);
775	else
776	earlyFactors= extEarlyFactorDetection (A, bufUniFactors,
777	newLiftBound, degs, earlySuccess,
778	info, eval, degree(A,y) + 1);
779	if (degs.getLength() > 1 && !earlySuccess)
780	{
781	if (newLiftBound > degree (A, y) + 1)
782	{
783	if (newLiftBound < newLiftBound)
784	liftBound= newLiftBound;
785	bufUniFactors.insert (LC(A, x));
786	henselLiftResume12 (A, bufUniFactors, degree (A, y) + 1, liftBound,
787	Pi, diophant, M);
788	}
789	}
790	else if (earlySuccess)
791	liftBound= newLiftBound;
792	}
793	else if (earlySuccess)
794	liftBound= newLiftBound;
795	}
796	if (newLiftBound > 0)
797	liftBound= newLiftBound;
798	return bufUniFactors;
799	}
800
801	CFList
802	extBiFactorize (const CanonicalForm& F, const ExtensionInfo& info);
803
804	/// bivariate factorization over finite fields as decribed in "Factoring
805	/// multivariate polynomials over a finite field" by L Bernardin.
806	CFList
807	biFactorize (const CanonicalForm& F, const ExtensionInfo& info)
808	{
809	if (F.inCoeffDomain())
810	return CFList(F);
811
812	CanonicalForm A= F;
813	bool GF= (CFFactory::gettype() == GaloisFieldDomain);
814
815	Variable alpha= info.getAlpha();
816	Variable beta= info.getBeta();
817	CanonicalForm gamma= info.getGamma();
818	CanonicalForm delta= info.getDelta();
819	int k= info.getGFDegree();
820	bool extension= info.isInExtension();
821	if (A.isUnivariate())
822	{
823	if (extension == false)
824	return uniFactorizer (F, alpha, GF);
825	else
826	{
827	CFList source, dest;
828	A= mapDown (A, info, source, dest);
829	return uniFactorizer (A, beta, GF);
830	}
831	}
832
833	CFMap N;
834	A= compress (A, N);
835	Variable y= A.mvar();
836
837	if (y.level() > 2) return CFList (F);
838	Variable x= Variable (1);
839
840	//remove and factorize content
841	CanonicalForm contentAx= content (A, x);
842	CanonicalForm contentAy= content (A);
843
844	A= A/(contentAx*contentAy);
845	CFList contentAxFactors, contentAyFactors;
846
847	if (!extension)
848	{
849	contentAxFactors= uniFactorizer (contentAx, alpha, GF);
850	contentAyFactors= uniFactorizer (contentAy, alpha, GF);
851	}
852
853	//trivial case
854	CFList factors;
855	if (A.inCoeffDomain())
856	{
857	append (factors, contentAxFactors);
858	append (factors, contentAyFactors);
859	decompress (factors, N);
860	return factors;
861	}
862	else if (A.isUnivariate())
863	{
864	factors= uniFactorizer (A, alpha, GF);
865	append (factors, contentAxFactors);
866	append (factors, contentAyFactors);
867	decompress (factors, N);
868	return factors;
869	}
870
871	// check derivatives
872	bool derivXZero= false;
873	CanonicalForm derivX= deriv (A, x);
874	CanonicalForm gcdDerivX;
875	if (derivX.isZero())
876	derivXZero= true;
877	else
878	{
879	gcdDerivX= gcd (A, derivX);
880	if (degree (gcdDerivX) > 0)
881	{
882	CanonicalForm g= A/gcdDerivX;
883	CFList factorsG=
884	Union (biFactorize (g, info), biFactorize (gcdDerivX, info));
885	append (factorsG, contentAxFactors);
886	append (factorsG, contentAyFactors);
887	decompress (factorsG, N);
888	normalize (factors);
889	return factorsG;
890	}
891	}
892	bool derivYZero= false;
893	CanonicalForm derivY= deriv (A, y);
894	CanonicalForm gcdDerivY;
895	if (derivY.isZero())
896	derivYZero= true;
897	else
898	{
899	gcdDerivY= gcd (A, derivY);
900	if (degree (gcdDerivY) > 0)
901	{
902	CanonicalForm g= A/gcdDerivY;
903	CFList factorsG=
904	Union (biFactorize (g, info), biFactorize (gcdDerivY, info));
905	append (factorsG, contentAxFactors);
906	append (factorsG, contentAyFactors);
907	decompress (factorsG, N);
908	normalize (factors);
909	return factorsG;
910	}
911	}
912	//main variable is chosen s.t. the degree in x is minimal
913	bool swap= false;
914	if ((degree (A) > degree (A, x)) \|\| derivXZero)
915	{
916	if (!derivYZero)
917	{
918	A= swapvar (A, y, x);
919	swap= derivXZero;
920	derivXZero= derivYZero;
921	derivYZero= swap;
922	swap= true;
923	}
924	}
925
926	bool fail;
927	CanonicalForm Aeval, evaluation, bufAeval, bufEvaluation, buf;
928	CFList uniFactors, list, bufUniFactors;
929	DegreePattern degs;
930	DegreePattern bufDegs;
931
932	bool fail2= false;
933	CanonicalForm Aeval2, evaluation2, bufAeval2, bufEvaluation2;
934	CFList bufUniFactors2, list2, uniFactors2;
935	DegreePattern degs2;
936	DegreePattern bufDegs2;
937	bool swap2= false;
938
939	// several univariate factorizations to obtain more information about the
940	// degree pattern therefore usually less combinations have to be tried during
941	// the recombination process
942	int factorNums= 5;
943	double logarithm= (double) ilog2 (totaldegree (A));
944	logarithm /= log2exp;
945	logarithm= ceil (logarithm);
946	if (factorNums < (int) logarithm)
947	factorNums= (int) logarithm;
948	int subCheck1= substituteCheck (A, x);
949	int subCheck2= substituteCheck (A, y);
950	for (int i= 0; i < factorNums; i++)
951	{
952	bufAeval= A;
953	bufEvaluation= evalPoint (A, bufAeval, alpha, list, GF, fail);
954	if (!derivXZero && !fail2)
955	{
956	buf= swapvar (A, x, y);
957	bufAeval2= buf;
958	bufEvaluation2= evalPoint (buf, bufAeval2, alpha, list2, GF, fail2);
959	}
960	// first try to change main variable if there is no valid evaluation point
961	if (fail && (i == 0))
962	{
963	if (!derivXZero)
964	bufEvaluation= evalPoint (buf, bufAeval, alpha, list, GF, fail);
965	else
966	fail= true;
967
968	if (!fail)
969	{
970	int dummy= subCheck2;
971	subCheck2= subCheck1;
972	subCheck1= dummy;
973	A= buf;
974	swap2= true;
975	}
976	}
977
978	// if there is no valid evaluation point pass to a field extension
979	if (fail && (i == 0))
980	{
981	factors= extBiFactorize (A, info);
982	appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
983	swap, swap2, N);
984	normalize (factors);
985	return factors;
986	}
987
988	// there is at least one valid evaluation point
989	// but we do not compute more univariate factorization over an extension
990	if (fail && (i != 0))
991	break;
992
993	// univariate factorization
994	TIMING_START (fac_uni_factorizer);
995	bufUniFactors= uniFactorizer (bufAeval, alpha, GF);
996	TIMING_END_AND_PRINT (fac_uni_factorizer,
997	"time for univariate factorization: ");
998	DEBOUTLN (cerr, "Lc (bufAeval)*prod (bufUniFactors)== bufAeval " <<
999	prod (bufUniFactors)*Lc (bufAeval) == bufAeval);
1000
1001	if (!derivXZero && !fail2)
1002	{
1003	TIMING_START (fac_uni_factorizer);
1004	bufUniFactors2= uniFactorizer (bufAeval2, alpha, GF);
1005	TIMING_END_AND_PRINT (fac_uni_factorizer,
1006	"time for univariate factorization in y: ");
1007	DEBOUTLN (cerr, "Lc (Aeval2)*prod (uniFactors2)== Aeval2 " <<
1008	prod (bufUniFactors2)*Lc (bufAeval2) == bufAeval2);
1009	}
1010
1011	if (bufUniFactors.length() == 1 \|\|
1012	(!fail2 && !derivXZero && (bufUniFactors2.length() == 1)))
1013	{
1014	if (extension)
1015	{
1016	CFList source, dest;
1017	ExtensionInfo info2= ExtensionInfo (beta, alpha, delta, gamma, k,
1018	info.getGFName(), info.isInExtension());
1019	appendMapDown (factors, A, info2, source, dest);
1020	}
1021	else
1022	factors.append (A);
1023
1024	appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1025	swap, swap2, N);
1026
1027	normalize (factors);
1028	return factors;
1029	}
1030
1031	if (i == 0)
1032	{
1033	if (subCheck1 > 0)
1034	{
1035	int subCheck= substituteCheck (bufUniFactors);
1036
1037	if (subCheck > 1 && (subCheck1%subCheck == 0))
1038	{
1039	CanonicalForm bufA= A;
1040	subst (bufA, bufA, subCheck, x);
1041	factors= biFactorize (bufA, info);
1042	reverseSubst (factors, subCheck, x);
1043	appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1044	swap, swap2, N);
1045	normalize (factors);
1046	return factors;
1047	}
1048	}
1049
1050	if (!derivXZero && !fail2 && subCheck2 > 0)
1051	{
1052	int subCheck= substituteCheck (bufUniFactors2);
1053
1054	if (subCheck > 1 && (subCheck2%subCheck == 0))
1055	{
1056	CanonicalForm bufA= A;
1057	subst (bufA, bufA, subCheck, y);
1058	factors= biFactorize (bufA, info);
1059	reverseSubst (factors, subCheck, y);
1060	appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1061	swap, swap2, N);
1062	normalize (factors);
1063	return factors;
1064	}
1065	}
1066	}
1067
1068	// degree analysis
1069	bufDegs = DegreePattern (bufUniFactors);
1070	if (!derivXZero && !fail2)
1071	bufDegs2= DegreePattern (bufUniFactors2);
1072
1073	if (i == 0)
1074	{
1075	Aeval= bufAeval;
1076	evaluation= bufEvaluation;
1077	uniFactors= bufUniFactors;
1078	degs= bufDegs;
1079	if (!derivXZero && !fail2)
1080	{
1081	Aeval2= bufAeval2;
1082	evaluation2= bufEvaluation2;
1083	uniFactors2= bufUniFactors2;
1084	degs2= bufDegs2;
1085	}
1086	}
1087	else
1088	{
1089	degs.intersect (bufDegs);
1090	if (!derivXZero && !fail2)
1091	{
1092	degs2.intersect (bufDegs2);
1093	if (bufUniFactors2.length() < uniFactors2.length())
1094	{
1095	uniFactors2= bufUniFactors2;
1096	Aeval2= bufAeval2;
1097	evaluation2= bufEvaluation2;
1098	}
1099	}
1100	if (bufUniFactors.length() < uniFactors.length())
1101	{
1102	uniFactors= bufUniFactors;
1103	Aeval= bufAeval;
1104	evaluation= bufEvaluation;
1105	}
1106	}
1107	list.append (bufEvaluation);
1108	if (!derivXZero && !fail2)
1109	list2.append (bufEvaluation2);
1110	}
1111
1112	if (!derivXZero && !fail2)
1113	{
1114	if (uniFactors.length() > uniFactors2.length() \|\|
1115	(uniFactors.length() == uniFactors2.length()
1116	&& degs.getLength() > degs2.getLength()))
1117	{
1118	degs= degs2;
1119	uniFactors= uniFactors2;
1120	evaluation= evaluation2;
1121	Aeval= Aeval2;
1122	A= buf;
1123	swap2= true;
1124	}
1125	}
1126
1127	if (degs.getLength() == 1) // A is irreducible
1128	{
1129	if (extension)
1130	{
1131	CFList source, dest;
1132	ExtensionInfo info2= ExtensionInfo (beta, alpha, delta, gamma, k,
1133	info.getGFName(), info.isInExtension());
1134	appendMapDown (factors, A, info2, source, dest);
1135	}
1136	else
1137	factors.append (A);
1138	appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1139	swap, swap2, N);
1140	normalize (factors);
1141	return factors;
1142	}
1143
1144	A= A (y + evaluation, y);
1145
1146	int liftBound;
1147	if (degree (A, y) == degree (LC (A, x)))
1148	liftBound= degree (LC (A, x)) + 1;
1149	else
1150	liftBound= degree (A, y) + 1 + degree (LC(A, x));
1151
1152	DEBOUTLN (cerr, "uniFactors= " << uniFactors);
1153	bool earlySuccess= false;
1154	CFList earlyFactors;
1155	TIMING_START (fac_hensel_lift);
1156	uniFactors= henselLiftAndEarly
1157	(A, earlySuccess, earlyFactors, degs, liftBound,
1158	uniFactors, info, evaluation);
1159	TIMING_END_AND_PRINT (fac_hensel_lift, "time for hensel lifting: ");
1160	DEBOUTLN (cerr, "lifted factors= " << uniFactors);
1161
1162	CanonicalForm MODl= power (y, liftBound);
1163	TIMING_START (fac_factor_recombination);
1164	if (!extension)
1165	factors= factorRecombination (uniFactors, A, MODl, degs);
1166	else
1167	factors= extFactorRecombination (uniFactors, A, MODl, info, degs,
1168	evaluation);
1169	TIMING_END_AND_PRINT (fac_factor_recombination,
1170	"time for factor recombination: ");
1171	if (earlySuccess)
1172	factors= Union (earlyFactors, factors);
1173	else if (!earlySuccess && degs.getLength() == 1)
1174	factors= earlyFactors;
1175
1176	if (!extension)
1177	{
1178	for (CFListIterator i= factors; i.hasItem(); i++)
1179	i.getItem()= i.getItem() (y - evaluation, y);
1180	}
1181
1182	appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1183	swap, swap2, N);
1184	normalize (factors);
1185
1186	return factors;
1187	}
1188
1189	CFList
1190	extBiFactorize (const CanonicalForm& F, const ExtensionInfo& info)
1191	{
1192
1193	CanonicalForm A= F;
1194	Variable alpha= info.getAlpha();
1195	Variable beta= info.getBeta();
1196	int k= info.getGFDegree();
1197	char cGFName= info.getGFName();
1198	CanonicalForm delta= info.getDelta();
1199
1200	bool GF= (CFFactory::gettype() == GaloisFieldDomain);
1201	Variable x= Variable (1);
1202	CFList factors;
1203	if (!GF && alpha == x) // we are in F_p
1204	{
1205	bool extension= true;
1206	int p= getCharacteristic();
1207	if (p*p < (1<<16)) // pass to GF if possible
1208	{
1209	setCharacteristic (getCharacteristic(), 2, 'Z');
1210	A= A.mapinto();
1211	ExtensionInfo info= ExtensionInfo (extension);
1212	factors= biFactorize (A, info);
1213
1214	Variable vBuf= rootOf (gf_mipo);
1215	setCharacteristic (getCharacteristic());
1216	for (CFListIterator j= factors; j.hasItem(); j++)
1217	j.getItem()= GF2FalphaRep (j.getItem(), vBuf);
1218	}
1219	else // not able to pass to GF, pass to F_p(\alpha)
1220	{
1221	CanonicalForm mipo= randomIrredpoly (3, Variable (1));
1222	Variable v= rootOf (mipo);
1223	ExtensionInfo info= ExtensionInfo (v);
1224	factors= biFactorize (A, info);
1225	}
1226	return factors;
1227	}
1228	else if (!GF && (alpha != x)) // we are in F_p(\alpha)
1229	{
1230	if (k == 1) // need factorization over F_p
1231	{
1232	int extDeg= degree (getMipo (alpha));
1233	extDeg++;
1234	CanonicalForm mipo= randomIrredpoly (extDeg + 1, Variable (1));
1235	Variable v= rootOf (mipo);
1236	ExtensionInfo info= ExtensionInfo (v);
1237	factors= biFactorize (A, info);
1238	}
1239	else
1240	{
1241	if (beta == Variable (1))
1242	{
1243	Variable v= chooseExtension (A, alpha);
1244	CanonicalForm primElem, imPrimElem;
1245	bool primFail= false;
1246	Variable vBuf;
1247	primElem= primitiveElement (alpha, vBuf, primFail);
1248	ASSERT (!primFail, "failure in integer factorizer");
1249	if (primFail)
1250	; //ERROR
1251	else
1252	imPrimElem= mapPrimElem (primElem, alpha, v);
1253
1254	CFList source, dest;
1255	CanonicalForm bufA= mapUp (A, alpha, v, primElem, imPrimElem,
1256	source, dest);
1257	ExtensionInfo info= ExtensionInfo (v, alpha, imPrimElem, primElem);
1258	factors= biFactorize (bufA, info);
1259	}
1260	else
1261	{
1262	Variable v= chooseExtension (A, alpha);
1263	CanonicalForm primElem, imPrimElem;
1264	bool primFail= false;
1265	Variable vBuf;
1266	ASSERT (!primFail, "failure in integer factorizer");
1267	if (primFail)
1268	; //ERROR
1269	else
1270	imPrimElem= mapPrimElem (delta, beta, v);
1271
1272	CFList source, dest;
1273	CanonicalForm bufA= mapDown (A, info, source, dest);
1274	source= CFList();
1275	dest= CFList();
1276	bufA= mapUp (bufA, beta, v, delta, imPrimElem, source, dest);
1277	ExtensionInfo info= ExtensionInfo (v, beta, imPrimElem, delta);
1278	factors= biFactorize (bufA, info);
1279	}
1280	}
1281	return factors;
1282	}
1283	else // we are in GF (p^k)
1284	{
1285	int p= getCharacteristic();
1286	int extensionDeg= getGFDegree();
1287	bool extension= true;
1288	if (k == 1) // need factorization over F_p
1289	{
1290	extensionDeg++;
1291	if (ipower (p, extensionDeg) < (1<<16))
1292	// pass to GF(p^k+1)
1293	{
1294	setCharacteristic (p);
1295	Variable vBuf= rootOf (gf_mipo);
1296	A= GF2FalphaRep (A, vBuf);
1297	setCharacteristic (p, extensionDeg, 'Z');
1298	ExtensionInfo info= ExtensionInfo (extension);
1299	factors= biFactorize (A.mapinto(), info);
1300	}
1301	else // not able to pass to another GF, pass to F_p(\alpha)
1302	{
1303	setCharacteristic (p);
1304	Variable vBuf= rootOf (gf_mipo);
1305	A= GF2FalphaRep (A, vBuf);
1306	Variable v= chooseExtension (A, beta);
1307	ExtensionInfo info= ExtensionInfo (v, extension);
1308	factors= biFactorize (A, info);
1309	}
1310	}
1311	else // need factorization over GF (p^k)
1312	{
1313	if (ipower (p, 2*extensionDeg) < (1<<16))
1314	// pass to GF (p^2k)
1315	{
1316	setCharacteristic (p, 2*extensionDeg, 'Z');
1317	ExtensionInfo info= ExtensionInfo (k, cGFName, extension);
1318	factors= biFactorize (GFMapUp (A, extensionDeg), info);
1319	setCharacteristic (p, extensionDeg, cGFName);
1320	}
1321	else // not able to pass to GF (p^2k), pass to F_p (\alpha)
1322	{
1323	setCharacteristic (p);
1324	Variable v1= rootOf (gf_mipo);
1325	A= GF2FalphaRep (A, v1);
1326	Variable v2= chooseExtension (A, v1);
1327	CanonicalForm primElem, imPrimElem;
1328	bool primFail= false;
1329	Variable vBuf;
1330	primElem= primitiveElement (v1, vBuf, primFail);
1331	ASSERT (!primFail, "failure in integer factorizer");
1332	if (primFail)
1333	; //ERROR
1334	else
1335	imPrimElem= mapPrimElem (primElem, v1, v2);
1336
1337	CFList source, dest;
1338	CanonicalForm bufA= mapUp (A, v1, v2, primElem, imPrimElem,
1339	source, dest);
1340	ExtensionInfo info= ExtensionInfo (v2, v1, imPrimElem, primElem);
1341	factors= biFactorize (bufA, info);
1342	setCharacteristic (p, k, cGFName);
1343	for (CFListIterator i= factors; i.hasItem(); i++)
1344	i.getItem()= Falpha2GFRep (i.getItem());
1345	}
1346	}
1347	return factors;
1348	}
1349	}
1350
1351	#endif
1352	/* HAVE_NTL */
1353
1354

Note: See TracBrowser for help on using the repository browser.

Download in other formats: