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

spielwiese
Last change on this file since fc2746 was fc2746, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
fix the usage of iiWriteMatrix NOTE: trailing (int)0 have been treated wrongly as (ring)NULL...
  • 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
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->minideal!=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->minideal->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->minideal!=NULL)
163    {
164      //Variable X(i);
165      CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->minideal->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->minideal!=NULL)
309    {
310      CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->minideal->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->minideal!=NULL)
377    {
378      CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->minideal->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->minideal!=NULL)
533      {
534        CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->minideal->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->minideal!=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->minideal!=NULL)
736    {
737      CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->minideal->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        intvec *w=NULL;
796        if (v!=NULL) w=*v;
797        if (r->cf->extRing->minideal==NULL)
798        {
799          if(!count_Factors(res,w,j,ff,convFactoryPSingTrP( J.getItem().factor(),r ),r))
800          {
801            if (w!=NULL)
802              (*w)[j]=1;
803            res->m[j]=p_One(r);
804          }
805        }
806        else
807        {
808          if (!count_Factors(res,w,j,ff,convFactoryAPSingAP( J.getItem().factor(),r ),r))
809          {
810            if (w!=NULL)
811              (*w)[j]=1;
812            res->m[j]=p_One(r);
813          }
814        }
815      }
816    }
817    if (rField_is_Extension(r) && (!p_IsConstantPoly(ff,r)))
818    {
819      singclap_factorize_retry++;
820      if (singclap_factorize_retry<3)
821      {
822        int jj;
823        #ifdef FACTORIZE2_DEBUG
824        printf("factorize_retry\n");
825        #endif
826        intvec *ww=NULL;
827        id_Test(res,r);
828        ideal h=singclap_factorize ( ff, &ww , with_exps, r );
829        id_Test(h,r);
830        int l=(*v)->length();
831        (*v)->resize(l+ww->length());
832        for(jj=0;jj<ww->length();jj++)
833          (**v)[jj+l]=(*ww)[jj];
834        delete ww;
835        ideal hh=idInit(IDELEMS(res)+IDELEMS(h),1);
836        for(jj=IDELEMS(res)-1;jj>=0;jj--)
837        {
838          hh->m[jj]=res->m[jj];
839          res->m[jj]=NULL;
840        }
841        for(jj=IDELEMS(h)-1;jj>=0;jj--)
842        {
843          hh->m[jj+IDELEMS(res)]=h->m[jj];
844          h->m[jj]=NULL;
845        }
846        id_Delete(&res,r);
847        id_Delete(&h,r);
848        res=hh;
849        id_Test(res,r);
850        ff=NULL;
851      }
852      else
853      {
854        WarnS("problem with factorize");
855        #if 0
856        pWrite(ff);
857        idShow(res);
858        #endif
859        id_Delete(&res,r);
860        res=idInit(2,1);
861        res->m[0]=p_One(r);
862        res->m[1]=ff; ff=NULL;
863      }
864    }
865    p_Delete(&ff,r);
866    if (N!=NULL)
867    {
868      p_Mult_nn(res->m[0],N,r);
869      n_Delete(&N,r->cf);
870      N=NULL;
871    }
872    // delete constants
873    if (res!=NULL)
874    {
875      int i=IDELEMS(res)-1;
876      int j=0;
877      for(;i>=0;i--)
878      {
879        if ((res->m[i]!=NULL)
880        && (pNext(res->m[i])==NULL)
881        && (p_IsConstant(res->m[i],r)))
882        {
883          if (with_exps!=0)
884          {
885            p_Delete(&(res->m[i]),r);
886            if ((v!=NULL) && ((*v)!=NULL))
887              (**v)[i]=0;
888            j++;
889          }
890          else if (i!=0)
891          {
892            while ((v!=NULL) && ((*v)!=NULL) && ((**v)[i]>1))
893            {
894              res->m[0]=p_Mult_q(res->m[0],p_Copy(res->m[i],r),r);
895              (**v)[i]--;
896            }
897            res->m[0]=p_Mult_q(res->m[0],res->m[i],r);
898            res->m[i]=NULL;
899            if ((v!=NULL) && ((*v)!=NULL))
900              (**v)[i]=1;
901            j++;
902          }
903        }
904      }
905      if (j>0)
906      {
907        idSkipZeroes(res);
908        if ((v!=NULL) && ((*v)!=NULL))
909        {
910          intvec *w=*v;
911          int len=IDELEMS(res);
912          *v = new intvec( len );
913          for (i=0,j=0;i<si_min(w->length(),len);i++)
914          {
915            if((*w)[i]!=0)
916            {
917              (**v)[j]=(*w)[i]; j++;
918            }
919          }
920          delete w;
921        }
922      }
923      if (res->m[0]==NULL)
924      {
925        res->m[0]=p_One(r);
926      }
927    }
928  }
929  if (rField_is_Q_a(r) && (r->cf->extRing->minideal!=NULL))
930  {
931    int i=IDELEMS(res)-1;
932    int stop=1;
933    if (with_exps!=0) stop=0;
934    for(;i>=stop;i--)
935    {
936      p_Norm(res->m[i],r);
937    }
938    if (with_exps==0) p_SetCoeff(res->m[0],old_lead_coeff,r);
939    else n_Delete(&old_lead_coeff,r->cf);
940  }
941  else
942    n_Delete(&old_lead_coeff,r->cf);
943  errorreported=save_errorreported;
944notImpl:
945  if (res==NULL)
946    WerrorS( feNotImplemented );
947  if (NN!=NULL)
948  {
949    n_Delete(&NN,r->cf);
950  }
951  if (N!=NULL)
952  {
953    n_Delete(&N,r->cf);
954  }
955  if (f!=NULL) p_Delete(&f,r);
956  //PrintS("......S\n");
957  return res;
958}
959ideal singclap_sqrfree ( poly f, intvec ** v , int with_exps, const ring r)
960{
961  p_Test(f,r);
962#ifdef FACTORIZE2_DEBUG
963  printf("singclap_sqrfree, degree %d\n",pTotaldegree(f));
964#endif
965  // with_exps: 3,1 return only true factors, no exponents
966  //            2 return true factors and exponents
967  //            0 return coeff, factors and exponents
968  BOOLEAN save_errorreported=errorreported;
969
970  ideal res=NULL;
971
972  // handle factorize(0) =========================================
973  if (f==NULL)
974  {
975    res=idInit(1,1);
976    if (with_exps!=1 && with_exps!=3)
977    {
978      (*v)=new intvec(1);
979      (**v)[0]=1;
980    }
981    return res;
982  }
983  // handle factorize(mon) =========================================
984  if (pNext(f)==NULL)
985  {
986    int i=0;
987    int n=0;
988    int e;
989    for(i=rVar(r);i>0;i--) if(p_GetExp(f,i,r)!=0) n++;
990    if (with_exps==0 || with_exps==3) n++; // with coeff
991    res=idInit(si_max(n,1),1);
992    switch(with_exps)
993    {
994      case 0: // with coef & exp.
995        res->m[0]=p_NSet(n_Copy(pGetCoeff(f),r->cf),r);
996        // no break
997      case 3: // with coef & exp.
998        res->m[0]=p_NSet(n_Copy(pGetCoeff(f),r->cf),r);
999        // no break
1000      case 2: // with exp.
1001        (*v)=new intvec(si_max(1,n));
1002        (**v)[0]=1;
1003        // no break
1004      case 1: ;
1005      #ifdef TEST
1006      default: ;
1007      #endif
1008    }
1009    res->m[0]=p_NSet(n_Copy(pGetCoeff(f),r->cf),r);
1010    if (n==0)
1011    {
1012      res->m[0]=p_One(r);
1013      // (**v)[0]=1; is already done
1014    }
1015    else
1016    {
1017      for(i=rVar(r);i>0;i--)
1018      {
1019        e=p_GetExp(f,i,r);
1020        if(e!=0)
1021        {
1022          n--;
1023          poly p=p_One(r);
1024          p_SetExp(p,i,1,r);
1025          p_Setm(p,r);
1026          res->m[n]=p;
1027          if (with_exps!=1) (**v)[n]=e;
1028        }
1029      }
1030    }
1031    p_Delete(&f,r);
1032    return res;
1033  }
1034  //PrintS("S:");pWrite(f);PrintLn();
1035  // use factory/libfac in general ==============================
1036  Off(SW_RATIONAL);
1037  On(SW_SYMMETRIC_FF);
1038  #ifdef HAVE_NTL
1039  extern int prime_number;
1040  if(rField_is_Q(r)) prime_number=0;
1041  #endif
1042  CFFList L;
1043  number N=NULL;
1044  number NN=NULL;
1045  number old_lead_coeff=n_Copy(pGetCoeff(f), r->cf);
1046
1047  if (!rField_is_Zp(r) && !rField_is_Zp_a(r)) /* Q, Q(a) */
1048  {
1049    //if (f!=NULL) // already tested at start of routine
1050    number n0=n_Copy(pGetCoeff(f),r->cf);
1051    if (with_exps==0 || with_exps==3)
1052      N=n_Copy(n0,r->cf);
1053    p_Cleardenom(f, r);
1054    //after here f should not have a denominator!!
1055    //PrintS("S:");p_Write(f,r);PrintLn();
1056    NN=n_Div(n0,pGetCoeff(f),r->cf);
1057    n_Delete(&n0,r->cf);
1058    if (with_exps==0 || with_exps==3)
1059    {
1060      n_Delete(&N,r->cf);
1061      N=n_Copy(NN,r->cf);
1062    }
1063  }
1064  else if (rField_is_Zp_a(r))
1065  {
1066    //if (f!=NULL) // already tested at start of routine
1067    if (singclap_factorize_retry==0)
1068    {
1069      number n0=n_Copy(pGetCoeff(f),r->cf);
1070      if (with_exps==0 || with_exps==3)
1071        N=n_Copy(n0,r->cf);
1072      p_Norm(f,r);
1073      p_Cleardenom(f, r);
1074      NN=n_Div(n0,pGetCoeff(f),r->cf);
1075      n_Delete(&n0,r->cf);
1076      if (with_exps==0 || with_exps==3)
1077      {
1078        n_Delete(&N,r->cf);
1079        N=n_Copy(NN,r->cf);
1080      }
1081    }
1082  }
1083  if (rField_is_Q(r) || rField_is_Zp(r))
1084  {
1085    setCharacteristic( rChar(r) );
1086    CanonicalForm F( convSingPFactoryP( f,r ) );
1087    L = sqrFree( F );
1088  }
1089  else if (rField_is_Extension(r))
1090  {
1091    if (rField_is_Q_a (r)) setCharacteristic (0);
1092    else                   setCharacteristic( rChar(r) );
1093    if (r->cf->extRing->minideal!=NULL)
1094    {
1095      CanonicalForm mipo=convSingPFactoryP(r->cf->extRing->minideal->m[0],
1096                                           r->cf->extRing);
1097      Variable a=rootOf(mipo);
1098      CanonicalForm F( convSingAPFactoryAP( f, a, r ) );
1099      L= sqrFree (F);
1100    }
1101    else
1102    {
1103      CanonicalForm F( convSingTrPFactoryP( f,r ) );
1104      L = sqrFree( F );
1105    }
1106  }
1107  #if 0
1108  else if (rField_is_GF())
1109  {
1110    int c=rChar(r);
1111    setCharacteristic( c, primepower(c) );
1112    CanonicalForm F( convSingGFFactoryGF( f ) );
1113    if (F.isUnivariate())
1114    {
1115      L = factorize( F );
1116    }
1117    else
1118    {
1119      goto notImpl;
1120    }
1121  }
1122  #endif
1123  else
1124  {
1125    goto notImpl;
1126  }
1127  {
1128    // convert into ideal
1129    int n = L.length();
1130    if (n==0) n=1;
1131    CFFListIterator J=L;
1132    int j=0;
1133    if (with_exps!=1)
1134    {
1135      if ((with_exps==2)&&(n>1))
1136      {
1137        n--;
1138        J++;
1139      }
1140      *v = new intvec( n );
1141    }
1142    else if (L.getFirst().factor().inCoeffDomain() && with_exps!=3)
1143    {
1144      n--;
1145      J++;
1146    }
1147    res = idInit( n ,1);
1148    for ( ; J.hasItem(); J++, j++ )
1149    {
1150      if (with_exps!=1 && with_exps!=3) (**v)[j] = J.getItem().exp();
1151      if (rField_is_Zp(r) || rField_is_Q(r))
1152        res->m[j] = convFactoryPSingP( J.getItem().factor(),r );
1153      else if (rField_is_Extension(r))     /* Q(a), Fp(a) */
1154      {
1155        if (r->cf->extRing->minideal==NULL)
1156          res->m[j]=convFactoryPSingTrP( J.getItem().factor(),r );
1157        else
1158          res->m[j]=convFactoryAPSingAP( J.getItem().factor(),r );
1159      }
1160    }
1161    if (res->m[0]==NULL)
1162    {
1163      res->m[0]=p_One(r);
1164    }
1165    if (N!=NULL)
1166    {
1167      p_Mult_nn(res->m[0],N,r);
1168      n_Delete(&N,r->cf);
1169      N=NULL;
1170    }
1171  }
1172  if (rField_is_Q_a(r) && (r->cf->extRing->minideal!=NULL))
1173  {
1174    int i=IDELEMS(res)-1;
1175    int stop=1;
1176    if (with_exps!=0 || with_exps==3) stop=0;
1177    for(;i>=stop;i--)
1178    {
1179      p_Norm(res->m[i],r);
1180    }
1181    if (with_exps==0 || with_exps==3) p_SetCoeff(res->m[0],old_lead_coeff,r);
1182    else n_Delete(&old_lead_coeff,r->cf);
1183  }
1184  else
1185    n_Delete(&old_lead_coeff,r->cf);
1186  p_Delete(&f,r);
1187  errorreported=save_errorreported;
1188notImpl:
1189  if (res==NULL)
1190    WerrorS( feNotImplemented );
1191  if (NN!=NULL)
1192  {
1193    n_Delete(&NN,r->cf);
1194  }
1195  if (N!=NULL)
1196  {
1197    n_Delete(&N,r->cf);
1198  }
1199  return res;
1200}
1201
1202#ifdef HAVE_LIBFAC
1203matrix singclap_irrCharSeries ( ideal I, const ring r)
1204{
1205  if (idIs0(I)) return mpNew(1,1);
1206
1207  // for now there is only the possibility to handle polynomials over
1208  // Q and Fp ...
1209  matrix res=NULL;
1210  int i;
1211  Off(SW_RATIONAL);
1212  On(SW_SYMMETRIC_FF);
1213  CFList L;
1214  ListCFList LL;
1215  if (((rChar(r) == 0) || (rChar(r) > 1) )
1216  && (rPar(r)==0))
1217  {
1218    setCharacteristic( rChar(r) );
1219    for(i=0;i<IDELEMS(I);i++)
1220    {
1221      poly p=I->m[i];
1222      if (p!=NULL)
1223      {
1224        p=p_Copy(p,r);
1225        p_Cleardenom(p, r);
1226        L.append(convSingPFactoryP(p,r));
1227      }
1228    }
1229  }
1230  // and over Q(a) / Fp(a)
1231  else if (( rChar(r)==1 ) // Q(a)
1232  || (rChar(r) <-1))       // Fp(a)
1233  {
1234    if (rChar(r)==1) setCharacteristic( 0 );
1235    else               setCharacteristic( -rChar(r) );
1236    for(i=0;i<IDELEMS(I);i++)
1237    {
1238      poly p=I->m[i];
1239      if (p!=NULL)
1240      {
1241        p=p_Copy(p,r);
1242        p_Cleardenom(p, r);
1243        L.append(convSingTrPFactoryP(p,r));
1244      }
1245    }
1246  }
1247  else
1248  {
1249    WerrorS( feNotImplemented );
1250    return res;
1251  }
1252
1253  // a very bad work-around --- FIX IT in libfac
1254  // should be fixed as of 2001/6/27
1255  int tries=0;
1256  int m,n;
1257  ListIterator<CFList> LLi;
1258  loop
1259  {
1260    LL=IrrCharSeries(L);
1261    m= LL.length(); // Anzahl Zeilen
1262    n=0;
1263    for ( LLi = LL; LLi.hasItem(); LLi++ )
1264    {
1265      n = si_max(LLi.getItem().length(),n);
1266    }
1267    if ((m!=0) && (n!=0)) break;
1268    tries++;
1269    if (tries>=5) break;
1270  }
1271  if ((m==0) || (n==0))
1272  {
1273    Warn("char_series returns %d x %d matrix from %d input polys (%d)",
1274      m,n,IDELEMS(I)+1,LL.length());
1275    iiWriteMatrix((matrix)I,"I",2,r,0);
1276    m=si_max(m,1);
1277    n=si_max(n,1);
1278  }
1279  res=mpNew(m,n);
1280  CFListIterator Li;
1281  for ( m=1, LLi = LL; LLi.hasItem(); LLi++, m++ )
1282  {
1283    for (n=1, Li = LLi.getItem(); Li.hasItem(); Li++, n++)
1284    {
1285      if ( (rChar(r) == 0) || (rChar(r) > 1) )
1286        MATELEM(res,m,n)=convFactoryPSingP(Li.getItem(),r);
1287      else
1288        MATELEM(res,m,n)=convFactoryPSingTrP(Li.getItem(),r);
1289    }
1290  }
1291  Off(SW_RATIONAL);
1292  return res;
1293}
1294
1295char* singclap_neworder ( ideal I, const ring r)
1296{
1297  int i;
1298  Off(SW_RATIONAL);
1299  On(SW_SYMMETRIC_FF);
1300  CFList L;
1301  if (((rChar(r) == 0) || (rChar(r) > 1) )
1302  && (rPar(r)==0))
1303  {
1304    setCharacteristic( rChar(r) );
1305    for(i=0;i<IDELEMS(I);i++)
1306    {
1307      L.append(convSingPFactoryP(I->m[i],r));
1308    }
1309  }
1310  // and over Q(a) / Fp(a)
1311  else if (( rChar(r)==1 ) // Q(a)
1312  || (rChar(r) <-1))       // Fp(a)
1313  {
1314    if (rChar(r)==1) setCharacteristic( 0 );
1315    else               setCharacteristic( -rChar(r) );
1316    for(i=0;i<IDELEMS(I);i++)
1317    {
1318      L.append(convSingTrPFactoryP(I->m[i],r));
1319    }
1320  }
1321  else
1322  {
1323    WerrorS( feNotImplemented );
1324    return NULL;
1325  }
1326
1327  List<int> IL=neworderint(L);
1328  ListIterator<int> Li;
1329  StringSetS("");
1330  Li = IL;
1331  int offs=rPar(r);
1332  int* mark=(int*)omAlloc0((rVar(r)+offs)*sizeof(int));
1333  int cnt=rVar(r)+offs;
1334  loop
1335  {
1336    if(! Li.hasItem()) break;
1337    BOOLEAN done=TRUE;
1338    i=Li.getItem()-1;
1339    mark[i]=1;
1340    if (i<offs)
1341    {
1342      done=FALSE;
1343      //StringAppendS(r->parameter[i]);
1344    }
1345    else
1346    {
1347      StringAppendS(r->names[i-offs]);
1348    }
1349    Li++;
1350    cnt--;
1351    if(cnt==0) break;
1352    if (done) StringAppendS(",");
1353  }
1354  for(i=0;i<rVar(r)+offs;i++)
1355  {
1356    BOOLEAN done=TRUE;
1357    if(mark[i]==0)
1358    {
1359      if (i<offs)
1360      {
1361        done=FALSE;
1362        //StringAppendS(r->parameter[i]);
1363      }
1364      else
1365      {
1366        StringAppendS(r->names[i-offs]);
1367      }
1368      cnt--;
1369      if(cnt==0) break;
1370      if (done) StringAppendS(",");
1371    }
1372  }
1373  char * s=omStrDup(StringAppendS(""));
1374  if (s[strlen(s)-1]==',') s[strlen(s)-1]='\0';
1375  return s;
1376}
1377#endif /*HAVE_LIBFAC*/
1378
1379BOOLEAN singclap_isSqrFree(poly f, const ring r)
1380{
1381  BOOLEAN b=FALSE;
1382  CanonicalForm F( convSingPFactoryP( f,r ) );
1383  if((r->cf->type==n_Zp)&&(!F.isUnivariate()))
1384      goto err;
1385  b=(BOOLEAN)isSqrFree(F);
1386  Off(SW_RATIONAL);
1387  return b;
1388err:
1389  WerrorS( feNotImplemented );
1390  return 0;
1391}
1392
1393poly singclap_det( const matrix m, const ring s )
1394{
1395  int r=m->rows();
1396  if (r!=m->cols())
1397  {
1398    Werror("det of %d x %d matrix",r,m->cols());
1399    return NULL;
1400  }
1401  poly res=NULL;
1402  CFMatrix M(r,r);
1403  int i,j;
1404  for(i=r;i>0;i--)
1405  {
1406    for(j=r;j>0;j--)
1407    {
1408      M(i,j)=convSingPFactoryP(MATELEM(m,i,j),s);
1409    }
1410  }
1411  res= convFactoryPSingP( determinant(M,r),s ) ;
1412  Off(SW_RATIONAL);
1413  return res;
1414}
1415
1416int singclap_det_i( intvec * m, const ring r)
1417{
1418//  assume( r == currRing ); // Anything else is not guaranted to work!
1419   
1420  setCharacteristic( 0 ); // ?
1421  CFMatrix M(m->rows(),m->cols());
1422  int i,j;
1423  for(i=m->rows();i>0;i--)
1424  {
1425    for(j=m->cols();j>0;j--)
1426    {
1427      M(i,j)=IMATELEM(*m,i,j);
1428    }
1429  }
1430  int res= convFactoryISingI( determinant(M,m->rows()) ) ;
1431  Off(SW_RATIONAL); // ?
1432  return res;
1433}
1434
1435#ifdef HAVE_NTL
1436#if 1
1437matrix singntl_HNF(matrix  m, const ring s )
1438{
1439  int r=m->rows();
1440  if (r!=m->cols())
1441  {
1442    Werror("HNF of %d x %d matrix",r,m->cols());
1443    return NULL;
1444  }
1445   
1446  matrix res=mp_New(r,r);
1447   
1448  if (rField_is_Q(s))
1449  {
1450
1451    CFMatrix M(r,r);
1452    int i,j;
1453    for(i=r;i>0;i--)
1454    {
1455      for(j=r;j>0;j--)
1456      {
1457        M(i,j)=convSingPFactoryP(MATELEM(m,i,j),s );
1458      }
1459    }
1460    CFMatrix *MM=cf_HNF(M);
1461    for(i=r;i>0;i--)
1462    {
1463      for(j=r;j>0;j--)
1464      {
1465        MATELEM(res,i,j)=convFactoryPSingP((*MM)(i,j),s);
1466      }
1467    }
1468    delete MM;
1469  }
1470  return res;
1471}
1472
1473intvec* singntl_HNF(intvec*  m, const ring)
1474{
1475  int r=m->rows();
1476  if (r!=m->cols())
1477  {
1478    Werror("HNF of %d x %d matrix",r,m->cols());
1479    return NULL;
1480  }
1481  setCharacteristic( 0 );
1482  CFMatrix M(r,r);
1483  int i,j;
1484  for(i=r;i>0;i--)
1485  {
1486    for(j=r;j>0;j--)
1487    {
1488      M(i,j)=IMATELEM(*m,i,j);
1489    }
1490  }
1491  CFMatrix *MM=cf_HNF(M);
1492  intvec *mm=ivCopy(m);
1493  for(i=r;i>0;i--)
1494  {
1495    for(j=r;j>0;j--)
1496    {
1497      IMATELEM(*mm,i,j)=convFactoryISingI((*MM)(i,j));
1498    }
1499  }
1500  delete MM;
1501  return mm;
1502}
1503
1504matrix singntl_LLL(matrix  m, const ring s )
1505{
1506  int r=m->rows();
1507  int c=m->cols();
1508  matrix res=mp_New(r,c);
1509  if (rField_is_Q(s))
1510  {
1511    CFMatrix M(r,c);
1512    int i,j;
1513    for(i=r;i>0;i--)
1514    {
1515      for(j=c;j>0;j--)
1516      {
1517        M(i,j)=convSingPFactoryP(MATELEM(m,i,j),s);
1518      }
1519    }
1520    CFMatrix *MM=cf_LLL(M);
1521    for(i=r;i>0;i--)
1522    {
1523      for(j=c;j>0;j--)
1524      {
1525        MATELEM(res,i,j)=convFactoryPSingP((*MM)(i,j),s);
1526      }
1527    }
1528    delete MM;
1529  }
1530  return res;
1531}
1532
1533intvec* singntl_LLL(intvec*  m, const ring)
1534{
1535  int r=m->rows();
1536  int c=m->cols();
1537  setCharacteristic( 0 );
1538  CFMatrix M(r,c);
1539  int i,j;
1540  for(i=r;i>0;i--)
1541  {
1542    for(j=r;j>0;j--)
1543    {
1544      M(i,j)=IMATELEM(*m,i,j);
1545    }
1546  }
1547  CFMatrix *MM=cf_LLL(M);
1548  intvec *mm=ivCopy(m);
1549  for(i=r;i>0;i--)
1550  {
1551    for(j=c;j>0;j--)
1552    {
1553      IMATELEM(*mm,i,j)=convFactoryISingI((*MM)(i,j));
1554    }
1555  }
1556  delete MM;
1557  return mm;
1558}
1559#endif
1560#endif
1561
1562
1563#endif /* HAVE_FACTORY */
Note: See TracBrowser for help on using the repository browser.