source: git/kernel/clapsing.cc @ 61944d0

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