Context Navigation

source: git/factory/facFqBivar.cc @ 1372ae

spielwiese

Last change on this file since 1372ae was 1372ae, checked in by Hans Schoenemann <hannes@…>, 13 years ago
Variable(i) ->v.level() git-svn-id: file:///usr/local/Singular/svn/trunk@13666 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 36.6 KB

Line
1	/*****************************************************************************\
2	* Computer Algebra System SINGULAR
3	\*****************************************************************************/
4	/** @file facFqBivar.cc
5	*
6	* This file provides functions for factorizing a bivariate polynomial over
7	* \f$ F_{p} \f$ , \f$ F_{p}(\alpha ) \f$ or GF, based on "Modern Computer
8	* Algebra, Chapter 15" by J. von zur Gathen & J. Gerhard and "Factoring
9	* multivariate polynomials over a finite field" by L. Bernardin.
10	*
11	* ABSTRACT: In contrast to biFactorizer() in facFqFactorice.cc we evaluate and
12	* factorize the polynomial in both variables. So far factor recombination is
13	* done naive!
14	*
15	* @author Martin Lee
16	*
17	* @internal @version \$Id$
18	*
19	**/
20	/*****************************************************************************/
21
22	#include <config.h>
23
24	#include "assert.h"
25	#include "debug.h"
26	#include "timing.h"
27
28	#include "canonicalform.h"
29	#include "cf_defs.h"
30	#include "cf_map_ext.h"
31	#include "cf_random.h"
32	#include "facHensel.h"
33	#include "cf_map.h"
34	#include "cf_gcd_smallp.h"
35	#include "facFqBivar.h"
36
37
38	#ifdef HAVE_NTL
39	#include "NTLconvert.h"
40
41	TIMING_DEFINE_PRINT(fac_uni_factorizer);
42	TIMING_DEFINE_PRINT(fac_hensel_lift);
43	TIMING_DEFINE_PRINT(fac_factor_recombination);
44
45	CanonicalForm prodMod0 (const CFList& L, const CanonicalForm& M)
46	{
47	if (L.isEmpty())
48	return 1;
49	else if (L.length() == 1)
50	return mod (L.getFirst()(0, 1) , M);
51	else if (L.length() == 2)
52	return mod (L.getFirst()(0, 1)*L.getLast()(0, 1), M);
53	else
54	{
55	int l= L.length()/2;
56	CFListIterator i= L;
57	CFList tmp1, tmp2;
58	CanonicalForm buf1, buf2;
59	for (int j= 1; j <= l; j++, i++)
60	tmp1.append (i.getItem());
61	tmp2= Difference (L, tmp1);
62	buf1= prodMod0 (tmp1, M);
63	buf2= prodMod0 (tmp2, M);
64	return mod (buf1*buf2, M);
65	}
66	}
67
68	CanonicalForm evalPoint (const CanonicalForm& F, CanonicalForm & eval,
69	const Variable& alpha, CFList& list, const bool& GF,
70	bool& fail)
71	{
72	fail= false;
73	Variable x= Variable(2);
74	FFRandom genFF;
75	GFRandom genGF;
76	CanonicalForm random, mipo;
77	double bound;
78	int p= getCharacteristic ();
79	if (alpha.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	if (GF)
863	gcdDerivX= GCD_GF (A, derivX);
864	else if (alpha == x)
865	gcdDerivX= GCD_small_p (A, derivX);
866	else
867	gcdDerivX= GCD_Fp_extension (A, derivX, alpha);
868	if (degree (gcdDerivX) > 0)
869	{
870	CanonicalForm g= A/gcdDerivX;
871	CFList factorsG=
872	Union (biFactorize (g, info), biFactorize (gcdDerivX, info));
873	append (factorsG, contentAxFactors);
874	append (factorsG, contentAyFactors);
875	decompress (factorsG, N);
876	normalize (factors);
877	return factorsG;
878	}
879	}
880	bool derivYZero= false;
881	CanonicalForm derivY= deriv (A, y);
882	CanonicalForm gcdDerivY;
883	if (derivY.isZero())
884	derivYZero= true;
885	else
886	{
887	if (GF)
888	gcdDerivY= GCD_GF (A, derivY);
889	else if (alpha == x)
890	gcdDerivY= GCD_small_p (A, derivY);
891	else
892	gcdDerivY= GCD_Fp_extension (A, derivY, alpha);
893	if (degree (gcdDerivY) > 0)
894	{
895	CanonicalForm g= A/gcdDerivY;
896	CFList factorsG=
897	Union (biFactorize (g, info), biFactorize (gcdDerivY, info));
898	append (factorsG, contentAxFactors);
899	append (factorsG, contentAyFactors);
900	decompress (factorsG, N);
901	normalize (factors);
902	return factorsG;
903	}
904	}
905	//main variable is chosen s.t. the degree in x is minimal
906	bool swap= false;
907	if ((degree (A) > degree (A, x)) \|\| derivXZero)
908	{
909	if (!derivYZero)
910	{
911	A= swapvar (A, y, x);
912	swap= derivXZero;
913	derivXZero= derivYZero;
914	derivYZero= swap;
915	swap= true;
916	}
917	}
918
919	bool fail;
920	CanonicalForm Aeval, evaluation, bufAeval, bufEvaluation, buf;
921	CFList uniFactors, list, bufUniFactors;
922	DegreePattern degs;
923	DegreePattern bufDegs;
924
925	bool fail2= false;
926	CanonicalForm Aeval2, evaluation2, bufAeval2, bufEvaluation2;
927	CFList bufUniFactors2, list2, uniFactors2;
928	DegreePattern degs2;
929	DegreePattern bufDegs2;
930	bool swap2= false;
931
932	// several univariate factorizations to obtain more information about the
933	// degree pattern therefore usually less combinations have to be tried during
934	// the recombination process
935	int factorNums= 5;
936	double logarithm= (double) ilog2 (totaldegree (A));
937	logarithm /= log2exp;
938	logarithm= ceil (logarithm);
939	if (factorNums < (int) logarithm)
940	factorNums= (int) logarithm;
941	for (int i= 0; i < factorNums; i++)
942	{
943	bufAeval= A;
944	bufEvaluation= evalPoint (A, bufAeval, alpha, list, GF, fail);
945	if (!derivXZero && !fail2)
946	{
947	buf= swapvar (A, x, y);
948	bufAeval2= buf;
949	bufEvaluation2= evalPoint (buf, bufAeval2, alpha, list2, GF, fail2);
950	}
951	// first try to change main variable if there is no valid evaluation point
952	if (fail && (i == 0))
953	{
954	if (!derivXZero)
955	bufEvaluation= evalPoint (buf, bufAeval, alpha, list, GF, fail);
956	else
957	fail= true;
958
959	if (!fail)
960	{
961	A= buf;
962	swap2= true;
963	}
964	}
965
966	// if there is no valid evaluation point pass to a field extension
967	if (fail && (i == 0))
968	{
969	factors= extBiFactorize (A, info);
970	appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
971	swap, swap2, N);
972	normalize (factors);
973	return factors;
974	}
975
976	// there is at least one valid evaluation point
977	// but we do not compute more univariate factorization over an extension
978	if (fail && (i != 0))
979	break;
980
981	// univariate factorization
982	TIMING_START (fac_uni_factorizer);
983	bufUniFactors= uniFactorizer (bufAeval, alpha, GF);
984	TIMING_END_AND_PRINT (fac_uni_factorizer,
985	"time for univariate factorization: ");
986	DEBOUTLN (cerr, "Lc (bufAeval)*prod (bufUniFactors)== bufAeval " <<
987	prod (bufUniFactors)*Lc (bufAeval) == bufAeval);
988
989	if (!derivXZero && !fail2)
990	{
991	TIMING_START (fac_uni_factorizer);
992	bufUniFactors2= uniFactorizer (bufAeval2, alpha, GF);
993	TIMING_END_AND_PRINT (fac_uni_factorizer,
994	"time for univariate factorization in y: ");
995	DEBOUTLN (cerr, "Lc (Aeval2)*prod (uniFactors2)== Aeval2 " <<
996	prod (bufUniFactors2)*Lc (bufAeval2) == bufAeval2);
997	}
998
999	if (bufUniFactors.length() == 1 \|\|
1000	(!fail2 && !derivXZero && (bufUniFactors2.length() == 1)))
1001	{
1002	if (extension)
1003	{
1004	CFList source, dest;
1005	ExtensionInfo info2= ExtensionInfo (beta, alpha, delta, gamma, k,
1006	info.getGFName(), info.isInExtension());
1007	appendMapDown (factors, A, info2, source, dest);
1008	}
1009	else
1010	factors.append (A);
1011
1012	appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1013	swap, swap2, N);
1014
1015	normalize (factors);
1016	return factors;
1017	}
1018
1019	// degree analysis
1020	bufDegs = DegreePattern (bufUniFactors);
1021	if (!derivXZero && !fail2)
1022	bufDegs2= DegreePattern (bufUniFactors2);
1023
1024	if (i == 0)
1025	{
1026	Aeval= bufAeval;
1027	evaluation= bufEvaluation;
1028	uniFactors= bufUniFactors;
1029	degs= bufDegs;
1030	if (!derivXZero && !fail2)
1031	{
1032	Aeval2= bufAeval2;
1033	evaluation2= bufEvaluation2;
1034	uniFactors2= bufUniFactors2;
1035	degs2= bufDegs2;
1036	}
1037	}
1038	else
1039	{
1040	degs.intersect (bufDegs);
1041	if (!derivXZero && !fail2)
1042	{
1043	degs2.intersect (bufDegs2);
1044	if (bufUniFactors2.length() < uniFactors2.length())
1045	{
1046	uniFactors2= bufUniFactors2;
1047	Aeval2= bufAeval2;
1048	evaluation2= bufEvaluation2;
1049	}
1050	}
1051	if (bufUniFactors.length() < uniFactors.length())
1052	{
1053	uniFactors= bufUniFactors;
1054	Aeval= bufAeval;
1055	evaluation= bufEvaluation;
1056	}
1057	}
1058	list.append (bufEvaluation);
1059	if (!derivXZero && !fail2)
1060	list2.append (bufEvaluation2);
1061	}
1062
1063	if (!derivXZero && !fail2)
1064	{
1065	if (uniFactors.length() > uniFactors2.length() \|\|
1066	(uniFactors.length() == uniFactors2.length()
1067	&& degs.getLength() > degs2.getLength()))
1068	{
1069	degs= degs2;
1070	uniFactors= uniFactors2;
1071	evaluation= evaluation2;
1072	Aeval= Aeval2;
1073	A= buf;
1074	swap2= true;
1075	}
1076	}
1077
1078	if (degs.getLength() == 1) // A is irreducible
1079	{
1080	if (extension)
1081	{
1082	CFList source, dest;
1083	ExtensionInfo info2= ExtensionInfo (beta, alpha, delta, gamma, k,
1084	info.getGFName(), info.isInExtension());
1085	appendMapDown (factors, A, info2, source, dest);
1086	}
1087	else
1088	factors.append (A);
1089	appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1090	swap, swap2, N);
1091	normalize (factors);
1092	return factors;
1093	}
1094
1095	A= A (y + evaluation, y);
1096
1097	int liftBound;
1098	if (degree (A, y) == degree (LC (A, x)))
1099	liftBound= degree (LC (A, x)) + 1;
1100	else
1101	liftBound= degree (A, y) + 1 + degree (LC(A, x));
1102
1103	DEBOUTLN (cerr, "uniFactors= " << uniFactors);
1104	bool earlySuccess= false;
1105	CFList earlyFactors;
1106	TIMING_START (fac_hensel_lift);
1107	uniFactors= henselLiftAndEarly
1108	(A, earlySuccess, earlyFactors, degs, liftBound,
1109	uniFactors, info, evaluation);
1110	TIMING_END_AND_PRINT (fac_hensel_lift, "time for hensel lifting: ");
1111	DEBOUTLN (cerr, "lifted factors= " << uniFactors);
1112
1113	CanonicalForm MODl= power (y, liftBound);
1114	TIMING_START (fac_factor_recombination);
1115	if (!extension)
1116	factors= factorRecombination (uniFactors, A, MODl, degs);
1117	else
1118	factors= extFactorRecombination (uniFactors, A, MODl, info, degs,
1119	evaluation);
1120	TIMING_END_AND_PRINT (fac_factor_recombination,
1121	"time for factor recombination: ");
1122	if (earlySuccess)
1123	factors= Union (earlyFactors, factors);
1124	else if (!earlySuccess && degs.getLength() == 1)
1125	factors= earlyFactors;
1126
1127	if (!extension)
1128	{
1129	for (CFListIterator i= factors; i.hasItem(); i++)
1130	i.getItem()= i.getItem() (y - evaluation, y);
1131	}
1132
1133	appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1134	swap, swap2, N);
1135	normalize (factors);
1136
1137	return factors;
1138	}
1139
1140	CFList
1141	extBiFactorize (const CanonicalForm& F, const ExtensionInfo& info)
1142	{
1143
1144	CanonicalForm A= F;
1145	Variable alpha= info.getAlpha();
1146	Variable beta= info.getBeta();
1147	int k= info.getGFDegree();
1148	char cGFName= info.getGFName();
1149	CanonicalForm delta= info.getDelta();
1150
1151	bool GF= (CFFactory::gettype() == GaloisFieldDomain);
1152	Variable x= Variable (1);
1153	CFList factors;
1154	if (!GF && alpha == x) // we are in F_p
1155	{
1156	bool extension= true;
1157	int p= getCharacteristic();
1158	if (p*p < (1<<16)) // pass to GF if possible
1159	{
1160	setCharacteristic (getCharacteristic(), 2, 'Z');
1161	A= A.mapinto();
1162	ExtensionInfo info= ExtensionInfo (extension);
1163	factors= biFactorize (A, info);
1164
1165	Variable vBuf= rootOf (gf_mipo);
1166	setCharacteristic (getCharacteristic());
1167	for (CFListIterator j= factors; j.hasItem(); j++)
1168	j.getItem()= GF2FalphaRep (j.getItem(), vBuf);
1169	}
1170	else // not able to pass to GF, pass to F_p(\alpha)
1171	{
1172	CanonicalForm mipo= randomIrredpoly (3, Variable (1));
1173	Variable v= rootOf (mipo);
1174	ExtensionInfo info= ExtensionInfo (v);
1175	factors= biFactorize (A, info);
1176	}
1177	return factors;
1178	}
1179	else if (!GF && (alpha != x)) // we are in F_p(\alpha)
1180	{
1181	if (k == 1) // need factorization over F_p
1182	{
1183	int extDeg= degree (getMipo (alpha));
1184	extDeg++;
1185	CanonicalForm mipo= randomIrredpoly (extDeg + 1, Variable (1));
1186	Variable v= rootOf (mipo);
1187	ExtensionInfo info= ExtensionInfo (v);
1188	factors= biFactorize (A, info);
1189	}
1190	else
1191	{
1192	if (beta == Variable (1))
1193	{
1194	Variable v= chooseExtension (A, alpha);
1195	CanonicalForm primElem, imPrimElem;
1196	bool primFail= false;
1197	Variable vBuf;
1198	primElem= primitiveElement (alpha, vBuf, primFail);
1199	ASSERT (!primFail, "failure in integer factorizer");
1200	if (primFail)
1201	; //ERROR
1202	else
1203	imPrimElem= mapPrimElem (primElem, vBuf, v);
1204
1205	CFList source, dest;
1206	CanonicalForm bufA= mapUp (A, alpha, v, primElem, imPrimElem,
1207	source, dest);
1208	ExtensionInfo info= ExtensionInfo (v, alpha, imPrimElem, primElem);
1209	factors= biFactorize (bufA, info);
1210	}
1211	else
1212	{
1213	Variable v= chooseExtension (A, alpha);
1214	CanonicalForm primElem, imPrimElem;
1215	bool primFail= false;
1216	Variable vBuf;
1217	ASSERT (!primFail, "failure in integer factorizer");
1218	if (primFail)
1219	; //ERROR
1220	else
1221	imPrimElem= mapPrimElem (delta, beta, v); //oder mapPrimElem (primElem, vBuf, v);
1222
1223	CFList source, dest;
1224	CanonicalForm bufA= mapDown (A, info, source, dest);
1225	source= CFList();
1226	dest= CFList();
1227	bufA= mapUp (bufA, beta, v, delta, imPrimElem, source, dest);
1228	ExtensionInfo info= ExtensionInfo (v, beta, imPrimElem, delta);
1229	factors= biFactorize (bufA, info);
1230	}
1231	}
1232	return factors;
1233	}
1234	else // we are in GF (p^k)
1235	{
1236	int p= getCharacteristic();
1237	int extensionDeg= getGFDegree();
1238	bool extension= true;
1239	if (k == 1) // need factorization over F_p
1240	{
1241	extensionDeg++;
1242	if (ipower (p, extensionDeg) < (1<<16))
1243	// pass to GF(p^k+1)
1244	{
1245	setCharacteristic (p);
1246	Variable vBuf= rootOf (gf_mipo);
1247	A= GF2FalphaRep (A, vBuf);
1248	setCharacteristic (p, extensionDeg, 'Z');
1249	ExtensionInfo info= ExtensionInfo (extension);
1250	factors= biFactorize (A.mapinto(), info);
1251	}
1252	else // not able to pass to another GF, pass to F_p(\alpha)
1253	{
1254	setCharacteristic (p);
1255	Variable vBuf= rootOf (gf_mipo);
1256	A= GF2FalphaRep (A, vBuf);
1257	Variable v= chooseExtension (A, beta);
1258	ExtensionInfo info= ExtensionInfo (v, extension);
1259	factors= biFactorize (A, info);
1260	}
1261	}
1262	else // need factorization over GF (p^k)
1263	{
1264	if (ipower (p, 2*extensionDeg) < (1<<16))
1265	// pass to GF (p^2k)
1266	{
1267	setCharacteristic (p, 2*extensionDeg, 'Z');
1268	ExtensionInfo info= ExtensionInfo (k, cGFName, extension);
1269	factors= biFactorize (GFMapUp (A, extensionDeg), info);
1270	setCharacteristic (p, extensionDeg, cGFName);
1271	}
1272	else // not able to pass to GF (p^2k), pass to F_p (\alpha)
1273	{
1274	setCharacteristic (p);
1275	Variable v1= rootOf (gf_mipo);
1276	A= GF2FalphaRep (A, v1);
1277	Variable v2= chooseExtension (A, v1);
1278	CanonicalForm primElem, imPrimElem;
1279	bool primFail= false;
1280	Variable vBuf;
1281	primElem= primitiveElement (v1, vBuf, primFail);
1282	ASSERT (!primFail, "failure in integer factorizer");
1283	if (primFail)
1284	; //ERROR
1285	else
1286	imPrimElem= mapPrimElem (primElem, vBuf, v2);
1287
1288	CFList source, dest;
1289	CanonicalForm bufA= mapUp (A, v1, v2, primElem, imPrimElem,
1290	source, dest);
1291	ExtensionInfo info= ExtensionInfo (v2, v1, imPrimElem, primElem);
1292	factors= biFactorize (bufA, info);
1293	setCharacteristic (p, k, cGFName);
1294	for (CFListIterator i= factors; i.hasItem(); i++)
1295	i.getItem()= Falpha2GFRep (i.getItem());
1296	}
1297	}
1298	return factors;
1299	}
1300	}
1301
1302	#endif
1303	/* HAVE_NTL */
1304
1305

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

Download in other formats: