source: git/kernel/clapsing.cc @ 02a069e

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