source: git/Singular/clapsing.cc @ f64447

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