source: git/libpolys/polys/clapsing.cc @ 255be23

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