source: git/kernel/clapsing.cc @ 1534d9

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