Context Navigation

source: git/factory/facFqBivar.cc @ 55608a7

fieker-DuValspielwiese

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

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

Download in other formats: