source: git/kernel/clapsing.cc @ cafd4ff

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