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