source: git/kernel/clapsing.cc @ d3e630

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