source: git/kernel/clapsing.cc @ 171950

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