source: git/kernel/clapsing.cc @ 3ad53dd

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