source: git/kernel/clapsing.cc @ 8fc55e2

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