source: git/Singular/clapsing.cc @ 07f20e1

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