source: git/kernel/clapsing.cc @ 147167f

fieker-DuValspielwiese
Last change on this file since 147167f was b07ba77, checked in by Hans Schönemann <hannes@…>, 14 years ago
clean up structs.h, part ... git-svn-id: file:///usr/local/Singular/svn/trunk@12437 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{
726  pTest(f);
727#ifdef FACTORIZE2_DEBUG
728  printf("singclap_factorize, degree %d\n",pTotaldegree(f));
729#endif
730  // with_exps: 3,1 return only true factors, no exponents
731  //            2 return true factors and exponents
732  //            0 return coeff, factors and exponents
733  BOOLEAN save_errorreported=errorreported;
734
735  ideal res=NULL;
736
737  // handle factorize(0) =========================================
738  if (f==NULL)
739  {
740    res=idInit(1,1);
741    if (with_exps!=1)
742    {
743      (*v)=new intvec(1);
744      (**v)[0]=1;
745    }
746    return res;
747  }
748  // handle factorize(mon) =========================================
749  if (pNext(f)==NULL)
750  {
751    int i=0;
752    int n=0;
753    int e;
754    for(i=pVariables;i>0;i--) if(pGetExp(f,i)!=0) n++;
755    if (with_exps==0) n++; // with coeff
756    res=idInit(si_max(n,1),1);
757    switch(with_exps)
758    {
759      case 0: // with coef & exp.
760        res->m[0]=pNSet(nCopy(pGetCoeff(f)));
761        // no break
762      case 2: // with exp.
763        (*v)=new intvec(si_max(1,n));
764        (**v)[0]=1;
765        // no break
766      case 1: ;
767      #ifdef TEST
768      default: ;
769      #endif
770    }
771    if (n==0)
772    {
773      res->m[0]=pOne();
774      // (**v)[0]=1; is already done
775      return res;
776    }
777    for(i=pVariables;i>0;i--)
778    {
779      e=pGetExp(f,i);
780      if(e!=0)
781      {
782        n--;
783        poly p=pOne();
784        pSetExp(p,i,1);
785        pSetm(p);
786        res->m[n]=p;
787        if (with_exps!=1) (**v)[n]=e;
788      }
789    }
790    return res;
791  }
792  //PrintS("S:");pWrite(f);PrintLn();
793  // use factory/libfac in general ==============================
794  Off(SW_RATIONAL);
795  On(SW_SYMMETRIC_FF);
796  #ifdef HAVE_NTL
797  extern int prime_number;
798  if(rField_is_Q()) prime_number=0;
799  #endif
800  CFFList L;
801  number N=NULL;
802  number NN=NULL;
803  number old_lead_coeff=nCopy(pGetCoeff(f));
804
805  if (!rField_is_Zp() && !rField_is_Zp_a()) /* Q, Q(a) */
806  {
807    //if (f!=NULL) // already tested at start of routine
808    {
809      number n0=nCopy(pGetCoeff(f));
810      if (with_exps==0)
811        N=nCopy(n0);
812      pCleardenom(f);
813      NN=nDiv(n0,pGetCoeff(f));
814      nDelete(&n0);
815      if (with_exps==0)
816      {
817        nDelete(&N);
818        N=nCopy(NN);
819      }
820    }
821  }
822  else if (rField_is_Zp_a())
823  {
824    //if (f!=NULL) // already tested at start of routine
825    if (singclap_factorize_retry==0)
826    {
827      number n0=nCopy(pGetCoeff(f));
828      if (with_exps==0)
829        N=nCopy(n0);
830      pNorm(f);
831      pCleardenom(f);
832      NN=nDiv(n0,pGetCoeff(f));
833      nDelete(&n0);
834      if (with_exps==0)
835      {
836        nDelete(&N);
837        N=nCopy(NN);
838      }
839    }
840  }
841  if (rField_is_Q() || rField_is_Zp())
842  {
843    setCharacteristic( nGetChar() );
844    CanonicalForm F( convSingPFactoryP( f ) );
845    if (rField_is_Q())
846    {
847      L = factorize( F );
848    }
849    else /* Fp */
850    {
851      do
852      {
853        libfac_interruptflag=0;
854        L = Factorize( F );
855      }
856      while ((libfac_interruptflag!=0) ||(L.isEmpty()));
857    }
858  }
859  #if 0
860  else if (rField_is_GF())
861  {
862    int c=rChar(currRing);
863    setCharacteristic( c, primepower(c) );
864    CanonicalForm F( convSingGFFactoryGF( f ) );
865    if (F.isUnivariate())
866    {
867      L = factorize( F );
868    }
869    else
870    {
871      goto notImpl;
872    }
873  }
874  #endif
875  // and over Q(a) / Fp(a)
876  else if (rField_is_Extension())
877  {
878    if (rField_is_Q_a()) setCharacteristic( 0 );
879    else                 setCharacteristic( -nGetChar() );
880    if (currRing->minpoly!=NULL)
881    {
882      CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
883                   currRing->algring);
884      Variable a=rootOf(mipo);
885      CanonicalForm F( convSingAPFactoryAP( f,a,currRing ) );
886      if (rField_is_Zp_a() && F.isUnivariate())
887      {
888        L = factorize( F, a );
889      }
890      else
891      {
892        CanonicalForm G( convSingTrPFactoryP( f ) );
893        //  over Q(a) / multivariate over Fp(a)
894        do
895        {
896          libfac_interruptflag=0;
897          L=Factorize2(G, mipo);
898        }
899        while ((libfac_interruptflag!=0) ||(L.isEmpty()));
900        #ifdef FACTORIZE2_DEBUG
901        printf("while okay\n");
902        #endif
903        libfac_interruptflag=0;
904      }
905    }
906    else
907    {
908      CanonicalForm F( convSingTrPFactoryP( f ) );
909      if (rField_is_Q_a())
910      {
911        L = factorize( F );
912      }
913      else /* Fp(a) */
914      {
915        L = Factorize( F );
916      }
917    }
918  }
919  else
920  {
921    goto notImpl;
922  }
923  {
924    poly ff=pCopy(f); // a copy for the retry stuff
925    // the first factor should be a constant
926    if ( ! L.getFirst().factor().inCoeffDomain() )
927      L.insert(CFFactor(1,1));
928    // convert into ideal
929    int n = L.length();
930    if (n==0) n=1;
931    CFFListIterator J=L;
932    int j=0;
933    if (with_exps!=1)
934    {
935      if ((with_exps==2)&&(n>1))
936      {
937        n--;
938        J++;
939      }
940      *v = new intvec( n );
941    }
942    res = idInit( n ,1);
943    for ( ; J.hasItem(); J++, j++ )
944    {
945      poly p;
946      if (with_exps!=1) (**v)[j] = J.getItem().exp();
947      if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
948      {
949        //count_Factors(res,*v,f, j, convFactoryPSingP( J.getItem().factor() );
950        res->m[j] = convFactoryPSingP( J.getItem().factor() );
951      }
952      #if 0
953      else if (rField_is_GF())
954        res->m[j] = convFactoryGFSingGF( J.getItem().factor() );
955      #endif
956      else if (rField_is_Extension())     /* Q(a), Fp(a) */
957      {
958        intvec *w=NULL;
959        if (v!=NULL) w=*v;
960        if (currRing->minpoly==NULL)
961        {
962          if(!count_Factors(res,w,j,ff,convFactoryPSingTrP( J.getItem().factor() )))
963          {
964            if (w!=NULL)
965              (*w)[j]=1;
966            res->m[j]=pOne();
967          }
968        }
969        else
970        {
971          if (!count_Factors(res,w,j,ff,convFactoryAPSingAP( J.getItem().factor(),currRing )))
972          {
973            if (w!=NULL)
974              (*w)[j]=1;
975            res->m[j]=pOne();
976          }
977        }
978      }
979    }
980    if (rField_is_Extension() && (!pIsConstantPoly(ff)))
981    {
982      singclap_factorize_retry++;
983      if (singclap_factorize_retry<3)
984      {
985        int jj;
986        #ifdef FACTORIZE2_DEBUG
987        printf("factorize_retry\n");
988        #endif
989        intvec *ww=NULL;
990        idTest(res);
991        ideal h=singclap_factorize ( ff, &ww , with_exps);
992        idTest(h);
993        int l=(*v)->length();
994        (*v)->resize(l+ww->length());
995        for(jj=0;jj<ww->length();jj++)
996          (**v)[jj+l]=(*ww)[jj];
997        delete ww;
998        ideal hh=idInit(IDELEMS(res)+IDELEMS(h),1);
999        for(jj=IDELEMS(res)-1;jj>=0;jj--)
1000        {
1001          hh->m[jj]=res->m[jj];
1002          res->m[jj]=NULL;
1003        }
1004        for(jj=IDELEMS(h)-1;jj>=0;jj--)
1005        {
1006          hh->m[jj+IDELEMS(res)]=h->m[jj];
1007          h->m[jj]=NULL;
1008        }
1009        idDelete(&res);
1010        idDelete(&h);
1011        res=hh;
1012        idTest(res);
1013        ff=NULL;
1014      }
1015      else
1016      {
1017        WarnS("problem with factorize");
1018        #if 0
1019        pWrite(ff);
1020        idShow(res);
1021        #endif
1022        idDelete(&res);
1023        res=idInit(2,1);
1024        res->m[0]=pOne();
1025        res->m[1]=ff; ff=NULL;
1026      }
1027    }
1028    pDelete(&ff);
1029    if (N!=NULL)
1030    {
1031      pMult_nn(res->m[0],N);
1032      nDelete(&N);
1033      N=NULL;
1034    }
1035    // delete constants
1036    if (res!=NULL)
1037    {
1038      int i=IDELEMS(res)-1;
1039      int j=0;
1040      for(;i>=0;i--)
1041      {
1042        if ((res->m[i]!=NULL)
1043        && (pNext(res->m[i])==NULL)
1044        && (pIsConstant(res->m[i])))
1045        {
1046          if (with_exps!=0)
1047          {
1048            pDelete(&(res->m[i]));
1049            if ((v!=NULL) && ((*v)!=NULL))
1050              (**v)[i]=0;
1051            j++;
1052          }
1053          else if (i!=0)
1054          {
1055            while ((v!=NULL) && ((*v)!=NULL) && ((**v)[i]>1))
1056            {
1057              res->m[0]=pMult(res->m[0],pCopy(res->m[i]));
1058              (**v)[i]--;
1059            }
1060            res->m[0]=pMult(res->m[0],res->m[i]);
1061            res->m[i]=NULL;
1062            if ((v!=NULL) && ((*v)!=NULL))
1063              (**v)[i]=1;
1064            j++;
1065          }
1066        }
1067      }
1068      if (j>0)
1069      {
1070        idSkipZeroes(res);
1071        if ((v!=NULL) && ((*v)!=NULL))
1072        {
1073          intvec *w=*v;
1074          int len=IDELEMS(res);
1075          *v = new intvec( len );
1076          for (i=0,j=0;i<si_min(w->length(),len);i++)
1077          {
1078            if((*w)[i]!=0)
1079            {
1080              (**v)[j]=(*w)[i]; j++;
1081            }
1082          }
1083          delete w;
1084        }
1085      }
1086      if (res->m[0]==NULL)
1087      {
1088        res->m[0]=pOne();
1089      }
1090    }
1091  }
1092  if (rField_is_Q_a() && (currRing->minpoly!=NULL))
1093  {
1094    int i=IDELEMS(res)-1;
1095    int stop=1;
1096    if (with_exps!=0) stop=0;
1097    for(;i>=stop;i--)
1098    {
1099      pNorm(res->m[i]);
1100    }
1101    if (with_exps==0) pSetCoeff(res->m[0],old_lead_coeff);
1102    else nDelete(&old_lead_coeff);
1103  }
1104  else
1105    nDelete(&old_lead_coeff);
1106  errorreported=save_errorreported;
1107notImpl:
1108  if (res==NULL)
1109    WerrorS( feNotImplemented );
1110  if (NN!=NULL)
1111  {
1112    nDelete(&NN);
1113  }
1114  if (N!=NULL)
1115  {
1116    nDelete(&N);
1117  }
1118  //if (f!=NULL) pDelete(&f);
1119  //PrintS("......S\n");
1120  return res;
1121}
1122ideal singclap_sqrfree ( poly f)
1123{
1124  pTest(f);
1125#ifdef FACTORIZE2_DEBUG
1126  printf("singclap_sqrfree, degree %d\n",pTotaldegree(f));
1127#endif
1128  // with_exps: 3,1 return only true factors, no exponents
1129  //            2 return true factors and exponents
1130  //            0 return coeff, factors and exponents
1131  BOOLEAN save_errorreported=errorreported;
1132
1133  ideal res=NULL;
1134
1135  // handle factorize(0) =========================================
1136  if (f==NULL)
1137  {
1138    res=idInit(1,1);
1139    return res;
1140  }
1141  // handle factorize(mon) =========================================
1142  if (pNext(f)==NULL)
1143  {
1144    int i=0;
1145    int n=0;
1146    int e;
1147    for(i=pVariables;i>0;i--) if(pGetExp(f,i)!=0) n++;
1148    n++; // with coeff
1149    res=idInit(si_max(n,1),1);
1150    res->m[0]=pNSet(nCopy(pGetCoeff(f)));
1151    if (n==0)
1152    {
1153      res->m[0]=pOne();
1154      // (**v)[0]=1; is already done
1155      return res;
1156    }
1157    for(i=pVariables;i>0;i--)
1158    {
1159      e=pGetExp(f,i);
1160      if(e!=0)
1161      {
1162        n--;
1163        poly p=pOne();
1164        pSetExp(p,i,1);
1165        pSetm(p);
1166        res->m[n]=p;
1167      }
1168    }
1169    return res;
1170  }
1171  //PrintS("S:");pWrite(f);PrintLn();
1172  // use factory/libfac in general ==============================
1173  Off(SW_RATIONAL);
1174  On(SW_SYMMETRIC_FF);
1175  #ifdef HAVE_NTL
1176  extern int prime_number;
1177  if(rField_is_Q()) prime_number=0;
1178  #endif
1179  CFFList L;
1180
1181  if (!rField_is_Zp() && !rField_is_Zp_a()) /* Q, Q(a) */
1182  {
1183    //if (f!=NULL) // already tested at start of routine
1184    {
1185      pCleardenom(f);
1186    }
1187  }
1188  else if (rField_is_Zp_a())
1189  {
1190    //if (f!=NULL) // already tested at start of routine
1191    if (singclap_factorize_retry==0)
1192    {
1193      pNorm(f);
1194      pCleardenom(f);
1195    }
1196  }
1197  if (rField_is_Q() || rField_is_Zp())
1198  {
1199    setCharacteristic( nGetChar() );
1200    CanonicalForm F( convSingPFactoryP( f ) );
1201    L = sqrFree( F );
1202  }
1203  #if 0
1204  else if (rField_is_GF())
1205  {
1206    int c=rChar(currRing);
1207    setCharacteristic( c, primepower(c) );
1208    CanonicalForm F( convSingGFFactoryGF( f ) );
1209    if (F.isUnivariate())
1210    {
1211      L = factorize( F );
1212    }
1213    else
1214    {
1215      goto notImpl;
1216    }
1217  }
1218  #endif
1219  // and over Q(a) / Fp(a)
1220  else if (rField_is_Extension())
1221  {
1222    if (rField_is_Q_a()) setCharacteristic( 0 );
1223    else                 setCharacteristic( -nGetChar() );
1224    if (currRing->minpoly!=NULL)
1225    {
1226      CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
1227                    currRing->algring);
1228      Variable a=rootOf(mipo);
1229      CanonicalForm F( convSingAPFactoryAP( f,a,currRing ) );
1230      CFFList SqrFreeMV( const CanonicalForm & f , const CanonicalForm & mipo=0) ;
1231
1232      L = SqrFreeMV( F,mipo );
1233      //WarnS("L = sqrFree( F,mipo );");
1234      //L = sqrFree( F );
1235    }
1236    else
1237    {
1238      CanonicalForm F( convSingTrPFactoryP( f ) );
1239      L = sqrFree( F );
1240    }
1241  }
1242  else
1243  {
1244    goto notImpl;
1245  }
1246  {
1247    poly ff=pCopy(f); // a copy for the retry stuff
1248    // convert into ideal
1249    int n = L.length();
1250    if (n==0) n=1;
1251    CFFListIterator J=L;
1252    int j=0;
1253    res = idInit( n ,1);
1254    for ( ; J.hasItem(); J++, j++ )
1255    {
1256      poly p;
1257      if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
1258        //count_Factors(res,*v,f, j, convFactoryPSingP( J.getItem().factor() );
1259        res->m[j] = convFactoryPSingP( J.getItem().factor() );
1260      #if 0
1261      else if (rField_is_GF())
1262        res->m[j] = convFactoryGFSingGF( J.getItem().factor() );
1263      #endif
1264      else if (rField_is_Extension())     /* Q(a), Fp(a) */
1265      {
1266        if (currRing->minpoly==NULL)
1267          res->m[j]=convFactoryPSingTrP( J.getItem().factor() );
1268        else
1269          res->m[j]=convFactoryAPSingAP( J.getItem().factor(),currRing );
1270      }
1271    }
1272    if (res->m[0]==NULL)
1273    {
1274      res->m[0]=pOne();
1275    }
1276  }
1277  errorreported=save_errorreported;
1278notImpl:
1279  if (res==NULL)
1280    WerrorS( feNotImplemented );
1281  return res;
1282}
1283matrix singclap_irrCharSeries ( ideal I)
1284{
1285  if (idIs0(I)) return mpNew(1,1);
1286
1287  // for now there is only the possibility to handle polynomials over
1288  // Q and Fp ...
1289  matrix res=NULL;
1290  int i;
1291  Off(SW_RATIONAL);
1292  On(SW_SYMMETRIC_FF);
1293  CFList L;
1294  ListCFList LL;
1295  if (((nGetChar() == 0) || (nGetChar() > 1) )
1296  && (currRing->parameter==NULL))
1297  {
1298    setCharacteristic( nGetChar() );
1299    for(i=0;i<IDELEMS(I);i++)
1300    {
1301      poly p=I->m[i];
1302      if (p!=NULL)
1303      {
1304        p=pCopy(p);
1305        pCleardenom(p);
1306        L.append(convSingPFactoryP(p));
1307      }
1308    }
1309  }
1310  // and over Q(a) / Fp(a)
1311  else if (( nGetChar()==1 ) /* Q(a) */
1312  || (nGetChar() <-1))       /* Fp(a) */
1313  {
1314    if (nGetChar()==1) setCharacteristic( 0 );
1315    else               setCharacteristic( -nGetChar() );
1316    for(i=0;i<IDELEMS(I);i++)
1317    {
1318      poly p=I->m[i];
1319      if (p!=NULL)
1320      {
1321        p=pCopy(p);
1322        pCleardenom(p);
1323        L.append(convSingTrPFactoryP(p));
1324      }
1325    }
1326  }
1327  else
1328  {
1329    WerrorS( feNotImplemented );
1330    return res;
1331  }
1332
1333  // a very bad work-around --- FIX IT in libfac
1334  // should be fixed as of 2001/6/27
1335  int tries=0;
1336  int m,n;
1337  ListIterator<CFList> LLi;
1338  loop
1339  {
1340    LL=IrrCharSeries(L);
1341    m= LL.length(); // Anzahl Zeilen
1342    n=0;
1343    for ( LLi = LL; LLi.hasItem(); LLi++ )
1344    {
1345      n = si_max(LLi.getItem().length(),n);
1346    }
1347    if ((m!=0) && (n!=0)) break;
1348    tries++;
1349    if (tries>=5) break;
1350  }
1351  if ((m==0) || (n==0))
1352  {
1353    Warn("char_series returns %d x %d matrix from %d input polys (%d)",
1354      m,n,IDELEMS(I)+1,LL.length());
1355    iiWriteMatrix((matrix)I,"I",2,0);
1356    m=si_max(m,1);
1357    n=si_max(n,1);
1358  }
1359  res=mpNew(m,n);
1360  CFListIterator Li;
1361  for ( m=1, LLi = LL; LLi.hasItem(); LLi++, m++ )
1362  {
1363    for (n=1, Li = LLi.getItem(); Li.hasItem(); Li++, n++)
1364    {
1365      if ( (nGetChar() == 0) || (nGetChar() > 1) )
1366        MATELEM(res,m,n)=convFactoryPSingP(Li.getItem());
1367      else
1368        MATELEM(res,m,n)=convFactoryPSingTrP(Li.getItem());
1369    }
1370  }
1371  Off(SW_RATIONAL);
1372  return res;
1373}
1374
1375char* singclap_neworder ( ideal I)
1376{
1377  int i;
1378  Off(SW_RATIONAL);
1379  On(SW_SYMMETRIC_FF);
1380  CFList L;
1381  if (((nGetChar() == 0) || (nGetChar() > 1) )
1382  && (currRing->parameter==NULL))
1383  {
1384    setCharacteristic( nGetChar() );
1385    for(i=0;i<IDELEMS(I);i++)
1386    {
1387      L.append(convSingPFactoryP(I->m[i]));
1388    }
1389  }
1390  // and over Q(a) / Fp(a)
1391  else if (( nGetChar()==1 ) /* Q(a) */
1392  || (nGetChar() <-1))       /* Fp(a) */
1393  {
1394    if (nGetChar()==1) setCharacteristic( 0 );
1395    else               setCharacteristic( -nGetChar() );
1396    for(i=0;i<IDELEMS(I);i++)
1397    {
1398      L.append(convSingTrPFactoryP(I->m[i]));
1399    }
1400  }
1401  else
1402  {
1403    WerrorS( feNotImplemented );
1404    return NULL;
1405  }
1406
1407  List<int> IL=neworderint(L);
1408  ListIterator<int> Li;
1409  StringSetS("");
1410  Li = IL;
1411  int offs=rPar(currRing);
1412  int* mark=(int*)omAlloc0((pVariables+offs)*sizeof(int));
1413  int cnt=pVariables+offs;
1414  loop
1415  {
1416    if(! Li.hasItem()) break;
1417    BOOLEAN done=TRUE;
1418    i=Li.getItem()-1;
1419    mark[i]=1;
1420    if (i<offs)
1421    {
1422      done=FALSE;
1423      //StringAppendS(currRing->parameter[i]);
1424    }
1425    else
1426    {
1427      StringAppendS(currRing->names[i-offs]);
1428    }
1429    Li++;
1430    cnt--;
1431    if(cnt==0) break;
1432    if (done) StringAppendS(",");
1433  }
1434  for(i=0;i<pVariables+offs;i++)
1435  {
1436    BOOLEAN done=TRUE;
1437    if(mark[i]==0)
1438    {
1439      if (i<offs)
1440      {
1441        done=FALSE;
1442        //StringAppendS(currRing->parameter[i]);
1443      }
1444      else
1445      {
1446        StringAppendS(currRing->names[i-offs]);
1447      }
1448      cnt--;
1449      if(cnt==0) break;
1450      if (done) StringAppendS(",");
1451    }
1452  }
1453  char * s=omStrDup(StringAppendS(""));
1454  if (s[strlen(s)-1]==',') s[strlen(s)-1]='\0';
1455  return s;
1456}
1457
1458BOOLEAN singclap_isSqrFree(poly f)
1459{
1460  BOOLEAN b=FALSE;
1461  Off(SW_RATIONAL);
1462  //  Q / Fp
1463  if (((nGetChar() == 0) || (nGetChar() > 1) )
1464  &&(currRing->parameter==NULL))
1465  {
1466    setCharacteristic( nGetChar() );
1467    CanonicalForm F( convSingPFactoryP( f ) );
1468    if((nGetChar()>1)&&(!F.isUnivariate()))
1469      goto err;
1470    b=(BOOLEAN)isSqrFree(F);
1471  }
1472  // and over Q(a) / Fp(a)
1473  else if (( nGetChar()==1 ) /* Q(a) */
1474  || (nGetChar() <-1))       /* Fp(a) */
1475  {
1476    if (nGetChar()==1) setCharacteristic( 0 );
1477    else               setCharacteristic( -nGetChar() );
1478    //if (currRing->minpoly!=NULL)
1479    //{
1480    //  CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
1481    //                                             currRing->algring);
1482    //  Variable a=rootOf(mipo);
1483    //  CanonicalForm F( convSingAPFactoryAP( f,a ) );
1484    //  ...
1485    //}
1486    //else
1487    {
1488      CanonicalForm F( convSingTrPFactoryP( f ) );
1489      b=(BOOLEAN)isSqrFree(F);
1490    }
1491    Off(SW_RATIONAL);
1492  }
1493  else
1494  {
1495err:
1496    WerrorS( feNotImplemented );
1497  }
1498  return b;
1499}
1500
1501poly singclap_det( const matrix m )
1502{
1503  int r=m->rows();
1504  if (r!=m->cols())
1505  {
1506    Werror("det of %d x %d matrix",r,m->cols());
1507    return NULL;
1508  }
1509  poly res=NULL;
1510  if (( nGetChar() == 0 || nGetChar() > 1 )
1511  && (currRing->parameter==NULL))
1512  {
1513    setCharacteristic( nGetChar() );
1514    CFMatrix M(r,r);
1515    int i,j;
1516    for(i=r;i>0;i--)
1517    {
1518      for(j=r;j>0;j--)
1519      {
1520        M(i,j)=convSingPFactoryP(MATELEM(m,i,j));
1521      }
1522    }
1523    res= convFactoryPSingP( determinant(M,r) ) ;
1524  }
1525  // and over Q(a) / Fp(a)
1526  else if (( nGetChar()==1 ) /* Q(a) */
1527  || (nGetChar() <-1))       /* Fp(a) */
1528  {
1529    if (nGetChar()==1) setCharacteristic( 0 );
1530    else               setCharacteristic( -nGetChar() );
1531    CFMatrix M(r,r);
1532    poly res;
1533    if (currRing->minpoly!=NULL)
1534    {
1535      CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
1536                                            currRing->algring);
1537      Variable a=rootOf(mipo);
1538      int i,j;
1539      for(i=r;i>0;i--)
1540      {
1541        for(j=r;j>0;j--)
1542        {
1543          M(i,j)=convSingAPFactoryAP(MATELEM(m,i,j),a,currRing);
1544        }
1545      }
1546      res= convFactoryAPSingAP( determinant(M,r),currRing ) ;
1547    }
1548    else
1549    {
1550      int i,j;
1551      for(i=r;i>0;i--)
1552      {
1553        for(j=r;j>0;j--)
1554        {
1555          M(i,j)=convSingTrPFactoryP(MATELEM(m,i,j));
1556        }
1557      }
1558      res= convFactoryPSingTrP( determinant(M,r) );
1559    }
1560  }
1561  else
1562    WerrorS( feNotImplemented );
1563  Off(SW_RATIONAL);
1564  return res;
1565}
1566
1567int singclap_det_i( intvec * m )
1568{
1569  setCharacteristic( 0 );
1570  CFMatrix M(m->rows(),m->cols());
1571  int i,j;
1572  for(i=m->rows();i>0;i--)
1573  {
1574    for(j=m->cols();j>0;j--)
1575    {
1576      M(i,j)=IMATELEM(*m,i,j);
1577    }
1578  }
1579  int res= convFactoryISingI( determinant(M,m->rows())) ;
1580  Off(SW_RATIONAL);
1581  return res;
1582}
1583napoly singclap_alglcm ( napoly f, napoly g )
1584{
1585
1586 // over Q(a) / Fp(a)
1587 if (nGetChar()==1) setCharacteristic( 0 );
1588 else               setCharacteristic( -nGetChar() );
1589 napoly res;
1590
1591 if (currRing->minpoly!=NULL)
1592 {
1593   CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
1594                                         currRing->algring);
1595   Variable a=rootOf(mipo);
1596   CanonicalForm F( convSingAFactoryA( f,a, currRing ) ),
1597                 G( convSingAFactoryA( g,a, currRing ) );
1598   CanonicalForm GCD;
1599
1600   // calculate gcd
1601   GCD = gcd( F, G );
1602
1603   // calculate lcm
1604   res= convFactoryASingA( (F/GCD)*G,currRing );
1605 }
1606 else
1607 {
1608   CanonicalForm F( convSingPFactoryP( f,currRing->algring ) ),
1609                 G( convSingPFactoryP( g,currRing->algring ) );
1610   CanonicalForm GCD;
1611   // calculate gcd
1612   GCD = gcd( F, G );
1613
1614   // calculate lcm
1615   res= convFactoryPSingP( (F/GCD)*G, currRing->algring );
1616 }
1617
1618 Off(SW_RATIONAL);
1619 return res;
1620}
1621
1622void singclap_algdividecontent ( napoly f, napoly g, napoly &ff, napoly &gg )
1623{
1624 // over Q(a) / Fp(a)
1625 if (nGetChar()==1) setCharacteristic( 0 );
1626 else               setCharacteristic( -nGetChar() );
1627 ff=gg=NULL;
1628 On(SW_RATIONAL);
1629
1630 if (currRing->minpoly!=NULL)
1631 {
1632   CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
1633                         currRing->algring);
1634   Variable a=rootOf(mipo);
1635   CanonicalForm F( convSingAFactoryA( f,a, currRing ) ),
1636                 G( convSingAFactoryA( g,a, currRing ) );
1637   CanonicalForm GCD;
1638
1639   GCD=gcd( F, G );
1640
1641   if ((GCD!=1) && (GCD!=0))
1642   {
1643     ff= convFactoryASingA( F/ GCD, currRing );
1644     gg= convFactoryASingA( G/ GCD, currRing );
1645   }
1646 }
1647 else
1648 {
1649   CanonicalForm F( convSingPFactoryP( f,currRing->algring ) ),
1650                 G( convSingPFactoryP( g,currRing->algring ) );
1651   CanonicalForm GCD;
1652
1653   GCD=gcd( F, G );
1654
1655   if ((GCD!=1) && (GCD!=0))
1656   {
1657     ff= convFactoryPSingP( F/ GCD, currRing->algring );
1658     gg= convFactoryPSingP( G/ GCD, currRing->algring );
1659   }
1660 }
1661
1662 Off(SW_RATIONAL);
1663}
1664
1665#if 0
1666lists singclap_chineseRemainder(lists x, lists q)
1667{
1668  //assume(x->nr == q->nr);
1669  //assume(x->nr >= 0);
1670  int n=x->nr+1;
1671  if ((x->nr<0) || (x->nr!=q->nr))
1672  {
1673    WerrorS("list are empty or not of equal length");
1674    return NULL;
1675  }
1676  lists res=(lists)omAlloc0Bin(slists_bin);
1677  CFArray X(1,n), Q(1,n);
1678  int i;
1679  for(i=0; i<n; i++)
1680  {
1681    if (x->m[i-1].Typ()==INT_CMD)
1682    {
1683      X[i]=(int)x->m[i-1].Data();
1684    }
1685    else if (x->m[i-1].Typ()==NUMBER_CMD)
1686    {
1687      number N=(number)x->m[i-1].Data();
1688      X[i]=convSingNFactoryN(N);
1689    }
1690    else
1691    {
1692      WerrorS("illegal type in chineseRemainder");
1693      omFreeBin(res,slists_bin);
1694      return NULL;
1695    }
1696    if (q->m[i-1].Typ()==INT_CMD)
1697    {
1698      Q[i]=(int)q->m[i-1].Data();
1699    }
1700    else if (q->m[i-1].Typ()==NUMBER_CMD)
1701    {
1702      number N=(number)x->m[i-1].Data();
1703      Q[i]=convSingNFactoryN(N);
1704    }
1705    else
1706    {
1707      WerrorS("illegal type in chineseRemainder");
1708      omFreeBin(res,slists_bin);
1709      return NULL;
1710    }
1711  }
1712  CanonicalForm r, prod;
1713  chineseRemainder( X, Q, r, prod );
1714  res->Init(2);
1715  res->m[0].rtyp=NUMBER_CMD;
1716  res->m[1].rtyp=NUMBER_CMD;
1717  res->m[0].data=(char *)convFactoryNSingN( r );
1718  res->m[1].data=(char *)convFactoryNSingN( prod );
1719  return res;
1720}
1721#endif
1722#endif
Note: See TracBrowser for help on using the repository browser.