source: git/kernel/clapsing.cc @ 59a7ca1

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