source: git/kernel/clapsing.cc @ 278418

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