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

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