source: git/kernel/clapsing.cc @ ca371d

spielwiese
Last change on this file since ca371d was 226dd7, checked in by Hans Schönemann <hannes@…>, 15 years ago
*hannes: better default in problem case git-svn-id: file:///usr/local/Singular/svn/trunk@11568 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 43.8 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2/****************************************
3*  Computer Algebra System SINGULAR     *
4****************************************/
5// $Id: clapsing.cc,v 1.36 2009-03-17 16:51:31 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            if (w!=NULL)
1049              (*w)[j]=1;
1050            res->m[j]=pOne();
1051          }
1052        }
1053        else
1054        {
1055          if (!count_Factors(res,w,j,ff,convFactoryAPSingAP( J.getItem().factor() )))
1056          {
1057            if (w!=NULL)
1058              (*w)[j]=1;
1059            res->m[j]=pOne();
1060          }
1061        }
1062      }
1063    }
1064    if (rField_is_Extension() && (!pIsConstantPoly(ff)))
1065    {
1066      singclap_factorize_retry++;
1067      if (singclap_factorize_retry<3)
1068      {
1069        int jj;
1070        #ifdef FACTORIZE2_DEBUG
1071        printf("factorize_retry\n");
1072        #endif
1073        intvec *ww=NULL;
1074        idTest(res);
1075        ideal h=singclap_factorize ( ff, &ww , with_exps);
1076        idTest(h);
1077        int l=(*v)->length();
1078        (*v)->resize(l+ww->length());
1079        for(jj=0;jj<ww->length();jj++)
1080          (**v)[jj+l]=(*ww)[jj];
1081        delete ww;
1082        ideal hh=idInit(IDELEMS(res)+IDELEMS(h),1);
1083        for(jj=IDELEMS(res)-1;jj>=0;jj--)
1084        {
1085          hh->m[jj]=res->m[jj];
1086          res->m[jj]=NULL;
1087        }
1088        for(jj=IDELEMS(h)-1;jj>=0;jj--)
1089        {
1090          hh->m[jj+IDELEMS(res)]=h->m[jj];
1091          h->m[jj]=NULL;
1092        }
1093        idDelete(&res);
1094        idDelete(&h);
1095        res=hh;
1096        idTest(res);
1097        ff=NULL;
1098      }
1099      else
1100      {
1101        WarnS("problem with factorize");
1102        #if 0
1103        pWrite(ff);
1104        idShow(res);
1105        #endif
1106        idDelete(&res);
1107        res=idInit(2,1);
1108        res->m[0]=pOne();
1109        res->m[1]=ff; ff=NULL;
1110      }
1111    }
1112    pDelete(&ff);
1113    if (N!=NULL)
1114    {
1115      pMult_nn(res->m[0],N);
1116      nDelete(&N);
1117      N=NULL;
1118    }
1119    // delete constants
1120    if (res!=NULL)
1121    {
1122      int i=IDELEMS(res)-1;
1123      int j=0;
1124      for(;i>=0;i--)
1125      {
1126        if ((res->m[i]!=NULL)
1127        && (pNext(res->m[i])==NULL)
1128        && (pIsConstant(res->m[i])))
1129        {
1130          if (with_exps!=0)
1131          {
1132            pDelete(&(res->m[i]));
1133            if ((v!=NULL) && ((*v)!=NULL))
1134              (**v)[i]=0;
1135            j++;
1136          }
1137          else if (i!=0)
1138          {
1139            while ((v!=NULL) && ((*v)!=NULL) && ((**v)[i]>1))
1140            {
1141              res->m[0]=pMult(res->m[0],pCopy(res->m[i]));
1142              (**v)[i]--;
1143            }
1144            res->m[0]=pMult(res->m[0],res->m[i]);
1145            res->m[i]=NULL;
1146            if ((v!=NULL) && ((*v)!=NULL))
1147              (**v)[i]=1;
1148            j++;
1149          }
1150        }
1151      }
1152      if (j>0)
1153      {
1154        idSkipZeroes(res);
1155        if ((v!=NULL) && ((*v)!=NULL))
1156        {
1157          intvec *w=*v;
1158          int len=IDELEMS(res);
1159          *v = new intvec( len );
1160          for (i=0,j=0;i<si_min(w->length(),len);i++)
1161          {
1162            if((*w)[i]!=0)
1163            {
1164              (**v)[j]=(*w)[i]; j++;
1165            }
1166          }
1167          delete w;
1168        }
1169      }
1170      if (res->m[0]==NULL)
1171      {
1172        res->m[0]=pOne();
1173      }
1174    }
1175  }
1176  if (rField_is_Q_a() && (currRing->minpoly!=NULL))
1177  {
1178    int i=IDELEMS(res)-1;
1179    int stop=1;
1180    if (with_exps!=0) stop=0;
1181    for(;i>=stop;i--)
1182    {
1183      pNorm(res->m[i]);
1184    }
1185    if (with_exps==0) pSetCoeff(res->m[0],old_lead_coeff);
1186    else nDelete(&old_lead_coeff);
1187  }
1188  else
1189    nDelete(&old_lead_coeff);
1190  errorreported=save_errorreported;
1191notImpl:
1192  if (res==NULL)
1193    WerrorS( feNotImplemented );
1194  if (NN!=NULL)
1195  {
1196    nDelete(&NN);
1197  }
1198  if (N!=NULL)
1199  {
1200    nDelete(&N);
1201  }
1202  //if (f!=NULL) pDelete(&f);
1203  //PrintS("......S\n");
1204  return res;
1205}
1206ideal singclap_sqrfree ( poly f)
1207{
1208  pTest(f);
1209#ifdef FACTORIZE2_DEBUG
1210  printf("singclap_sqrfree, degree %d\n",pTotaldegree(f));
1211#endif
1212  // with_exps: 3,1 return only true factors, no exponents
1213  //            2 return true factors and exponents
1214  //            0 return coeff, factors and exponents
1215  BOOLEAN save_errorreported=errorreported;
1216
1217  ideal res=NULL;
1218
1219  // handle factorize(0) =========================================
1220  if (f==NULL)
1221  {
1222    res=idInit(1,1);
1223    return res;
1224  }
1225  // handle factorize(mon) =========================================
1226  if (pNext(f)==NULL)
1227  {
1228    int i=0;
1229    int n=0;
1230    int e;
1231    for(i=pVariables;i>0;i--) if(pGetExp(f,i)!=0) n++;
1232    n++; // with coeff
1233    res=idInit(si_max(n,1),1);
1234    res->m[0]=pOne();
1235    pSetCoeff(res->m[0],nCopy(pGetCoeff(f)));
1236    if (n==0)
1237    {
1238      res->m[0]=pOne();
1239      // (**v)[0]=1; is already done
1240      return res;
1241    }
1242    for(i=pVariables;i>0;i--)
1243    {
1244      e=pGetExp(f,i);
1245      if(e!=0)
1246      {
1247        n--;
1248        poly p=pOne();
1249        pSetExp(p,i,1);
1250        pSetm(p);
1251        res->m[n]=p;
1252      }
1253    }
1254    return res;
1255  }
1256  //PrintS("S:");pWrite(f);PrintLn();
1257  // use factory/libfac in general ==============================
1258  Off(SW_RATIONAL);
1259  On(SW_SYMMETRIC_FF);
1260  #ifdef HAVE_NTL
1261  extern int prime_number;
1262  if(rField_is_Q()) prime_number=0;
1263  #endif
1264  CFFList L;
1265
1266  if (!rField_is_Zp() && !rField_is_Zp_a()) /* Q, Q(a) */
1267  {
1268    //if (f!=NULL) // already tested at start of routine
1269    {
1270      pCleardenom(f);
1271    }
1272  }
1273  else if (rField_is_Zp_a())
1274  {
1275    //if (f!=NULL) // already tested at start of routine
1276    if (singclap_factorize_retry==0)
1277    {
1278      pNorm(f);
1279      pCleardenom(f);
1280    }
1281  }
1282  if (rField_is_Q() || rField_is_Zp())
1283  {
1284    setCharacteristic( nGetChar() );
1285    CanonicalForm F( convSingPFactoryP( f ) );
1286    L = sqrFree( F );
1287  }
1288  #if 0
1289  else if (rField_is_GF())
1290  {
1291    int c=rChar(currRing);
1292    setCharacteristic( c, primepower(c) );
1293    CanonicalForm F( convSingGFFactoryGF( f ) );
1294    if (F.isUnivariate())
1295    {
1296      L = factorize( F );
1297    }
1298    else
1299    {
1300      goto notImpl;
1301    }
1302  }
1303  #endif
1304  // and over Q(a) / Fp(a)
1305  else if (rField_is_Extension())
1306  {
1307    if (rField_is_Q_a()) setCharacteristic( 0 );
1308    else                 setCharacteristic( -nGetChar() );
1309    if (currRing->minpoly!=NULL)
1310    {
1311      CanonicalForm mipo=convSingTrFactoryP(((lnumber)currRing->minpoly)->z);
1312      Variable a=rootOf(mipo);
1313      CanonicalForm F( convSingAPFactoryAP( f,a ) );
1314      CFFList SqrFreeMV( const CanonicalForm & f , const CanonicalForm & mipo=0) ;
1315
1316      L = SqrFreeMV( F,mipo );
1317      //WarnS("L = sqrFree( F,mipo );");
1318      //L = sqrFree( F );
1319    }
1320    else
1321    {
1322      CanonicalForm F( convSingTrPFactoryP( f ) );
1323      L = sqrFree( F );
1324    }
1325  }
1326  else
1327  {
1328    goto notImpl;
1329  }
1330  {
1331    poly ff=pCopy(f); // a copy for the retry stuff
1332    // convert into ideal
1333    int n = L.length();
1334    if (n==0) n=1;
1335    CFFListIterator J=L;
1336    int j=0;
1337    res = idInit( n ,1);
1338    for ( ; J.hasItem(); J++, j++ )
1339    {
1340      poly p;
1341      if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
1342        //count_Factors(res,*v,f, j, convFactoryPSingP( J.getItem().factor() );
1343        res->m[j] = convFactoryPSingP( J.getItem().factor() );
1344      #if 0
1345      else if (rField_is_GF())
1346        res->m[j] = convFactoryGFSingGF( J.getItem().factor() );
1347      #endif
1348      else if (rField_is_Extension())     /* Q(a), Fp(a) */
1349      {
1350        if (currRing->minpoly==NULL)
1351          res->m[j]=convFactoryPSingTrP( J.getItem().factor() );
1352        else
1353          res->m[j]=convFactoryAPSingAP( J.getItem().factor() );
1354      }
1355    }
1356    if (res->m[0]==NULL)
1357    {
1358      res->m[0]=pOne();
1359    }
1360  }
1361  errorreported=save_errorreported;
1362notImpl:
1363  if (res==NULL)
1364    WerrorS( feNotImplemented );
1365  return res;
1366}
1367matrix singclap_irrCharSeries ( ideal I)
1368{
1369#ifdef HAVE_LIBFAC_P
1370  if (idIs0(I)) return mpNew(1,1);
1371
1372  // for now there is only the possibility to handle polynomials over
1373  // Q and Fp ...
1374  matrix res=NULL;
1375  int i;
1376  Off(SW_RATIONAL);
1377  On(SW_SYMMETRIC_FF);
1378  CFList L;
1379  ListCFList LL;
1380  if (((nGetChar() == 0) || (nGetChar() > 1) )
1381  && (currRing->parameter==NULL))
1382  {
1383    setCharacteristic( nGetChar() );
1384    for(i=0;i<IDELEMS(I);i++)
1385    {
1386      poly p=I->m[i];
1387      if (p!=NULL)
1388      {
1389        p=pCopy(p);
1390        pCleardenom(p);
1391        L.append(convSingPFactoryP(p));
1392      }
1393    }
1394  }
1395  // and over Q(a) / Fp(a)
1396  else if (( nGetChar()==1 ) /* Q(a) */
1397  || (nGetChar() <-1))       /* Fp(a) */
1398  {
1399    if (nGetChar()==1) setCharacteristic( 0 );
1400    else               setCharacteristic( -nGetChar() );
1401    for(i=0;i<IDELEMS(I);i++)
1402    {
1403      poly p=I->m[i];
1404      if (p!=NULL)
1405      {
1406        p=pCopy(p);
1407        pCleardenom(p);
1408        L.append(convSingTrPFactoryP(p));
1409      }
1410    }
1411  }
1412  else
1413  {
1414    WerrorS( feNotImplemented );
1415    return res;
1416  }
1417
1418  // a very bad work-around --- FIX IT in libfac
1419  // should be fixed as of 2001/6/27
1420  int tries=0;
1421  int m,n;
1422  ListIterator<CFList> LLi;
1423  loop
1424  {
1425    LL=IrrCharSeries(L);
1426    m= LL.length(); // Anzahl Zeilen
1427    n=0;
1428    for ( LLi = LL; LLi.hasItem(); LLi++ )
1429    {
1430      n = si_max(LLi.getItem().length(),n);
1431    }
1432    if ((m!=0) && (n!=0)) break;
1433    tries++;
1434    if (tries>=5) break;
1435  }
1436  if ((m==0) || (n==0))
1437  {
1438    Warn("char_series returns %d x %d matrix from %d input polys (%d)",
1439      m,n,IDELEMS(I)+1,LL.length());
1440    iiWriteMatrix((matrix)I,"I",2,0);
1441    m=si_max(m,1);
1442    n=si_max(n,1);
1443  }
1444  res=mpNew(m,n);
1445  CFListIterator Li;
1446  for ( m=1, LLi = LL; LLi.hasItem(); LLi++, m++ )
1447  {
1448    for (n=1, Li = LLi.getItem(); Li.hasItem(); Li++, n++)
1449    {
1450      if ( (nGetChar() == 0) || (nGetChar() > 1) )
1451        MATELEM(res,m,n)=convFactoryPSingP(Li.getItem());
1452      else
1453        MATELEM(res,m,n)=convFactoryPSingTrP(Li.getItem());
1454    }
1455  }
1456  Off(SW_RATIONAL);
1457  return res;
1458#else
1459  return NULL;
1460#endif
1461}
1462
1463char* singclap_neworder ( ideal I)
1464{
1465#ifdef HAVE_LIBFAC_P
1466  int i;
1467  Off(SW_RATIONAL);
1468  On(SW_SYMMETRIC_FF);
1469  CFList L;
1470  if (((nGetChar() == 0) || (nGetChar() > 1) )
1471  && (currRing->parameter==NULL))
1472  {
1473    setCharacteristic( nGetChar() );
1474    for(i=0;i<IDELEMS(I);i++)
1475    {
1476      L.append(convSingPFactoryP(I->m[i]));
1477    }
1478  }
1479  // and over Q(a) / Fp(a)
1480  else if (( nGetChar()==1 ) /* Q(a) */
1481  || (nGetChar() <-1))       /* Fp(a) */
1482  {
1483    if (nGetChar()==1) setCharacteristic( 0 );
1484    else               setCharacteristic( -nGetChar() );
1485    for(i=0;i<IDELEMS(I);i++)
1486    {
1487      L.append(convSingTrPFactoryP(I->m[i]));
1488    }
1489  }
1490  else
1491  {
1492    WerrorS( feNotImplemented );
1493    return NULL;
1494  }
1495
1496  List<int> IL=neworderint(L);
1497  ListIterator<int> Li;
1498  StringSetS("");
1499  Li = IL;
1500  int offs=rPar(currRing);
1501  int* mark=(int*)omAlloc0((pVariables+offs)*sizeof(int));
1502  int cnt=pVariables+offs;
1503  loop
1504  {
1505    if(! Li.hasItem()) break;
1506    BOOLEAN done=TRUE;
1507    i=Li.getItem()-1;
1508    mark[i]=1;
1509    if (i<offs)
1510    {
1511      done=FALSE;
1512      //StringAppendS(currRing->parameter[i]);
1513    }
1514    else
1515    {
1516      StringAppendS(currRing->names[i-offs]);
1517    }
1518    Li++;
1519    cnt--;
1520    if(cnt==0) break;
1521    if (done) StringAppendS(",");
1522  }
1523  for(i=0;i<pVariables+offs;i++)
1524  {
1525    BOOLEAN done=TRUE;
1526    if(mark[i]==0)
1527    {
1528      if (i<offs)
1529      {
1530        done=FALSE;
1531        //StringAppendS(currRing->parameter[i]);
1532      }
1533      else
1534      {
1535        StringAppendS(currRing->names[i-offs]);
1536      }
1537      cnt--;
1538      if(cnt==0) break;
1539      if (done) StringAppendS(",");
1540    }
1541  }
1542  char * s=omStrDup(StringAppendS(""));
1543  if (s[strlen(s)-1]==',') s[strlen(s)-1]='\0';
1544  return s;
1545#else
1546  return NULL;
1547#endif
1548}
1549
1550BOOLEAN singclap_isSqrFree(poly f)
1551{
1552  BOOLEAN b=FALSE;
1553  Off(SW_RATIONAL);
1554  //  Q / Fp
1555  if (((nGetChar() == 0) || (nGetChar() > 1) )
1556  &&(currRing->parameter==NULL))
1557  {
1558    setCharacteristic( nGetChar() );
1559    CanonicalForm F( convSingPFactoryP( f ) );
1560    if((nGetChar()>1)&&(!F.isUnivariate()))
1561      goto err;
1562    b=(BOOLEAN)isSqrFree(F);
1563  }
1564  // and over Q(a) / Fp(a)
1565  else if (( nGetChar()==1 ) /* Q(a) */
1566  || (nGetChar() <-1))       /* Fp(a) */
1567  {
1568    if (nGetChar()==1) setCharacteristic( 0 );
1569    else               setCharacteristic( -nGetChar() );
1570    //if (currRing->minpoly!=NULL)
1571    //{
1572    //  CanonicalForm mipo=convSingTrFactoryP(((lnumber)currRing->minpoly)->z);
1573    //  Variable a=rootOf(mipo);
1574    //  CanonicalForm F( convSingAPFactoryAP( f,a ) );
1575    //  ...
1576    //}
1577    //else
1578    {
1579      CanonicalForm F( convSingTrPFactoryP( f ) );
1580      b=(BOOLEAN)isSqrFree(F);
1581    }
1582    Off(SW_RATIONAL);
1583  }
1584  else
1585  {
1586err:
1587    WerrorS( feNotImplemented );
1588  }
1589  return b;
1590}
1591
1592poly singclap_det( const matrix m )
1593{
1594  int r=m->rows();
1595  if (r!=m->cols())
1596  {
1597    Werror("det of %d x %d matrix",r,m->cols());
1598    return NULL;
1599  }
1600  poly res=NULL;
1601  if (( nGetChar() == 0 || nGetChar() > 1 )
1602  && (currRing->parameter==NULL))
1603  {
1604    setCharacteristic( nGetChar() );
1605    CFMatrix M(r,r);
1606    int i,j;
1607    for(i=r;i>0;i--)
1608    {
1609      for(j=r;j>0;j--)
1610      {
1611        M(i,j)=convSingPFactoryP(MATELEM(m,i,j));
1612      }
1613    }
1614    res= convFactoryPSingP( determinant(M,r) ) ;
1615  }
1616  // and over Q(a) / Fp(a)
1617  else if (( nGetChar()==1 ) /* Q(a) */
1618  || (nGetChar() <-1))       /* Fp(a) */
1619  {
1620    if (nGetChar()==1) setCharacteristic( 0 );
1621    else               setCharacteristic( -nGetChar() );
1622    CFMatrix M(r,r);
1623    poly res;
1624    if (currRing->minpoly!=NULL)
1625    {
1626      CanonicalForm mipo=convSingTrFactoryP(((lnumber)currRing->minpoly)->z);
1627      Variable a=rootOf(mipo);
1628      int i,j;
1629      for(i=r;i>0;i--)
1630      {
1631        for(j=r;j>0;j--)
1632        {
1633          M(i,j)=convSingAPFactoryAP(MATELEM(m,i,j),a);
1634        }
1635      }
1636      res= convFactoryAPSingAP( determinant(M,r) ) ;
1637    }
1638    else
1639    {
1640      int i,j;
1641      for(i=r;i>0;i--)
1642      {
1643        for(j=r;j>0;j--)
1644        {
1645          M(i,j)=convSingTrPFactoryP(MATELEM(m,i,j));
1646        }
1647      }
1648      res= convFactoryPSingTrP( determinant(M,r) );
1649    }
1650  }
1651  else
1652    WerrorS( feNotImplemented );
1653  Off(SW_RATIONAL);
1654  return res;
1655}
1656
1657int singclap_det_i( intvec * m )
1658{
1659  setCharacteristic( 0 );
1660  CFMatrix M(m->rows(),m->cols());
1661  int i,j;
1662  for(i=m->rows();i>0;i--)
1663  {
1664    for(j=m->cols();j>0;j--)
1665    {
1666      M(i,j)=IMATELEM(*m,i,j);
1667    }
1668  }
1669  int res= convFactoryISingI( determinant(M,m->rows())) ;
1670  Off(SW_RATIONAL);
1671  return res;
1672}
1673napoly singclap_alglcm ( napoly f, napoly g )
1674{
1675 FACTORY_ALGOUT( "f", f );
1676 FACTORY_ALGOUT( "g", g );
1677
1678 // over Q(a) / Fp(a)
1679 if (nGetChar()==1) setCharacteristic( 0 );
1680 else               setCharacteristic( -nGetChar() );
1681 napoly res;
1682
1683 if (currRing->minpoly!=NULL)
1684 {
1685   CanonicalForm mipo=convSingTrFactoryP(((lnumber)currRing->minpoly)->z);
1686   Variable a=rootOf(mipo);
1687   CanonicalForm F( convSingAFactoryA( f,a ) ), G( convSingAFactoryA( g,a ) );
1688   CanonicalForm GCD;
1689
1690   TIMING_START( algLcmTimer );
1691   // calculate gcd
1692#ifdef FACTORY_GCD_TEST
1693   GCD = CFPrimitiveGcdUtil::gcd( F, G );
1694#else
1695   GCD = gcd( F, G );
1696#endif
1697   TIMING_END( algLcmTimer );
1698
1699   FACTORY_CFAOUT( "d", GCD );
1700   FACTORY_GCDSTAT( "alcm:", F, G, GCD );
1701
1702   // calculate lcm
1703   res= convFactoryASingA( (F/GCD)*G );
1704 }
1705 else
1706 {
1707   CanonicalForm F( convSingTrFactoryP( f ) ), G( convSingTrFactoryP( g ) );
1708   CanonicalForm GCD;
1709   TIMING_START( algLcmTimer );
1710   // calculate gcd
1711#ifdef FACTORY_GCD_TEST
1712   GCD = CFPrimitiveGcdUtil::gcd( F, G );
1713#else
1714   GCD = gcd( F, G );
1715#endif
1716   TIMING_END( algLcmTimer );
1717
1718   FACTORY_CFTROUT( "d", GCD );
1719   FACTORY_GCDSTAT( "alcm:", F, G, GCD );
1720
1721   // calculate lcm
1722   res= convFactoryPSingTr( (F/GCD)*G );
1723 }
1724
1725 Off(SW_RATIONAL);
1726 return res;
1727}
1728
1729void singclap_algdividecontent ( napoly f, napoly g, napoly &ff, napoly &gg )
1730{
1731 FACTORY_ALGOUT( "f", f );
1732 FACTORY_ALGOUT( "g", g );
1733
1734 // over Q(a) / Fp(a)
1735 if (nGetChar()==1) setCharacteristic( 0 );
1736 else               setCharacteristic( -nGetChar() );
1737 ff=gg=NULL;
1738 On(SW_RATIONAL);
1739
1740 if (currRing->minpoly!=NULL)
1741 {
1742   CanonicalForm mipo=convSingTrFactoryP(((lnumber)currRing->minpoly)->z);
1743   Variable a=rootOf(mipo);
1744   CanonicalForm F( convSingAFactoryA( f,a ) ), G( convSingAFactoryA( g,a ) );
1745   CanonicalForm GCD;
1746
1747   TIMING_START( algContentTimer );
1748#ifdef FACTORY_GCD_TEST
1749   GCD=CFPrimitiveGcdUtil::gcd( F, G );
1750#else
1751   GCD=gcd( F, G );
1752#endif
1753   TIMING_END( algContentTimer );
1754
1755   FACTORY_CFAOUT( "d", GCD );
1756   FACTORY_GCDSTAT( "acnt:", F, G, GCD );
1757
1758   if ((GCD!=1) && (GCD!=0))
1759   {
1760     ff= convFactoryASingA( F/ GCD );
1761     gg= convFactoryASingA( G/ GCD );
1762   }
1763 }
1764 else
1765 {
1766   CanonicalForm F( convSingTrFactoryP( f ) ), G( convSingTrFactoryP( g ) );
1767   CanonicalForm GCD;
1768
1769   TIMING_START( algContentTimer );
1770#ifdef FACTORY_GCD_TEST
1771   GCD=CFPrimitiveGcdUtil::gcd( F, G );
1772#else
1773   GCD=gcd( F, G );
1774#endif
1775   TIMING_END( algContentTimer );
1776
1777   FACTORY_CFTROUT( "d", GCD );
1778   FACTORY_GCDSTAT( "acnt:", F, G, GCD );
1779
1780   if ((GCD!=1) && (GCD!=0))
1781   {
1782     ff= convFactoryPSingTr( F/ GCD );
1783     gg= convFactoryPSingTr( G/ GCD );
1784   }
1785 }
1786
1787 Off(SW_RATIONAL);
1788}
1789
1790#if 0
1791lists singclap_chineseRemainder(lists x, lists q)
1792{
1793  //assume(x->nr == q->nr);
1794  //assume(x->nr >= 0);
1795  int n=x->nr+1;
1796  if ((x->nr<0) || (x->nr!=q->nr))
1797  {
1798    WerrorS("list are empty or not of equal length");
1799    return NULL;
1800  }
1801  lists res=(lists)omAlloc0Bin(slists_bin);
1802  CFArray X(1,n), Q(1,n);
1803  int i;
1804  for(i=0; i<n; i++)
1805  {
1806    if (x->m[i-1].Typ()==INT_CMD)
1807    {
1808      X[i]=(int)x->m[i-1].Data();
1809    }
1810    else if (x->m[i-1].Typ()==NUMBER_CMD)
1811    {
1812      number N=(number)x->m[i-1].Data();
1813      X[i]=convSingNFactoryN(N);
1814    }
1815    else
1816    {
1817      WerrorS("illegal type in chineseRemainder");
1818      omFreeBin(res,slists_bin);
1819      return NULL;
1820    }
1821    if (q->m[i-1].Typ()==INT_CMD)
1822    {
1823      Q[i]=(int)q->m[i-1].Data();
1824    }
1825    else if (q->m[i-1].Typ()==NUMBER_CMD)
1826    {
1827      number N=(number)x->m[i-1].Data();
1828      Q[i]=convSingNFactoryN(N);
1829    }
1830    else
1831    {
1832      WerrorS("illegal type in chineseRemainder");
1833      omFreeBin(res,slists_bin);
1834      return NULL;
1835    }
1836  }
1837  CanonicalForm r, prod;
1838  chineseRemainder( X, Q, r, prod );
1839  res->Init(2);
1840  res->m[0].rtyp=NUMBER_CMD;
1841  res->m[1].rtyp=NUMBER_CMD;
1842  res->m[0].data=(char *)convFactoryNSingN( r );
1843  res->m[1].data=(char *)convFactoryNSingN( prod );
1844  return res;
1845}
1846#endif
1847#endif
Note: See TracBrowser for help on using the repository browser.