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

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