Context Navigation

source: git/factory/facFqBivar.cc @ 93134c

spielwiese

Last change on this file since 93134c was 93134c, checked in by Martin Lee <martinlee84@…>, 13 years ago
bug fixes for factor recombination in case an extension is needed git-svn-id: file:///usr/local/Singular/svn/trunk@13801 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 37.7 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	for (int i= 0; i < factorNums; i++)
949	{
950	bufAeval= A;
951	bufEvaluation= evalPoint (A, bufAeval, alpha, list, GF, fail);
952	if (!derivXZero && !fail2)
953	{
954	buf= swapvar (A, x, y);
955	bufAeval2= buf;
956	bufEvaluation2= evalPoint (buf, bufAeval2, alpha, list2, GF, fail2);
957	}
958	// first try to change main variable if there is no valid evaluation point
959	if (fail && (i == 0))
960	{
961	if (!derivXZero)
962	bufEvaluation= evalPoint (buf, bufAeval, alpha, list, GF, fail);
963	else
964	fail= true;
965
966	if (!fail)
967	{
968	A= buf;
969	swap2= true;
970	}
971	}
972
973	// if there is no valid evaluation point pass to a field extension
974	if (fail && (i == 0))
975	{
976	factors= extBiFactorize (A, info);
977	appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
978	swap, swap2, N);
979	normalize (factors);
980	return factors;
981	}
982
983	// there is at least one valid evaluation point
984	// but we do not compute more univariate factorization over an extension
985	if (fail && (i != 0))
986	break;
987
988	// univariate factorization
989	TIMING_START (fac_uni_factorizer);
990	bufUniFactors= uniFactorizer (bufAeval, alpha, GF);
991	TIMING_END_AND_PRINT (fac_uni_factorizer,
992	"time for univariate factorization: ");
993	DEBOUTLN (cerr, "Lc (bufAeval)*prod (bufUniFactors)== bufAeval " <<
994	prod (bufUniFactors)*Lc (bufAeval) == bufAeval);
995
996	if (!derivXZero && !fail2)
997	{
998	TIMING_START (fac_uni_factorizer);
999	bufUniFactors2= uniFactorizer (bufAeval2, alpha, GF);
1000	TIMING_END_AND_PRINT (fac_uni_factorizer,
1001	"time for univariate factorization in y: ");
1002	DEBOUTLN (cerr, "Lc (Aeval2)*prod (uniFactors2)== Aeval2 " <<
1003	prod (bufUniFactors2)*Lc (bufAeval2) == bufAeval2);
1004	}
1005
1006	if (bufUniFactors.length() == 1 \|\|
1007	(!fail2 && !derivXZero && (bufUniFactors2.length() == 1)))
1008	{
1009	if (extension)
1010	{
1011	CFList source, dest;
1012	ExtensionInfo info2= ExtensionInfo (beta, alpha, delta, gamma, k,
1013	info.getGFName(), info.isInExtension());
1014	appendMapDown (factors, A, info2, source, dest);
1015	}
1016	else
1017	factors.append (A);
1018
1019	appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1020	swap, swap2, N);
1021
1022	normalize (factors);
1023	return factors;
1024	}
1025
1026	if (i == 0)
1027	{
1028	int subCheck= substituteCheck (bufUniFactors);
1029
1030	if (subCheck > 1)
1031	{
1032	CanonicalForm bufA= A;
1033	subst (bufA, bufA, subCheck, x);
1034	factors= biFactorize (bufA, info);
1035	reverseSubst (factors, subCheck, x);
1036	appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1037	swap, swap2, N);
1038	normalize (factors);
1039	return factors;
1040	}
1041
1042	if (!derivXZero && !fail2)
1043	{
1044	subCheck= substituteCheck (bufUniFactors2);
1045	if (subCheck > 1)
1046	{
1047	CanonicalForm bufA= A;
1048	subst (bufA, bufA, subCheck, y);
1049	factors= biFactorize (bufA, info);
1050	reverseSubst (factors, subCheck, x);
1051	appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1052	swap, swap2, N);
1053	normalize (factors);
1054	return factors;
1055	}
1056	}
1057	}
1058
1059	// degree analysis
1060	bufDegs = DegreePattern (bufUniFactors);
1061	if (!derivXZero && !fail2)
1062	bufDegs2= DegreePattern (bufUniFactors2);
1063
1064	if (i == 0)
1065	{
1066	Aeval= bufAeval;
1067	evaluation= bufEvaluation;
1068	uniFactors= bufUniFactors;
1069	degs= bufDegs;
1070	if (!derivXZero && !fail2)
1071	{
1072	Aeval2= bufAeval2;
1073	evaluation2= bufEvaluation2;
1074	uniFactors2= bufUniFactors2;
1075	degs2= bufDegs2;
1076	}
1077	}
1078	else
1079	{
1080	degs.intersect (bufDegs);
1081	if (!derivXZero && !fail2)
1082	{
1083	degs2.intersect (bufDegs2);
1084	if (bufUniFactors2.length() < uniFactors2.length())
1085	{
1086	uniFactors2= bufUniFactors2;
1087	Aeval2= bufAeval2;
1088	evaluation2= bufEvaluation2;
1089	}
1090	}
1091	if (bufUniFactors.length() < uniFactors.length())
1092	{
1093	uniFactors= bufUniFactors;
1094	Aeval= bufAeval;
1095	evaluation= bufEvaluation;
1096	}
1097	}
1098	list.append (bufEvaluation);
1099	if (!derivXZero && !fail2)
1100	list2.append (bufEvaluation2);
1101	}
1102
1103	if (!derivXZero && !fail2)
1104	{
1105	if (uniFactors.length() > uniFactors2.length() \|\|
1106	(uniFactors.length() == uniFactors2.length()
1107	&& degs.getLength() > degs2.getLength()))
1108	{
1109	degs= degs2;
1110	uniFactors= uniFactors2;
1111	evaluation= evaluation2;
1112	Aeval= Aeval2;
1113	A= buf;
1114	swap2= true;
1115	}
1116	}
1117
1118	if (degs.getLength() == 1) // A is irreducible
1119	{
1120	if (extension)
1121	{
1122	CFList source, dest;
1123	ExtensionInfo info2= ExtensionInfo (beta, alpha, delta, gamma, k,
1124	info.getGFName(), info.isInExtension());
1125	appendMapDown (factors, A, info2, source, dest);
1126	}
1127	else
1128	factors.append (A);
1129	appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1130	swap, swap2, N);
1131	normalize (factors);
1132	return factors;
1133	}
1134
1135	A= A (y + evaluation, y);
1136
1137	int liftBound;
1138	if (degree (A, y) == degree (LC (A, x)))
1139	liftBound= degree (LC (A, x)) + 1;
1140	else
1141	liftBound= degree (A, y) + 1 + degree (LC(A, x));
1142
1143	DEBOUTLN (cerr, "uniFactors= " << uniFactors);
1144	bool earlySuccess= false;
1145	CFList earlyFactors;
1146	TIMING_START (fac_hensel_lift);
1147	uniFactors= henselLiftAndEarly
1148	(A, earlySuccess, earlyFactors, degs, liftBound,
1149	uniFactors, info, evaluation);
1150	TIMING_END_AND_PRINT (fac_hensel_lift, "time for hensel lifting: ");
1151	DEBOUTLN (cerr, "lifted factors= " << uniFactors);
1152
1153	CanonicalForm MODl= power (y, liftBound);
1154	TIMING_START (fac_factor_recombination);
1155	if (!extension)
1156	factors= factorRecombination (uniFactors, A, MODl, degs);
1157	else
1158	factors= extFactorRecombination (uniFactors, A, MODl, info, degs,
1159	evaluation);
1160	TIMING_END_AND_PRINT (fac_factor_recombination,
1161	"time for factor recombination: ");
1162	if (earlySuccess)
1163	factors= Union (earlyFactors, factors);
1164	else if (!earlySuccess && degs.getLength() == 1)
1165	factors= earlyFactors;
1166
1167	if (!extension)
1168	{
1169	for (CFListIterator i= factors; i.hasItem(); i++)
1170	i.getItem()= i.getItem() (y - evaluation, y);
1171	}
1172
1173	appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1174	swap, swap2, N);
1175	normalize (factors);
1176
1177	return factors;
1178	}
1179
1180	CFList
1181	extBiFactorize (const CanonicalForm& F, const ExtensionInfo& info)
1182	{
1183
1184	CanonicalForm A= F;
1185	Variable alpha= info.getAlpha();
1186	Variable beta= info.getBeta();
1187	int k= info.getGFDegree();
1188	char cGFName= info.getGFName();
1189	CanonicalForm delta= info.getDelta();
1190
1191	bool GF= (CFFactory::gettype() == GaloisFieldDomain);
1192	Variable x= Variable (1);
1193	CFList factors;
1194	if (!GF && alpha == x) // we are in F_p
1195	{
1196	bool extension= true;
1197	int p= getCharacteristic();
1198	if (p*p < (1<<16)) // pass to GF if possible
1199	{
1200	setCharacteristic (getCharacteristic(), 2, 'Z');
1201	A= A.mapinto();
1202	ExtensionInfo info= ExtensionInfo (extension);
1203	factors= biFactorize (A, info);
1204
1205	Variable vBuf= rootOf (gf_mipo);
1206	setCharacteristic (getCharacteristic());
1207	for (CFListIterator j= factors; j.hasItem(); j++)
1208	j.getItem()= GF2FalphaRep (j.getItem(), vBuf);
1209	}
1210	else // not able to pass to GF, pass to F_p(\alpha)
1211	{
1212	CanonicalForm mipo= randomIrredpoly (3, Variable (1));
1213	Variable v= rootOf (mipo);
1214	ExtensionInfo info= ExtensionInfo (v);
1215	factors= biFactorize (A, info);
1216	}
1217	return factors;
1218	}
1219	else if (!GF && (alpha != x)) // we are in F_p(\alpha)
1220	{
1221	if (k == 1) // need factorization over F_p
1222	{
1223	int extDeg= degree (getMipo (alpha));
1224	extDeg++;
1225	CanonicalForm mipo= randomIrredpoly (extDeg + 1, Variable (1));
1226	Variable v= rootOf (mipo);
1227	ExtensionInfo info= ExtensionInfo (v);
1228	factors= biFactorize (A, info);
1229	}
1230	else
1231	{
1232	if (beta == Variable (1))
1233	{
1234	Variable v= chooseExtension (A, alpha);
1235	CanonicalForm primElem, imPrimElem;
1236	bool primFail= false;
1237	Variable vBuf;
1238	primElem= primitiveElement (alpha, vBuf, primFail);
1239	ASSERT (!primFail, "failure in integer factorizer");
1240	if (primFail)
1241	; //ERROR
1242	else
1243	imPrimElem= mapPrimElem (primElem, vBuf, v);
1244
1245	CFList source, dest;
1246	CanonicalForm bufA= mapUp (A, alpha, v, primElem, imPrimElem,
1247	source, dest);
1248	ExtensionInfo info= ExtensionInfo (v, alpha, imPrimElem, primElem);
1249	factors= biFactorize (bufA, info);
1250	}
1251	else
1252	{
1253	Variable v= chooseExtension (A, alpha);
1254	CanonicalForm primElem, imPrimElem;
1255	bool primFail= false;
1256	Variable vBuf;
1257	ASSERT (!primFail, "failure in integer factorizer");
1258	if (primFail)
1259	; //ERROR
1260	else
1261	imPrimElem= mapPrimElem (delta, beta, v); //oder mapPrimElem (primElem, vBuf, v);
1262
1263	CFList source, dest;
1264	CanonicalForm bufA= mapDown (A, info, source, dest);
1265	source= CFList();
1266	dest= CFList();
1267	bufA= mapUp (bufA, beta, v, delta, imPrimElem, source, dest);
1268	ExtensionInfo info= ExtensionInfo (v, beta, imPrimElem, delta);
1269	factors= biFactorize (bufA, info);
1270	}
1271	}
1272	return factors;
1273	}
1274	else // we are in GF (p^k)
1275	{
1276	int p= getCharacteristic();
1277	int extensionDeg= getGFDegree();
1278	bool extension= true;
1279	if (k == 1) // need factorization over F_p
1280	{
1281	extensionDeg++;
1282	if (ipower (p, extensionDeg) < (1<<16))
1283	// pass to GF(p^k+1)
1284	{
1285	setCharacteristic (p);
1286	Variable vBuf= rootOf (gf_mipo);
1287	A= GF2FalphaRep (A, vBuf);
1288	setCharacteristic (p, extensionDeg, 'Z');
1289	ExtensionInfo info= ExtensionInfo (extension);
1290	factors= biFactorize (A.mapinto(), info);
1291	}
1292	else // not able to pass to another GF, pass to F_p(\alpha)
1293	{
1294	setCharacteristic (p);
1295	Variable vBuf= rootOf (gf_mipo);
1296	A= GF2FalphaRep (A, vBuf);
1297	Variable v= chooseExtension (A, beta);
1298	ExtensionInfo info= ExtensionInfo (v, extension);
1299	factors= biFactorize (A, info);
1300	}
1301	}
1302	else // need factorization over GF (p^k)
1303	{
1304	if (ipower (p, 2*extensionDeg) < (1<<16))
1305	// pass to GF (p^2k)
1306	{
1307	setCharacteristic (p, 2*extensionDeg, 'Z');
1308	ExtensionInfo info= ExtensionInfo (k, cGFName, extension);
1309	factors= biFactorize (GFMapUp (A, extensionDeg), info);
1310	setCharacteristic (p, extensionDeg, cGFName);
1311	}
1312	else // not able to pass to GF (p^2k), pass to F_p (\alpha)
1313	{
1314	setCharacteristic (p);
1315	Variable v1= rootOf (gf_mipo);
1316	A= GF2FalphaRep (A, v1);
1317	Variable v2= chooseExtension (A, v1);
1318	CanonicalForm primElem, imPrimElem;
1319	bool primFail= false;
1320	Variable vBuf;
1321	primElem= primitiveElement (v1, vBuf, primFail);
1322	ASSERT (!primFail, "failure in integer factorizer");
1323	if (primFail)
1324	; //ERROR
1325	else
1326	imPrimElem= mapPrimElem (primElem, vBuf, v2);
1327
1328	CFList source, dest;
1329	CanonicalForm bufA= mapUp (A, v1, v2, primElem, imPrimElem,
1330	source, dest);
1331	ExtensionInfo info= ExtensionInfo (v2, v1, imPrimElem, primElem);
1332	factors= biFactorize (bufA, info);
1333	setCharacteristic (p, k, cGFName);
1334	for (CFListIterator i= factors; i.hasItem(); i++)
1335	i.getItem()= Falpha2GFRep (i.getItem());
1336	}
1337	}
1338	return factors;
1339	}
1340	}
1341
1342	#endif
1343	/* HAVE_NTL */
1344
1345

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

Download in other formats: