1 // emacs edit mode for this file is -*- C++ -*-
2 /****************************************
3 *  Computer Algebra System SINGULAR     *
4 ****************************************/
5 /*
6 * ABSTRACT: interface between Singular and factory
7 */
8 
9 //#define FACTORIZE2_DEBUG
10 
11 #include "misc/auxiliary.h"
12 #include "clapsing.h"
13 
14 #include "factory/factory.h"
15 #include "factory/cf_roots.h"
16 
17 #include "coeffs/numbers.h"
18 #include "coeffs/coeffs.h"
19 #include "coeffs/bigintmat.h"
20 
21 #include "monomials/ring.h"
22 #include "simpleideals.h"
23 #include "polys/flintconv.h"
24 #include "polys/flint_mpoly.h"
25 
26 
27 //#include "polys.h"
28 #define TRANSEXT_PRIVATES
29 
30 #include "ext_fields/transext.h"
31 
32 
33 #include "clapconv.h"
34 
35 #include "polys/monomials/p_polys.h"
36 #include "polys/monomials/ring.h"
37 #include "polys/simpleideals.h"
38 #include "misc/intvec.h"
39 #include "polys/matpol.h"
40 #include "coeffs/bigintmat.h"
41 
42 
43 void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
44 
singclap_gcd_r(poly f,poly g,const ring r)45 poly singclap_gcd_r ( poly f, poly g, const ring r )
46 {
47   poly res=NULL;
48 
49   assume(f!=NULL);
50   assume(g!=NULL);
51 
52   if(pNext(f)==NULL)
53   {
54     return p_GcdMon(f,g,r);
55   }
56   else if(pNext(g)==NULL)
57   {
58     return p_GcdMon(g,f,r);
59   }
60   #ifdef HAVE_FLINT
61   #if __FLINT_RELEASE >= 20503
62   if (rField_is_Zp(r) && (r->cf->ch>10))
63   {
64     nmod_mpoly_ctx_t ctx;
65     if (!convSingRFlintR(ctx,r))
66     {
67       // leading coef. 1
68       return Flint_GCD_MP(f,pLength(f),g,pLength(g),ctx,r);
69     }
70   }
71   else
72   if (rField_is_Q(r))
73   {
74     fmpq_mpoly_ctx_t ctx;
75     if (!convSingRFlintR(ctx,r))
76     {
77       // leading coef. positive, all coeffs in Z
78       poly res=Flint_GCD_MP(f,pLength(f),g,pLength(g),ctx,r);
79       res=p_Cleardenom(res,r);
80       return res;
81     }
82   }
83   #endif
84   #endif
85   Off(SW_RATIONAL);
86   if (rField_is_Q(r) || rField_is_Zp(r) || rField_is_Z(r)
87   || (rField_is_Zn(r)&&(r->cf->convSingNFactoryN!=ndConvSingNFactoryN)))
88   {
89     setCharacteristic( rChar(r) );
90     CanonicalForm F( convSingPFactoryP( f,r ) ), G( convSingPFactoryP( g, r ) );
91     res=convFactoryPSingP( gcd( F, G ) , r);
92     if ( rField_is_Zp(r))
93       p_Norm(res,r); // leading coef. 1
94     else if (rField_is_Q(r) && (!n_GreaterZero(pGetCoeff(res),r->cf)))
95       res = p_Neg(res,r); // leading coef. positive, all coeffs in Z
96   }
97   // and over Q(a) / Fp(a)
98   else if ( r->cf->extRing!=NULL )
99   {
100     if ( rField_is_Q_a(r)) setCharacteristic( 0 );
101     else                   setCharacteristic( rChar(r) );
102     if (r->cf->extRing->qideal!=NULL)
103     {
104       bool b1=isOn(SW_USE_QGCD);
105       if ( rField_is_Q_a(r) ) On(SW_USE_QGCD);
106       CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->qideal->m[0],
107                                            r->cf->extRing);
108       Variable a=rootOf(mipo);
109       CanonicalForm F( convSingAPFactoryAP( f,a,r ) ),
110                     G( convSingAPFactoryAP( g,a,r ) );
111       res= convFactoryAPSingAP( gcd( F, G ),r );
112       prune (a);
113       if (!b1) Off(SW_USE_QGCD);
114       if ( rField_is_Zp_a(r)) p_Norm(res,r); // leading coef. 1
115     }
116     else
117     {
118       convSingTrP(f,r);
119       convSingTrP(g,r);
120       CanonicalForm F( convSingTrPFactoryP( f,r ) ), G( convSingTrPFactoryP( g,r ) );
121       res= convFactoryPSingTrP( gcd( F, G ),r );
122     }
123   }
124   else if (r->cf->convSingNFactoryN==ndConvSingNFactoryN)
125     WerrorS( feNotImplemented );
126   else
127   { // handle user type coeffs:
128     setCharacteristic( rChar(r) );
129     CanonicalForm F( convSingPFactoryP( f,r ) ), G( convSingPFactoryP( g, r ) );
130     res=convFactoryPSingP( gcd( F, G ) , r);
131   }
132   Off(SW_RATIONAL);
133   return res;
134 }
135 
singclap_gcd_and_divide(poly & f,poly & g,const ring r)136 poly singclap_gcd_and_divide ( poly& f, poly& g, const ring r)
137 {
138   poly res=NULL;
139 
140   if (g == NULL)
141   {
142     res= f;
143     f=p_One (r);
144     return res;
145   }
146   if (f==NULL)
147   {
148     res= g;
149     g=p_One (r);
150     return res;
151   }
152   if (pNext(g)==NULL)
153   {
154     poly G=p_GcdMon(g,f,r);
155     if (!n_IsOne(pGetCoeff(G),r->cf)
156     || (!p_IsConstant(G,r)))
157     {
158       f=p_Div_mm(f,G,r);
159       g=p_Div_mm(g,G,r);
160     }
161     return G;
162   }
163   else if (pNext(f)==NULL)
164   {
165     poly G=p_GcdMon(f,g,r);
166     if (!n_IsOne(pGetCoeff(G),r->cf)
167     || (!p_IsConstant(G,r)))
168     {
169       f=p_Div_mm(f,G,r);
170       g=p_Div_mm(g,G,r);
171     }
172     return G;
173   }
174 
175   Off(SW_RATIONAL);
176   CanonicalForm F,G,GCD;
177   if (rField_is_Q(r) || (rField_is_Zp(r))
178   || (rField_is_Zn(r)&&(r->cf->convSingNFactoryN!=ndConvSingNFactoryN)))
179   {
180     bool b1=isOn(SW_USE_EZGCD_P);
181     setCharacteristic( rChar(r) );
182     F=convSingPFactoryP( f,r );
183     G=convSingPFactoryP( g,r );
184     GCD=gcd(F,G);
185     if (!GCD.isOne())
186     {
187       p_Delete(&f,r);
188       p_Delete(&g,r);
189       if (getCharacteristic() == 0)
190         On (SW_RATIONAL);
191       F /= GCD;
192       G /= GCD;
193       if (getCharacteristic() == 0)
194       {
195         CanonicalForm denF= bCommonDen (F);
196         CanonicalForm denG= bCommonDen (G);
197         G *= denG;
198         F *= denF;
199         Off (SW_RATIONAL);
200         CanonicalForm gcddenFdenG= gcd (denG, denF);
201         denG /= gcddenFdenG;
202         denF /= gcddenFdenG;
203         On (SW_RATIONAL);
204         G *= denF;
205         F *= denG;
206       }
207       f=convFactoryPSingP( F, r);
208       g=convFactoryPSingP( G, r);
209     }
210     res=convFactoryPSingP( GCD , r);
211     if (!b1) Off (SW_USE_EZGCD_P);
212   }
213   // and over Q(a) / Fp(a)
214   else if ( r->cf->extRing )
215   {
216     if ( rField_is_Q_a(r)) setCharacteristic( 0 );
217     else                   setCharacteristic( rChar(r) );
218     if (r->cf->extRing->qideal!=NULL)
219     {
220       bool b1=isOn(SW_USE_QGCD);
221       if ( rField_is_Q_a(r) ) On(SW_USE_QGCD);
222       CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->qideal->m[0],
223                                            r->cf->extRing);
224       Variable a=rootOf(mipo);
225       F=( convSingAPFactoryAP( f,a,r ) );
226       G=( convSingAPFactoryAP( g,a,r ) );
227       GCD=gcd(F,G);
228       if (!GCD.isOne())
229       {
230         p_Delete(&f,r);
231         p_Delete(&g,r);
232         if (getCharacteristic() == 0)
233           On (SW_RATIONAL);
234         F /= GCD;
235         G /= GCD;
236         if (getCharacteristic() == 0)
237         {
238           CanonicalForm denF= bCommonDen (F);
239           CanonicalForm denG= bCommonDen (G);
240           G *= denG;
241           F *= denF;
242           Off (SW_RATIONAL);
243           CanonicalForm gcddenFdenG= gcd (denG, denF);
244           denG /= gcddenFdenG;
245           denF /= gcddenFdenG;
246           On (SW_RATIONAL);
247           G *= denF;
248           F *= denG;
249         }
250         f= convFactoryAPSingAP( F,r );
251         g= convFactoryAPSingAP( G,r );
252       }
253       res= convFactoryAPSingAP( GCD,r );
254       prune (a);
255       if (!b1) Off(SW_USE_QGCD);
256     }
257     else
258     {
259       F=( convSingTrPFactoryP( f,r ) );
260       G=( convSingTrPFactoryP( g,r ) );
261       GCD=gcd(F,G);
262       if (!GCD.isOne())
263       {
264         p_Delete(&f,r);
265         p_Delete(&g,r);
266         if (getCharacteristic() == 0)
267           On (SW_RATIONAL);
268         F /= GCD;
269         G /= GCD;
270         if (getCharacteristic() == 0)
271         {
272           CanonicalForm denF= bCommonDen (F);
273           CanonicalForm denG= bCommonDen (G);
274           G *= denG;
275           F *= denF;
276           Off (SW_RATIONAL);
277           CanonicalForm gcddenFdenG= gcd (denG, denF);
278           denG /= gcddenFdenG;
279           denF /= gcddenFdenG;
280           On (SW_RATIONAL);
281           G *= denF;
282           F *= denG;
283         }
284         f= convFactoryPSingTrP( F,r );
285         g= convFactoryPSingTrP( G,r );
286       }
287       res= convFactoryPSingTrP( GCD,r );
288     }
289   }
290   else
291     WerrorS( feNotImplemented );
292   Off(SW_RATIONAL);
293   return res;
294 }
295 
296 /*2 find the maximal exponent of var(i) in poly p*/
pGetExp_Var(poly p,int i,const ring r)297 int pGetExp_Var(poly p, int i, const ring r)
298 {
299   int m=0;
300   int mm;
301   while (p!=NULL)
302   {
303     mm=p_GetExp(p,i,r);
304     if (mm>m) m=mm;
305     pIter(p);
306   }
307   return m;
308 }
309 
310 // destroys f,g,x
singclap_resultant(poly f,poly g,poly x,const ring r)311 poly singclap_resultant ( poly f, poly g , poly x, const ring r)
312 {
313   poly res=NULL;
314   int i=p_IsPurePower(x, r);
315   if (i==0)
316   {
317     WerrorS("3rd argument must be a ring variable");
318     goto resultant_returns_res;
319   }
320   if ((f==NULL) || (g==NULL))
321     goto resultant_returns_res;
322   // for now there is only the possibility to handle polynomials over
323   // Q and Fp ...
324   if (rField_is_Zp(r) || rField_is_Q(r)
325   || (rField_is_Zn(r)&&(r->cf->convSingNFactoryN!=ndConvSingNFactoryN)))
326   {
327     Variable X(i);
328     setCharacteristic( rChar(r) );
329     CanonicalForm F( convSingPFactoryP( f,r ) ), G( convSingPFactoryP( g,r ) );
330     res=convFactoryPSingP( resultant( F, G, X),r );
331     Off(SW_RATIONAL);
332     goto resultant_returns_res;
333   }
334   // and over Q(a) / Fp(a)
335   else if (r->cf->extRing!=NULL)
336   {
337     if (rField_is_Q_a(r)) setCharacteristic( 0 );
338     else               setCharacteristic( rChar(r) );
339     Variable X(i+rPar(r));
340     if (r->cf->extRing->qideal!=NULL)
341     {
342       //Variable X(i);
343       CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->qideal->m[0],
344                                            r->cf->extRing);
345       Variable a=rootOf(mipo);
346       CanonicalForm F( convSingAPFactoryAP( f,a,r ) ),
347                     G( convSingAPFactoryAP( g,a,r ) );
348       res= convFactoryAPSingAP( resultant( F, G, X ),r );
349       prune (a);
350     }
351     else
352     {
353       //Variable X(i+rPar(currRing));
354       number nf,ng;
355       p_Cleardenom_n(f,r,nf);p_Cleardenom_n(g,r,ng);
356       int ef,eg;
357       ef=pGetExp_Var(f,i,r);
358       eg=pGetExp_Var(g,i,r);
359       CanonicalForm F( convSingTrPFactoryP( f,r ) ), G( convSingTrPFactoryP( g,r ) );
360       res= convFactoryPSingTrP( resultant( F, G, X ),r );
361       if ((nf!=NULL)&&(!n_IsOne(nf,r->cf)))
362       {
363         number n=n_Invers(nf,r->cf);
364         while(eg>0)
365         {
366           res=__p_Mult_nn(res,n,r);
367           eg--;
368         }
369         n_Delete(&n,r->cf);
370       }
371       n_Delete(&nf,r->cf);
372       if ((ng!=NULL)&&(!n_IsOne(ng,r->cf)))
373       {
374         number n=n_Invers(ng,r->cf);
375         while(ef>0)
376         {
377           res=__p_Mult_nn(res,n,r);
378           ef--;
379         }
380         n_Delete(&n,r->cf);
381       }
382       n_Delete(&ng,r->cf);
383     }
384     Off(SW_RATIONAL);
385     goto resultant_returns_res;
386   }
387   else
388     WerrorS( feNotImplemented );
389 resultant_returns_res:
390   p_Delete(&f,r);
391   p_Delete(&g,r);
392   p_Delete(&x,r);
393   return res;
394 }
395 //poly singclap_resultant ( poly f, poly g , poly x)
396 //{
397 //  int i=pVar(x);
398 //  if (i==0)
399 //  {
400 //    WerrorS("ringvar expected");
401 //    return NULL;
402 //  }
403 //  ideal I=idInit(1,1);
404 //
405 //  // get the coeffs von f wrt. x:
406 //  I->m[0]=pCopy(f);
407 //  matrix ffi=mpCoeffs(I,i);
408 //  ffi->rank=1;
409 //  ffi->ncols=ffi->nrows;
410 //  ffi->nrows=1;
411 //  ideal fi=(ideal)ffi;
412 //
413 //  // get the coeffs von g wrt. x:
414 //  I->m[0]=pCopy(g);
415 //  matrix ggi=mpCoeffs(I,i);
416 //  ggi->rank=1;
417 //  ggi->ncols=ggi->nrows;
418 //  ggi->nrows=1;
419 //  ideal gi=(ideal)ggi;
420 //
421 //  // contruct the matrix:
422 //  int fn=IDELEMS(fi); //= deg(f,x)+1
423 //  int gn=IDELEMS(gi); //= deg(g,x)+1
424 //  matrix m=mpNew(fn+gn-2,fn+gn-2);
425 //  if(m==NULL)
426 //  {
427 //    return NULL;
428 //  }
429 //
430 //  // enter the coeffs into m:
431 //  int j;
432 //  for(i=0;i<gn-1;i++)
433 //  {
434 //    for(j=0;j<fn;j++)
435 //    {
436 //      MATELEM(m,i+1,fn-j+i)=pCopy(fi->m[j]);
437 //    }
438 //  }
439 //  for(i=0;i<fn-1;i++)
440 //  {
441 //    for(j=0;j<gn;j++)
442 //    {
443 //      MATELEM(m,gn+i,gn-j+i)=pCopy(gi->m[j]);
444 //    }
445 //  }
446 //
447 //  poly r=mpDet(m);
448 //
449 //  idDelete(&fi);
450 //  idDelete(&gi);
451 //  idDelete((ideal *)&m);
452 //  return r;
453 //}
454 
singclap_extgcd(poly f,poly g,poly & res,poly & pa,poly & pb,const ring r)455 BOOLEAN singclap_extgcd ( poly f, poly g, poly &res, poly &pa, poly &pb , const ring r)
456 {
457   // for now there is only the possibility to handle univariate
458   // polynomials over
459   // Q and Fp ...
460   res=NULL;pa=NULL;pb=NULL;
461   On(SW_SYMMETRIC_FF);
462   if ( rField_is_Q(r) || rField_is_Zp(r)
463   || (rField_is_Zn(r)&&(r->cf->convSingNFactoryN!=ndConvSingNFactoryN)))
464   {
465     setCharacteristic( rChar(r) );
466     CanonicalForm F( convSingPFactoryP( f,r ) ), G( convSingPFactoryP( g,r) );
467     CanonicalForm FpG=F+G;
468     if (!(FpG.isUnivariate()|| FpG.inCoeffDomain()))
469     //if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
470     {
471       Off(SW_RATIONAL);
472       WerrorS("not univariate");
473       return TRUE;
474     }
475     CanonicalForm Fa,Gb;
476     On(SW_RATIONAL);
477     res=convFactoryPSingP( extgcd( F, G, Fa, Gb ),r );
478     pa=convFactoryPSingP(Fa,r);
479     pb=convFactoryPSingP(Gb,r);
480     Off(SW_RATIONAL);
481   }
482   // and over Q(a) / Fp(a)
483   else if ( r->cf->extRing!=NULL )
484   {
485     if (rField_is_Q_a(r)) setCharacteristic( 0 );
486     else                 setCharacteristic( rChar(r) );
487     CanonicalForm Fa,Gb;
488     if (r->cf->extRing->qideal!=NULL)
489     {
490       CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->qideal->m[0],
491                                            r->cf->extRing);
492       Variable a=rootOf(mipo);
493       CanonicalForm F( convSingAPFactoryAP( f,a,r ) ),
494                     G( convSingAPFactoryAP( g,a,r ) );
495       CanonicalForm FpG=F+G;
496       if (!(FpG.isUnivariate()|| FpG.inCoeffDomain()))
497       //if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
498       {
499         WerrorS("not univariate");
500         return TRUE;
501       }
502       res= convFactoryAPSingAP( extgcd( F, G, Fa, Gb ),r );
503       pa=convFactoryAPSingAP(Fa,r);
504       pb=convFactoryAPSingAP(Gb,r);
505       prune (a);
506     }
507     else
508     {
509       CanonicalForm F( convSingTrPFactoryP( f, r ) ), G( convSingTrPFactoryP( g, r ) );
510       CanonicalForm FpG=F+G;
511       if (!(FpG.isUnivariate()|| FpG.inCoeffDomain()))
512       //if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
513       {
514         Off(SW_RATIONAL);
515         WerrorS("not univariate");
516         return TRUE;
517       }
518       res= convFactoryPSingTrP( extgcd( F, G, Fa, Gb ), r );
519       pa=convFactoryPSingTrP(Fa, r);
520       pb=convFactoryPSingTrP(Gb, r);
521     }
522     Off(SW_RATIONAL);
523   }
524   else
525   {
526     WerrorS( feNotImplemented );
527     return TRUE;
528   }
529 #ifndef SING_NDEBUG
530   // checking the result of extgcd:
531   poly dummy;
532   dummy=p_Sub(p_Add_q(pp_Mult_qq(f,pa,r),pp_Mult_qq(g,pb,r),r),p_Copy(res,r),r);
533   if (dummy!=NULL)
534   {
535     PrintS("extgcd( ");p_Write(f,r);p_Write0(g,r);PrintS(" )\n");
536     PrintS("extgcd( ");p_Write(f,r);p_Write0(g,r);PrintS(" )\n");
537     p_Delete(&dummy,r);
538   }
539 #endif
540   return FALSE;
541 }
542 
singclap_pmult(poly f,poly g,const ring r)543 poly singclap_pmult ( poly f, poly g, const ring r )
544 {
545   poly res=NULL;
546   On(SW_RATIONAL);
547   if (rField_is_Zp(r) || rField_is_Q(r) || rField_is_Z(r)
548   || (rField_is_Zn(r)&&(r->cf->convSingNFactoryN!=ndConvSingNFactoryN)))
549   {
550     if (rField_is_Z(r)) Off(SW_RATIONAL);
551     setCharacteristic( rChar(r) );
552     CanonicalForm F( convSingPFactoryP( f,r ) ), G( convSingPFactoryP( g,r ) );
553     res = convFactoryPSingP( F * G,r );
554   }
555   else if (r->cf->extRing!=NULL)
556   {
557     if (rField_is_Q_a(r)) setCharacteristic( 0 );
558     else               setCharacteristic( rChar(r) );
559     if (r->cf->extRing->qideal!=NULL)
560     {
561       CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->qideal->m[0],
562                                                  r->cf->extRing);
563       Variable a=rootOf(mipo);
564       CanonicalForm F( convSingAPFactoryAP( f,a,r ) ),
565                     G( convSingAPFactoryAP( g,a,r ) );
566       res= convFactoryAPSingAP(  F * G, r  );
567       prune (a);
568     }
569     else
570     {
571       CanonicalForm F( convSingTrPFactoryP( f,r ) ), G( convSingTrPFactoryP( g,r ) );
572       res= convFactoryPSingTrP(  F * G,r  );
573     }
574   }
575 #if 0 // not yet working
576   else if (rField_is_GF())
577   {
578     //Print("GF(%d^%d)\n",nfCharP,nfMinPoly[0]);
579     setCharacteristic( nfCharP,nfMinPoly[0], currRing->parameter[0][0] );
580     CanonicalForm F( convSingGFFactoryGF( f ) ), G( convSingGFFactoryGF( g ) );
581     res = convFactoryGFSingGF( F * G );
582   }
583 #endif
584   else
585     WerrorS( feNotImplemented );
586   Off(SW_RATIONAL);
587   return res;
588 }
589 
singclap_pdivide(poly f,poly g,const ring r)590 poly singclap_pdivide ( poly f, poly g, const ring r )
591 {
592   poly res=NULL;
593 
594   #ifdef HAVE_FLINT
595   #if __FLINT_RELEASE >= 20503
596   /*
597     If the division is not exact, control will pass to factory where the
598     polynomials can be divided using the ordering that factory chooses.
599   */
600   if (rField_is_Zp(r))
601   {
602     nmod_mpoly_ctx_t ctx;
603     if (!convSingRFlintR(ctx,r))
604     {
605       res = Flint_Divide_MP(f,0,g,0,ctx,r);
606       if (res != NULL)
607         return res;
608     }
609   }
610   else
611   if (rField_is_Q(r))
612   {
613     fmpq_mpoly_ctx_t ctx;
614     if (!convSingRFlintR(ctx,r))
615     {
616       res = Flint_Divide_MP(f,0,g,0,ctx,r);
617       if (res != NULL)
618         return res;
619     }
620   }
621   #endif
622   #endif
623 
624   On(SW_RATIONAL);
625   if (rField_is_Zp(r) || rField_is_Q(r)
626   || (rField_is_Zn(r)&&(r->cf->convSingNFactoryN!=ndConvSingNFactoryN)))
627   {
628     setCharacteristic( rChar(r) );
629     CanonicalForm F( convSingPFactoryP( f,r ) ), G( convSingPFactoryP( g,r ) );
630     res = convFactoryPSingP( F / G,r );
631   }
632   // div is not implemented for ZZ coeffs in factory
633   else if (r->cf->extRing!=NULL)
634   {
635     if (rField_is_Q_a(r)) setCharacteristic( 0 );
636     else               setCharacteristic( rChar(r) );
637     if (r->cf->extRing->qideal!=NULL)
638     {
639       CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->qideal->m[0],
640                                                  r->cf->extRing);
641       Variable a=rootOf(mipo);
642       CanonicalForm F( convSingAPFactoryAP( f,a,r ) ),
643                     G( convSingAPFactoryAP( g,a,r ) );
644       res= convFactoryAPSingAP(  F / G, r  );
645       prune (a);
646     }
647     else
648     {
649       CanonicalForm F( convSingTrPFactoryP( f,r ) ), G( convSingTrPFactoryP( g,r ) );
650       res= convFactoryPSingTrP(  F / G,r  );
651     }
652   }
653 #if 0 // not yet working
654   else if (rField_is_GF())
655   {
656     //Print("GF(%d^%d)\n",nfCharP,nfMinPoly[0]);
657     setCharacteristic( nfCharP,nfMinPoly[0], currRing->parameter[0][0] );
658     CanonicalForm F( convSingGFFactoryGF( f ) ), G( convSingGFFactoryGF( g ) );
659     res = convFactoryGFSingGF( F / G );
660   }
661 #endif
662   else
663     WerrorS( feNotImplemented );
664   Off(SW_RATIONAL);
665   return res;
666 }
667 
singclap_pmod(poly f,poly g,const ring r)668 poly singclap_pmod ( poly f, poly g, const ring r )
669 {
670   poly res=NULL;
671   On(SW_RATIONAL);
672   if (rField_is_Zp(r) || rField_is_Q(r)
673   || (rField_is_Zn(r)&&(r->cf->convSingNFactoryN!=ndConvSingNFactoryN)))
674   {
675     setCharacteristic( rChar(r) );
676     CanonicalForm F( convSingPFactoryP( f,r ) ), G( convSingPFactoryP( g,r ) );
677     CanonicalForm Q,R;
678     divrem(F,G,Q,R);
679     res = convFactoryPSingP(R,r);
680     //res = convFactoryPSingP( F-(F/G)*G,r );
681   }
682   // mod is not implemented for ZZ coeffs in factory
683   else if (r->cf->extRing!=NULL)
684   {
685     if (rField_is_Q_a(r)) setCharacteristic( 0 );
686     else               setCharacteristic( rChar(r) );
687     if (r->cf->extRing->qideal!=NULL)
688     {
689       CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->qideal->m[0],
690                                                  r->cf->extRing);
691       Variable a=rootOf(mipo);
692       CanonicalForm F( convSingAPFactoryAP( f,a,r ) ),
693                     G( convSingAPFactoryAP( g,a,r ) );
694       CanonicalForm Q,R;
695       divrem(F,G,Q,R);
696       res = convFactoryAPSingAP(R,r);
697       //res= convFactoryAPSingAP( F-(F/G)*G, r  );
698       prune (a);
699     }
700     else
701     {
702       CanonicalForm F( convSingTrPFactoryP( f,r ) ), G( convSingTrPFactoryP( g,r ) );
703       CanonicalForm Q,R;
704       divrem(F,G,Q,R);
705       res = convFactoryPSingTrP(R,r);
706       //res= convFactoryPSingTrP(  F-(F/G)*G,r  );
707     }
708   }
709 #if 0 // not yet working
710   else if (rField_is_GF())
711   {
712     //Print("GF(%d^%d)\n",nfCharP,nfMinPoly[0]);
713     setCharacteristic( nfCharP,nfMinPoly[0], currRing->parameter[0][0] );
714     CanonicalForm F( convSingGFFactoryGF( f ) ), G( convSingGFFactoryGF( g ) );
715     res = convFactoryGFSingGF( F / G );
716   }
717 #endif
718   else
719     WerrorS( feNotImplemented );
720   Off(SW_RATIONAL);
721   return res;
722 }
723 
724 #if 0
725 // unused
726 void singclap_divide_content ( poly f, const ring r )
727 {
728   if ( f==NULL )
729   {
730     return;
731   }
732   else  if ( pNext( f ) == NULL )
733   {
734     p_SetCoeff( f, n_Init( 1, r->cf ), r );
735     return;
736   }
737   else
738   {
739     if ( rField_is_Q_a(r) )
740       setCharacteristic( 0 );
741     else  if ( rField_is_Zp_a(r) )
742       setCharacteristic( -rChar(r) );
743     else
744       return; /* not implemented*/
745 
746     CFList L;
747     CanonicalForm g, h;
748     poly p = pNext(f);
749 
750     // first attemp: find 2 smallest g:
751 
752     number g1=pGetCoeff(f);
753     number g2=pGetCoeff(p); // p==pNext(f);
754     pIter(p);
755     int sz1=n_Size(g1, r->cf);
756     int sz2=n_Size(g2, r->cf);
757     if (sz1>sz2)
758     {
759       number gg=g1;
760       g1=g2; g2=gg;
761       int sz=sz1;
762       sz1=sz2; sz2=sz;
763     }
764     while (p!=NULL)
765     {
766       int n_sz=n_Size(pGetCoeff(p),r->cf);
767       if (n_sz<sz1)
768       {
769         sz2=sz1;
770         g2=g1;
771         g1=pGetCoeff(p);
772         sz1=n_sz;
773         if (sz1<=3) break;
774       }
775       else if(n_sz<sz2)
776       {
777         sz2=n_sz;
778         g2=pGetCoeff(p);
779         sz2=n_sz;
780       }
781       pIter(p);
782     }
783     g = convSingPFactoryP( NUM(((fraction)g1)), r->cf->extRing );
784     g = gcd( g, convSingPFactoryP( NUM(((fraction)g2)) , r->cf->extRing));
785 
786     // second run: gcd's
787 
788     p = f;
789     while ( (p != NULL) && (g != 1)  && ( g != 0))
790     {
791       h = convSingPFactoryP( NUM(((fraction)pGetCoeff(p))), r->cf->extRing );
792       pIter( p );
793 
794       g = gcd( g, h );
795 
796       L.append( h );
797     }
798     if (( g == 1 ) || (g == 0))
799     {
800       // pTest(f);
801       return;
802     }
803     else
804     {
805       CFListIterator i;
806       for ( i = L, p = f; i.hasItem(); i++, p=pNext(p) )
807       {
808         fraction c=(fraction)pGetCoeff(p);
809         p_Delete(&(NUM(c)),r->cf->extRing); // 2nd arg used to be nacRing
810         NUM(c)=convFactoryPSingP( i.getItem() / g, r->cf->extRing );
811         //nTest((number)c);
812         //#ifdef LDEBUG
813         //number cn=(number)c;
814         //StringSetS(""); nWrite(nt); StringAppend(" ==> ");
815         //nWrite(cn);PrintS(StringEndS("\n")); // NOTE/TODO: use StringAppendS("\n"); omFree(s);
816         //#endif
817       }
818     }
819     // pTest(f);
820   }
821 }
822 #endif
823 
count_Factors(ideal I,intvec * v,int j,poly & f,poly fac,const ring r)824 BOOLEAN count_Factors(ideal I, intvec *v,int j, poly &f, poly fac, const ring r)
825 {
826   p_Test(f,r);
827   p_Test(fac,r);
828   int e=0;
829   if (!p_IsConstant(fac,r))
830   {
831 #ifdef FACTORIZE2_DEBUG
832     printf("start count_Factors(%d), Fdeg=%d, factor deg=%d\n",j,p_Totaldegree(f,r),p_Totaldegree(fac,r));
833     p_wrp(fac,r);PrintLn();
834 #endif
835     On(SW_RATIONAL);
836     CanonicalForm F, FAC,Q,R;
837     Variable a;
838     if (rField_is_Zp(r) || rField_is_Q(r)
839     || (rField_is_Zn(r)&&(r->cf->convSingNFactoryN!=ndConvSingNFactoryN)))
840     {
841       F=convSingPFactoryP( f,r );
842       FAC=convSingPFactoryP( fac,r );
843     }
844     else if (r->cf->extRing!=NULL)
845     {
846       if (r->cf->extRing->qideal!=NULL)
847       {
848         CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->qideal->m[0],
849                                     r->cf->extRing);
850         a=rootOf(mipo);
851         F=convSingAPFactoryAP( f,a,r );
852         FAC=convSingAPFactoryAP( fac,a,r );
853       }
854       else
855       {
856         F=convSingTrPFactoryP( f,r );
857         FAC=convSingTrPFactoryP( fac,r );
858       }
859     }
860     else
861       WerrorS( feNotImplemented );
862 
863     poly q;
864     loop
865     {
866       Q=F;
867       Q/=FAC;
868       R=Q;
869       R*=FAC;
870       R-=F;
871       if (R.isZero())
872       {
873         if (rField_is_Zp(r) || rField_is_Q(r)
874         || (rField_is_Zn(r)&&(r->cf->convSingNFactoryN!=ndConvSingNFactoryN)))
875         {
876           q = convFactoryPSingP( Q,r );
877         }
878         else if (r->cf->extRing!=NULL)
879         {
880           if (r->cf->extRing->qideal!=NULL)
881           {
882             q= convFactoryAPSingAP( Q,r );
883           }
884           else
885           {
886             q= convFactoryPSingTrP( Q,r );
887           }
888         }
889         e++; p_Delete(&f,r); f=q; q=NULL; F=Q;
890       }
891       else
892       {
893         break;
894       }
895     }
896     if (r->cf->extRing!=NULL)
897       if (r->cf->extRing->qideal!=NULL)
898         prune (a);
899     if (e==0)
900     {
901       Off(SW_RATIONAL);
902       return FALSE;
903     }
904   }
905   else e=1;
906   I->m[j]=fac;
907   if (v!=NULL) (*v)[j]=e;
908   Off(SW_RATIONAL);
909   return TRUE;
910 }
911 
912 VAR int singclap_factorize_retry;
913 
singclap_factorize(poly f,intvec ** v,int with_exps,const ring r)914 ideal singclap_factorize ( poly f, intvec ** v , int with_exps, const ring r)
915 /* destroys f, sets *v */
916 {
917   p_Test(f,r);
918 #ifdef FACTORIZE2_DEBUG
919   printf("singclap_factorize, degree %ld\n",p_Totaldegree(f,r));
920 #endif
921   // with_exps: 3,1 return only true factors, no exponents
922   //            2 return true factors and exponents
923   //            0 return coeff, factors and exponents
924   BOOLEAN save_errorreported=errorreported;
925 
926   ideal res=NULL;
927 
928   // handle factorize(0) =========================================
929   if (f==NULL)
930   {
931     res=idInit(1,1);
932     if (with_exps!=1)
933     {
934       (*v)=new intvec(1);
935       (**v)[0]=1;
936     }
937     return res;
938   }
939   // handle factorize(mon) =========================================
940   if (pNext(f)==NULL)
941   {
942     int i=0;
943     int n=0;
944     int e;
945     for(i=rVar(r);i>0;i--) if(p_GetExp(f,i,r)!=0) n++;
946     if (with_exps==0) n++; // with coeff
947     res=idInit(si_max(n,1),1);
948     switch(with_exps)
949     {
950       case 0: // with coef & exp.
951         res->m[0]=p_NSet(n_Copy(pGetCoeff(f),r->cf),r);
952         // no break
953       case 2: // with exp.
954         (*v)=new intvec(si_max(1,n));
955         (**v)[0]=1;
956         // no break
957       case 1: ;
958 #ifdef TEST
959       default: ;
960 #endif
961     }
962     if (n==0)
963     {
964       if (res->m[0]==NULL) res->m[0]=p_One(r);
965       // (**v)[0]=1; is already done
966     }
967     else
968     {
969       for(i=rVar(r);i>0;i--)
970       {
971         e=p_GetExp(f,i,r);
972         if(e!=0)
973         {
974           n--;
975           poly p=p_One(r);
976           p_SetExp(p,i,1,r);
977           p_Setm(p,r);
978           res->m[n]=p;
979           if (with_exps!=1) (**v)[n]=e;
980         }
981       }
982     }
983     p_Delete(&f,r);
984     return res;
985   }
986   //PrintS("S:");p_Write(f,r);PrintLn();
987   // use factory/libfac in general ==============================
988   Variable dummy(-1); prune(dummy); // remove all (tmp.) extensions
989   Off(SW_RATIONAL);
990   On(SW_SYMMETRIC_FF);
991   CFFList L;
992   number N=NULL;
993   number NN=NULL;
994   number old_lead_coeff=n_Copy(pGetCoeff(f), r->cf);
995 
996   Variable a;
997   if (r->cf->convSingNFactoryN!=ndConvSingNFactoryN)
998   {
999     if (rField_is_Q(r) || rField_is_Q_a(r) || rField_is_Z(r)) /* Q, Q(a), Z */
1000     {
1001       //if (f!=NULL) // already tested at start of routine
1002       {
1003         number n0=n_Copy(pGetCoeff(f),r->cf);
1004         if (with_exps==0)
1005           N=n_Copy(n0,r->cf);
1006         if (rField_is_Z(r)) p_Content(f, r);
1007         else                p_Cleardenom(f, r);
1008         //after here f should not have a denominator!! and no content
1009         //PrintS("S:");p_Write(f,r);PrintLn();
1010         NN=n_Div(n0,pGetCoeff(f),r->cf);
1011         n_Delete(&n0,r->cf);
1012         if (with_exps==0)
1013         {
1014           n_Delete(&N,r->cf);
1015           N=n_Copy(NN,r->cf);
1016         }
1017       }
1018     }
1019     else if (rField_is_Zp_a(r))
1020     {
1021       //if (f!=NULL) // already tested at start of routine
1022       if (singclap_factorize_retry==0)
1023       {
1024         number n0=n_Copy(pGetCoeff(f),r->cf);
1025         if (with_exps==0)
1026           N=n_Copy(n0,r->cf);
1027         p_Norm(f,r);
1028         p_Cleardenom(f, r);
1029         NN=n_Div(n0,pGetCoeff(f),r->cf);
1030         n_Delete(&n0,r->cf);
1031         if (with_exps==0)
1032         {
1033           n_Delete(&N,r->cf);
1034           N=n_Copy(NN,r->cf);
1035         }
1036       }
1037     }
1038     if (rField_is_Q(r) || rField_is_Zp(r) || rField_is_Z(r) || rField_is_Zn(r))
1039     {
1040       setCharacteristic( rChar(r) );
1041       if (errorreported) goto notImpl; // char too large
1042       CanonicalForm F( convSingPFactoryP( f,r ) );
1043       L = factorize( F );
1044     }
1045     // and over Q(a) / Fp(a)
1046     else if (r->cf->extRing!=NULL)
1047     {
1048       if (rField_is_Q_a (r)) setCharacteristic (0);
1049       else                   setCharacteristic( rChar(r) );
1050       if (errorreported) goto notImpl; // char too large
1051       if (r->cf->extRing->qideal!=NULL) /*algebraic extension */
1052       {
1053         CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->qideal->m[0],
1054                                              r->cf->extRing);
1055         a=rootOf(mipo);
1056         CanonicalForm F( convSingAPFactoryAP( f, a, r ) );
1057         L = factorize( F, a );
1058         prune(a);
1059       }
1060       else /* rational functions */
1061       {
1062         CanonicalForm F( convSingTrPFactoryP( f,r ) );
1063         L = factorize( F );
1064       }
1065     }
1066     else
1067     {
1068       goto notImpl;
1069     }
1070   }
1071   else
1072   {
1073     goto notImpl;
1074   }
1075   if (errorreported)
1076   {
1077     errorreported=FALSE;
1078   }
1079   {
1080     poly ff=p_Copy(f,r); // a copy for the retry stuff
1081     // the first factor should be a constant
1082     if ( ! L.getFirst().factor().inCoeffDomain() )
1083       L.insert(CFFactor(1,1));
1084     // convert into ideal
1085     int n = L.length();
1086     if (n==0) n=1;
1087     CFFListIterator J=L;
1088     int j=0;
1089     if (with_exps!=1)
1090     {
1091       if ((with_exps==2)&&(n>1))
1092       {
1093         n--;
1094         J++;
1095       }
1096       *v = new intvec( n );
1097     }
1098     res = idInit( n ,1);
1099     for ( ; J.hasItem(); J++, j++ )
1100     {
1101       if (with_exps!=1) (**v)[j] = J.getItem().exp();
1102       if (rField_is_Zp(r) || rField_is_Q(r)||  rField_is_Z(r)
1103       || rField_is_Zn(r))           /* Q, Fp, Z */
1104       {
1105         //count_Factors(res,*v,f, j, convFactoryPSingP( J.getItem().factor() );
1106         res->m[j] = convFactoryPSingP( J.getItem().factor(),r );
1107       }
1108 #if 0
1109       else if (rField_is_GF())
1110         res->m[j] = convFactoryGFSingGF( J.getItem().factor() );
1111 #endif
1112       else if (r->cf->extRing!=NULL)     /* Q(a), Fp(a) */
1113       {
1114 #ifndef SING_NDEBUG
1115         intvec *w=NULL;
1116         if (v!=NULL) w=*v;
1117 #endif
1118         if (r->cf->extRing->qideal==NULL)
1119         {
1120 #ifdef SING_NDEBUG
1121           res->m[j]= convFactoryPSingTrP( J.getItem().factor(),r );
1122 #else
1123           if(!count_Factors(res,w,j,ff,convFactoryPSingTrP( J.getItem().factor(),r ),r))
1124           {
1125             if (w!=NULL)
1126               (*w)[j]=1;
1127             res->m[j]=p_One(r);
1128           }
1129 #endif
1130         }
1131         else
1132         {
1133 #ifdef SING_NDEBUG
1134           res->m[j]= convFactoryAPSingAP( J.getItem().factor(),r );
1135 #else
1136           if (!count_Factors(res,w,j,ff,convFactoryAPSingAP( J.getItem().factor(),r ),r))
1137           {
1138             if (w!=NULL)
1139               (*w)[j]=1;
1140             res->m[j]=p_One(r);
1141           }
1142 #endif
1143         }
1144       }
1145     }
1146     if (r->cf->extRing!=NULL)
1147       if (r->cf->extRing->qideal!=NULL)
1148         prune (a);
1149 #ifndef SING_NDEBUG
1150     if ((r->cf->extRing!=NULL) && (!p_IsConstant(ff,r)))
1151     {
1152       singclap_factorize_retry++;
1153       if (singclap_factorize_retry<3)
1154       {
1155         int jj;
1156 #ifdef FACTORIZE2_DEBUG
1157         printf("factorize_retry\n");
1158 #endif
1159         intvec *ww=NULL;
1160         id_Test(res,r);
1161         ideal h=singclap_factorize ( ff, &ww , with_exps, r );
1162         id_Test(h,r);
1163         int l=(*v)->length();
1164         (*v)->resize(l+ww->length());
1165         for(jj=0;jj<ww->length();jj++)
1166           (**v)[jj+l]=(*ww)[jj];
1167         delete ww;
1168         ideal hh=idInit(IDELEMS(res)+IDELEMS(h),1);
1169         for(jj=IDELEMS(res)-1;jj>=0;jj--)
1170         {
1171           hh->m[jj]=res->m[jj];
1172           res->m[jj]=NULL;
1173         }
1174         for(jj=IDELEMS(h)-1;jj>=0;jj--)
1175         {
1176           hh->m[jj+IDELEMS(res)]=h->m[jj];
1177           h->m[jj]=NULL;
1178         }
1179         id_Delete(&res,r);
1180         id_Delete(&h,r);
1181         res=hh;
1182         id_Test(res,r);
1183         ff=NULL;
1184       }
1185       else
1186       {
1187         WarnS("problem with factorize");
1188 #if 0
1189         pWrite(ff);
1190         idShow(res);
1191 #endif
1192         id_Delete(&res,r);
1193         res=idInit(2,1);
1194         res->m[0]=p_One(r);
1195         res->m[1]=ff; ff=NULL;
1196       }
1197     }
1198 #endif
1199     p_Delete(&ff,r);
1200     if (N!=NULL)
1201     {
1202       __p_Mult_nn(res->m[0],N,r);
1203       n_Delete(&N,r->cf);
1204       N=NULL;
1205     }
1206     // delete constants
1207     if (res!=NULL)
1208     {
1209       int i=IDELEMS(res)-1;
1210       int j=0;
1211       for(;i>=0;i--)
1212       {
1213         if ((res->m[i]!=NULL)
1214         && (pNext(res->m[i])==NULL)
1215         && (p_IsConstant(res->m[i],r)))
1216         {
1217           if (with_exps!=0)
1218           {
1219             p_Delete(&(res->m[i]),r);
1220             if ((v!=NULL) && ((*v)!=NULL))
1221               (**v)[i]=0;
1222             j++;
1223           }
1224           else if (i!=0)
1225           {
1226             while ((v!=NULL) && ((*v)!=NULL) && ((**v)[i]>1))
1227             {
1228               res->m[0]=p_Mult_q(res->m[0],p_Copy(res->m[i],r),r);
1229               (**v)[i]--;
1230             }
1231             res->m[0]=p_Mult_q(res->m[0],res->m[i],r);
1232             res->m[i]=NULL;
1233             if ((v!=NULL) && ((*v)!=NULL))
1234               (**v)[i]=1;
1235             j++;
1236           }
1237         }
1238       }
1239       if (j>0)
1240       {
1241         idSkipZeroes(res);
1242         if ((v!=NULL) && ((*v)!=NULL))
1243         {
1244           intvec *w=*v;
1245           int len=IDELEMS(res);
1246           *v = new intvec( len );
1247           for (i=0,j=0;i<si_min(w->length(),len);i++)
1248           {
1249             if((*w)[i]!=0)
1250             {
1251               (**v)[j]=(*w)[i]; j++;
1252             }
1253           }
1254           delete w;
1255         }
1256       }
1257       if (res->m[0]==NULL)
1258       {
1259         res->m[0]=p_One(r);
1260       }
1261     }
1262   }
1263   if (rField_is_Q_a(r) && (r->cf->extRing->qideal!=NULL))
1264   {
1265     int i=IDELEMS(res)-1;
1266     int stop=1;
1267     if (with_exps!=0) stop=0;
1268     for(;i>=stop;i--)
1269     {
1270       p_Norm(res->m[i],r);
1271     }
1272     if (with_exps==0) p_SetCoeff(res->m[0],old_lead_coeff,r);
1273     else n_Delete(&old_lead_coeff,r->cf);
1274   }
1275   else
1276     n_Delete(&old_lead_coeff,r->cf);
1277   errorreported=save_errorreported;
1278 notImpl:
1279   prune(a);
1280   if (res==NULL)
1281   {
1282     WerrorS( feNotImplemented );
1283     if ((v!=NULL) && ((*v)!=NULL) &&(with_exps==2))
1284     {
1285        *v = new intvec( 1 );
1286        (*v)[0]=1;
1287     }
1288     res=idInit(2,1);
1289     res->m[0]=p_One(r);
1290     res->m[1]=f;
1291   }
1292   else p_Delete(&f,r);
1293   if (NN!=NULL)
1294   {
1295     n_Delete(&NN,r->cf);
1296   }
1297   if (N!=NULL)
1298   {
1299     n_Delete(&N,r->cf);
1300   }
1301   return res;
1302 }
1303 
singclap_sqrfree(poly f,intvec ** v,int with_exps,const ring r)1304 ideal singclap_sqrfree ( poly f, intvec ** v , int with_exps, const ring r)
1305 {
1306   p_Test(f,r);
1307 #ifdef FACTORIZE2_DEBUG
1308   printf("singclap_sqrfree, degree %d\n",pTotaldegree(f));
1309 #endif
1310   // with_exps: 3,1 return only true factors, no exponents
1311   //            2 return true factors and exponents
1312   //            0 return coeff, factors and exponents
1313   BOOLEAN save_errorreported=errorreported;
1314 
1315   ideal res=NULL;
1316 
1317   // handle factorize(0) =========================================
1318   if (f==NULL)
1319   {
1320     res=idInit(1,1);
1321     if (with_exps!=1 && with_exps!=3)
1322     {
1323       (*v)=new intvec(1);
1324       (**v)[0]=1;
1325     }
1326     return res;
1327   }
1328   // handle factorize(mon) =========================================
1329   if (pNext(f)==NULL)
1330   {
1331     int i=0;
1332     int n=0;
1333     int e;
1334     for(i=rVar(r);i>0;i--) if(p_GetExp(f,i,r)!=0) n++;
1335     if (with_exps==0 || with_exps==3) n++; // with coeff
1336     res=idInit(si_max(n,1),1);
1337     if(with_exps!=1)
1338     {
1339         (*v)=new intvec(si_max(1,n));
1340         (**v)[0]=1;
1341     }
1342     res->m[0]=p_NSet(n_Copy(pGetCoeff(f),r->cf),r);
1343     if (n==0)
1344     {
1345       res->m[0]=p_One(r);
1346       // (**v)[0]=1; is already done
1347     }
1348     else
1349     {
1350       for(i=rVar(r);i>0;i--)
1351       {
1352         e=p_GetExp(f,i,r);
1353         if(e!=0)
1354         {
1355           n--;
1356           poly p=p_One(r);
1357           p_SetExp(p,i,1,r);
1358           p_Setm(p,r);
1359           res->m[n]=p;
1360           if (with_exps!=1) (**v)[n]=e;
1361         }
1362       }
1363     }
1364     p_Delete(&f,r);
1365     return res;
1366   }
1367   //PrintS("S:");pWrite(f);PrintLn();
1368   // use factory/libfac in general ==============================
1369   Off(SW_RATIONAL);
1370   On(SW_SYMMETRIC_FF);
1371   CFFList L;
1372   number N=NULL;
1373   number NN=NULL;
1374   number old_lead_coeff=n_Copy(pGetCoeff(f), r->cf);
1375   Variable a;
1376 
1377   if (!rField_is_Zp(r) && !rField_is_Zp_a(r)) /* Q, Q(a) */
1378   {
1379     //if (f!=NULL) // already tested at start of routine
1380     number n0=n_Copy(old_lead_coeff,r->cf);
1381     if (with_exps==0 || with_exps==3)
1382       N=n_Copy(n0,r->cf);
1383     p_Cleardenom(f, r);
1384     //after here f should not have a denominator!!
1385     //PrintS("S:");p_Write(f,r);PrintLn();
1386     NN=n_Div(n0,pGetCoeff(f),r->cf);
1387     n_Delete(&n0,r->cf);
1388     if (with_exps==0 || with_exps==3)
1389     {
1390       n_Delete(&N,r->cf);
1391       N=n_Copy(NN,r->cf);
1392     }
1393   }
1394   else if (rField_is_Zp_a(r))
1395   {
1396     //if (f!=NULL) // already tested at start of routine
1397     if (singclap_factorize_retry==0)
1398     {
1399       number n0=n_Copy(old_lead_coeff,r->cf);
1400       if (with_exps==0 || with_exps==3)
1401         N=n_Copy(n0,r->cf);
1402       p_Norm(f,r);
1403       p_Cleardenom(f, r);
1404       NN=n_Div(n0,pGetCoeff(f),r->cf);
1405       n_Delete(&n0,r->cf);
1406       if (with_exps==0 || with_exps==3)
1407       {
1408         n_Delete(&N,r->cf);
1409         N=n_Copy(NN,r->cf);
1410       }
1411     }
1412   }
1413   if (rField_is_Q(r) || rField_is_Zp(r)
1414   || (rField_is_Zn(r)&&(r->cf->convSingNFactoryN!=ndConvSingNFactoryN)))
1415   {
1416     setCharacteristic( rChar(r) );
1417     CanonicalForm F( convSingPFactoryP( f,r ) );
1418     L = sqrFree( F );
1419   }
1420   else if (r->cf->extRing!=NULL)
1421   {
1422     if (rField_is_Q_a (r)) setCharacteristic (0);
1423     else                   setCharacteristic( rChar(r) );
1424     if (r->cf->extRing->qideal!=NULL)
1425     {
1426       CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->qideal->m[0],
1427                                            r->cf->extRing);
1428       a=rootOf(mipo);
1429       CanonicalForm F( convSingAPFactoryAP( f, a, r ) );
1430       L= sqrFree (F);
1431     }
1432     else
1433     {
1434       CanonicalForm F( convSingTrPFactoryP( f,r ) );
1435       L = sqrFree( F );
1436     }
1437   }
1438   #if 0
1439   else if (rField_is_GF())
1440   {
1441     int c=rChar(r);
1442     setCharacteristic( c, primepower(c) );
1443     CanonicalForm F( convSingGFFactoryGF( f ) );
1444     if (F.isUnivariate())
1445     {
1446       L = factorize( F );
1447     }
1448     else
1449     {
1450       goto notImpl;
1451     }
1452   }
1453   #endif
1454   else
1455   {
1456     goto notImpl;
1457   }
1458   {
1459     // convert into ideal
1460     int n = L.length();
1461     if (n==0) n=1;
1462     CFFListIterator J=L;
1463     int j=0;
1464     if (with_exps!=1)
1465     {
1466       if ((with_exps==2)&&(n>1))
1467       {
1468         n--;
1469         J++;
1470       }
1471       *v = new intvec( n );
1472     }
1473     else if (L.getFirst().factor().inCoeffDomain() && with_exps!=3)
1474     {
1475       n--;
1476       J++;
1477     }
1478     res = idInit( n ,1);
1479     for ( ; J.hasItem(); J++, j++ )
1480     {
1481       if (with_exps!=1 && with_exps!=3) (**v)[j] = J.getItem().exp();
1482       if (rField_is_Zp(r) || rField_is_Q(r)
1483       || (rField_is_Zn(r)&&(r->cf->convSingNFactoryN!=ndConvSingNFactoryN)))
1484         res->m[j] = convFactoryPSingP( J.getItem().factor(),r );
1485       else if (r->cf->extRing!=NULL)     /* Q(a), Fp(a) */
1486       {
1487         if (r->cf->extRing->qideal==NULL)
1488           res->m[j]=convFactoryPSingTrP( J.getItem().factor(),r );
1489         else
1490           res->m[j]=convFactoryAPSingAP( J.getItem().factor(),r );
1491       }
1492     }
1493     if (res->m[0]==NULL)
1494     {
1495       res->m[0]=p_One(r);
1496     }
1497     if (N!=NULL)
1498     {
1499       __p_Mult_nn(res->m[0],N,r);
1500       n_Delete(&N,r->cf);
1501       N=NULL;
1502     }
1503   }
1504   if (r->cf->extRing!=NULL)
1505     if (r->cf->extRing->qideal!=NULL)
1506       prune (a);
1507   if (rField_is_Q_a(r) && (r->cf->extRing->qideal!=NULL))
1508   {
1509     int i=IDELEMS(res)-1;
1510     int stop=1;
1511     if (with_exps!=0 || with_exps==3) stop=0;
1512     for(;i>=stop;i--)
1513     {
1514       p_Norm(res->m[i],r);
1515     }
1516     if (with_exps==0 || with_exps==3) p_SetCoeff(res->m[0],old_lead_coeff,r);
1517     else n_Delete(&old_lead_coeff,r->cf);
1518   }
1519   else
1520     n_Delete(&old_lead_coeff,r->cf);
1521   p_Delete(&f,r);
1522   errorreported=save_errorreported;
1523 notImpl:
1524   if (res==NULL)
1525     WerrorS( feNotImplemented );
1526   if (NN!=NULL)
1527   {
1528     n_Delete(&NN,r->cf);
1529   }
1530   if (N!=NULL)
1531   {
1532     n_Delete(&N,r->cf);
1533   }
1534   return res;
1535 }
1536 
singclap_irrCharSeries(ideal I,const ring r)1537 matrix singclap_irrCharSeries ( ideal I, const ring r)
1538 {
1539   if (idIs0(I)) return mpNew(1,1);
1540 
1541   // for now there is only the possibility to handle polynomials over
1542   // Q and Fp ...
1543   matrix res=NULL;
1544   int i;
1545   Off(SW_RATIONAL);
1546   On(SW_SYMMETRIC_FF);
1547   CFList L;
1548   ListCFList LL;
1549   if (rField_is_Q(r) || rField_is_Zp(r)
1550   || (rField_is_Zn(r)&&(r->cf->convSingNFactoryN!=ndConvSingNFactoryN)))
1551   {
1552     setCharacteristic( rChar(r) );
1553     for(i=0;i<IDELEMS(I);i++)
1554     {
1555       poly p=I->m[i];
1556       if (p!=NULL)
1557       {
1558         p=p_Copy(p,r);
1559         p_Cleardenom(p, r);
1560         L.append(convSingPFactoryP(p,r));
1561         p_Delete(&p,r);
1562       }
1563     }
1564   }
1565   // and over Q(a) / Fp(a)
1566   else if (nCoeff_is_transExt (r->cf))
1567   {
1568     setCharacteristic( rChar(r) );
1569     for(i=0;i<IDELEMS(I);i++)
1570     {
1571       poly p=I->m[i];
1572       if (p!=NULL)
1573       {
1574         p=p_Copy(p,r);
1575         p_Cleardenom(p, r);
1576         L.append(convSingTrPFactoryP(p,r));
1577         p_Delete(&p,r);
1578       }
1579     }
1580   }
1581   else
1582   {
1583     WerrorS( feNotImplemented );
1584     return res;
1585   }
1586 
1587   // a very bad work-around --- FIX IT in libfac
1588   // should be fixed as of 2001/6/27
1589   int tries=0;
1590   int m,n;
1591   ListIterator<CFList> LLi;
1592   loop
1593   {
1594     LL=irrCharSeries(L);
1595     m= LL.length(); // Anzahl Zeilen
1596     n=0;
1597     for ( LLi = LL; LLi.hasItem(); LLi++ )
1598     {
1599       n = si_max(LLi.getItem().length(),n);
1600     }
1601     if ((m!=0) && (n!=0)) break;
1602     tries++;
1603     if (tries>=5) break;
1604   }
1605   if ((m==0) || (n==0))
1606   {
1607     Warn("char_series returns %d x %d matrix from %d input polys (%d)",
1608       m,n,IDELEMS(I)+1,LL.length());
1609     iiWriteMatrix((matrix)I,"I",2,r,0);
1610     m=si_max(m,1);
1611     n=si_max(n,1);
1612   }
1613   res=mpNew(m,n);
1614   CFListIterator Li;
1615   for ( m=1, LLi = LL; LLi.hasItem(); LLi++, m++ )
1616   {
1617     for (n=1, Li = LLi.getItem(); Li.hasItem(); Li++, n++)
1618     {
1619       if (rField_is_Q(r) || rField_is_Zp(r)
1620       || (rField_is_Zn(r)&&(r->cf->convSingNFactoryN!=ndConvSingNFactoryN)))
1621         MATELEM(res,m,n)=convFactoryPSingP(Li.getItem(),r);
1622       else
1623         MATELEM(res,m,n)=convFactoryPSingTrP(Li.getItem(),r);
1624     }
1625   }
1626   Off(SW_RATIONAL);
1627   return res;
1628 }
1629 
singclap_neworder(ideal I,const ring r)1630 char* singclap_neworder ( ideal I, const ring r)
1631 {
1632   int i;
1633   Off(SW_RATIONAL);
1634   On(SW_SYMMETRIC_FF);
1635   CFList L;
1636   if (rField_is_Q(r) || rField_is_Zp(r)
1637   || (rField_is_Zn(r)&&(r->cf->convSingNFactoryN!=ndConvSingNFactoryN)))
1638   {
1639     setCharacteristic( rChar(r) );
1640     for(i=0;i<IDELEMS(I);i++)
1641     {
1642       poly p=I->m[i];
1643       if (p!=NULL)
1644       {
1645         p=p_Copy(p,r);
1646         p_Cleardenom(p, r);
1647         L.append(convSingPFactoryP(p,r));
1648       }
1649     }
1650   }
1651   // and over Q(a) / Fp(a)
1652   else if (nCoeff_is_transExt (r->cf))
1653   {
1654     setCharacteristic( rChar(r) );
1655     for(i=0;i<IDELEMS(I);i++)
1656     {
1657       poly p=I->m[i];
1658       if (p!=NULL)
1659       {
1660         p=p_Copy(p,r);
1661         p_Cleardenom(p, r);
1662         L.append(convSingTrPFactoryP(p,r));
1663       }
1664     }
1665   }
1666   else
1667   {
1668     WerrorS( feNotImplemented );
1669     return NULL;
1670   }
1671 
1672   List<int> IL=neworderint(L);
1673   ListIterator<int> Li;
1674   StringSetS("");
1675   Li = IL;
1676   int offs=rPar(r);
1677   int* mark=(int*)omAlloc0((rVar(r)+offs)*sizeof(int));
1678   int cnt=rVar(r)+offs;
1679   loop
1680   {
1681     if(! Li.hasItem()) break;
1682     BOOLEAN done=TRUE;
1683     i=Li.getItem()-1;
1684     mark[i]=1;
1685     if (i<offs)
1686     {
1687       done=FALSE;
1688       //StringAppendS(r->parameter[i]);
1689     }
1690     else
1691     {
1692       StringAppendS(r->names[i-offs]);
1693     }
1694     Li++;
1695     cnt--;
1696     if(cnt==0) break;
1697     if (done) StringAppendS(",");
1698   }
1699   for(i=0;i<rVar(r)+offs;i++)
1700   {
1701     BOOLEAN done=TRUE;
1702     if(mark[i]==0)
1703     {
1704       if (i<offs)
1705       {
1706         done=FALSE;
1707         //StringAppendS(r->parameter[i]);
1708       }
1709       else
1710       {
1711         StringAppendS(r->names[i-offs]);
1712       }
1713       cnt--;
1714       if(cnt==0) break;
1715       if (done) StringAppendS(",");
1716     }
1717   }
1718   char * s=StringEndS();
1719   if (s[strlen(s)-1]==',') s[strlen(s)-1]='\0';
1720   return s;
1721 }
1722 
singclap_det(const matrix m,const ring s)1723 poly singclap_det( const matrix m, const ring s )
1724 {
1725   int r=m->rows();
1726   if (r!=m->cols())
1727   {
1728     Werror("det of %d x %d matrix",r,m->cols());
1729     return NULL;
1730   }
1731   poly res=NULL;
1732   CFMatrix M(r,r);
1733   int i,j;
1734   for(i=r;i>0;i--)
1735   {
1736     for(j=r;j>0;j--)
1737     {
1738       M(i,j)=convSingPFactoryP(MATELEM(m,i,j),s);
1739     }
1740   }
1741   res= convFactoryPSingP( determinant(M,r),s ) ;
1742   Off(SW_RATIONAL);
1743   return res;
1744 }
1745 
singclap_det_i(intvec * m,const ring)1746 int singclap_det_i( intvec * m, const ring /*r*/)
1747 {
1748 //  assume( r == currRing ); // Anything else is not guaranted to work!
1749 
1750   setCharacteristic( 0 ); // ?
1751   CFMatrix M(m->rows(),m->cols());
1752   int i,j;
1753   for(i=m->rows();i>0;i--)
1754   {
1755     for(j=m->cols();j>0;j--)
1756     {
1757       M(i,j)=IMATELEM(*m,i,j);
1758     }
1759   }
1760   int res= convFactoryISingI( determinant(M,m->rows()) ) ;
1761   return res;
1762 }
1763 
singclap_det_bi(bigintmat * m,const coeffs cf)1764 number singclap_det_bi( bigintmat * m, const coeffs cf)
1765 {
1766   assume(m->basecoeffs()==cf);
1767   CFMatrix M(m->rows(),m->cols());
1768   int i,j;
1769   BOOLEAN setchar=TRUE;
1770   for(i=m->rows();i>0;i--)
1771   {
1772     for(j=m->cols();j>0;j--)
1773     {
1774       M(i,j)=n_convSingNFactoryN(BIMATELEM(*m,i,j),setchar,cf);
1775       setchar=FALSE;
1776     }
1777   }
1778   number res=n_convFactoryNSingN( determinant(M,m->rows()),cf ) ;
1779   return res;
1780 }
1781 
1782 #if defined(HAVE_NTL) || defined(AHVE_FLINT)
singntl_HNF(matrix m,const ring s)1783 matrix singntl_HNF(matrix  m, const ring s )
1784 {
1785   int r=m->rows();
1786   if (r!=m->cols())
1787   {
1788     Werror("HNF of %d x %d matrix",r,m->cols());
1789     return NULL;
1790   }
1791 
1792   matrix res=mpNew(r,r);
1793 
1794   if (rField_is_Q(s))
1795   {
1796 
1797     CFMatrix M(r,r);
1798     int i,j;
1799     for(i=r;i>0;i--)
1800     {
1801       for(j=r;j>0;j--)
1802       {
1803         M(i,j)=convSingPFactoryP(MATELEM(m,i,j),s );
1804       }
1805     }
1806     CFMatrix *MM=cf_HNF(M);
1807     for(i=r;i>0;i--)
1808     {
1809       for(j=r;j>0;j--)
1810       {
1811         MATELEM(res,i,j)=convFactoryPSingP((*MM)(i,j),s);
1812       }
1813     }
1814     delete MM;
1815   }
1816   return res;
1817 }
1818 
singntl_HNF(intvec * m)1819 intvec* singntl_HNF(intvec*  m)
1820 {
1821   int r=m->rows();
1822   if (r!=m->cols())
1823   {
1824     Werror("HNF of %d x %d matrix",r,m->cols());
1825     return NULL;
1826   }
1827   setCharacteristic( 0 );
1828   CFMatrix M(r,r);
1829   int i,j;
1830   for(i=r;i>0;i--)
1831   {
1832     for(j=r;j>0;j--)
1833     {
1834       M(i,j)=IMATELEM(*m,i,j);
1835     }
1836   }
1837   CFMatrix *MM=cf_HNF(M);
1838   intvec *mm=ivCopy(m);
1839   for(i=r;i>0;i--)
1840   {
1841     for(j=r;j>0;j--)
1842     {
1843       IMATELEM(*mm,i,j)=convFactoryISingI((*MM)(i,j));
1844     }
1845   }
1846   delete MM;
1847   return mm;
1848 }
1849 
singntl_HNF(bigintmat * b)1850 bigintmat* singntl_HNF(bigintmat*  b)
1851 {
1852   int r=b->rows();
1853   if (r!=b->cols())
1854   {
1855     Werror("HNF of %d x %d matrix",r,b->cols());
1856     return NULL;
1857   }
1858   setCharacteristic( 0 );
1859   CFMatrix M(r,r);
1860   int i,j;
1861   for(i=r;i>0;i--)
1862   {
1863     for(j=r;j>0;j--)
1864     {
1865       M(i,j)=n_convSingNFactoryN(BIMATELEM(*b,i,j),FALSE,b->basecoeffs());
1866     }
1867   }
1868   CFMatrix *MM=cf_HNF(M);
1869   bigintmat *mm=bimCopy(b);
1870   for(i=r;i>0;i--)
1871   {
1872     for(j=r;j>0;j--)
1873     {
1874       BIMATELEM(*mm,i,j)=n_convFactoryNSingN((*MM)(i,j),b->basecoeffs());
1875     }
1876   }
1877   delete MM;
1878   return mm;
1879 }
1880 
singntl_LLL(matrix m,const ring s)1881 matrix singntl_LLL(matrix  m, const ring s )
1882 {
1883   int r=m->rows();
1884   int c=m->cols();
1885   matrix res=mpNew(r,c);
1886   if (rField_is_Q(s))
1887   {
1888     CFMatrix M(r,c);
1889     int i,j;
1890     for(i=r;i>0;i--)
1891     {
1892       for(j=c;j>0;j--)
1893       {
1894         M(i,j)=convSingPFactoryP(MATELEM(m,i,j),s);
1895       }
1896     }
1897     CFMatrix *MM=cf_LLL(M);
1898     for(i=r;i>0;i--)
1899     {
1900       for(j=c;j>0;j--)
1901       {
1902         MATELEM(res,i,j)=convFactoryPSingP((*MM)(i,j),s);
1903       }
1904     }
1905     delete MM;
1906   }
1907   return res;
1908 }
1909 
singntl_LLL(intvec * m)1910 intvec* singntl_LLL(intvec*  m)
1911 {
1912   int r=m->rows();
1913   int c=m->cols();
1914   setCharacteristic( 0 );
1915   CFMatrix M(r,c);
1916   int i,j;
1917   for(i=r;i>0;i--)
1918   {
1919     for(j=c;j>0;j--)
1920     {
1921       M(i,j)=IMATELEM(*m,i,j);
1922     }
1923   }
1924   CFMatrix *MM=cf_LLL(M);
1925   intvec *mm=ivCopy(m);
1926   for(i=r;i>0;i--)
1927   {
1928     for(j=c;j>0;j--)
1929     {
1930       IMATELEM(*mm,i,j)=convFactoryISingI((*MM)(i,j));
1931     }
1932   }
1933   delete MM;
1934   return mm;
1935 }
1936 #else
singntl_HNF(matrix m,const ring s)1937 matrix singntl_HNF(matrix  m, const ring s )
1938 {
1939   WerrorS("NTL/FLINT missing");
1940   return NULL;
1941 }
1942 
singntl_HNF(intvec * m)1943 intvec* singntl_HNF(intvec*  m)
1944 {
1945   WerrorS("NTL/FLINT missing");
1946   return NULL;
1947 }
1948 
singntl_LLL(matrix m,const ring s)1949 matrix singntl_LLL(matrix  m, const ring s )
1950 {
1951   WerrorS("NTL/FLINT missing");
1952   return NULL;
1953 }
1954 
singntl_LLL(intvec * m)1955 intvec* singntl_LLL(intvec*  m)
1956 {
1957   WerrorS("NTL/FLINT missing");
1958   return NULL;
1959 }
1960 #endif
1961 
1962 #if defined(HAVE_NTL) || defined(HAVE_FLINT)
singclap_absFactorize(poly f,ideal & mipos,intvec ** exps,int & numFactors,const ring r)1963 ideal singclap_absFactorize ( poly f, ideal & mipos, intvec ** exps, int & numFactors, const ring r)
1964 {
1965   p_Test(f, r);
1966 
1967   ideal res=NULL;
1968 
1969   int offs = rPar(r);
1970   if (f==NULL)
1971   {
1972     res= idInit (1, 1);
1973     mipos= idInit (1, 1);
1974     mipos->m[0]= convFactoryPSingTrP (Variable (offs), r); //overkill
1975     (*exps)=new intvec (1);
1976     (**exps)[0]= 1;
1977     numFactors= 0;
1978     return res;
1979   }
1980   CanonicalForm F( convSingTrPFactoryP( f, r) );
1981 
1982   bool isRat= isOn (SW_RATIONAL);
1983   if (!isRat)
1984     On (SW_RATIONAL);
1985 
1986   CFAFList absFactors= absFactorize (F);
1987 
1988   int n= absFactors.length();
1989   *exps=new intvec (n);
1990 
1991   res= idInit (n, 1);
1992 
1993   mipos= idInit (n, 1);
1994 
1995   Variable x= Variable (offs);
1996   Variable alpha;
1997   int i= 0;
1998   numFactors= 0;
1999   int count;
2000   CFAFListIterator iter= absFactors;
2001   CanonicalForm lead= iter.getItem().factor();
2002   if (iter.getItem().factor().inCoeffDomain())
2003   {
2004     i++;
2005     iter++;
2006   }
2007   for (; iter.hasItem(); iter++, i++)
2008   {
2009     (**exps)[i]= iter.getItem().exp();
2010     alpha= iter.getItem().minpoly().mvar();
2011     if (iter.getItem().minpoly().isOne())
2012       lead /= power (bCommonDen (iter.getItem().factor()), iter.getItem().exp());
2013     else
2014       lead /= power (power (bCommonDen (iter.getItem().factor()), degree (iter.getItem().minpoly())), iter.getItem().exp());
2015     res->m[i]= convFactoryPSingTrP (replacevar (iter.getItem().factor()*bCommonDen (iter.getItem().factor()), alpha, x), r);
2016     if (iter.getItem().minpoly().isOne())
2017     {
2018       count= iter.getItem().exp();
2019       mipos->m[i]= convFactoryPSingTrP (x,r);
2020     }
2021     else
2022     {
2023       count= iter.getItem().exp()*degree (iter.getItem().minpoly());
2024       mipos->m[i]= convFactoryPSingTrP (replacevar (iter.getItem().minpoly(), alpha, x), r);
2025     }
2026     if (!iter.getItem().minpoly().isOne())
2027       prune (alpha);
2028     numFactors += count;
2029   }
2030   if (!isRat)
2031     Off (SW_RATIONAL);
2032 
2033   (**exps)[0]= 1;
2034   res->m[0]= convFactoryPSingTrP (lead, r);
2035   mipos->m[0]= convFactoryPSingTrP (x, r);
2036   return res;
2037 }
2038 
2039 #else
singclap_absFactorize(poly f,ideal & mipos,intvec ** exps,int & numFactors,const ring r)2040 ideal singclap_absFactorize ( poly f, ideal & mipos, intvec ** exps, int & numFactors, const ring r)
2041 {
2042   WerrorS("NTL/FLINT missing: absFactorize");
2043   return NULL;
2044 }
2045 
2046 #endif /* HAVE_NTL */
2047 
Zp_roots(poly p,const ring r)2048 int * Zp_roots(poly p, const ring r)
2049 {
2050   CanonicalForm pp=convSingPFactoryP(p,r);
2051   return Zp_roots(pp);
2052 }
2053