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

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