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

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