source: git/libpolys/polys/clapsing.cc @ 5c0bf0

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