source: git/kernel/clapsing.cc @ 68349d

spielwiese
Last change on this file since 68349d was d8932ed, checked in by Hans Schönemann <hannes@…>, 20 years ago
*hannes: v2-0 git-svn-id: file:///usr/local/Singular/svn/trunk@7045 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 33.4 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2/****************************************
3*  Computer Algebra System SINGULAR     *
4****************************************/
5// $Id: clapsing.cc,v 1.2 2004-02-17 17:40:33 Singular Exp $
6/*
7* ABSTRACT: interface between Singular and factory
8*/
9
10#include "mod2.h"
11#include "omalloc.h"
12#ifdef HAVE_FACTORY
13#define SI_DONT_HAVE_GLOBAL_VARS
14#include "structs.h"
15#include "clapsing.h"
16#include "numbers.h"
17#include "ring.h"
18#include <factory.h>
19#include "clapconv.h"
20#ifdef HAVE_LIBFAC_P
21#include <factor.h>
22//CanonicalForm algcd(const CanonicalForm & F, const CanonicalForm & g, const CFList & as, const Varlist & order);
23CanonicalForm alg_gcd(const CanonicalForm &, const CanonicalForm &, const CFList &);
24#endif
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=convClapPSingTr( 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=convClapASingA( 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=convSingTrClapP(((lnumber)currRing->minpoly)->z); \
166    Variable a=rootOf(mipo); \
167    printPolyPattern( tag, convSingAClapA( f,a ), rPar( currRing ) ); \
168  } \
169  else \
170  { \
171    printPolyPattern( tag, convSingTrClapP( 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
195
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 pCopy(f); // g==0 => gcd=f (but do a pCleardenom)
205  if (f==NULL) return pCopy(g); // f==0 => gcd=g (but do a pCleardenom)
206
207  if (pIsConstant(f) || pIsConstant(g)) return pOne();
208
209  // for now there is only the possibility to handle polynomials over
210  // Q and Fp ...
211  if (( nGetChar() == 0 || nGetChar() > 1 )
212  && (currRing->parameter==NULL))
213  {
214    setCharacteristic( nGetChar() );
215    CanonicalForm F( convSingPClapP( f ) ), G( convSingPClapP( g ) );
216    res=convClapPSingP( gcd( F, G ) );
217    Off(SW_RATIONAL);
218  }
219  // and over Q(a) / Fp(a)
220  else if (( nGetChar()==1 ) /* Q(a) */
221  || (nGetChar() <-1))       /* Fp(a) */
222  {
223    if (nGetChar()==1) setCharacteristic( 0 );
224    else               setCharacteristic( -nGetChar() );
225    if (currRing->minpoly!=NULL)
226    {
227      if ( nGetChar()==1 ) /* Q(a) */
228      {
229      //  WerrorS( feNotImplemented );
230        CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
231        //Varlist ord;
232        //ord.append(mipo.mvar());
233        CFList as(mipo);
234        Variable a=rootOf(mipo);
235        //CanonicalForm F( convSingAPClapAP( f,a ) ), G( convSingAPClapAP( g,a ) );
236        CanonicalForm F( convSingTrPClapP(f) ), G( convSingTrPClapP(g) );
237        //res= convClapAPSingAP( algcd( F, G, as, ord) );
238        //res= convClapAPSingAP( alg_gcd( F, G, as) );
239        res= convClapAPSingAP( alg_gcd( F, G, as) );
240      }
241      else
242      {
243        CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
244        Variable a=rootOf(mipo);
245        CanonicalForm F( convSingAPClapAP( f,a ) ), G( convSingAPClapAP( g,a ) );
246        res= convClapAPSingAP( gcd( F, G ) );
247      }
248    }
249    else
250    {
251      CanonicalForm F( convSingTrPClapP( f ) ), G( convSingTrPClapP( g ) );
252      res= convClapPSingTrP( gcd( F, G ) );
253    }
254    Off(SW_RATIONAL);
255  }
256  #if 0
257  else if (( nGetChar()>1 )&&(currRing->parameter!=NULL)) /* GF(q) */
258  {
259    int p=rChar(currRing);
260    int n=2;
261    int t=p*p;
262    while (t!=nChar) { t*=p;n++; }
263    setCharacteristic(p,n,'a');
264    CanonicalForm F( convSingGFClapGF( f ) ), G( convSingGFClapGF( g ) );
265    res= convClapGFSingGF( gcd( F, G ) );
266  }
267  #endif
268  else
269    WerrorS( feNotImplemented );
270
271  pDelete(&f);
272  pDelete(&g);
273  pTest(res);
274  return res;
275}
276
277poly singclap_resultant ( poly f, poly g , poly x)
278{
279  int i=pIsPurePower(x);
280  if (i==0)
281  {
282    WerrorS("3rd argument must be a ring variable");
283    return NULL;
284  }
285  if ((f==NULL) || (g==NULL))
286    return NULL;
287  // for now there is only the possibility to handle polynomials over
288  // Q and Fp ...
289  if (( nGetChar() == 0 || nGetChar() > 1 )
290  && (currRing->parameter==NULL))
291  {
292    Variable X(i);
293    setCharacteristic( nGetChar() );
294    CanonicalForm F( convSingPClapP( f ) ), G( convSingPClapP( g ) );
295    poly res=convClapPSingP( resultant( F, G, X ) );
296    Off(SW_RATIONAL);
297    return res;
298  }
299  // and over Q(a) / Fp(a)
300  else if (( nGetChar()==1 ) /* Q(a) */
301  || (nGetChar() <-1))       /* Fp(a) */
302  {
303    if (nGetChar()==1) setCharacteristic( 0 );
304    else               setCharacteristic( -nGetChar() );
305    poly res;
306    if (currRing->minpoly!=NULL)
307    {
308      Variable X(i);
309      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
310      Variable a=rootOf(mipo);
311      CanonicalForm F( convSingAPClapAP( f,a ) ), G( convSingAPClapAP( g,a ) );
312      res= convClapAPSingAP( resultant( F, G, X ) );
313    }
314    else
315    {
316      Variable X(i+rPar(currRing));
317      CanonicalForm F( convSingTrPClapP( f ) ), G( convSingTrPClapP( g ) );
318      res= convClapPSingTrP( resultant( F, G, X ) );
319    }
320    Off(SW_RATIONAL);
321    return res;
322  }
323  else
324    WerrorS( feNotImplemented );
325  return NULL;
326}
327//poly singclap_resultant ( poly f, poly g , poly x)
328//{
329//  int i=pVar(x);
330//  if (i==0)
331//  {
332//    WerrorS("ringvar expected");
333//    return NULL;
334//  }
335//  ideal I=idInit(1,1);
336//
337//  // get the coeffs von f wrt. x:
338//  I->m[0]=pCopy(f);
339//  matrix ffi=mpCoeffs(I,i);
340//  ffi->rank=1;
341//  ffi->ncols=ffi->nrows;
342//  ffi->nrows=1;
343//  ideal fi=(ideal)ffi;
344//
345//  // get the coeffs von g wrt. x:
346//  I->m[0]=pCopy(g);
347//  matrix ggi=mpCoeffs(I,i);
348//  ggi->rank=1;
349//  ggi->ncols=ggi->nrows;
350//  ggi->nrows=1;
351//  ideal gi=(ideal)ggi;
352//
353//  // contruct the matrix:
354//  int fn=IDELEMS(fi); //= deg(f,x)+1
355//  int gn=IDELEMS(gi); //= deg(g,x)+1
356//  matrix m=mpNew(fn+gn-2,fn+gn-2);
357//  if(m==NULL)
358//  {
359//    return NULL;
360//  }
361//
362//  // enter the coeffs into m:
363//  int j;
364//  for(i=0;i<gn-1;i++)
365//  {
366//    for(j=0;j<fn;j++)
367//    {
368//      MATELEM(m,i+1,fn-j+i)=pCopy(fi->m[j]);
369//    }
370//  }
371//  for(i=0;i<fn-1;i++)
372//  {
373//    for(j=0;j<gn;j++)
374//    {
375//      MATELEM(m,gn+i,gn-j+i)=pCopy(gi->m[j]);
376//    }
377//  }
378//
379//  poly r=mpDet(m);
380//
381//  idDelete(&fi);
382//  idDelete(&gi);
383//  idDelete((ideal *)&m);
384//  return r;
385//}
386
387BOOLEAN singclap_extgcd ( poly f, poly g, poly &res, poly &pa, poly &pb )
388{
389  // for now there is only the possibility to handle univariate
390  // polynomials over
391  // Q and Fp ...
392  res=NULL;pa=NULL;pb=NULL;
393  On(SW_SYMMETRIC_FF);
394  if (( nGetChar() == 0 || nGetChar() > 1 )
395  && (currRing->parameter==NULL))
396  {
397    setCharacteristic( nGetChar() );
398    CanonicalForm F( convSingPClapP( f ) ), G( convSingPClapP( g ) );
399    if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
400    {
401      Off(SW_RATIONAL);
402      WerrorS("not univariate");
403      return TRUE;
404    }
405    CanonicalForm Fa,Gb;
406    On(SW_RATIONAL);
407    res=convClapPSingP( extgcd( F, G, Fa, Gb ) );
408    pa=convClapPSingP(Fa);
409    pb=convClapPSingP(Gb);
410    Off(SW_RATIONAL);
411  }
412  // and over Q(a) / Fp(a)
413  else if (( nGetChar()==1 ) /* Q(a) */
414  || (nGetChar() <-1))       /* Fp(a) */
415  {
416    if (nGetChar()==1) setCharacteristic( 0 );
417    else               setCharacteristic( -nGetChar() );
418    CanonicalForm Fa,Gb;
419    if (currRing->minpoly!=NULL)
420    {
421      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
422      Variable a=rootOf(mipo);
423      CanonicalForm F( convSingAPClapAP( f,a ) ), G( convSingAPClapAP( g,a ) );
424      if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
425      {
426        WerrorS("not univariate");
427        return TRUE;
428      }
429      res= convClapAPSingAP( extgcd( F, G, Fa, Gb ) );
430      pa=convClapAPSingAP(Fa);
431      pb=convClapAPSingAP(Gb);
432    }
433    else
434    {
435      CanonicalForm F( convSingTrPClapP( f ) ), G( convSingTrPClapP( g ) );
436      if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
437      {
438        Off(SW_RATIONAL);
439        WerrorS("not univariate");
440        return TRUE;
441      }
442      res= convClapPSingTrP( extgcd( F, G, Fa, Gb ) );
443      pa=convClapPSingTrP(Fa);
444      pb=convClapPSingTrP(Gb);
445    }
446    Off(SW_RATIONAL);
447  }
448  else
449  {
450    WerrorS( feNotImplemented );
451    return TRUE;
452  }
453  return FALSE;
454}
455
456poly singclap_pdivide ( poly f, poly g )
457{
458  // for now there is only the possibility to handle polynomials over
459  // Q and Fp ...
460  poly res=NULL;
461  On(SW_RATIONAL);
462  if (( nGetChar() == 0 || nGetChar() > 1 )
463  && (currRing->parameter==NULL))
464  {
465    setCharacteristic( nGetChar() );
466    CanonicalForm F( convSingPClapP( f ) ), G( convSingPClapP( g ) );
467    res = convClapPSingP( F / G );
468  }
469  // and over Q(a) / Fp(a)
470  else if (( nGetChar()==1 ) /* Q(a) */
471  || (nGetChar() <-1))       /* Fp(a) */
472  {
473    if (nGetChar()==1) setCharacteristic( 0 );
474    else               setCharacteristic( -nGetChar() );
475    if (currRing->minpoly!=NULL)
476    {
477      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
478      Variable a=rootOf(mipo);
479      CanonicalForm F( convSingAPClapAP( f,a ) ), G( convSingAPClapAP( g,a ) );
480      res= convClapAPSingAP(  F / G  );
481    }
482    else
483    {
484      CanonicalForm F( convSingTrPClapP( f ) ), G( convSingTrPClapP( g ) );
485      res= convClapPSingTrP(  F / G  );
486    }
487  }
488  else
489    WerrorS( feNotImplemented );
490  Off(SW_RATIONAL);
491  return res;
492}
493
494void singclap_divide_content ( poly f )
495{
496  if ( f==NULL )
497  {
498    return;
499  }
500  else  if ( pNext( f ) == NULL )
501  {
502    pSetCoeff( f, nInit( 1 ) );
503    return;
504  }
505  else
506  {
507    if ( nGetChar() == 1 )
508      setCharacteristic( 0 );
509    else  if ( nGetChar() == -1 )
510      return; /* not implemented for R */
511    else  if ( nGetChar() < 0 )
512      setCharacteristic( -nGetChar() );
513    else if (currRing->parameter==NULL) /* not GF(q) */
514      setCharacteristic( nGetChar() );
515    else
516      return; /* not implemented*/
517
518    CFList L;
519    CanonicalForm g, h;
520    poly p = pNext(f);
521
522    // first attemp: find 2 smallest g:
523
524    number g1=pGetCoeff(f);
525    number g2=pGetCoeff(p); // p==pNext(f);
526    pIter(p);
527    int sz1=nSize(g1);
528    int sz2=nSize(g2);
529    if (sz1>sz2)
530    {
531      number gg=g1;
532      g1=g2; g2=gg;
533      int sz=sz1;
534      sz1=sz2; sz2=sz;
535    }
536    while (p!=NULL)
537    {
538      int n_sz=nSize(pGetCoeff(p));
539      if (n_sz<sz1)
540      {
541        sz2=sz1;
542        g2=g1;
543        g1=pGetCoeff(p);
544        sz1=n_sz;
545        if (sz1<=3) break;
546      }
547      else if(n_sz<sz2)
548      {
549        sz2=n_sz;
550        g2=pGetCoeff(p);
551        sz2=n_sz;
552      }
553      pIter(p);
554    }
555    FACTORY_ALGOUT( "G", ((lnumber)g1)->z );
556    g = convSingTrClapP( ((lnumber)g1)->z );
557    g = gcd( g, convSingTrClapP( ((lnumber)g2)->z ));
558
559    // second run: gcd's
560
561    p = f;
562    TIMING_START( contentTimer );
563    while ( (p != NULL) && (g != 1)  && ( g != 0))
564    {
565      FACTORY_ALGOUT( "h", (((lnumber)pGetCoeff(p))->z) );
566      h = convSingTrClapP( ((lnumber)pGetCoeff(p))->z );
567      pIter( p );
568#ifdef FACTORY_GCD_STAT
569      // save g
570      CanonicalForm gOld = g;
571#endif
572
573#ifdef FACTORY_GCD_TEST
574      g = CFPrimitiveGcdUtil::gcd( g, h );
575#else
576      g = gcd( g, h );
577#endif
578
579      FACTORY_GCDSTAT( "gcnt:", gOld, h, g );
580      FACTORY_CFTROUT( "g", g );
581      L.append( h );
582    }
583    TIMING_END( contentTimer );
584    FACTORY_CONTSTAT( "cont:", g );
585    if (( g == 1 ) || (g == 0))
586    {
587      // pTest(f);
588      return;
589    }
590    else
591    {
592      CFListIterator i;
593      for ( i = L, p = f; i.hasItem(); i++, p=pNext(p) )
594      {
595        lnumber c=(lnumber)pGetCoeff(p);
596        napDelete(&c->z);
597        c->z=convClapPSingTr( i.getItem() / g );
598        //nTest((number)c);
599        //#ifdef LDEBUG
600        //number cn=(number)c;
601        //StringSetS(""); nWrite(nt); StringAppend(" ==> ");
602        //nWrite(cn);PrintS(StringAppend("\n"));
603        //#endif
604      }
605    }
606    // pTest(f);
607  }
608}
609
610static int primepower(int c)
611{
612  int p=1;
613  int cc=c;
614  while(cc!= rInternalChar(currRing)) { cc*=c; p++; }
615  return p;
616}
617
618ideal singclap_factorize ( poly f, intvec ** v , int with_exps)
619{
620  // with_exps: 3,1 return only true factors, no exponents
621  //            2 return true factors and exponents
622  //            0 return coeff, factors and exponents
623
624
625  ideal res=NULL;
626
627  // handle factorize(0) =========================================
628  if (f==NULL)
629  {
630    res=idInit(1,1);
631    if (with_exps!=1)
632    {
633      (*v)=new intvec(1);
634      (**v)[0]=1;
635    }
636    return res;
637  }
638  // handle factorize(mon) =========================================
639  if (pNext(f)==NULL)
640  {
641    int i=0;
642    int n=0;
643    int e;
644    for(i=pVariables;i>0;i--) if(pGetExp(f,i)!=0) n++;
645    if (with_exps==0) n++; // with coeff
646    res=idInit(si_max(n,1),1);
647    switch(with_exps)
648    {
649      case 0: // with coef & exp.
650        res->m[0]=pOne();
651        pSetCoeff(res->m[0],nCopy(pGetCoeff(f)));
652        // no break
653      case 2: // with exp.
654        (*v)=new intvec(si_max(n,1));
655        (**v)[0]=1;
656        // no break
657      case 1: ;
658      #ifdef TEST
659      default: ;
660      #endif
661    }
662    if (n==0)
663    {
664      res->m[0]=pOne();
665      // (**v)[0]=1; is already done
666      return res;
667    }
668    for(i=pVariables;i>0;i--)
669    {
670      e=pGetExp(f,i);
671      if(e!=0)
672      {
673        n--;
674        poly p=pOne();
675        pSetExp(p,i,1);
676        pSetm(p);
677        res->m[n]=p;
678        if (with_exps!=1) (**v)[n]=e;
679      }
680    }
681    return res;
682  }
683  //PrintS("S:");pWrite(f);PrintLn();
684  // use factory/liffac in general ==============================
685  Off(SW_RATIONAL);
686  On(SW_SYMMETRIC_FF);
687  On(SW_USE_NTL);
688  #ifdef HAVE_NTL
689  extern int prime_number;
690  if(rField_is_Q()) prime_number=0;
691  #endif
692  CFFList L;
693  number N=NULL;
694  number NN=NULL;
695  number old_lead_coeff=nCopy(pGetCoeff(f));
696
697  if (rField_is_Q() || rField_is_Zp())
698  {
699    setCharacteristic( nGetChar() );
700    if (nGetChar()==0) /* Q */
701    {
702      //if (f!=NULL) // already tested at start of routine
703      {
704        number n0=nCopy(pGetCoeff(f));
705        if (with_exps==0)
706          N=nCopy(n0);
707        pCleardenom(f);
708        NN=nDiv(n0,pGetCoeff(f));
709        nDelete(&n0);
710        if (with_exps==0)
711        {
712          nDelete(&N);
713          N=nCopy(NN);
714        }
715      }
716    }
717    CanonicalForm F( convSingPClapP( f ) );
718    if (nGetChar()==0) /* Q */
719    {
720      L = factorize( F );
721    }
722    else /* Fp */
723    {
724#ifdef HAVE_LIBFAC_P
725      L = Factorize( F );
726#else
727      goto notImpl;
728#endif
729    }
730  }
731  #if 0
732  else if (rField_is_GF())
733  {
734    int c=rChar(currRing);
735    setCharacteristic( c, primepower(c) );
736    CanonicalForm F( convSingGFClapGF( f ) );
737    if (F.isUnivariate())
738    {
739      L = factorize( F );
740    }
741    else
742    {
743      goto notImpl;
744    }
745  }
746  #endif
747  // and over Q(a) / Fp(a)
748  else if (rField_is_Extension())
749  {
750    if (rField_is_Q_a()) setCharacteristic( 0 );
751    else                 setCharacteristic( -nGetChar() );
752    if (currRing->minpoly!=NULL)
753    {
754      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
755      Variable a=rootOf(mipo);
756      CanonicalForm F( convSingAPClapAP( f,a ) );
757      L.insert(F);
758      if (rField_is_Zp_a() && F.isUnivariate())
759      {
760        L = factorize( F, a );
761      }
762      else
763      {
764        CanonicalForm G( convSingTrPClapP( f ) );
765#ifdef HAVE_LIBFAC_P
766        if (rField_is_Q_a())
767        {
768          CFList as(mipo);
769          L = newfactoras( G, as, 1);
770        }
771        else
772        {
773          L=Factorize(G, mipo);
774        }
775#else
776        WarnS("complete factorization only for univariate polynomials");
777        if (rField_is_Q_a()||(!F.isUnivariate()) /* Q(a) */
778        {
779          L = factorize( G );
780        }
781        else
782        {
783          L = factorize( G, a );
784        }
785#endif
786      }
787    }
788    else
789    {
790      CanonicalForm F( convSingTrPClapP( f ) );
791      if (rField_is_Q_a())
792      {
793        L = factorize( F );
794      }
795      else /* Fp(a) */
796      {
797#ifdef HAVE_LIBFAC_P
798        L = Factorize( F );
799        //L = factorize( F );
800#else
801        goto notImpl;
802#endif
803      }
804    }
805  }
806  else
807  {
808    goto notImpl;
809  }
810  {
811    // the first factor should be a constant
812    if ( ! L.getFirst().factor().inCoeffDomain() )
813      L.insert(CFFactor(1,1));
814    // convert into ideal
815    int n = L.length();
816    CFFListIterator J=L;
817    int j=0;
818    if (with_exps!=1)
819    {
820      if ((with_exps==2)&&(n>1))
821      {
822        n--;
823        J++;
824      }
825      *v = new intvec( n );
826    }
827    res = idInit( n ,1);
828    for ( ; J.hasItem(); J++, j++ )
829    {
830      if (with_exps!=1) (**v)[j] = J.getItem().exp();
831      if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
832        res->m[j] = convClapPSingP( J.getItem().factor() );
833      #if 0
834      else if (rField_is_GF())
835        res->m[j] = convClapGFSingGF( J.getItem().factor() );
836      #endif
837      else if (rField_is_Extension())     /* Q(a), Fp(a) */
838      {
839        if (currRing->minpoly==NULL)
840          res->m[j] = convClapPSingTrP( J.getItem().factor() );
841        else
842          res->m[j] = convClapAPSingAP( J.getItem().factor() );
843      }
844    }
845    if (N!=NULL)
846    {
847      pMult_nn(res->m[0],N);
848      nDelete(&N);
849      N=NULL;
850    }
851    // delete constants
852    if (res!=NULL)
853    {
854      int i=IDELEMS(res)-1;
855      int j=0;
856      for(;i>=0;i--)
857      {
858        if ((res->m[i]!=NULL)
859        && (pNext(res->m[i])==NULL)
860        && (pIsConstant(res->m[i])))
861        {
862          if (with_exps!=0)
863          {
864            pDelete(&(res->m[i]));
865            if ((v!=NULL) && ((*v)!=NULL))
866              (**v)[i]=0;
867            j++;
868          }
869          else if (i!=0)
870          {
871            while ((v!=NULL) && ((*v)!=NULL) && ((**v)[i]>1))
872            {
873              res->m[0]=pMult(res->m[0],pCopy(res->m[i]));
874              (**v)[i]--;
875            }
876            res->m[0]=pMult(res->m[0],res->m[i]);
877            res->m[i]=NULL;
878            if ((v!=NULL) && ((*v)!=NULL))
879              (**v)[i]=0;
880            j++;
881          }
882        }
883      }
884      if (j>0)
885      {
886        idSkipZeroes(res);
887        if ((v!=NULL) && ((*v)!=NULL))
888        {
889          intvec *w=*v;
890          *v = new intvec( si_max(n-j,1) );
891          for (i=0,j=0;i<w->length();i++)
892          {
893            if((*w)[i]!=0)
894            {
895              (**v)[j]=(*w)[i]; j++;
896            }
897          }
898          delete w;
899        }
900      }
901      if (res->m[0]==NULL)
902      {
903        res->m[0]=pOne();
904      }
905    }
906  }
907  if (rField_is_Q_a() && (currRing->minpoly!=NULL))
908  {
909    int i=IDELEMS(res)-1;
910    int stop=1;
911    if (with_exps!=0) stop=0;
912    for(;i>=stop;i--)
913    {
914      pNorm(res->m[i]);
915    }
916    if (with_exps==0) pSetCoeff(res->m[0],old_lead_coeff);
917    else nDelete(&old_lead_coeff);
918  }
919  else
920    nDelete(&old_lead_coeff);
921notImpl:
922  if (res==NULL)
923    WerrorS( feNotImplemented );
924  if (NN!=NULL)
925  {
926    pMult_nn(f,NN);
927    nDelete(&NN);
928  }
929  if (N!=NULL)
930  {
931    nDelete(&N);
932  }
933  //PrintS("......S\n");
934  return res;
935}
936
937matrix singclap_irrCharSeries ( ideal I)
938{
939#ifdef HAVE_LIBFAC_P
940  if (idIs0(I)) return mpNew(1,1);
941
942  // for now there is only the possibility to handle polynomials over
943  // Q and Fp ...
944  matrix res=NULL;
945  int i;
946  Off(SW_RATIONAL);
947  On(SW_SYMMETRIC_FF);
948  On(SW_USE_NTL);
949  //Off(SW_USE_NTL);
950  CFList L;
951  ListCFList LL;
952  if (((nGetChar() == 0) || (nGetChar() > 1) )
953  && (currRing->parameter==NULL))
954  {
955    setCharacteristic( nGetChar() );
956    for(i=0;i<IDELEMS(I);i++)
957    {
958      poly p=I->m[i];
959      if (p!=NULL)
960      {
961        p=pCopy(p);
962        pCleardenom(p);
963        L.append(convSingPClapP(p));
964      }
965    }
966  }
967  // and over Q(a) / Fp(a)
968  else if (( nGetChar()==1 ) /* Q(a) */
969  || (nGetChar() <-1))       /* Fp(a) */
970  {
971    if (nGetChar()==1) setCharacteristic( 0 );
972    else               setCharacteristic( -nGetChar() );
973    for(i=0;i<IDELEMS(I);i++)
974    {
975      poly p=I->m[i];
976      if (p!=NULL)
977      {
978        p=pCopy(p);
979        pCleardenom(p);
980        L.append(convSingTrPClapP(p));
981      }
982    }
983  }
984  else
985  {
986    WerrorS( feNotImplemented );
987    return res;
988  }
989
990  // a very bad work-around --- FIX IT in libfac
991  // should be fixed as of 2001/6/27
992  int tries=0;
993  int m,n;
994  ListIterator<CFList> LLi;
995  loop
996  {
997    LL=IrrCharSeries(L);
998    m= LL.length(); // Anzahl Zeilen
999    n=0;
1000    for ( LLi = LL; LLi.hasItem(); LLi++ )
1001    {
1002      n = si_max(LLi.getItem().length(),n);
1003    }
1004    if ((m!=0) && (n!=0)) break;
1005    tries++;
1006    if (tries>=5) break;
1007  }
1008  if ((m==0) || (n==0))
1009  {
1010    Warn("char_series returns %d x %d matrix from %d input polys (%d)",
1011      m,n,IDELEMS(I)+1,LL.length());
1012    iiWriteMatrix((matrix)I,"I",2,0);
1013    m=si_max(m,1);
1014    n=si_max(n,1);
1015  }
1016  res=mpNew(m,n);
1017  CFListIterator Li;
1018  for ( m=1, LLi = LL; LLi.hasItem(); LLi++, m++ )
1019  {
1020    for (n=1, Li = LLi.getItem(); Li.hasItem(); Li++, n++)
1021    {
1022      if ( (nGetChar() == 0) || (nGetChar() > 1) )
1023        MATELEM(res,m,n)=convClapPSingP(Li.getItem());
1024      else
1025        MATELEM(res,m,n)=convClapPSingTrP(Li.getItem());
1026    }
1027  }
1028  Off(SW_RATIONAL);
1029  return res;
1030#else
1031  return NULL;
1032#endif
1033}
1034
1035char* singclap_neworder ( ideal I)
1036{
1037#ifdef HAVE_LIBFAC_P
1038  int i;
1039  Off(SW_RATIONAL);
1040  On(SW_SYMMETRIC_FF);
1041  CFList L;
1042  if (((nGetChar() == 0) || (nGetChar() > 1) )
1043  && (currRing->parameter==NULL))
1044  {
1045    setCharacteristic( nGetChar() );
1046    for(i=0;i<IDELEMS(I);i++)
1047    {
1048      L.append(convSingPClapP(I->m[i]));
1049    }
1050  }
1051  // and over Q(a) / Fp(a)
1052  else if (( nGetChar()==1 ) /* Q(a) */
1053  || (nGetChar() <-1))       /* Fp(a) */
1054  {
1055    if (nGetChar()==1) setCharacteristic( 0 );
1056    else               setCharacteristic( -nGetChar() );
1057    for(i=0;i<IDELEMS(I);i++)
1058    {
1059      L.append(convSingTrPClapP(I->m[i]));
1060    }
1061  }
1062  else
1063  {
1064    WerrorS( feNotImplemented );
1065    return NULL;
1066  }
1067
1068  List<int> IL=neworderint(L);
1069  ListIterator<int> Li;
1070  StringSetS("");
1071  Li = IL;
1072  int offs=rPar(currRing);
1073  int* mark=(int*)omAlloc0((pVariables+offs)*sizeof(int));
1074  int cnt=pVariables+offs;
1075  loop
1076  {
1077    if(! Li.hasItem()) break;
1078    BOOLEAN done=TRUE;
1079    i=Li.getItem()-1;
1080    mark[i]=1;
1081    if (i<offs)
1082    {
1083      done=FALSE;
1084      //StringAppendS(currRing->parameter[i]);
1085    }
1086    else
1087    {
1088      StringAppendS(currRing->names[i-offs]);
1089    }
1090    Li++;
1091    cnt--;
1092    if(cnt==0) break;
1093    if (done) StringAppendS(",");
1094  }
1095  for(i=0;i<pVariables+offs;i++)
1096  {
1097    BOOLEAN done=TRUE;
1098    if(mark[i]==0)
1099    {
1100      if (i<offs)
1101      {
1102        done=FALSE;
1103        //StringAppendS(currRing->parameter[i]);
1104      }
1105      else
1106      {
1107        StringAppendS(currRing->names[i-offs]);
1108      }
1109      cnt--;
1110      if(cnt==0) break;
1111      if (done) StringAppendS(",");
1112    }
1113  }
1114  char * s=omStrDup(StringAppendS(""));
1115  if (s[strlen(s)-1]==',') s[strlen(s)-1]='\0';
1116  return s;
1117#else
1118  return NULL;
1119#endif
1120}
1121
1122BOOLEAN singclap_isSqrFree(poly f)
1123{
1124  BOOLEAN b=FALSE;
1125  Off(SW_RATIONAL);
1126  //  Q / Fp
1127  if (((nGetChar() == 0) || (nGetChar() > 1) )
1128  &&(currRing->parameter==NULL))
1129  {
1130    setCharacteristic( nGetChar() );
1131    CanonicalForm F( convSingPClapP( f ) );
1132    if((nGetChar()>1)&&(!F.isUnivariate()))
1133      goto err;
1134    b=(BOOLEAN)isSqrFree(F);
1135  }
1136  // and over Q(a) / Fp(a)
1137  else if (( nGetChar()==1 ) /* Q(a) */
1138  || (nGetChar() <-1))       /* Fp(a) */
1139  {
1140    if (nGetChar()==1) setCharacteristic( 0 );
1141    else               setCharacteristic( -nGetChar() );
1142    //if (currRing->minpoly!=NULL)
1143    //{
1144    //  CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1145    //  Variable a=rootOf(mipo);
1146    //  CanonicalForm F( convSingAPClapAP( f,a ) );
1147    //  ...
1148    //}
1149    //else
1150    {
1151      CanonicalForm F( convSingTrPClapP( f ) );
1152      b=(BOOLEAN)isSqrFree(F);
1153    }
1154    Off(SW_RATIONAL);
1155  }
1156  else
1157  {
1158err:
1159    WerrorS( feNotImplemented );
1160  }
1161  return b;
1162}
1163
1164poly singclap_det( const matrix m )
1165{
1166  int r=m->rows();
1167  if (r!=m->cols())
1168  {
1169    Werror("det of %d x %d matrix",r,m->cols());
1170    return NULL;
1171  }
1172  poly res=NULL;
1173  if (( nGetChar() == 0 || nGetChar() > 1 )
1174  && (currRing->parameter==NULL))
1175  {
1176    setCharacteristic( nGetChar() );
1177    CFMatrix M(r,r);
1178    int i,j;
1179    for(i=r;i>0;i--)
1180    {
1181      for(j=r;j>0;j--)
1182      {
1183        M(i,j)=convSingPClapP(MATELEM(m,i,j));
1184      }
1185    }
1186    res= convClapPSingP( determinant(M,r) ) ;
1187  }
1188  // and over Q(a) / Fp(a)
1189  else if (( nGetChar()==1 ) /* Q(a) */
1190  || (nGetChar() <-1))       /* Fp(a) */
1191  {
1192    if (nGetChar()==1) setCharacteristic( 0 );
1193    else               setCharacteristic( -nGetChar() );
1194    CFMatrix M(r,r);
1195    poly res;
1196    if (currRing->minpoly!=NULL)
1197    {
1198      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1199      Variable a=rootOf(mipo);
1200      int i,j;
1201      for(i=r;i>0;i--)
1202      {
1203        for(j=r;j>0;j--)
1204        {
1205          M(i,j)=convSingAPClapAP(MATELEM(m,i,j),a);
1206        }
1207      }
1208      res= convClapAPSingAP( determinant(M,r) ) ;
1209    }
1210    else
1211    {
1212      int i,j;
1213      for(i=r;i>0;i--)
1214      {
1215        for(j=r;j>0;j--)
1216        {
1217          M(i,j)=convSingTrPClapP(MATELEM(m,i,j));
1218        }
1219      }
1220      res= convClapPSingTrP( determinant(M,r) );
1221    }
1222  }
1223  else
1224    WerrorS( feNotImplemented );
1225  Off(SW_RATIONAL);
1226  return res;
1227}
1228
1229int singclap_det_i( intvec * m )
1230{
1231  setCharacteristic( 0 );
1232  CFMatrix M(m->rows(),m->cols());
1233  int i,j;
1234  for(i=m->rows();i>0;i--)
1235  {
1236    for(j=m->cols();j>0;j--)
1237    {
1238      M(i,j)=IMATELEM(*m,i,j);
1239    }
1240  }
1241  int res= convClapISingI( determinant(M,m->rows())) ;
1242  Off(SW_RATIONAL);
1243  return res;
1244}
1245napoly singclap_alglcm ( napoly f, napoly g )
1246{
1247 FACTORY_ALGOUT( "f", f );
1248 FACTORY_ALGOUT( "g", g );
1249
1250 // over Q(a) / Fp(a)
1251 if (nGetChar()==1) setCharacteristic( 0 );
1252 else               setCharacteristic( -nGetChar() );
1253 napoly res;
1254
1255 if (currRing->minpoly!=NULL)
1256 {
1257   CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1258   Variable a=rootOf(mipo);
1259   CanonicalForm F( convSingAClapA( f,a ) ), G( convSingAClapA( g,a ) );
1260   CanonicalForm GCD;
1261
1262   TIMING_START( algLcmTimer );
1263   // calculate gcd
1264#ifdef FACTORY_GCD_TEST
1265   GCD = CFPrimitiveGcdUtil::gcd( F, G );
1266#else
1267   GCD = gcd( F, G );
1268#endif
1269   TIMING_END( algLcmTimer );
1270
1271   FACTORY_CFAOUT( "d", GCD );
1272   FACTORY_GCDSTAT( "alcm:", F, G, GCD );
1273
1274   // calculate lcm
1275   res= convClapASingA( (F/GCD)*G );
1276 }
1277 else
1278 {
1279   CanonicalForm F( convSingTrClapP( f ) ), G( convSingTrClapP( g ) );
1280   CanonicalForm GCD;
1281   TIMING_START( algLcmTimer );
1282   // calculate gcd
1283#ifdef FACTORY_GCD_TEST
1284   GCD = CFPrimitiveGcdUtil::gcd( F, G );
1285#else
1286   GCD = gcd( F, G );
1287#endif
1288   TIMING_END( algLcmTimer );
1289
1290   FACTORY_CFTROUT( "d", GCD );
1291   FACTORY_GCDSTAT( "alcm:", F, G, GCD );
1292
1293   // calculate lcm
1294   res= convClapPSingTr( (F/GCD)*G );
1295 }
1296
1297 Off(SW_RATIONAL);
1298 return res;
1299}
1300
1301void singclap_algdividecontent ( napoly f, napoly g, napoly &ff, napoly &gg )
1302{
1303 FACTORY_ALGOUT( "f", f );
1304 FACTORY_ALGOUT( "g", g );
1305
1306 // over Q(a) / Fp(a)
1307 if (nGetChar()==1) setCharacteristic( 0 );
1308 else               setCharacteristic( -nGetChar() );
1309 ff=gg=NULL;
1310 On(SW_RATIONAL);
1311
1312 if (currRing->minpoly!=NULL)
1313 {
1314   CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1315   Variable a=rootOf(mipo);
1316   CanonicalForm F( convSingAClapA( f,a ) ), G( convSingAClapA( g,a ) );
1317   CanonicalForm GCD;
1318
1319   TIMING_START( algContentTimer );
1320#ifdef FACTORY_GCD_TEST
1321   GCD=CFPrimitiveGcdUtil::gcd( F, G );
1322#else
1323   GCD=gcd( F, G );
1324#endif
1325   TIMING_END( algContentTimer );
1326
1327   FACTORY_CFAOUT( "d", GCD );
1328   FACTORY_GCDSTAT( "acnt:", F, G, GCD );
1329
1330   if ((GCD!=1) && (GCD!=0))
1331   {
1332     ff= convClapASingA( F/ GCD );
1333     gg= convClapASingA( G/ GCD );
1334   }
1335 }
1336 else
1337 {
1338   CanonicalForm F( convSingTrClapP( f ) ), G( convSingTrClapP( g ) );
1339   CanonicalForm GCD;
1340
1341   TIMING_START( algContentTimer );
1342#ifdef FACTORY_GCD_TEST
1343   GCD=CFPrimitiveGcdUtil::gcd( F, G );
1344#else
1345   GCD=gcd( F, G );
1346#endif
1347   TIMING_END( algContentTimer );
1348
1349   FACTORY_CFTROUT( "d", GCD );
1350   FACTORY_GCDSTAT( "acnt:", F, G, GCD );
1351
1352   if ((GCD!=1) && (GCD!=0))
1353   {
1354     ff= convClapPSingTr( F/ GCD );
1355     gg= convClapPSingTr( G/ GCD );
1356   }
1357 }
1358
1359 Off(SW_RATIONAL);
1360}
1361
1362#if 0
1363lists singclap_chineseRemainder(lists x, lists q)
1364{
1365  //assume(x->nr == q->nr);
1366  //assume(x->nr >= 0);
1367  int n=x->nr+1;
1368  if ((x->nr<0) || (x->nr!=q->nr))
1369  {
1370    WerrorS("list are empty or not of equal length");
1371    return NULL;
1372  }
1373  lists res=(lists)omAlloc0Bin(slists_bin);
1374  CFArray X(1,n), Q(1,n);
1375  int i;
1376  for(i=0; i<n; i++)
1377  {
1378    if (x->m[i-1].Typ()==INT_CMD)
1379    {
1380      X[i]=(int)x->m[i-1].Data();
1381    }
1382    else if (x->m[i-1].Typ()==NUMBER_CMD)
1383    {
1384      number N=(number)x->m[i-1].Data();
1385      X[i]=convSingNClapN(N);
1386    }
1387    else
1388    {
1389      WerrorS("illegal type in chineseRemainder");
1390      omFreeBin(res,slists_bin);
1391      return NULL;
1392    }
1393    if (q->m[i-1].Typ()==INT_CMD)
1394    {
1395      Q[i]=(int)q->m[i-1].Data();
1396    }
1397    else if (q->m[i-1].Typ()==NUMBER_CMD)
1398    {
1399      number N=(number)x->m[i-1].Data();
1400      Q[i]=convSingNClapN(N);
1401    }
1402    else
1403    {
1404      WerrorS("illegal type in chineseRemainder");
1405      omFreeBin(res,slists_bin);
1406      return NULL;
1407    }
1408  }
1409  CanonicalForm r, prod;
1410  chineseRemainder( X, Q, r, prod );
1411  res->Init(2);
1412  res->m[0].rtyp=NUMBER_CMD;
1413  res->m[1].rtyp=NUMBER_CMD;
1414  res->m[0].data=(char *)convClapNSingN( r );
1415  res->m[1].data=(char *)convClapNSingN( prod );
1416  return res;
1417}
1418#endif
1419#endif
Note: See TracBrowser for help on using the repository browser.