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

spielwiese
Last change on this file since ea1d44c was ea1d44c, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: count factors only in debug version
  • Property mode set to 100644
File size: 35.3 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
503static int primepower(int c, const ring r)
504{
505  int p=1;
506  int cc=c;
507  while(cc!= rChar(r)) { cc*=c; p++; }
508  return p;
509}
510
511BOOLEAN count_Factors(ideal I, intvec *v,int j, poly &f, poly fac, const ring r)
512{
513  p_Test(f,r);
514  p_Test(fac,r);
515  int e=0;
516  if (!p_IsConstantPoly(fac,r))
517  {
518#ifdef FACTORIZE2_DEBUG
519    printf("start count_Factors(%d), Fdeg=%d, factor deg=%d\n",j,p_Totaldegree(f,r),p_Totaldegree(fac,r));
520    p_wrp(fac,r);PrintLn();
521#endif
522    On(SW_RATIONAL);
523    CanonicalForm F, FAC,Q,R;
524    Variable a;
525    if (rField_is_Zp(r) || rField_is_Q(r))
526    {
527      F=convSingPFactoryP( f,r );
528      FAC=convSingPFactoryP( fac,r );
529    }
530    else if (rField_is_Extension(r))
531    {
532      if (r->cf->extRing->qideal!=NULL)
533      {
534        CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->qideal->m[0],
535                                    r->cf->extRing);
536        a=rootOf(mipo);
537        F=convSingAPFactoryAP( f,a,r );
538        FAC=convSingAPFactoryAP( fac,a,r );
539      }
540      else
541      {
542        F=convSingTrPFactoryP( f,r );
543        FAC=convSingTrPFactoryP( fac,r );
544      }
545    }
546    else
547      WerrorS( feNotImplemented );
548
549    poly q;
550    loop
551    {
552      Q=F;
553      Q/=FAC;
554      R=Q;
555      R*=FAC;
556      R-=F;
557      if (R.isZero())
558      {
559        if (rField_is_Zp(r) || rField_is_Q(r))
560        {
561          q = convFactoryPSingP( Q,r );
562        }
563        else if (rField_is_Extension(r))
564        {
565          if (r->cf->extRing->qideal!=NULL)
566          {
567            q= convFactoryAPSingAP( Q,r );
568          }
569          else
570          {
571            q= convFactoryPSingTrP( Q,r );
572          }
573        }
574        e++; p_Delete(&f,r); f=q; q=NULL; F=Q;
575      }
576      else
577      {
578        break;
579      }
580    }
581    if (e==0)
582    {
583      Off(SW_RATIONAL);
584      return FALSE;
585    }
586  }
587  else e=1;
588  I->m[j]=fac;
589  if (v!=NULL) (*v)[j]=e;
590  Off(SW_RATIONAL);
591  return TRUE;
592}
593
594#ifdef HAVE_FACTORY
595int singclap_factorize_retry;
596# ifdef HAVE_LIBFAC
597  extern int libfac_interruptflag;
598# endif
599#endif
600
601ideal singclap_factorize ( poly f, intvec ** v , int with_exps, const ring r)
602/* destroys f, sets *v */
603{
604  p_Test(f,r);
605#ifdef FACTORIZE2_DEBUG
606  printf("singclap_factorize, degree %ld\n",p_Totaldegree(f,r));
607#endif
608  // with_exps: 3,1 return only true factors, no exponents
609  //            2 return true factors and exponents
610  //            0 return coeff, factors and exponents
611  BOOLEAN save_errorreported=errorreported;
612
613  ideal res=NULL;
614
615  // handle factorize(0) =========================================
616  if (f==NULL)
617  {
618    res=idInit(1,1);
619    if (with_exps!=1)
620    {
621      (*v)=new intvec(1);
622      (**v)[0]=1;
623    }
624    return res;
625  }
626  // handle factorize(mon) =========================================
627  if (pNext(f)==NULL)
628  {
629    int i=0;
630    int n=0;
631    int e;
632    for(i=rVar(r);i>0;i--) if(p_GetExp(f,i,r)!=0) n++;
633    if (with_exps==0) n++; // with coeff
634    res=idInit(si_max(n,1),1);
635    switch(with_exps)
636    {
637      case 0: // with coef & exp.
638        res->m[0]=p_NSet(n_Copy(pGetCoeff(f),r->cf),r);
639        // no break
640      case 2: // with exp.
641        (*v)=new intvec(si_max(1,n));
642        (**v)[0]=1;
643        // no break
644      case 1: ;
645      #ifdef TEST
646      default: ;
647      #endif
648    }
649    if (n==0)
650    {
651      res->m[0]=p_One(r);
652      // (**v)[0]=1; is already done
653    }
654    else
655    {
656      for(i=rVar(r);i>0;i--)
657      {
658        e=p_GetExp(f,i,r);
659        if(e!=0)
660        {
661          n--;
662          poly p=p_One(r);
663          p_SetExp(p,i,1,r);
664          p_Setm(p,r);
665          res->m[n]=p;
666          if (with_exps!=1) (**v)[n]=e;
667        }
668      }
669    }
670    p_Delete(&f,r);
671    return res;
672  }
673  //PrintS("S:");p_Write(f,r);PrintLn();
674  // use factory/libfac in general ==============================
675  Off(SW_RATIONAL);
676  On(SW_SYMMETRIC_FF);
677  #ifdef HAVE_NTL
678  extern int prime_number;
679  if(rField_is_Q(r)) prime_number=0;
680  #endif
681  CFFList L;
682  number N=NULL;
683  number NN=NULL;
684  number old_lead_coeff=n_Copy(pGetCoeff(f), r->cf);
685
686  if (!rField_is_Zp(r) && !rField_is_Zp_a(r)) /* Q, Q(a) */
687  {
688    //if (f!=NULL) // already tested at start of routine
689    {
690      number n0=n_Copy(pGetCoeff(f),r->cf);
691      if (with_exps==0)
692        N=n_Copy(n0,r->cf);
693      p_Cleardenom(f, r);
694      //after here f should not have a denominator!!
695      //PrintS("S:");p_Write(f,r);PrintLn();
696      NN=n_Div(n0,pGetCoeff(f),r->cf);
697      n_Delete(&n0,r->cf);
698      if (with_exps==0)
699      {
700        n_Delete(&N,r->cf);
701        N=n_Copy(NN,r->cf);
702      }
703    }
704  }
705  else if (rField_is_Zp_a(r))
706  {
707    //if (f!=NULL) // already tested at start of routine
708    if (singclap_factorize_retry==0)
709    {
710      number n0=n_Copy(pGetCoeff(f),r->cf);
711      if (with_exps==0)
712        N=n_Copy(n0,r->cf);
713      p_Norm(f,r);
714      p_Cleardenom(f, r);
715      NN=n_Div(n0,pGetCoeff(f),r->cf);
716      n_Delete(&n0,r->cf);
717      if (with_exps==0)
718      {
719        n_Delete(&N,r->cf);
720        N=n_Copy(NN,r->cf);
721      }
722    }
723  }
724  if (rField_is_Q(r) || rField_is_Zp(r))
725  {
726    setCharacteristic( rChar(r) );
727    CanonicalForm F( convSingPFactoryP( f,r ) );
728    L = factorize( F );
729  }
730  // and over Q(a) / Fp(a)
731  else if (rField_is_Extension(r))
732  {
733    if (rField_is_Q_a (r)) setCharacteristic (0);
734    else                   setCharacteristic( rChar(r) );
735    if (r->cf->extRing->qideal!=NULL)
736    {
737      CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->qideal->m[0],
738                                           r->cf->extRing);
739      Variable a=rootOf(mipo);
740      CanonicalForm F( convSingAPFactoryAP( f, a, r ) );
741      if (rField_is_Zp_a(r))
742      {
743        L = factorize( F, a );
744      }
745      else
746      {
747        //  over Q(a)
748        L= factorize (F, a);
749      }
750    }
751    else
752    {
753      CanonicalForm F( convSingTrPFactoryP( f,r ) );
754      L = factorize( F );
755    }
756  }
757  else
758  {
759    goto notImpl;
760  }
761  {
762    poly ff=p_Copy(f,r); // a copy for the retry stuff
763    // the first factor should be a constant
764    if ( ! L.getFirst().factor().inCoeffDomain() )
765      L.insert(CFFactor(1,1));
766    // convert into ideal
767    int n = L.length();
768    if (n==0) n=1;
769    CFFListIterator J=L;
770    int j=0;
771    if (with_exps!=1)
772    {
773      if ((with_exps==2)&&(n>1))
774      {
775        n--;
776        J++;
777      }
778      *v = new intvec( n );
779    }
780    res = idInit( n ,1);
781    for ( ; J.hasItem(); J++, j++ )
782    {
783      if (with_exps!=1) (**v)[j] = J.getItem().exp();
784      if (rField_is_Zp(r) || rField_is_Q(r))           /* Q, Fp */
785      {
786        //count_Factors(res,*v,f, j, convFactoryPSingP( J.getItem().factor() );
787        res->m[j] = convFactoryPSingP( J.getItem().factor(),r );
788      }
789      #if 0
790      else if (rField_is_GF())
791        res->m[j] = convFactoryGFSingGF( J.getItem().factor() );
792      #endif
793      else if (rField_is_Extension(r))     /* Q(a), Fp(a) */
794      {
795#ifndef NDEBUG
796        intvec *w=NULL;
797        if (v!=NULL) w=*v;
798#endif
799        if (r->cf->extRing->qideal==NULL)
800        {
801#ifdef NDEBUG
802          res->m[j]= convFactoryPSingTrP( J.getItem().factor(),r );
803#else
804          if(!count_Factors(res,w,j,ff,convFactoryPSingTrP( J.getItem().factor(),r ),r))
805          {
806            if (w!=NULL)
807              (*w)[j]=1;
808            res->m[j]=p_One(r);
809          }
810#endif
811        }
812        else
813        {
814#ifdef NDEBUG
815          res->m[j]= convFactoryAPSingAP( J.getItem().factor(),r );
816#else
817          if (!count_Factors(res,w,j,ff,convFactoryAPSingAP( J.getItem().factor(),r ),r))
818          {
819            if (w!=NULL)
820              (*w)[j]=1;
821            res->m[j]=p_One(r);
822          }
823#endif
824        }
825      }
826    }
827#ifndef NDEBUG
828    if (rField_is_Extension(r) && (!p_IsConstantPoly(ff,r)))
829    {
830      singclap_factorize_retry++;
831      if (singclap_factorize_retry<3)
832      {
833        int jj;
834        #ifdef FACTORIZE2_DEBUG
835        printf("factorize_retry\n");
836        #endif
837        intvec *ww=NULL;
838        id_Test(res,r);
839        ideal h=singclap_factorize ( ff, &ww , with_exps, r );
840        id_Test(h,r);
841        int l=(*v)->length();
842        (*v)->resize(l+ww->length());
843        for(jj=0;jj<ww->length();jj++)
844          (**v)[jj+l]=(*ww)[jj];
845        delete ww;
846        ideal hh=idInit(IDELEMS(res)+IDELEMS(h),1);
847        for(jj=IDELEMS(res)-1;jj>=0;jj--)
848        {
849          hh->m[jj]=res->m[jj];
850          res->m[jj]=NULL;
851        }
852        for(jj=IDELEMS(h)-1;jj>=0;jj--)
853        {
854          hh->m[jj+IDELEMS(res)]=h->m[jj];
855          h->m[jj]=NULL;
856        }
857        id_Delete(&res,r);
858        id_Delete(&h,r);
859        res=hh;
860        id_Test(res,r);
861        ff=NULL;
862      }
863      else
864      {
865        WarnS("problem with factorize");
866        #if 0
867        pWrite(ff);
868        idShow(res);
869        #endif
870        id_Delete(&res,r);
871        res=idInit(2,1);
872        res->m[0]=p_One(r);
873        res->m[1]=ff; ff=NULL;
874      }
875    }
876#endif
877    p_Delete(&ff,r);
878    if (N!=NULL)
879    {
880      p_Mult_nn(res->m[0],N,r);
881      n_Delete(&N,r->cf);
882      N=NULL;
883    }
884    // delete constants
885    if (res!=NULL)
886    {
887      int i=IDELEMS(res)-1;
888      int j=0;
889      for(;i>=0;i--)
890      {
891        if ((res->m[i]!=NULL)
892        && (pNext(res->m[i])==NULL)
893        && (p_IsConstant(res->m[i],r)))
894        {
895          if (with_exps!=0)
896          {
897            p_Delete(&(res->m[i]),r);
898            if ((v!=NULL) && ((*v)!=NULL))
899              (**v)[i]=0;
900            j++;
901          }
902          else if (i!=0)
903          {
904            while ((v!=NULL) && ((*v)!=NULL) && ((**v)[i]>1))
905            {
906              res->m[0]=p_Mult_q(res->m[0],p_Copy(res->m[i],r),r);
907              (**v)[i]--;
908            }
909            res->m[0]=p_Mult_q(res->m[0],res->m[i],r);
910            res->m[i]=NULL;
911            if ((v!=NULL) && ((*v)!=NULL))
912              (**v)[i]=1;
913            j++;
914          }
915        }
916      }
917      if (j>0)
918      {
919        idSkipZeroes(res);
920        if ((v!=NULL) && ((*v)!=NULL))
921        {
922          intvec *w=*v;
923          int len=IDELEMS(res);
924          *v = new intvec( len );
925          for (i=0,j=0;i<si_min(w->length(),len);i++)
926          {
927            if((*w)[i]!=0)
928            {
929              (**v)[j]=(*w)[i]; j++;
930            }
931          }
932          delete w;
933        }
934      }
935      if (res->m[0]==NULL)
936      {
937        res->m[0]=p_One(r);
938      }
939    }
940  }
941  if (rField_is_Q_a(r) && (r->cf->extRing->qideal!=NULL))
942  {
943    int i=IDELEMS(res)-1;
944    int stop=1;
945    if (with_exps!=0) stop=0;
946    for(;i>=stop;i--)
947    {
948      p_Norm(res->m[i],r);
949    }
950    if (with_exps==0) p_SetCoeff(res->m[0],old_lead_coeff,r);
951    else n_Delete(&old_lead_coeff,r->cf);
952  }
953  else
954    n_Delete(&old_lead_coeff,r->cf);
955  errorreported=save_errorreported;
956notImpl:
957  if (res==NULL)
958    WerrorS( feNotImplemented );
959  if (NN!=NULL)
960  {
961    n_Delete(&NN,r->cf);
962  }
963  if (N!=NULL)
964  {
965    n_Delete(&N,r->cf);
966  }
967  if (f!=NULL) p_Delete(&f,r);
968  //PrintS("......S\n");
969  return res;
970}
971ideal singclap_sqrfree ( poly f, intvec ** v , int with_exps, const ring r)
972{
973  p_Test(f,r);
974#ifdef FACTORIZE2_DEBUG
975  printf("singclap_sqrfree, degree %d\n",pTotaldegree(f));
976#endif
977  // with_exps: 3,1 return only true factors, no exponents
978  //            2 return true factors and exponents
979  //            0 return coeff, factors and exponents
980  BOOLEAN save_errorreported=errorreported;
981
982  ideal res=NULL;
983
984  // handle factorize(0) =========================================
985  if (f==NULL)
986  {
987    res=idInit(1,1);
988    if (with_exps!=1 && with_exps!=3)
989    {
990      (*v)=new intvec(1);
991      (**v)[0]=1;
992    }
993    return res;
994  }
995  // handle factorize(mon) =========================================
996  if (pNext(f)==NULL)
997  {
998    int i=0;
999    int n=0;
1000    int e;
1001    for(i=rVar(r);i>0;i--) if(p_GetExp(f,i,r)!=0) n++;
1002    if (with_exps==0 || with_exps==3) n++; // with coeff
1003    res=idInit(si_max(n,1),1);
1004    switch(with_exps)
1005    {
1006      case 0: // with coef & exp.
1007        res->m[0]=p_NSet(n_Copy(pGetCoeff(f),r->cf),r);
1008        // no break
1009      case 3: // with coef & exp.
1010        res->m[0]=p_NSet(n_Copy(pGetCoeff(f),r->cf),r);
1011        // no break
1012      case 2: // with exp.
1013        (*v)=new intvec(si_max(1,n));
1014        (**v)[0]=1;
1015        // no break
1016      case 1: ;
1017      #ifdef TEST
1018      default: ;
1019      #endif
1020    }
1021    res->m[0]=p_NSet(n_Copy(pGetCoeff(f),r->cf),r);
1022    if (n==0)
1023    {
1024      res->m[0]=p_One(r);
1025      // (**v)[0]=1; is already done
1026    }
1027    else
1028    {
1029      for(i=rVar(r);i>0;i--)
1030      {
1031        e=p_GetExp(f,i,r);
1032        if(e!=0)
1033        {
1034          n--;
1035          poly p=p_One(r);
1036          p_SetExp(p,i,1,r);
1037          p_Setm(p,r);
1038          res->m[n]=p;
1039          if (with_exps!=1) (**v)[n]=e;
1040        }
1041      }
1042    }
1043    p_Delete(&f,r);
1044    return res;
1045  }
1046  //PrintS("S:");pWrite(f);PrintLn();
1047  // use factory/libfac in general ==============================
1048  Off(SW_RATIONAL);
1049  On(SW_SYMMETRIC_FF);
1050  #ifdef HAVE_NTL
1051  extern int prime_number;
1052  if(rField_is_Q(r)) prime_number=0;
1053  #endif
1054  CFFList L;
1055  number N=NULL;
1056  number NN=NULL;
1057  number old_lead_coeff=n_Copy(pGetCoeff(f), r->cf);
1058
1059  if (!rField_is_Zp(r) && !rField_is_Zp_a(r)) /* Q, Q(a) */
1060  {
1061    //if (f!=NULL) // already tested at start of routine
1062    number n0=n_Copy(pGetCoeff(f),r->cf);
1063    if (with_exps==0 || with_exps==3)
1064      N=n_Copy(n0,r->cf);
1065    p_Cleardenom(f, r);
1066    //after here f should not have a denominator!!
1067    //PrintS("S:");p_Write(f,r);PrintLn();
1068    NN=n_Div(n0,pGetCoeff(f),r->cf);
1069    n_Delete(&n0,r->cf);
1070    if (with_exps==0 || with_exps==3)
1071    {
1072      n_Delete(&N,r->cf);
1073      N=n_Copy(NN,r->cf);
1074    }
1075  }
1076  else if (rField_is_Zp_a(r))
1077  {
1078    //if (f!=NULL) // already tested at start of routine
1079    if (singclap_factorize_retry==0)
1080    {
1081      number n0=n_Copy(pGetCoeff(f),r->cf);
1082      if (with_exps==0 || with_exps==3)
1083        N=n_Copy(n0,r->cf);
1084      p_Norm(f,r);
1085      p_Cleardenom(f, r);
1086      NN=n_Div(n0,pGetCoeff(f),r->cf);
1087      n_Delete(&n0,r->cf);
1088      if (with_exps==0 || with_exps==3)
1089      {
1090        n_Delete(&N,r->cf);
1091        N=n_Copy(NN,r->cf);
1092      }
1093    }
1094  }
1095  if (rField_is_Q(r) || rField_is_Zp(r))
1096  {
1097    setCharacteristic( rChar(r) );
1098    CanonicalForm F( convSingPFactoryP( f,r ) );
1099    L = sqrFree( F );
1100  }
1101  else if (rField_is_Extension(r))
1102  {
1103    if (rField_is_Q_a (r)) setCharacteristic (0);
1104    else                   setCharacteristic( rChar(r) );
1105    if (r->cf->extRing->qideal!=NULL)
1106    {
1107      CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->qideal->m[0],
1108                                           r->cf->extRing);
1109      Variable a=rootOf(mipo);
1110      CanonicalForm F( convSingAPFactoryAP( f, a, r ) );
1111      L= sqrFree (F);
1112    }
1113    else
1114    {
1115      CanonicalForm F( convSingTrPFactoryP( f,r ) );
1116      L = sqrFree( F );
1117    }
1118  }
1119  #if 0
1120  else if (rField_is_GF())
1121  {
1122    int c=rChar(r);
1123    setCharacteristic( c, primepower(c) );
1124    CanonicalForm F( convSingGFFactoryGF( f ) );
1125    if (F.isUnivariate())
1126    {
1127      L = factorize( F );
1128    }
1129    else
1130    {
1131      goto notImpl;
1132    }
1133  }
1134  #endif
1135  else
1136  {
1137    goto notImpl;
1138  }
1139  {
1140    // convert into ideal
1141    int n = L.length();
1142    if (n==0) n=1;
1143    CFFListIterator J=L;
1144    int j=0;
1145    if (with_exps!=1)
1146    {
1147      if ((with_exps==2)&&(n>1))
1148      {
1149        n--;
1150        J++;
1151      }
1152      *v = new intvec( n );
1153    }
1154    else if (L.getFirst().factor().inCoeffDomain() && with_exps!=3)
1155    {
1156      n--;
1157      J++;
1158    }
1159    res = idInit( n ,1);
1160    for ( ; J.hasItem(); J++, j++ )
1161    {
1162      if (with_exps!=1 && with_exps!=3) (**v)[j] = J.getItem().exp();
1163      if (rField_is_Zp(r) || rField_is_Q(r))
1164        res->m[j] = convFactoryPSingP( J.getItem().factor(),r );
1165      else if (rField_is_Extension(r))     /* Q(a), Fp(a) */
1166      {
1167        if (r->cf->extRing->qideal==NULL)
1168          res->m[j]=convFactoryPSingTrP( J.getItem().factor(),r );
1169        else
1170          res->m[j]=convFactoryAPSingAP( J.getItem().factor(),r );
1171      }
1172    }
1173    if (res->m[0]==NULL)
1174    {
1175      res->m[0]=p_One(r);
1176    }
1177    if (N!=NULL)
1178    {
1179      p_Mult_nn(res->m[0],N,r);
1180      n_Delete(&N,r->cf);
1181      N=NULL;
1182    }
1183  }
1184  if (rField_is_Q_a(r) && (r->cf->extRing->qideal!=NULL))
1185  {
1186    int i=IDELEMS(res)-1;
1187    int stop=1;
1188    if (with_exps!=0 || with_exps==3) stop=0;
1189    for(;i>=stop;i--)
1190    {
1191      p_Norm(res->m[i],r);
1192    }
1193    if (with_exps==0 || with_exps==3) p_SetCoeff(res->m[0],old_lead_coeff,r);
1194    else n_Delete(&old_lead_coeff,r->cf);
1195  }
1196  else
1197    n_Delete(&old_lead_coeff,r->cf);
1198  p_Delete(&f,r);
1199  errorreported=save_errorreported;
1200notImpl:
1201  if (res==NULL)
1202    WerrorS( feNotImplemented );
1203  if (NN!=NULL)
1204  {
1205    n_Delete(&NN,r->cf);
1206  }
1207  if (N!=NULL)
1208  {
1209    n_Delete(&N,r->cf);
1210  }
1211  return res;
1212}
1213
1214#ifdef HAVE_LIBFAC
1215matrix singclap_irrCharSeries ( ideal I, const ring r)
1216{
1217  if (idIs0(I)) return mpNew(1,1);
1218
1219  // for now there is only the possibility to handle polynomials over
1220  // Q and Fp ...
1221  matrix res=NULL;
1222  int i;
1223  Off(SW_RATIONAL);
1224  On(SW_SYMMETRIC_FF);
1225  CFList L;
1226  ListCFList LL;
1227  if (((rChar(r) == 0) || (rChar(r) > 1) )
1228  && (rPar(r)==0))
1229  {
1230    setCharacteristic( rChar(r) );
1231    for(i=0;i<IDELEMS(I);i++)
1232    {
1233      poly p=I->m[i];
1234      if (p!=NULL)
1235      {
1236        p=p_Copy(p,r);
1237        p_Cleardenom(p, r);
1238        L.append(convSingPFactoryP(p,r));
1239      }
1240    }
1241  }
1242  // and over Q(a) / Fp(a)
1243  else if (( rChar(r)==1 ) // Q(a)
1244  || (rChar(r) <-1))       // Fp(a)
1245  {
1246    if (rChar(r)==1) setCharacteristic( 0 );
1247    else               setCharacteristic( -rChar(r) );
1248    for(i=0;i<IDELEMS(I);i++)
1249    {
1250      poly p=I->m[i];
1251      if (p!=NULL)
1252      {
1253        p=p_Copy(p,r);
1254        p_Cleardenom(p, r);
1255        L.append(convSingTrPFactoryP(p,r));
1256      }
1257    }
1258  }
1259  else
1260  {
1261    WerrorS( feNotImplemented );
1262    return res;
1263  }
1264
1265  // a very bad work-around --- FIX IT in libfac
1266  // should be fixed as of 2001/6/27
1267  int tries=0;
1268  int m,n;
1269  ListIterator<CFList> LLi;
1270  loop
1271  {
1272    LL=IrrCharSeries(L);
1273    m= LL.length(); // Anzahl Zeilen
1274    n=0;
1275    for ( LLi = LL; LLi.hasItem(); LLi++ )
1276    {
1277      n = si_max(LLi.getItem().length(),n);
1278    }
1279    if ((m!=0) && (n!=0)) break;
1280    tries++;
1281    if (tries>=5) break;
1282  }
1283  if ((m==0) || (n==0))
1284  {
1285    Warn("char_series returns %d x %d matrix from %d input polys (%d)",
1286      m,n,IDELEMS(I)+1,LL.length());
1287    iiWriteMatrix((matrix)I,"I",2,r,0);
1288    m=si_max(m,1);
1289    n=si_max(n,1);
1290  }
1291  res=mpNew(m,n);
1292  CFListIterator Li;
1293  for ( m=1, LLi = LL; LLi.hasItem(); LLi++, m++ )
1294  {
1295    for (n=1, Li = LLi.getItem(); Li.hasItem(); Li++, n++)
1296    {
1297      if ( (rChar(r) == 0) || (rChar(r) > 1) )
1298        MATELEM(res,m,n)=convFactoryPSingP(Li.getItem(),r);
1299      else
1300        MATELEM(res,m,n)=convFactoryPSingTrP(Li.getItem(),r);
1301    }
1302  }
1303  Off(SW_RATIONAL);
1304  return res;
1305}
1306
1307char* singclap_neworder ( ideal I, const ring r)
1308{
1309  int i;
1310  Off(SW_RATIONAL);
1311  On(SW_SYMMETRIC_FF);
1312  CFList L;
1313  if (((rChar(r) == 0) || (rChar(r) > 1) )
1314  && (rPar(r)==0))
1315  {
1316    setCharacteristic( rChar(r) );
1317    for(i=0;i<IDELEMS(I);i++)
1318    {
1319      L.append(convSingPFactoryP(I->m[i],r));
1320    }
1321  }
1322  // and over Q(a) / Fp(a)
1323  else if (( rChar(r)==1 ) // Q(a)
1324  || (rChar(r) <-1))       // Fp(a)
1325  {
1326    if (rChar(r)==1) setCharacteristic( 0 );
1327    else               setCharacteristic( -rChar(r) );
1328    for(i=0;i<IDELEMS(I);i++)
1329    {
1330      L.append(convSingTrPFactoryP(I->m[i],r));
1331    }
1332  }
1333  else
1334  {
1335    WerrorS( feNotImplemented );
1336    return NULL;
1337  }
1338
1339  List<int> IL=neworderint(L);
1340  ListIterator<int> Li;
1341  StringSetS("");
1342  Li = IL;
1343  int offs=rPar(r);
1344  int* mark=(int*)omAlloc0((rVar(r)+offs)*sizeof(int));
1345  int cnt=rVar(r)+offs;
1346  loop
1347  {
1348    if(! Li.hasItem()) break;
1349    BOOLEAN done=TRUE;
1350    i=Li.getItem()-1;
1351    mark[i]=1;
1352    if (i<offs)
1353    {
1354      done=FALSE;
1355      //StringAppendS(r->parameter[i]);
1356    }
1357    else
1358    {
1359      StringAppendS(r->names[i-offs]);
1360    }
1361    Li++;
1362    cnt--;
1363    if(cnt==0) break;
1364    if (done) StringAppendS(",");
1365  }
1366  for(i=0;i<rVar(r)+offs;i++)
1367  {
1368    BOOLEAN done=TRUE;
1369    if(mark[i]==0)
1370    {
1371      if (i<offs)
1372      {
1373        done=FALSE;
1374        //StringAppendS(r->parameter[i]);
1375      }
1376      else
1377      {
1378        StringAppendS(r->names[i-offs]);
1379      }
1380      cnt--;
1381      if(cnt==0) break;
1382      if (done) StringAppendS(",");
1383    }
1384  }
1385  char * s=omStrDup(StringAppendS(""));
1386  if (s[strlen(s)-1]==',') s[strlen(s)-1]='\0';
1387  return s;
1388}
1389#endif /*HAVE_LIBFAC*/
1390
1391BOOLEAN singclap_isSqrFree(poly f, const ring r)
1392{
1393  BOOLEAN b=FALSE;
1394  CanonicalForm F( convSingPFactoryP( f,r ) );
1395  if((r->cf->type==n_Zp)&&(!F.isUnivariate()))
1396      goto err;
1397  b=(BOOLEAN)isSqrFree(F);
1398  Off(SW_RATIONAL);
1399  return b;
1400err:
1401  WerrorS( feNotImplemented );
1402  return 0;
1403}
1404
1405poly singclap_det( const matrix m, const ring s )
1406{
1407  int r=m->rows();
1408  if (r!=m->cols())
1409  {
1410    Werror("det of %d x %d matrix",r,m->cols());
1411    return NULL;
1412  }
1413  poly res=NULL;
1414  CFMatrix M(r,r);
1415  int i,j;
1416  for(i=r;i>0;i--)
1417  {
1418    for(j=r;j>0;j--)
1419    {
1420      M(i,j)=convSingPFactoryP(MATELEM(m,i,j),s);
1421    }
1422  }
1423  res= convFactoryPSingP( determinant(M,r),s ) ;
1424  Off(SW_RATIONAL);
1425  return res;
1426}
1427
1428int singclap_det_i( intvec * m, const ring r)
1429{
1430//  assume( r == currRing ); // Anything else is not guaranted to work!
1431   
1432  setCharacteristic( 0 ); // ?
1433  CFMatrix M(m->rows(),m->cols());
1434  int i,j;
1435  for(i=m->rows();i>0;i--)
1436  {
1437    for(j=m->cols();j>0;j--)
1438    {
1439      M(i,j)=IMATELEM(*m,i,j);
1440    }
1441  }
1442  int res= convFactoryISingI( determinant(M,m->rows()) ) ;
1443  Off(SW_RATIONAL); // ?
1444  return res;
1445}
1446
1447#ifdef HAVE_NTL
1448#if 1
1449matrix singntl_HNF(matrix  m, const ring s )
1450{
1451  int r=m->rows();
1452  if (r!=m->cols())
1453  {
1454    Werror("HNF of %d x %d matrix",r,m->cols());
1455    return NULL;
1456  }
1457   
1458  matrix res=mp_New(r,r);
1459   
1460  if (rField_is_Q(s))
1461  {
1462
1463    CFMatrix M(r,r);
1464    int i,j;
1465    for(i=r;i>0;i--)
1466    {
1467      for(j=r;j>0;j--)
1468      {
1469        M(i,j)=convSingPFactoryP(MATELEM(m,i,j),s );
1470      }
1471    }
1472    CFMatrix *MM=cf_HNF(M);
1473    for(i=r;i>0;i--)
1474    {
1475      for(j=r;j>0;j--)
1476      {
1477        MATELEM(res,i,j)=convFactoryPSingP((*MM)(i,j),s);
1478      }
1479    }
1480    delete MM;
1481  }
1482  return res;
1483}
1484
1485intvec* singntl_HNF(intvec*  m, const ring)
1486{
1487  int r=m->rows();
1488  if (r!=m->cols())
1489  {
1490    Werror("HNF of %d x %d matrix",r,m->cols());
1491    return NULL;
1492  }
1493  setCharacteristic( 0 );
1494  CFMatrix M(r,r);
1495  int i,j;
1496  for(i=r;i>0;i--)
1497  {
1498    for(j=r;j>0;j--)
1499    {
1500      M(i,j)=IMATELEM(*m,i,j);
1501    }
1502  }
1503  CFMatrix *MM=cf_HNF(M);
1504  intvec *mm=ivCopy(m);
1505  for(i=r;i>0;i--)
1506  {
1507    for(j=r;j>0;j--)
1508    {
1509      IMATELEM(*mm,i,j)=convFactoryISingI((*MM)(i,j));
1510    }
1511  }
1512  delete MM;
1513  return mm;
1514}
1515
1516matrix singntl_LLL(matrix  m, const ring s )
1517{
1518  int r=m->rows();
1519  int c=m->cols();
1520  matrix res=mp_New(r,c);
1521  if (rField_is_Q(s))
1522  {
1523    CFMatrix M(r,c);
1524    int i,j;
1525    for(i=r;i>0;i--)
1526    {
1527      for(j=c;j>0;j--)
1528      {
1529        M(i,j)=convSingPFactoryP(MATELEM(m,i,j),s);
1530      }
1531    }
1532    CFMatrix *MM=cf_LLL(M);
1533    for(i=r;i>0;i--)
1534    {
1535      for(j=c;j>0;j--)
1536      {
1537        MATELEM(res,i,j)=convFactoryPSingP((*MM)(i,j),s);
1538      }
1539    }
1540    delete MM;
1541  }
1542  return res;
1543}
1544
1545intvec* singntl_LLL(intvec*  m, const ring)
1546{
1547  int r=m->rows();
1548  int c=m->cols();
1549  setCharacteristic( 0 );
1550  CFMatrix M(r,c);
1551  int i,j;
1552  for(i=r;i>0;i--)
1553  {
1554    for(j=r;j>0;j--)
1555    {
1556      M(i,j)=IMATELEM(*m,i,j);
1557    }
1558  }
1559  CFMatrix *MM=cf_LLL(M);
1560  intvec *mm=ivCopy(m);
1561  for(i=r;i>0;i--)
1562  {
1563    for(j=c;j>0;j--)
1564    {
1565      IMATELEM(*mm,i,j)=convFactoryISingI((*MM)(i,j));
1566    }
1567  }
1568  delete MM;
1569  return mm;
1570}
1571#endif
1572#endif
1573
1574
1575#endif /* HAVE_FACTORY */
Note: See TracBrowser for help on using the repository browser.