source: git/kernel/clapsing.cc @ bafaec0

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