source: git/kernel/clapsing.cc @ 56fba1

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