source: git/kernel/clapsing.cc @ cf74cd6

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