source: git/kernel/clapsing.cc @ 35aab3

spielwiese
Last change on this file since 35aab3 was 35aab3, checked in by Hans Schönemann <hannes@…>, 20 years ago
This commit was generated by cvs2svn to compensate for changes in r6879, which included commits to RCS files with non-trunk default branches. git-svn-id: file:///usr/local/Singular/svn/trunk@6880 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 33.3 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2/****************************************
3*  Computer Algebra System SINGULAR     *
4****************************************/
5// $Id: clapsing.cc,v 1.1.1.1 2003-10-06 12:15:50 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( F, a );
784        }
785#endif
786      }
787    }
788    else
789    {
790      CanonicalForm F( convSingTrPClapP( f ) );
791      if ((rField_is_Q_a())&&(currRing->minpoly!=NULL))
792      {
793        WarnS("factorization may be incomplete");
794        L = factorize( F );
795      }
796      else /* Fp(a) */
797      {
798#ifdef HAVE_LIBFAC_P
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            res->m[0]=pMult(res->m[0],res->m[i]);
872            res->m[i]=NULL;
873            if ((v!=NULL) && ((*v)!=NULL))
874              (**v)[i]=0;
875            j++;
876          }
877        }
878      }
879      if (j>0)
880      {
881        idSkipZeroes(res);
882        if ((v!=NULL) && ((*v)!=NULL))
883        {
884          intvec *w=*v;
885          *v = new intvec( si_max(n-j,1) );
886          for (i=0,j=0;i<w->length();i++)
887          {
888            if((*w)[i]!=0)
889            {
890              (**v)[j]=(*w)[i]; j++;
891            }
892          }
893          delete w;
894        }
895      }
896      if (res->m[0]==NULL)
897      {
898        res->m[0]=pOne();
899      }
900    }
901  }
902  if (rField_is_Q_a() && (currRing->minpoly!=NULL))
903  {
904    int i=IDELEMS(res)-1;
905    int stop=1;
906    if (with_exps!=0) stop=0;
907    for(;i>=stop;i--)
908    {
909      pNorm(res->m[i]);
910    }
911    if (with_exps==0) pSetCoeff(res->m[0],old_lead_coeff);
912    else nDelete(&old_lead_coeff);
913  }
914  else
915    nDelete(&old_lead_coeff);
916notImpl:
917  if (res==NULL)
918    WerrorS( feNotImplemented );
919  if (NN!=NULL)
920  {
921    pMult_nn(f,NN);
922    nDelete(&NN);
923  }
924  if (N!=NULL)
925  {
926    nDelete(&N);
927  }
928  //PrintS("......S\n");
929  return res;
930}
931
932matrix singclap_irrCharSeries ( ideal I)
933{
934#ifdef HAVE_LIBFAC_P
935  if (idIs0(I)) return mpNew(1,1);
936
937  // for now there is only the possibility to handle polynomials over
938  // Q and Fp ...
939  matrix res=NULL;
940  int i;
941  Off(SW_RATIONAL);
942  On(SW_SYMMETRIC_FF);
943  On(SW_USE_NTL);
944  //Off(SW_USE_NTL);
945  CFList L;
946  ListCFList LL;
947  if (((nGetChar() == 0) || (nGetChar() > 1) )
948  && (currRing->parameter==NULL))
949  {
950    setCharacteristic( nGetChar() );
951    for(i=0;i<IDELEMS(I);i++)
952    {
953      poly p=I->m[i];
954      if (p!=NULL)
955      {
956        p=pCopy(p);
957        pCleardenom(p);
958        L.append(convSingPClapP(p));
959      }
960    }
961  }
962  // and over Q(a) / Fp(a)
963  else if (( nGetChar()==1 ) /* Q(a) */
964  || (nGetChar() <-1))       /* Fp(a) */
965  {
966    if (nGetChar()==1) setCharacteristic( 0 );
967    else               setCharacteristic( -nGetChar() );
968    for(i=0;i<IDELEMS(I);i++)
969    {
970      poly p=I->m[i];
971      if (p!=NULL)
972      {
973        p=pCopy(p);
974        pCleardenom(p);
975        L.append(convSingTrPClapP(p));
976      }
977    }
978  }
979  else
980  {
981    WerrorS( feNotImplemented );
982    return res;
983  }
984
985  // a very bad work-around --- FIX IT in libfac
986  // should be fixed as of 2001/6/27
987  int tries=0;
988  int m,n;
989  ListIterator<CFList> LLi;
990  loop
991  {
992    LL=IrrCharSeries(L);
993    m= LL.length(); // Anzahl Zeilen
994    n=0;
995    for ( LLi = LL; LLi.hasItem(); LLi++ )
996    {
997      n = si_max(LLi.getItem().length(),n);
998    }
999    if ((m!=0) && (n!=0)) break;
1000    tries++;
1001    if (tries>=5) break;
1002  }
1003  if ((m==0) || (n==0))
1004  {
1005    Warn("char_series returns %d x %d matrix from %d input polys (%d)",
1006      m,n,IDELEMS(I)+1,LL.length());
1007    iiWriteMatrix((matrix)I,"I",2,0);
1008    m=si_max(m,1);
1009    n=si_max(n,1);
1010  }
1011  res=mpNew(m,n);
1012  CFListIterator Li;
1013  for ( m=1, LLi = LL; LLi.hasItem(); LLi++, m++ )
1014  {
1015    for (n=1, Li = LLi.getItem(); Li.hasItem(); Li++, n++)
1016    {
1017      if ( (nGetChar() == 0) || (nGetChar() > 1) )
1018        MATELEM(res,m,n)=convClapPSingP(Li.getItem());
1019      else
1020        MATELEM(res,m,n)=convClapPSingTrP(Li.getItem());
1021    }
1022  }
1023  Off(SW_RATIONAL);
1024  return res;
1025#else
1026  return NULL;
1027#endif
1028}
1029
1030char* singclap_neworder ( ideal I)
1031{
1032#ifdef HAVE_LIBFAC_P
1033  int i;
1034  Off(SW_RATIONAL);
1035  On(SW_SYMMETRIC_FF);
1036  CFList L;
1037  if (((nGetChar() == 0) || (nGetChar() > 1) )
1038  && (currRing->parameter==NULL))
1039  {
1040    setCharacteristic( nGetChar() );
1041    for(i=0;i<IDELEMS(I);i++)
1042    {
1043      L.append(convSingPClapP(I->m[i]));
1044    }
1045  }
1046  // and over Q(a) / Fp(a)
1047  else if (( nGetChar()==1 ) /* Q(a) */
1048  || (nGetChar() <-1))       /* Fp(a) */
1049  {
1050    if (nGetChar()==1) setCharacteristic( 0 );
1051    else               setCharacteristic( -nGetChar() );
1052    for(i=0;i<IDELEMS(I);i++)
1053    {
1054      L.append(convSingTrPClapP(I->m[i]));
1055    }
1056  }
1057  else
1058  {
1059    WerrorS( feNotImplemented );
1060    return NULL;
1061  }
1062
1063  List<int> IL=neworderint(L);
1064  ListIterator<int> Li;
1065  StringSetS("");
1066  Li = IL;
1067  int offs=rPar(currRing);
1068  int* mark=(int*)omAlloc0((pVariables+offs)*sizeof(int));
1069  int cnt=pVariables+offs;
1070  loop
1071  {
1072    if(! Li.hasItem()) break;
1073    BOOLEAN done=TRUE;
1074    i=Li.getItem()-1;
1075    mark[i]=1;
1076    if (i<offs)
1077    {
1078      done=FALSE;
1079      //StringAppendS(currRing->parameter[i]);
1080    }
1081    else
1082    {
1083      StringAppendS(currRing->names[i-offs]);
1084    }
1085    Li++;
1086    cnt--;
1087    if(cnt==0) break;
1088    if (done) StringAppendS(",");
1089  }
1090  for(i=0;i<pVariables+offs;i++)
1091  {
1092    BOOLEAN done=TRUE;
1093    if(mark[i]==0)
1094    {
1095      if (i<offs)
1096      {
1097        done=FALSE;
1098        //StringAppendS(currRing->parameter[i]);
1099      }
1100      else
1101      {
1102        StringAppendS(currRing->names[i-offs]);
1103      }
1104      cnt--;
1105      if(cnt==0) break;
1106      if (done) StringAppendS(",");
1107    }
1108  }
1109  char * s=omStrDup(StringAppendS(""));
1110  if (s[strlen(s)-1]==',') s[strlen(s)-1]='\0';
1111  return s;
1112#else
1113  return NULL;
1114#endif
1115}
1116
1117BOOLEAN singclap_isSqrFree(poly f)
1118{
1119  BOOLEAN b=FALSE;
1120  Off(SW_RATIONAL);
1121  //  Q / Fp
1122  if (((nGetChar() == 0) || (nGetChar() > 1) )
1123  &&(currRing->parameter==NULL))
1124  {
1125    setCharacteristic( nGetChar() );
1126    CanonicalForm F( convSingPClapP( f ) );
1127    if((nGetChar()>1)&&(!F.isUnivariate()))
1128      goto err;
1129    b=(BOOLEAN)isSqrFree(F);
1130  }
1131  // and over Q(a) / Fp(a)
1132  else if (( nGetChar()==1 ) /* Q(a) */
1133  || (nGetChar() <-1))       /* Fp(a) */
1134  {
1135    if (nGetChar()==1) setCharacteristic( 0 );
1136    else               setCharacteristic( -nGetChar() );
1137    //if (currRing->minpoly!=NULL)
1138    //{
1139    //  CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1140    //  Variable a=rootOf(mipo);
1141    //  CanonicalForm F( convSingAPClapAP( f,a ) );
1142    //  ...
1143    //}
1144    //else
1145    {
1146      CanonicalForm F( convSingTrPClapP( f ) );
1147      b=(BOOLEAN)isSqrFree(F);
1148    }
1149    Off(SW_RATIONAL);
1150  }
1151  else
1152  {
1153err:
1154    WerrorS( feNotImplemented );
1155  }
1156  return b;
1157}
1158
1159poly singclap_det( const matrix m )
1160{
1161  int r=m->rows();
1162  if (r!=m->cols())
1163  {
1164    Werror("det of %d x %d matrix",r,m->cols());
1165    return NULL;
1166  }
1167  poly res=NULL;
1168  if (( nGetChar() == 0 || nGetChar() > 1 )
1169  && (currRing->parameter==NULL))
1170  {
1171    setCharacteristic( nGetChar() );
1172    CFMatrix M(r,r);
1173    int i,j;
1174    for(i=r;i>0;i--)
1175    {
1176      for(j=r;j>0;j--)
1177      {
1178        M(i,j)=convSingPClapP(MATELEM(m,i,j));
1179      }
1180    }
1181    res= convClapPSingP( determinant(M,r) ) ;
1182  }
1183  // and over Q(a) / Fp(a)
1184  else if (( nGetChar()==1 ) /* Q(a) */
1185  || (nGetChar() <-1))       /* Fp(a) */
1186  {
1187    if (nGetChar()==1) setCharacteristic( 0 );
1188    else               setCharacteristic( -nGetChar() );
1189    CFMatrix M(r,r);
1190    poly res;
1191    if (currRing->minpoly!=NULL)
1192    {
1193      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1194      Variable a=rootOf(mipo);
1195      int i,j;
1196      for(i=r;i>0;i--)
1197      {
1198        for(j=r;j>0;j--)
1199        {
1200          M(i,j)=convSingAPClapAP(MATELEM(m,i,j),a);
1201        }
1202      }
1203      res= convClapAPSingAP( determinant(M,r) ) ;
1204    }
1205    else
1206    {
1207      int i,j;
1208      for(i=r;i>0;i--)
1209      {
1210        for(j=r;j>0;j--)
1211        {
1212          M(i,j)=convSingTrPClapP(MATELEM(m,i,j));
1213        }
1214      }
1215      res= convClapPSingTrP( determinant(M,r) );
1216    }
1217  }
1218  else
1219    WerrorS( feNotImplemented );
1220  Off(SW_RATIONAL);
1221  return res;
1222}
1223
1224int singclap_det_i( intvec * m )
1225{
1226  setCharacteristic( 0 );
1227  CFMatrix M(m->rows(),m->cols());
1228  int i,j;
1229  for(i=m->rows();i>0;i--)
1230  {
1231    for(j=m->cols();j>0;j--)
1232    {
1233      M(i,j)=IMATELEM(*m,i,j);
1234    }
1235  }
1236  int res= convClapISingI( determinant(M,m->rows())) ;
1237  Off(SW_RATIONAL);
1238  return res;
1239}
1240napoly singclap_alglcm ( napoly f, napoly g )
1241{
1242 FACTORY_ALGOUT( "f", f );
1243 FACTORY_ALGOUT( "g", g );
1244
1245 // over Q(a) / Fp(a)
1246 if (nGetChar()==1) setCharacteristic( 0 );
1247 else               setCharacteristic( -nGetChar() );
1248 napoly res;
1249
1250 if (currRing->minpoly!=NULL)
1251 {
1252   CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1253   Variable a=rootOf(mipo);
1254   CanonicalForm F( convSingAClapA( f,a ) ), G( convSingAClapA( g,a ) );
1255   CanonicalForm GCD;
1256
1257   TIMING_START( algLcmTimer );
1258   // calculate gcd
1259#ifdef FACTORY_GCD_TEST
1260   GCD = CFPrimitiveGcdUtil::gcd( F, G );
1261#else
1262   GCD = gcd( F, G );
1263#endif
1264   TIMING_END( algLcmTimer );
1265
1266   FACTORY_CFAOUT( "d", GCD );
1267   FACTORY_GCDSTAT( "alcm:", F, G, GCD );
1268
1269   // calculate lcm
1270   res= convClapASingA( (F/GCD)*G );
1271 }
1272 else
1273 {
1274   CanonicalForm F( convSingTrClapP( f ) ), G( convSingTrClapP( g ) );
1275   CanonicalForm GCD;
1276   TIMING_START( algLcmTimer );
1277   // calculate gcd
1278#ifdef FACTORY_GCD_TEST
1279   GCD = CFPrimitiveGcdUtil::gcd( F, G );
1280#else
1281   GCD = gcd( F, G );
1282#endif
1283   TIMING_END( algLcmTimer );
1284
1285   FACTORY_CFTROUT( "d", GCD );
1286   FACTORY_GCDSTAT( "alcm:", F, G, GCD );
1287
1288   // calculate lcm
1289   res= convClapPSingTr( (F/GCD)*G );
1290 }
1291
1292 Off(SW_RATIONAL);
1293 return res;
1294}
1295
1296void singclap_algdividecontent ( napoly f, napoly g, napoly &ff, napoly &gg )
1297{
1298 FACTORY_ALGOUT( "f", f );
1299 FACTORY_ALGOUT( "g", g );
1300
1301 // over Q(a) / Fp(a)
1302 if (nGetChar()==1) setCharacteristic( 0 );
1303 else               setCharacteristic( -nGetChar() );
1304 ff=gg=NULL;
1305 On(SW_RATIONAL);
1306
1307 if (currRing->minpoly!=NULL)
1308 {
1309   CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1310   Variable a=rootOf(mipo);
1311   CanonicalForm F( convSingAClapA( f,a ) ), G( convSingAClapA( g,a ) );
1312   CanonicalForm GCD;
1313
1314   TIMING_START( algContentTimer );
1315#ifdef FACTORY_GCD_TEST
1316   GCD=CFPrimitiveGcdUtil::gcd( F, G );
1317#else
1318   GCD=gcd( F, G );
1319#endif
1320   TIMING_END( algContentTimer );
1321
1322   FACTORY_CFAOUT( "d", GCD );
1323   FACTORY_GCDSTAT( "acnt:", F, G, GCD );
1324
1325   if ((GCD!=1) && (GCD!=0))
1326   {
1327     ff= convClapASingA( F/ GCD );
1328     gg= convClapASingA( G/ GCD );
1329   }
1330 }
1331 else
1332 {
1333   CanonicalForm F( convSingTrClapP( f ) ), G( convSingTrClapP( g ) );
1334   CanonicalForm GCD;
1335
1336   TIMING_START( algContentTimer );
1337#ifdef FACTORY_GCD_TEST
1338   GCD=CFPrimitiveGcdUtil::gcd( F, G );
1339#else
1340   GCD=gcd( F, G );
1341#endif
1342   TIMING_END( algContentTimer );
1343
1344   FACTORY_CFTROUT( "d", GCD );
1345   FACTORY_GCDSTAT( "acnt:", F, G, GCD );
1346
1347   if ((GCD!=1) && (GCD!=0))
1348   {
1349     ff= convClapPSingTr( F/ GCD );
1350     gg= convClapPSingTr( G/ GCD );
1351   }
1352 }
1353
1354 Off(SW_RATIONAL);
1355}
1356
1357#if 0
1358lists singclap_chineseRemainder(lists x, lists q)
1359{
1360  //assume(x->nr == q->nr);
1361  //assume(x->nr >= 0);
1362  int n=x->nr+1;
1363  if ((x->nr<0) || (x->nr!=q->nr))
1364  {
1365    WerrorS("list are empty or not of equal length");
1366    return NULL;
1367  }
1368  lists res=(lists)omAlloc0Bin(slists_bin);
1369  CFArray X(1,n), Q(1,n);
1370  int i;
1371  for(i=0; i<n; i++)
1372  {
1373    if (x->m[i-1].Typ()==INT_CMD)
1374    {
1375      X[i]=(int)x->m[i-1].Data();
1376    }
1377    else if (x->m[i-1].Typ()==NUMBER_CMD)
1378    {
1379      number N=(number)x->m[i-1].Data();
1380      X[i]=convSingNClapN(N);
1381    }
1382    else
1383    {
1384      WerrorS("illegal type in chineseRemainder");
1385      omFreeBin(res,slists_bin);
1386      return NULL;
1387    }
1388    if (q->m[i-1].Typ()==INT_CMD)
1389    {
1390      Q[i]=(int)q->m[i-1].Data();
1391    }
1392    else if (q->m[i-1].Typ()==NUMBER_CMD)
1393    {
1394      number N=(number)x->m[i-1].Data();
1395      Q[i]=convSingNClapN(N);
1396    }
1397    else
1398    {
1399      WerrorS("illegal type in chineseRemainder");
1400      omFreeBin(res,slists_bin);
1401      return NULL;
1402    }
1403  }
1404  CanonicalForm r, prod;
1405  chineseRemainder( X, Q, r, prod );
1406  res->Init(2);
1407  res->m[0].rtyp=NUMBER_CMD;
1408  res->m[1].rtyp=NUMBER_CMD;
1409  res->m[0].data=(char *)convClapNSingN( r );
1410  res->m[1].data=(char *)convClapNSingN( prod );
1411  return res;
1412}
1413#endif
1414#endif
Note: See TracBrowser for help on using the repository browser.