Context Navigation

source: git/factory/facFqBivar.cc @ 04cdf06

spielwiese

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

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

Download in other formats: