source: git/kernel/clapsing.cc @ f60485d

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