Context Navigation

source: git/factory/facFqBivar.cc @ e558c41

spielwiese

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

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

Download in other formats: