source: git/kernel/clapsing.cc @ fcb8022

fieker-DuValspielwiese
Last change on this file since fcb8022 was 5d68a6, checked in by Hans Schönemann <hannes@…>, 15 years ago
*hannes: factory conversion:GF git-svn-id: file:///usr/local/Singular/svn/trunk@11253 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 43.4 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2/****************************************
3*  Computer Algebra System SINGULAR     *
4****************************************/
5// $Id: clapsing.cc,v 1.32 2008-12-17 15:08:50 Singular Exp $
6/*
7* ABSTRACT: interface between Singular and factory
8*/
9
10//#define FACTORIZE2_DEBUG
11#include "mod2.h"
12#include "omalloc.h"
13#ifdef HAVE_FACTORY
14#define SI_DONT_HAVE_GLOBAL_VARS
15#include "structs.h"
16#include "clapsing.h"
17#include "numbers.h"
18#include "ring.h"
19#include "ffields.h"
20#include <factory.h>
21#include "clapconv.h"
22#ifdef HAVE_LIBFAC_P
23#include <factor.h>
24//CanonicalForm algcd(const CanonicalForm & F, const CanonicalForm & g, const CFList & as, const Varlist & order);
25CanonicalForm alg_gcd(const CanonicalForm &, const CanonicalForm &, const CFList &);
26#endif
27#include "ring.h"
28
29//
30// FACTORY_GCD_TEST: use new gcd instead of old one.  Does not work
31//   without new gcd-implementation which is not publicly available.
32//
33// FACTORY_GCD_STAT: print statistics on polynomials.  Works only
34//   with the file `gcd_stat.cc' and `gcd_stat.h which may be found
35//   in the repository, module `factory-devel'.
36//   Overall statistics may printed using `system("gcdstat");'.
37//
38// FACTORY_GCD_TIMING: accumulate time used for gcd calculations.
39//   Time may be printed (and reset) with `system("gcdtime");'.
40//   For this define, `timing.h' from the factory source directory
41//   has to be copied to the Singular source directory.
42//   Note: for better readability, the macros `TIMING_START()' and
43//   `TIMING_END()' are used in any case.  However, they expand to
44//   nothing if `FACTORY_GCD_TIMING' is off.
45//
46// FACTORY_GCD_DEBOUT: print polynomials involved in gcd calculations.
47//   The polynomials are printed by means of the macros
48//   `FACTORY_*OUT_POLY' which are defined to be empty if
49//   `FACTORY_GCD_DEBOUT' is off.
50//
51// FACTORY_GCD_DEBOUT_PATTERN: print degree patterns of polynomials
52//   involved in gcd calculations.
53//   The patterns are printed by means of the macros
54//   `FACTORY_*OUT_PAT' which are defined to be empty if
55//   `FACTORY_GCD_DEBOUT_PATTERN' is off.
56//
57//   A degree pattern looks like this:
58//
59//   totDeg  size    deg(v1) deg(v2) ...
60//
61//   where "totDeg" means total degree, "size" the number of terms,
62//   and "deg(vi)" is the degree with respect to variable i.
63//   In univariate case, the "deg(vi)" are missing.  For this feature
64//   you need the files `gcd_stat.cc' and `gcd_stat.h'.
65//
66//
67// And here is what the functions print if `FACTORY_GCD_DEBOUT' (1),
68// `FACTORY_GCD_STAT' (2), or `FACTORY_GCD_DEBOUT_PATTERN' (3) is on:
69//
70// sinclap_divide_content:
71// (1) G = <firstCoeff>
72// (3) G#= <firstCoeff, pattern>
73// (1) h = <nextCoeff>
74// (3) h#= <nextCoeff, pattern>
75// (2) gcnt: <statistics on gcd as explained above>
76// (1) g = <intermediateResult>
77// (3) g#= <intermediateResult, pattern>
78// (1) h = <nextCoeff>
79// (3) h#= <nextCoeff, pattern>
80// (2) gcnt: <statistics on gcd as explained above>
81//  ...
82// (1) h = <lastCoeff>
83// (3) h#= <lastCoeff, pattern>
84// (1) g = <finalResult>
85// (3) g#= <finalResult, pattern>
86// (2) gcnt: <statistics on gcd as explained above>
87// (2) cont: <statistics on content as explained above>
88//
89// singclap_alglcm:
90// (1) f = <inputPolyF>
91// (3) f#= <inputPolyF, pattern>
92// (1) g = <inputPolyG>
93// (3) g#= <inputPolyG, pattern>
94// (1) d = <its gcd>
95// (3) d#= <its gcd, pattern>
96// (2) alcm: <statistics as explained above>
97//
98// singclap_algdividecontent:
99// (1) f = <inputPolyF>
100// (3) f#= <inputPolyF, pattern>
101// (1) g = <inputPolyG>
102// (3) g#= <inputPolyG, pattern>
103// (1) d = <its gcd>
104// (3) d#= <its gcd, pattern>
105// (2) acnt: <statistics as explained above>
106//
107
108#ifdef FACTORY_GCD_STAT
109#include "gcd_stat.h"
110#define FACTORY_GCDSTAT( tag, f, g, d ) \
111  printGcdStat( tag, f, g, d )
112#define FACTORY_CONTSTAT( tag, f ) \
113  printContStat( tag, f )
114#else
115#define FACTORY_GCDSTAT( tag, f, g, d )
116#define FACTORY_CONTSTAT( tag, f )
117#endif
118
119#ifdef FACTORY_GCD_TIMING
120#define TIMING
121#include "timing.h"
122TIMING_DEFINE_PRINT( contentTimer );
123TIMING_DEFINE_PRINT( algContentTimer );
124TIMING_DEFINE_PRINT( algLcmTimer );
125#else
126#define TIMING_START( timer )
127#define TIMING_END( timer )
128#endif
129
130#ifdef FACTORY_GCD_DEBOUT
131#include "longalg.h"
132#include "febase.h"
133// napoly f
134#define FACTORY_ALGOUT_POLY( tag, f ) \
135  StringSetS( tag ); \
136  napWrite( f ); \
137  pRINtS(StringAppendS("\n"));
138// CanonicalForm f, represents transcendent extension
139#define FACTORY_CFTROUT_POLY( tag, f ) \
140  { \
141    napoly F=convFactoryPSingTr( f ); \
142    StringSetS( tag ); \
143    napWrite( F ); \
144    PrintS(StringAppendS("\n")); \
145    napDelete(&F); \
146  }
147// CanonicalForm f, represents algebraic extension
148#define FACTORY_CFAOUT_POLY( tag, f ) \
149  { \
150    napoly F=convFactoryASingA( f ); \
151    StringSetS( tag ); \
152    napWrite( F ); \
153    PrintS(StringAppendS("\n")); \
154    napDelete(&F); \
155  }
156#else /* ! FACTORY_GCD_DEBOUT */
157#define FACTORY_ALGOUT_POLY( tag, f )
158#define FACTORY_CFTROUT_POLY( tag, f )
159#define FACTORY_CFAOUT_POLY( tag, f )
160#endif /* ! FACTORY_GCD_DEBOUT */
161
162#ifdef FACTORY_GCD_DEBOUT_PATTERN
163// napoly f
164#define FACTORY_ALGOUT_PAT( tag, f ) \
165  if (currRing->minpoly!=NULL) \
166  { \
167    CanonicalForm mipo=convSingTrFactoryP(((lnumber)currRing->minpoly)->z); \
168    Variable a=rootOf(mipo); \
169    printPolyPattern( tag, convSingAFactoryA( f,a ), rPar( currRing ) ); \
170  } \
171  else \
172  { \
173    printPolyPattern( tag, convSingTrFactoryP( f ), rPar( currRing ) ); \
174  }
175// CanonicalForm f, represents transcendent extension
176#define FACTORY_CFTROUT_PAT( tag, f ) printPolyPattern( tag, f, rPar( currRing ) )
177// CanonicalForm f, represents algebraic extension
178#define FACTORY_CFAOUT_PAT( tag, f ) printPolyPattern( tag, f, rPar( currRing ) )
179#else /* ! FACTORY_GCD_DEBOUT_PATTERN */
180#define FACTORY_ALGOUT_PAT( tag, f )
181#define FACTORY_CFTROUT_PAT( tag, f )
182#define FACTORY_CFAOUT_PAT( tag, f )
183#endif /* ! FACTORY_GCD_DEBOUT_PATTERN */
184
185// these macors combine both print macros
186#define FACTORY_ALGOUT( tag, f ) \
187  FACTORY_ALGOUT_POLY( tag " = ", f ); \
188  FACTORY_ALGOUT_PAT( tag "#= ", f )
189#define FACTORY_CFTROUT( tag, f ) \
190  FACTORY_CFTROUT_POLY( tag " = ", f ); \
191  FACTORY_CFTROUT_PAT( tag "#= ", f )
192#define FACTORY_CFAOUT( tag, f ) \
193  FACTORY_CFAOUT_POLY( tag " = ", f ); \
194  FACTORY_CFAOUT_PAT( tag "#= ", f )
195
196
197
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        //  L=Factorize(G, mipo);
975#else
976        WarnS("complete factorization only for univariate polynomials");
977        if (rField_is_Q_a() ||(!F.isUnivariate())) /* Q(a) */
978        {
979          L = factorize( G );
980        }
981        else
982        {
983          L = factorize( G, a );
984        }
985#endif
986      }
987    }
988    else
989    {
990      CanonicalForm F( convSingTrPFactoryP( f ) );
991      if (rField_is_Q_a())
992      {
993        L = factorize( F );
994      }
995      else /* Fp(a) */
996      {
997#ifdef HAVE_LIBFAC_P
998        L = Factorize( F );
999#else
1000        goto notImpl;
1001#endif
1002      }
1003    }
1004  }
1005  else
1006  {
1007    goto notImpl;
1008  }
1009  {
1010    poly ff=pCopy(f); // a copy for the retry stuff
1011    // the first factor should be a constant
1012    if ( ! L.getFirst().factor().inCoeffDomain() )
1013      L.insert(CFFactor(1,1));
1014    // convert into ideal
1015    int n = L.length();
1016    if (n==0) n=1;
1017    CFFListIterator J=L;
1018    int j=0;
1019    if (with_exps!=1)
1020    {
1021      if ((with_exps==2)&&(n>1))
1022      {
1023        n--;
1024        J++;
1025      }
1026      *v = new intvec( n );
1027    }
1028    res = idInit( n ,1);
1029    for ( ; J.hasItem(); J++, j++ )
1030    {
1031      poly p;
1032      if (with_exps!=1) (**v)[j] = J.getItem().exp();
1033      if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
1034        //count_Factors(res,*v,f, j, convFactoryPSingP( J.getItem().factor() );
1035        res->m[j] = convFactoryPSingP( J.getItem().factor() );
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          count_Factors(res,w,j,ff,convFactoryPSingTrP( J.getItem().factor() ));
1046        else
1047          count_Factors(res,w,j,ff,convFactoryAPSingAP( J.getItem().factor() ));
1048      }
1049    }
1050    if (rField_is_Extension() && (!pIsConstantPoly(ff)))
1051    {
1052      singclap_factorize_retry++;
1053      if (singclap_factorize_retry<3)
1054      {
1055        int jj;
1056        #ifdef FACTORIZE2_DEBUG
1057        printf("factorize_retry\n");
1058        #endif
1059        intvec *ww=NULL;
1060        idTest(res);
1061        ideal h=singclap_factorize ( ff, &ww , with_exps);
1062        idTest(h);
1063        int l=(*v)->length();
1064        (*v)->resize(l+ww->length());
1065        for(jj=0;jj<ww->length();jj++)
1066          (**v)[jj+l]=(*ww)[jj];
1067        delete ww;
1068        ideal hh=idInit(IDELEMS(res)+IDELEMS(h),1);
1069        for(jj=IDELEMS(res)-1;jj>=0;jj--)
1070        {
1071          hh->m[jj]=res->m[jj];
1072          res->m[jj]=NULL;
1073        }
1074        for(jj=IDELEMS(h)-1;jj>=0;jj--)
1075        {
1076          hh->m[jj+IDELEMS(res)]=h->m[jj];
1077          h->m[jj]=NULL;
1078        }
1079        idDelete(&res);
1080        idDelete(&h);
1081        res=hh;
1082        idTest(res);
1083        ff=NULL;
1084      }
1085      else
1086      {
1087        WarnS("problem with factorize");
1088      }
1089    }
1090    pDelete(&ff);
1091    if (N!=NULL)
1092    {
1093      pMult_nn(res->m[0],N);
1094      nDelete(&N);
1095      N=NULL;
1096    }
1097    // delete constants
1098    if (res!=NULL)
1099    {
1100      int i=IDELEMS(res)-1;
1101      int j=0;
1102      for(;i>=0;i--)
1103      {
1104        if ((res->m[i]!=NULL)
1105        && (pNext(res->m[i])==NULL)
1106        && (pIsConstant(res->m[i])))
1107        {
1108          if (with_exps!=0)
1109          {
1110            pDelete(&(res->m[i]));
1111            if ((v!=NULL) && ((*v)!=NULL))
1112              (**v)[i]=0;
1113            j++;
1114          }
1115          else if (i!=0)
1116          {
1117            while ((v!=NULL) && ((*v)!=NULL) && ((**v)[i]>1))
1118            {
1119              res->m[0]=pMult(res->m[0],pCopy(res->m[i]));
1120              (**v)[i]--;
1121            }
1122            res->m[0]=pMult(res->m[0],res->m[i]);
1123            res->m[i]=NULL;
1124            if ((v!=NULL) && ((*v)!=NULL))
1125              (**v)[i]=0;
1126            j++;
1127          }
1128        }
1129      }
1130      if (j>0)
1131      {
1132        idSkipZeroes(res);
1133        if ((v!=NULL) && ((*v)!=NULL))
1134        {
1135          intvec *w=*v;
1136          int len=IDELEMS(res);
1137          *v = new intvec( len );
1138          for (i=0,j=0;i<si_min(w->length(),len);i++)
1139          {
1140            if((*w)[i]!=0)
1141            {
1142              (**v)[j]=(*w)[i]; j++;
1143            }
1144          }
1145          delete w;
1146        }
1147      }
1148      if (res->m[0]==NULL)
1149      {
1150        res->m[0]=pOne();
1151      }
1152    }
1153  }
1154  if (rField_is_Q_a() && (currRing->minpoly!=NULL))
1155  {
1156    int i=IDELEMS(res)-1;
1157    int stop=1;
1158    if (with_exps!=0) stop=0;
1159    for(;i>=stop;i--)
1160    {
1161      pNorm(res->m[i]);
1162    }
1163    if (with_exps==0) pSetCoeff(res->m[0],old_lead_coeff);
1164    else nDelete(&old_lead_coeff);
1165  }
1166  else
1167    nDelete(&old_lead_coeff);
1168  errorreported=save_errorreported;
1169notImpl:
1170  if (res==NULL)
1171    WerrorS( feNotImplemented );
1172  if (NN!=NULL)
1173  {
1174    nDelete(&NN);
1175  }
1176  if (N!=NULL)
1177  {
1178    nDelete(&N);
1179  }
1180  //if (f!=NULL) pDelete(&f);
1181  //PrintS("......S\n");
1182  return res;
1183}
1184ideal singclap_sqrfree ( poly f)
1185{
1186  pTest(f);
1187#ifdef FACTORIZE2_DEBUG
1188  printf("singclap_sqrfree, degree %d\n",pTotaldegree(f));
1189#endif
1190  // with_exps: 3,1 return only true factors, no exponents
1191  //            2 return true factors and exponents
1192  //            0 return coeff, factors and exponents
1193  BOOLEAN save_errorreported=errorreported;
1194
1195  ideal res=NULL;
1196
1197  // handle factorize(0) =========================================
1198  if (f==NULL)
1199  {
1200    res=idInit(1,1);
1201    return res;
1202  }
1203  // handle factorize(mon) =========================================
1204  if (pNext(f)==NULL)
1205  {
1206    int i=0;
1207    int n=0;
1208    int e;
1209    for(i=pVariables;i>0;i--) if(pGetExp(f,i)!=0) n++;
1210    n++; // with coeff
1211    res=idInit(si_max(n,1),1);
1212    res->m[0]=pOne();
1213    pSetCoeff(res->m[0],nCopy(pGetCoeff(f)));
1214    if (n==0)
1215    {
1216      res->m[0]=pOne();
1217      // (**v)[0]=1; is already done
1218      return res;
1219    }
1220    for(i=pVariables;i>0;i--)
1221    {
1222      e=pGetExp(f,i);
1223      if(e!=0)
1224      {
1225        n--;
1226        poly p=pOne();
1227        pSetExp(p,i,1);
1228        pSetm(p);
1229        res->m[n]=p;
1230      }
1231    }
1232    return res;
1233  }
1234  //PrintS("S:");pWrite(f);PrintLn();
1235  // use factory/libfac in general ==============================
1236  Off(SW_RATIONAL);
1237  On(SW_SYMMETRIC_FF);
1238  #ifdef HAVE_NTL
1239  extern int prime_number;
1240  if(rField_is_Q()) prime_number=0;
1241  #endif
1242  CFFList L;
1243
1244  if (!rField_is_Zp() && !rField_is_Zp_a()) /* Q, Q(a) */
1245  {
1246    //if (f!=NULL) // already tested at start of routine
1247    {
1248      pCleardenom(f);
1249    }
1250  }
1251  else if (rField_is_Zp_a())
1252  {
1253    //if (f!=NULL) // already tested at start of routine
1254    if (singclap_factorize_retry==0)
1255    {
1256      pNorm(f);
1257      pCleardenom(f);
1258    }
1259  }
1260  if (rField_is_Q() || rField_is_Zp())
1261  {
1262    setCharacteristic( nGetChar() );
1263    CanonicalForm F( convSingPFactoryP( f ) );
1264    L = sqrFree( F );
1265  }
1266  #if 0
1267  else if (rField_is_GF())
1268  {
1269    int c=rChar(currRing);
1270    setCharacteristic( c, primepower(c) );
1271    CanonicalForm F( convSingGFFactoryGF( f ) );
1272    if (F.isUnivariate())
1273    {
1274      L = factorize( F );
1275    }
1276    else
1277    {
1278      goto notImpl;
1279    }
1280  }
1281  #endif
1282  // and over Q(a) / Fp(a)
1283  else if (rField_is_Extension())
1284  {
1285    if (rField_is_Q_a()) setCharacteristic( 0 );
1286    else                 setCharacteristic( -nGetChar() );
1287    if (currRing->minpoly!=NULL)
1288    {
1289      CanonicalForm mipo=convSingTrFactoryP(((lnumber)currRing->minpoly)->z);
1290      Variable a=rootOf(mipo);
1291      CanonicalForm F( convSingAPFactoryAP( f,a ) );
1292      CFFList SqrFreeMV( const CanonicalForm & f , const CanonicalForm & mipo=0) ;
1293
1294      L = SqrFreeMV( F,mipo );
1295      //WarnS("L = sqrFree( F,mipo );");
1296      //L = sqrFree( F );
1297    }
1298    else
1299    {
1300      CanonicalForm F( convSingTrPFactoryP( f ) );
1301      L = sqrFree( F );
1302    }
1303  }
1304  else
1305  {
1306    goto notImpl;
1307  }
1308  {
1309    poly ff=pCopy(f); // a copy for the retry stuff
1310    // convert into ideal
1311    int n = L.length();
1312    if (n==0) n=1;
1313    CFFListIterator J=L;
1314    int j=0;
1315    res = idInit( n ,1);
1316    for ( ; J.hasItem(); J++, j++ )
1317    {
1318      poly p;
1319      if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
1320        //count_Factors(res,*v,f, j, convFactoryPSingP( J.getItem().factor() );
1321        res->m[j] = convFactoryPSingP( J.getItem().factor() );
1322      #if 0
1323      else if (rField_is_GF())
1324        res->m[j] = convFactoryGFSingGF( J.getItem().factor() );
1325      #endif
1326      else if (rField_is_Extension())     /* Q(a), Fp(a) */
1327      {
1328        if (currRing->minpoly==NULL)
1329          res->m[j]=convFactoryPSingTrP( J.getItem().factor() );
1330        else
1331          res->m[j]=convFactoryAPSingAP( J.getItem().factor() );
1332      }
1333    }
1334    if (res->m[0]==NULL)
1335    {
1336      res->m[0]=pOne();
1337    }
1338  }
1339  errorreported=save_errorreported;
1340notImpl:
1341  if (res==NULL)
1342    WerrorS( feNotImplemented );
1343  return res;
1344}
1345matrix singclap_irrCharSeries ( ideal I)
1346{
1347#ifdef HAVE_LIBFAC_P
1348  if (idIs0(I)) return mpNew(1,1);
1349
1350  // for now there is only the possibility to handle polynomials over
1351  // Q and Fp ...
1352  matrix res=NULL;
1353  int i;
1354  Off(SW_RATIONAL);
1355  On(SW_SYMMETRIC_FF);
1356  CFList L;
1357  ListCFList LL;
1358  if (((nGetChar() == 0) || (nGetChar() > 1) )
1359  && (currRing->parameter==NULL))
1360  {
1361    setCharacteristic( nGetChar() );
1362    for(i=0;i<IDELEMS(I);i++)
1363    {
1364      poly p=I->m[i];
1365      if (p!=NULL)
1366      {
1367        p=pCopy(p);
1368        pCleardenom(p);
1369        L.append(convSingPFactoryP(p));
1370      }
1371    }
1372  }
1373  // and over Q(a) / Fp(a)
1374  else if (( nGetChar()==1 ) /* Q(a) */
1375  || (nGetChar() <-1))       /* Fp(a) */
1376  {
1377    if (nGetChar()==1) setCharacteristic( 0 );
1378    else               setCharacteristic( -nGetChar() );
1379    for(i=0;i<IDELEMS(I);i++)
1380    {
1381      poly p=I->m[i];
1382      if (p!=NULL)
1383      {
1384        p=pCopy(p);
1385        pCleardenom(p);
1386        L.append(convSingTrPFactoryP(p));
1387      }
1388    }
1389  }
1390  else
1391  {
1392    WerrorS( feNotImplemented );
1393    return res;
1394  }
1395
1396  // a very bad work-around --- FIX IT in libfac
1397  // should be fixed as of 2001/6/27
1398  int tries=0;
1399  int m,n;
1400  ListIterator<CFList> LLi;
1401  loop
1402  {
1403    LL=IrrCharSeries(L);
1404    m= LL.length(); // Anzahl Zeilen
1405    n=0;
1406    for ( LLi = LL; LLi.hasItem(); LLi++ )
1407    {
1408      n = si_max(LLi.getItem().length(),n);
1409    }
1410    if ((m!=0) && (n!=0)) break;
1411    tries++;
1412    if (tries>=5) break;
1413  }
1414  if ((m==0) || (n==0))
1415  {
1416    Warn("char_series returns %d x %d matrix from %d input polys (%d)",
1417      m,n,IDELEMS(I)+1,LL.length());
1418    iiWriteMatrix((matrix)I,"I",2,0);
1419    m=si_max(m,1);
1420    n=si_max(n,1);
1421  }
1422  res=mpNew(m,n);
1423  CFListIterator Li;
1424  for ( m=1, LLi = LL; LLi.hasItem(); LLi++, m++ )
1425  {
1426    for (n=1, Li = LLi.getItem(); Li.hasItem(); Li++, n++)
1427    {
1428      if ( (nGetChar() == 0) || (nGetChar() > 1) )
1429        MATELEM(res,m,n)=convFactoryPSingP(Li.getItem());
1430      else
1431        MATELEM(res,m,n)=convFactoryPSingTrP(Li.getItem());
1432    }
1433  }
1434  Off(SW_RATIONAL);
1435  return res;
1436#else
1437  return NULL;
1438#endif
1439}
1440
1441char* singclap_neworder ( ideal I)
1442{
1443#ifdef HAVE_LIBFAC_P
1444  int i;
1445  Off(SW_RATIONAL);
1446  On(SW_SYMMETRIC_FF);
1447  CFList L;
1448  if (((nGetChar() == 0) || (nGetChar() > 1) )
1449  && (currRing->parameter==NULL))
1450  {
1451    setCharacteristic( nGetChar() );
1452    for(i=0;i<IDELEMS(I);i++)
1453    {
1454      L.append(convSingPFactoryP(I->m[i]));
1455    }
1456  }
1457  // and over Q(a) / Fp(a)
1458  else if (( nGetChar()==1 ) /* Q(a) */
1459  || (nGetChar() <-1))       /* Fp(a) */
1460  {
1461    if (nGetChar()==1) setCharacteristic( 0 );
1462    else               setCharacteristic( -nGetChar() );
1463    for(i=0;i<IDELEMS(I);i++)
1464    {
1465      L.append(convSingTrPFactoryP(I->m[i]));
1466    }
1467  }
1468  else
1469  {
1470    WerrorS( feNotImplemented );
1471    return NULL;
1472  }
1473
1474  List<int> IL=neworderint(L);
1475  ListIterator<int> Li;
1476  StringSetS("");
1477  Li = IL;
1478  int offs=rPar(currRing);
1479  int* mark=(int*)omAlloc0((pVariables+offs)*sizeof(int));
1480  int cnt=pVariables+offs;
1481  loop
1482  {
1483    if(! Li.hasItem()) break;
1484    BOOLEAN done=TRUE;
1485    i=Li.getItem()-1;
1486    mark[i]=1;
1487    if (i<offs)
1488    {
1489      done=FALSE;
1490      //StringAppendS(currRing->parameter[i]);
1491    }
1492    else
1493    {
1494      StringAppendS(currRing->names[i-offs]);
1495    }
1496    Li++;
1497    cnt--;
1498    if(cnt==0) break;
1499    if (done) StringAppendS(",");
1500  }
1501  for(i=0;i<pVariables+offs;i++)
1502  {
1503    BOOLEAN done=TRUE;
1504    if(mark[i]==0)
1505    {
1506      if (i<offs)
1507      {
1508        done=FALSE;
1509        //StringAppendS(currRing->parameter[i]);
1510      }
1511      else
1512      {
1513        StringAppendS(currRing->names[i-offs]);
1514      }
1515      cnt--;
1516      if(cnt==0) break;
1517      if (done) StringAppendS(",");
1518    }
1519  }
1520  char * s=omStrDup(StringAppendS(""));
1521  if (s[strlen(s)-1]==',') s[strlen(s)-1]='\0';
1522  return s;
1523#else
1524  return NULL;
1525#endif
1526}
1527
1528BOOLEAN singclap_isSqrFree(poly f)
1529{
1530  BOOLEAN b=FALSE;
1531  Off(SW_RATIONAL);
1532  //  Q / Fp
1533  if (((nGetChar() == 0) || (nGetChar() > 1) )
1534  &&(currRing->parameter==NULL))
1535  {
1536    setCharacteristic( nGetChar() );
1537    CanonicalForm F( convSingPFactoryP( f ) );
1538    if((nGetChar()>1)&&(!F.isUnivariate()))
1539      goto err;
1540    b=(BOOLEAN)isSqrFree(F);
1541  }
1542  // and over Q(a) / Fp(a)
1543  else if (( nGetChar()==1 ) /* Q(a) */
1544  || (nGetChar() <-1))       /* Fp(a) */
1545  {
1546    if (nGetChar()==1) setCharacteristic( 0 );
1547    else               setCharacteristic( -nGetChar() );
1548    //if (currRing->minpoly!=NULL)
1549    //{
1550    //  CanonicalForm mipo=convSingTrFactoryP(((lnumber)currRing->minpoly)->z);
1551    //  Variable a=rootOf(mipo);
1552    //  CanonicalForm F( convSingAPFactoryAP( f,a ) );
1553    //  ...
1554    //}
1555    //else
1556    {
1557      CanonicalForm F( convSingTrPFactoryP( f ) );
1558      b=(BOOLEAN)isSqrFree(F);
1559    }
1560    Off(SW_RATIONAL);
1561  }
1562  else
1563  {
1564err:
1565    WerrorS( feNotImplemented );
1566  }
1567  return b;
1568}
1569
1570poly singclap_det( const matrix m )
1571{
1572  int r=m->rows();
1573  if (r!=m->cols())
1574  {
1575    Werror("det of %d x %d matrix",r,m->cols());
1576    return NULL;
1577  }
1578  poly res=NULL;
1579  if (( nGetChar() == 0 || nGetChar() > 1 )
1580  && (currRing->parameter==NULL))
1581  {
1582    setCharacteristic( nGetChar() );
1583    CFMatrix M(r,r);
1584    int i,j;
1585    for(i=r;i>0;i--)
1586    {
1587      for(j=r;j>0;j--)
1588      {
1589        M(i,j)=convSingPFactoryP(MATELEM(m,i,j));
1590      }
1591    }
1592    res= convFactoryPSingP( determinant(M,r) ) ;
1593  }
1594  // and over Q(a) / Fp(a)
1595  else if (( nGetChar()==1 ) /* Q(a) */
1596  || (nGetChar() <-1))       /* Fp(a) */
1597  {
1598    if (nGetChar()==1) setCharacteristic( 0 );
1599    else               setCharacteristic( -nGetChar() );
1600    CFMatrix M(r,r);
1601    poly res;
1602    if (currRing->minpoly!=NULL)
1603    {
1604      CanonicalForm mipo=convSingTrFactoryP(((lnumber)currRing->minpoly)->z);
1605      Variable a=rootOf(mipo);
1606      int i,j;
1607      for(i=r;i>0;i--)
1608      {
1609        for(j=r;j>0;j--)
1610        {
1611          M(i,j)=convSingAPFactoryAP(MATELEM(m,i,j),a);
1612        }
1613      }
1614      res= convFactoryAPSingAP( determinant(M,r) ) ;
1615    }
1616    else
1617    {
1618      int i,j;
1619      for(i=r;i>0;i--)
1620      {
1621        for(j=r;j>0;j--)
1622        {
1623          M(i,j)=convSingTrPFactoryP(MATELEM(m,i,j));
1624        }
1625      }
1626      res= convFactoryPSingTrP( determinant(M,r) );
1627    }
1628  }
1629  else
1630    WerrorS( feNotImplemented );
1631  Off(SW_RATIONAL);
1632  return res;
1633}
1634
1635int singclap_det_i( intvec * m )
1636{
1637  setCharacteristic( 0 );
1638  CFMatrix M(m->rows(),m->cols());
1639  int i,j;
1640  for(i=m->rows();i>0;i--)
1641  {
1642    for(j=m->cols();j>0;j--)
1643    {
1644      M(i,j)=IMATELEM(*m,i,j);
1645    }
1646  }
1647  int res= convFactoryISingI( determinant(M,m->rows())) ;
1648  Off(SW_RATIONAL);
1649  return res;
1650}
1651napoly singclap_alglcm ( napoly f, napoly g )
1652{
1653 FACTORY_ALGOUT( "f", f );
1654 FACTORY_ALGOUT( "g", g );
1655
1656 // over Q(a) / Fp(a)
1657 if (nGetChar()==1) setCharacteristic( 0 );
1658 else               setCharacteristic( -nGetChar() );
1659 napoly res;
1660
1661 if (currRing->minpoly!=NULL)
1662 {
1663   CanonicalForm mipo=convSingTrFactoryP(((lnumber)currRing->minpoly)->z);
1664   Variable a=rootOf(mipo);
1665   CanonicalForm F( convSingAFactoryA( f,a ) ), G( convSingAFactoryA( g,a ) );
1666   CanonicalForm GCD;
1667
1668   TIMING_START( algLcmTimer );
1669   // calculate gcd
1670#ifdef FACTORY_GCD_TEST
1671   GCD = CFPrimitiveGcdUtil::gcd( F, G );
1672#else
1673   GCD = gcd( F, G );
1674#endif
1675   TIMING_END( algLcmTimer );
1676
1677   FACTORY_CFAOUT( "d", GCD );
1678   FACTORY_GCDSTAT( "alcm:", F, G, GCD );
1679
1680   // calculate lcm
1681   res= convFactoryASingA( (F/GCD)*G );
1682 }
1683 else
1684 {
1685   CanonicalForm F( convSingTrFactoryP( f ) ), G( convSingTrFactoryP( g ) );
1686   CanonicalForm GCD;
1687   TIMING_START( algLcmTimer );
1688   // calculate gcd
1689#ifdef FACTORY_GCD_TEST
1690   GCD = CFPrimitiveGcdUtil::gcd( F, G );
1691#else
1692   GCD = gcd( F, G );
1693#endif
1694   TIMING_END( algLcmTimer );
1695
1696   FACTORY_CFTROUT( "d", GCD );
1697   FACTORY_GCDSTAT( "alcm:", F, G, GCD );
1698
1699   // calculate lcm
1700   res= convFactoryPSingTr( (F/GCD)*G );
1701 }
1702
1703 Off(SW_RATIONAL);
1704 return res;
1705}
1706
1707void singclap_algdividecontent ( napoly f, napoly g, napoly &ff, napoly &gg )
1708{
1709 FACTORY_ALGOUT( "f", f );
1710 FACTORY_ALGOUT( "g", g );
1711
1712 // over Q(a) / Fp(a)
1713 if (nGetChar()==1) setCharacteristic( 0 );
1714 else               setCharacteristic( -nGetChar() );
1715 ff=gg=NULL;
1716 On(SW_RATIONAL);
1717
1718 if (currRing->minpoly!=NULL)
1719 {
1720   CanonicalForm mipo=convSingTrFactoryP(((lnumber)currRing->minpoly)->z);
1721   Variable a=rootOf(mipo);
1722   CanonicalForm F( convSingAFactoryA( f,a ) ), G( convSingAFactoryA( g,a ) );
1723   CanonicalForm GCD;
1724
1725   TIMING_START( algContentTimer );
1726#ifdef FACTORY_GCD_TEST
1727   GCD=CFPrimitiveGcdUtil::gcd( F, G );
1728#else
1729   GCD=gcd( F, G );
1730#endif
1731   TIMING_END( algContentTimer );
1732
1733   FACTORY_CFAOUT( "d", GCD );
1734   FACTORY_GCDSTAT( "acnt:", F, G, GCD );
1735
1736   if ((GCD!=1) && (GCD!=0))
1737   {
1738     ff= convFactoryASingA( F/ GCD );
1739     gg= convFactoryASingA( G/ GCD );
1740   }
1741 }
1742 else
1743 {
1744   CanonicalForm F( convSingTrFactoryP( f ) ), G( convSingTrFactoryP( g ) );
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_CFTROUT( "d", GCD );
1756   FACTORY_GCDSTAT( "acnt:", F, G, GCD );
1757
1758   if ((GCD!=1) && (GCD!=0))
1759   {
1760     ff= convFactoryPSingTr( F/ GCD );
1761     gg= convFactoryPSingTr( G/ GCD );
1762   }
1763 }
1764
1765 Off(SW_RATIONAL);
1766}
1767
1768#if 0
1769lists singclap_chineseRemainder(lists x, lists q)
1770{
1771  //assume(x->nr == q->nr);
1772  //assume(x->nr >= 0);
1773  int n=x->nr+1;
1774  if ((x->nr<0) || (x->nr!=q->nr))
1775  {
1776    WerrorS("list are empty or not of equal length");
1777    return NULL;
1778  }
1779  lists res=(lists)omAlloc0Bin(slists_bin);
1780  CFArray X(1,n), Q(1,n);
1781  int i;
1782  for(i=0; i<n; i++)
1783  {
1784    if (x->m[i-1].Typ()==INT_CMD)
1785    {
1786      X[i]=(int)x->m[i-1].Data();
1787    }
1788    else if (x->m[i-1].Typ()==NUMBER_CMD)
1789    {
1790      number N=(number)x->m[i-1].Data();
1791      X[i]=convSingNFactoryN(N);
1792    }
1793    else
1794    {
1795      WerrorS("illegal type in chineseRemainder");
1796      omFreeBin(res,slists_bin);
1797      return NULL;
1798    }
1799    if (q->m[i-1].Typ()==INT_CMD)
1800    {
1801      Q[i]=(int)q->m[i-1].Data();
1802    }
1803    else if (q->m[i-1].Typ()==NUMBER_CMD)
1804    {
1805      number N=(number)x->m[i-1].Data();
1806      Q[i]=convSingNFactoryN(N);
1807    }
1808    else
1809    {
1810      WerrorS("illegal type in chineseRemainder");
1811      omFreeBin(res,slists_bin);
1812      return NULL;
1813    }
1814  }
1815  CanonicalForm r, prod;
1816  chineseRemainder( X, Q, r, prod );
1817  res->Init(2);
1818  res->m[0].rtyp=NUMBER_CMD;
1819  res->m[1].rtyp=NUMBER_CMD;
1820  res->m[0].data=(char *)convFactoryNSingN( r );
1821  res->m[1].data=(char *)convFactoryNSingN( prod );
1822  return res;
1823}
1824#endif
1825#endif
Note: See TracBrowser for help on using the repository browser.