source: git/libpolys/polys/clapsing.cc @ cd964c

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