source: git/Singular/clapsing.cc @ 6227ad

spielwiese
Last change on this file since 6227ad was 6227ad, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes: index shift in lamat, minor fixes in fglm, clapsing(factor) TCL-fixes in grammar.y, febase.inc, misc.cc git-svn-id: file:///usr/local/Singular/svn/trunk@2388 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 25.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.39 1998-07-29 13:55:43 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// FACTORY_GCD_TEST: use new gcd instead of old one.
28// FACTORY_GCD_STAT: print statistics on polynomials.  Works only
29//   with the file `gcd_stat.cc' and `gcd_stat.h'.
30// FACTORY_GCD_TIMING: accumulate time used for gcd calculations.
31//   Time may be printed (and reset) with `system("gcdtime");'.
32//   For tis define, `timing.h' from the factory source directory
33//   has to be copied to the Singualr source directory.
34//   Note: for better readability, the macros `TIMING_START()' and
35//   `TIMING_END()' are used in any case.  However, they expand to
36//   nothing if `FACTORY_GCD_TIMING' is off.
37// FACTORY_GCD_DEBOUT: print polynomials involved in gcd calculations.
38//   The polynomials are printed by means of the macros
39//   `FACTORY_*OUT' which are defined to be empty if
40//   `FACTORY_GCD_DEBOUT' is off.
41
42#ifdef FACTORY_GCD_STAT
43#include "gcd_stat.h"
44#define FACTORY_GCDSTAT( tag, f, g, d ) \
45  printGcdStat( tag, f, g, d, rPar( currRing ) == 1 )
46#define FACTORY_CONTSTAT( tag, f ) \
47  printContStat( tag, f, rPar( currRing ) == 1 )
48#else
49#define FACTORY_GCDSTAT( tag, f, g, d )
50#define FACTORY_CONTSTAT( tag, f )
51#endif
52
53#ifdef FACTORY_GCD_TIMING
54#define TIMING
55#include "timing.h"
56TIMING_DEFINE_PRINT( contentTimer );
57TIMING_DEFINE_PRINT( algContentTimer );
58TIMING_DEFINE_PRINT( algLcmTimer );
59#else
60#define TIMING_START( timer )
61#define TIMING_END( timer )
62#endif
63
64#ifdef FACTORY_GCD_DEBOUT
65#include "longalg.h"
66#include "febase.h"
67// alg f
68#define FACTORY_ALGOUT( tag, f ) \
69  StringSetS( tag ); \
70  napWrite( f ); \
71  PrintS(StringAppend("\n"));
72// CanonicalForm f, represents transcendent extension
73#define FACTORY_CFTROUT( tag, f ) \
74  { \
75    alg F=convClapPSingTr( f ); \
76    StringSetS( tag ); \
77    napWrite( F ); \
78    PrintS(StringAppend("\n")); \
79    napDelete( &F ); \
80  }
81// CanonicalForm f, represents algebraic extension
82#define FACTORY_CFAOUT( tag, f ) \
83  { \
84    alg F=convClapASingA( f ); \
85    StringSetS( tag ); \
86    napWrite( F ); \
87    PrintS(StringAppend("\n")); \
88    napDelete( &F ); \
89  }
90#else
91#define FACTORY_ALGOUT( tag, f )
92#define FACTORY_CFTROUT( tag, f )
93#define FACTORY_CFAOUT( tag, f )
94#endif
95
96poly singclap_gcd ( poly f, poly g )
97{
98  poly res=NULL;
99
100  if (f) pCleardenom(f);
101  if (g) pCleardenom(g);
102
103  // for now there is only the possibility to handle polynomials over
104  // Q and Fp ...
105  if (( nGetChar() == 0 || nGetChar() > 1 )
106  && (currRing->parameter==NULL))
107  {
108    setCharacteristic( nGetChar() );
109    CanonicalForm F( convSingPClapP( f ) ), G( convSingPClapP( g ) );
110    res=convClapPSingP( gcd( F, G ) );
111    Off(SW_RATIONAL);
112  }
113  // and over Q(a) / Fp(a)
114  else if (( nGetChar()==1 ) /* Q(a) */
115  || (nGetChar() <-1))       /* Fp(a) */
116  {
117    if (nGetChar()==1) setCharacteristic( 0 );
118    else               setCharacteristic( -nGetChar() );
119    if (currRing->minpoly!=NULL)
120    {
121      if ( nGetChar()==1 ) /* Q(a) */
122      {
123        WerrorS( feNotImplemented );
124      }
125      else
126      {
127        CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
128        Variable a=rootOf(mipo);
129        CanonicalForm F( convSingAPClapAP( f,a ) ), G( convSingAPClapAP( g,a ) );
130        res= convClapAPSingAP( gcd( F, G ) );
131      }
132    }
133    else
134    {
135      CanonicalForm F( convSingTrPClapP( f ) ), G( convSingTrPClapP( g ) );
136      res= convClapPSingTrP( gcd( F, G ) );
137    }
138    Off(SW_RATIONAL);
139  }
140  #if 0
141  else if (( nGetChar()>1 )&&(currRing->parameter!=NULL)) /* GF(q) */
142  {
143    int p=rChar(currRing);
144    int n=2;
145    int t=p*p;
146    while (t!=nChar) { t*=p;n++; }
147    setCharacteristic(p,n,'a');
148    CanonicalForm F( convSingGFClapGF( f ) ), G( convSingGFClapGF( g ) );
149    res= convClapGFSingGF( gcd( F, G ) );
150  }
151  #endif
152  else
153    WerrorS( feNotImplemented );
154
155  pDelete(&f);
156  pDelete(&g);
157  return res;
158}
159
160poly singclap_resultant ( poly f, poly g , poly x)
161{
162  int i=pIsPurePower(x);
163  if (i==0)
164  {
165    WerrorS("3rd argument must be a ring variable");
166    return NULL;
167  }
168  // for now there is only the possibility to handle polynomials over
169  // Q and Fp ...
170  if (( nGetChar() == 0 || nGetChar() > 1 )
171  && (currRing->parameter==NULL))
172  {
173    Variable X(i);
174    setCharacteristic( nGetChar() );
175    CanonicalForm F( convSingPClapP( f ) ), G( convSingPClapP( g ) );
176    poly res=convClapPSingP( resultant( F, G, X ) );
177    Off(SW_RATIONAL);
178    return res;
179  }
180  // and over Q(a) / Fp(a)
181  else if (( nGetChar()==1 ) /* Q(a) */
182  || (nGetChar() <-1))       /* Fp(a) */
183  {
184    if (nGetChar()==1) setCharacteristic( 0 );
185    else               setCharacteristic( -nGetChar() );
186    poly res;
187    if (currRing->minpoly!=NULL)
188    {
189      Variable X(i);
190      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
191      Variable a=rootOf(mipo);
192      CanonicalForm F( convSingAPClapAP( f,a ) ), G( convSingAPClapAP( g,a ) );
193      res= convClapAPSingAP( resultant( F, G, X ) );
194    }
195    else
196    {
197      Variable X(i+rPar(currRing));
198      CanonicalForm F( convSingTrPClapP( f ) ), G( convSingTrPClapP( g ) );
199      res= convClapPSingTrP( resultant( F, G, X ) );
200    }
201    Off(SW_RATIONAL);
202    return res;
203  }
204  else
205    WerrorS( feNotImplemented );
206  return NULL;
207}
208//poly singclap_resultant ( poly f, poly g , poly x)
209//{
210//  int i=pVar(x);
211//  if (i==0)
212//  {
213//    WerrorS("ringvar expected");
214//    return NULL;
215//  }
216//  ideal I=idInit(1,1);
217//
218//  // get the coeffs von f wrt. x:
219//  I->m[0]=pCopy(f);
220//  matrix ffi=mpCoeffs(I,i);
221//  ffi->rank=1;
222//  ffi->ncols=ffi->nrows;
223//  ffi->nrows=1;
224//  ideal fi=(ideal)ffi;
225//
226//  // get the coeffs von g wrt. x:
227//  I->m[0]=pCopy(g);
228//  matrix ggi=mpCoeffs(I,i);
229//  ggi->rank=1;
230//  ggi->ncols=ggi->nrows;
231//  ggi->nrows=1;
232//  ideal gi=(ideal)ggi;
233//
234//  // contruct the matrix:
235//  int fn=IDELEMS(fi); //= deg(f,x)+1
236//  int gn=IDELEMS(gi); //= deg(g,x)+1
237//  matrix m=mpNew(fn+gn-2,fn+gn-2);
238//  if(m==NULL)
239//  {
240//    return NULL;
241//  }
242//
243//  // enter the coeffs into m:
244//  int j;
245//  for(i=0;i<gn-1;i++)
246//  {
247//    for(j=0;j<fn;j++)
248//    {
249//      MATELEM(m,i+1,fn-j+i)=pCopy(fi->m[j]);
250//    }
251//  }
252//  for(i=0;i<fn-1;i++)
253//  {
254//    for(j=0;j<gn;j++)
255//    {
256//      MATELEM(m,gn+i,gn-j+i)=pCopy(gi->m[j]);
257//    }
258//  }
259//
260//  poly r=mpDet(m);
261//
262//  idDelete(&fi);
263//  idDelete(&gi);
264//  idDelete((ideal *)&m);
265//  return r;
266//}
267
268lists singclap_extgcd ( poly f, poly g )
269{
270  // for now there is only the possibility to handle univariate
271  // polynomials over
272  // Q and Fp ...
273  poly res=NULL,pa=NULL,pb=NULL;
274  On(SW_SYMMETRIC_FF);
275  if (( nGetChar() == 0 || nGetChar() > 1 )
276  && (currRing->parameter==NULL))
277  {
278    setCharacteristic( nGetChar() );
279    CanonicalForm F( convSingPClapP( f ) ), G( convSingPClapP( g ) );
280    if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
281    {
282      Off(SW_RATIONAL);
283      WerrorS("not univariate");
284      return NULL;
285    }
286    CanonicalForm Fa,Gb;
287    On(SW_RATIONAL);
288    res=convClapPSingP( extgcd( F, G, Fa, Gb ) );
289    pa=convClapPSingP(Fa);
290    pb=convClapPSingP(Gb);
291    Off(SW_RATIONAL);
292  }
293  // and over Q(a) / Fp(a)
294  else if (( nGetChar()==1 ) /* Q(a) */
295  || (nGetChar() <-1))       /* Fp(a) */
296  {
297    if (nGetChar()==1) setCharacteristic( 0 );
298    else               setCharacteristic( -nGetChar() );
299    CanonicalForm Fa,Gb;
300    if (currRing->minpoly!=NULL)
301    {
302      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
303      Variable a=rootOf(mipo);
304      CanonicalForm F( convSingAPClapAP( f,a ) ), G( convSingAPClapAP( g,a ) );
305      if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
306      {
307        WerrorS("not univariate");
308        return NULL;
309      }
310      res= convClapAPSingAP( extgcd( F, G, Fa, Gb ) );
311      pa=convClapAPSingAP(Fa);
312      pb=convClapAPSingAP(Gb);
313    }
314    else
315    {
316      CanonicalForm F( convSingTrPClapP( f ) ), G( convSingTrPClapP( g ) );
317      if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
318      {
319        Off(SW_RATIONAL);
320        WerrorS("not univariate");
321        return NULL;
322      }
323      res= convClapPSingTrP( extgcd( F, G, Fa, Gb ) );
324      pa=convClapPSingTrP(Fa);
325      pb=convClapPSingTrP(Gb);
326    }
327    Off(SW_RATIONAL);
328  }
329  else
330  {
331    WerrorS( feNotImplemented );
332    return NULL;
333  }
334  lists L=(lists)Alloc(sizeof(slists));
335  L->Init(3);
336  L->m[0].rtyp=POLY_CMD;
337  L->m[0].data=(void *)res;
338  L->m[1].rtyp=POLY_CMD;
339  L->m[1].data=(void *)pa;
340  L->m[2].rtyp=POLY_CMD;
341  L->m[2].data=(void *)pb;
342  return L;
343}
344
345poly singclap_pdivide ( poly f, poly g )
346{
347  // for now there is only the possibility to handle polynomials over
348  // Q and Fp ...
349  poly res=NULL;
350  On(SW_RATIONAL);
351  if (( nGetChar() == 0 || nGetChar() > 1 )
352  && (currRing->parameter==NULL))
353  {
354    setCharacteristic( nGetChar() );
355    CanonicalForm F( convSingPClapP( f ) ), G( convSingPClapP( g ) );
356    res = convClapPSingP( F / G );
357  }
358  // and over Q(a) / Fp(a)
359  else if (( nGetChar()==1 ) /* Q(a) */
360  || (nGetChar() <-1))       /* Fp(a) */
361  {
362    if (nGetChar()==1) setCharacteristic( 0 );
363    else               setCharacteristic( -nGetChar() );
364    if (currRing->minpoly!=NULL)
365    {
366      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
367      Variable a=rootOf(mipo);
368      CanonicalForm F( convSingAPClapAP( f,a ) ), G( convSingAPClapAP( g,a ) );
369      res= convClapAPSingAP(  F / G  );
370    }
371    else
372    {
373      CanonicalForm F( convSingTrPClapP( f ) ), G( convSingTrPClapP( g ) );
374      res= convClapPSingTrP(  F / G  );
375    }
376  }
377  else
378    WerrorS( feNotImplemented );
379  Off(SW_RATIONAL);
380  return res;
381}
382
383void singclap_divide_content ( poly f )
384{
385  if ( nGetChar() == 1 )
386    setCharacteristic( 0 );
387  else  if ( nGetChar() == -1 )
388    return; /* not implemented for R */
389  else  if ( nGetChar() < 0 )
390    setCharacteristic( -nGetChar() );
391  else if (currRing->parameter==NULL) /* not GF(q) */
392    setCharacteristic( nGetChar() );
393  else
394    return; /* not implemented*/
395  if ( f==NULL )
396  {
397    return;
398  }
399  else  if ( pNext( f ) == NULL )
400  {
401    pSetCoeff( f, nInit( 1 ) );
402    return;
403  }
404  else
405  {
406    CFList L;
407    CanonicalForm g, h;
408    poly p = pNext(f);
409    nTest(pGetCoeff(f));
410    FACTORY_ALGOUT( "G = ", (((lnumber)pGetCoeff(f))->z) );
411    g = convSingTrClapP( ((lnumber)pGetCoeff(f))->z );
412    L.append( g );
413    TIMING_START( contentTimer );
414    while ( (p != NULL) && (g != 1) )
415    {
416      nTest(pGetCoeff(p));
417      FACTORY_ALGOUT( "h = ", (((lnumber)pGetCoeff(p))->z) );
418      h = convSingTrClapP( ((lnumber)pGetCoeff(p))->z );
419      p = pNext( p );
420#ifdef FACTORY_GCD_STAT
421      // save g
422      CanonicalForm gOld = g;
423#endif
424
425#ifdef FACTORY_GCD_TEST
426      g = CFPrimitiveGcdUtil::gcd( g, h );
427#else
428      g = gcd( g, h );
429#endif
430
431      FACTORY_GCDSTAT( "cont", gOld, h, g );
432      FACTORY_CFTROUT( "g = ", g );
433      L.append( h );
434    }
435    TIMING_END( contentTimer );
436    FACTORY_CFTROUT( "C = ", g );
437    FACTORY_CONTSTAT( "CONT", g );
438    if ( g == 1 )
439    {
440      pTest(f);
441      return;
442    }
443    #ifdef LDEBUG
444    else if ( g == 0 )
445    {
446      pTest(f);
447      pWrite(f);
448      PrintS("=> gcd 0 in divide_content\n");
449      return;
450    }
451    #endif
452    else
453    {
454      CFListIterator i;
455      for ( i = L, p = f; i.hasItem(); i++, p=pNext(p) )
456      {
457        lnumber c=(lnumber)pGetCoeff(p);
458        napDelete(&c->z);
459        #ifdef LDEBUG
460        number nt=(number)Alloc0(sizeof(rnumber));
461        lnumber nnt=(lnumber)nt;
462        nnt->z=convClapPSingTr( i.getItem());
463        nTest(nt);
464        #endif
465        c->z=convClapPSingTr( i.getItem() / g );
466        nTest((number)c);
467        //#ifdef LDEBUG
468        //number cn=(number)c;
469        //StringSet(""); nWrite(nt); StringAppend(" ==> ");
470        //nWrite(cn);PrintS(StringAppend("\n"));
471        //#endif
472      }
473    }
474    pTest(f);
475  }
476}
477
478ideal singclap_factorize ( poly f, intvec ** v , int with_exps)
479{
480  // with_exps: 1 return only true factors
481  //            2 return true factors and exponents
482  //            0 return factors and exponents
483
484  ideal res=NULL;
485  if (f==NULL)
486  {
487    res=idInit(1,1);
488    if (with_exps!=1)
489    {
490      (*v)=new intvec(1);
491      (*v)[1]=1;
492    }
493    return res;
494  }
495  Off(SW_RATIONAL);
496  On(SW_SYMMETRIC_FF);
497  CFFList L;
498  number N=NULL;
499
500  if (( (nGetChar() == 0) || (nGetChar() > 1) )
501  && (currRing->parameter==NULL))
502  {
503    setCharacteristic( nGetChar() );
504    if (nGetChar()==0) /* Q */
505    {
506      if (f!=NULL)
507      {
508        if (with_exps==0)
509          N=nCopy(pGetCoeff(f));
510        pCleardenom(f);
511        if (with_exps==0)
512        {
513          number nn=nDiv(N,pGetCoeff(f));
514          nDelete(&N);
515          N=nn;
516        }
517      }
518    }
519    CanonicalForm F( convSingPClapP( f ) );
520    if (nGetChar()==0) /* Q */
521    {
522      L = factorize( F );
523    }
524    else /* Fp */
525    {
526#ifdef HAVE_LIBFAC_P
527      L = Factorize( F );
528#else
529      goto notImpl;
530#endif
531    }
532  }
533  // and over Q(a) / Fp(a)
534  else if (( nGetChar()==1 ) /* Q(a) */
535  || (nGetChar() <-1))       /* Fp(a) */
536  {
537    if (nGetChar()==1) setCharacteristic( 0 );
538    else               setCharacteristic( -nGetChar() );
539    if ((currRing->minpoly!=NULL)
540    && (nGetChar()<(-1)))
541    {
542      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
543      Variable a=rootOf(mipo);
544      CanonicalForm F( convSingAPClapAP( f,a ) );
545      if (F.isUnivariate())
546      {
547        L = factorize( F, a );
548      }
549      else
550      {
551        CanonicalForm G( convSingTrPClapP( f ) );
552        L = factorize( G );
553      }
554    }
555    else
556    {
557      CanonicalForm F( convSingTrPClapP( f ) );
558      if (nGetChar()==1) /* Q(a) */
559      {
560        L = factorize( F );
561      }
562      else /* Fp(a) */
563      {
564#ifdef HAVE_LIBFAC_P
565        L = Factorize( F );
566#else
567        goto notImpl;
568#endif
569      }
570    }
571  }
572  else
573  {
574    goto notImpl;
575  }
576  {
577    // the first factor should be a constant
578    if ( ! L.getFirst().factor().inCoeffDomain() )
579      L.insert(CFFactor(1,1));
580    // convert into ideal
581    int n = L.length();
582    CFFListIterator J=L;
583    int j=0;
584    if (with_exps!=1)
585    {
586      if ((with_exps==2)&&(n>1))
587      {
588        n--;
589        J++;
590      }
591      *v = new intvec( n );
592    }
593    res = idInit( n ,1);
594    for ( ; J.hasItem(); J++, j++ )
595    {
596      if (with_exps!=1) (**v)[j] = J.getItem().exp();
597      if ((nGetChar()==0)||(nGetChar()>1))           /* Q, Fp */
598        res->m[j] = convClapPSingP( J.getItem().factor() );
599      else if ((nGetChar()==1)||(nGetChar()<-1))     /* Q(a), Fp(a) */
600      {
601        if (currRing->minpoly==NULL)
602          res->m[j] = convClapPSingTrP( J.getItem().factor() );
603        else
604          res->m[j] = convClapAPSingAP( J.getItem().factor() );
605      }
606    }
607    if (N!=NULL)
608    {
609      pMultN(res->m[0],N);
610      nDelete(&N);
611    }
612    // delete constants
613    if ((with_exps!=0) && (res!=NULL))
614    {
615      int i=IDELEMS(res)-1;
616      int j=0;
617      for(;i>=0;i--)
618      {
619        if (pIsConstant(res->m[i]))
620        {
621          pDelete(&(res->m[i]));
622          if ((v!=NULL) && ((*v)!=NULL))
623            (**v)[i]=0;
624          j++;
625        }
626      }
627      if (j>0)
628      {
629        idSkipZeroes(res);
630        if ((v!=NULL) && ((*v)!=NULL))
631        {
632          intvec *w=*v;
633          *v = new intvec( max(n-j,1) );
634          for (i=0,j=0;i<w->length();i++)
635          {
636            if((*w)[i]!=0)
637            {
638              (**v)[j]=(*w)[i]; j++;
639            }
640          }
641          delete w;
642        }
643      }
644      if (res->m[0]==NULL)
645      {
646        res->m[0]=pOne();
647      }
648    }
649  }
650  return res;
651notImpl:
652  WerrorS( feNotImplemented );
653  return NULL;
654}
655
656matrix singclap_irrCharSeries ( ideal I)
657{
658#ifdef HAVE_LIBFAC_P
659  // for now there is only the possibility to handle polynomials over
660  // Q and Fp ...
661  matrix res=NULL;
662  int i;
663  Off(SW_RATIONAL);
664  On(SW_SYMMETRIC_FF);
665  CFList L;
666  ListCFList LL;
667  if (((nGetChar() == 0) || (nGetChar() > 1) )
668  && (currRing->parameter==NULL))
669  {
670    setCharacteristic( nGetChar() );
671    for(i=0;i<IDELEMS(I);i++)
672    {
673      L.append(convSingPClapP(I->m[i]));
674    }
675  }
676  // and over Q(a) / Fp(a)
677  else if (( nGetChar()==1 ) /* Q(a) */
678  || (nGetChar() <-1))       /* Fp(a) */
679  {
680    if (nGetChar()==1) setCharacteristic( 0 );
681    else               setCharacteristic( -nGetChar() );
682    for(i=0;i<IDELEMS(I);i++)
683    {
684      L.append(convSingTrPClapP(I->m[i]));
685    }
686  }
687  else
688  {
689    WerrorS( feNotImplemented );
690    return res;
691  }
692
693  LL=IrrCharSeries(L);
694  int m= LL.length(); // Anzahl Zeilen
695  int n=0;
696  ListIterator<CFList> LLi;
697  CFListIterator Li;
698  for ( LLi = LL; LLi.hasItem(); LLi++ )
699  {
700    n = max(LLi.getItem().length(),n);
701  }
702  res=mpNew(m,n);
703  if ((m==0) || (n==0))
704  {
705    Warn("char_series returns %d x %d matrix from %d input polys (%d)\n",m,n,IDELEMS(I)+1,LL.length());
706    iiWriteMatrix((matrix)I,"I",2,0);
707  }
708  for ( m=1, LLi = LL; LLi.hasItem(); LLi++, m++ )
709  {
710    for (n=1, Li = LLi.getItem(); Li.hasItem(); Li++, n++)
711    {
712      if ( (nGetChar() == 0) || (nGetChar() > 1) )
713        MATELEM(res,m,n)=convClapPSingP(Li.getItem());
714      else
715        MATELEM(res,m,n)=convClapPSingTrP(Li.getItem());
716    }
717  }
718  Off(SW_RATIONAL);
719  return res;
720#else
721  return NULL;
722#endif
723}
724
725char* singclap_neworder ( ideal I)
726{
727#ifdef HAVE_LIBFAC_P
728  int i;
729  Off(SW_RATIONAL);
730  On(SW_SYMMETRIC_FF);
731  CFList L;
732  if (((nGetChar() == 0) || (nGetChar() > 1) )
733  && (currRing->parameter==NULL))
734  {
735    setCharacteristic( nGetChar() );
736    for(i=0;i<IDELEMS(I);i++)
737    {
738      L.append(convSingPClapP(I->m[i]));
739    }
740  }
741  // and over Q(a) / Fp(a)
742  else if (( nGetChar()==1 ) /* Q(a) */
743  || (nGetChar() <-1))       /* Fp(a) */
744  {
745    if (nGetChar()==1) setCharacteristic( 0 );
746    else               setCharacteristic( -nGetChar() );
747    for(i=0;i<IDELEMS(I);i++)
748    {
749      L.append(convSingTrPClapP(I->m[i]));
750    }
751  }
752  else
753  {
754    WerrorS( feNotImplemented );
755    return NULL;
756  }
757
758  List<int> IL=neworderint(L);
759  ListIterator<int> Li;
760  StringSet("");
761  Li = IL;
762  int* mark=(int*)Alloc0(pVariables*sizeof(int));
763  int cnt=pVariables;
764  loop
765  {
766    i=Li.getItem()-1;
767    mark[i]=1;
768    StringAppend(currRing->names[i]);
769    Li++;
770    cnt--;
771    if(cnt==0) break;
772    StringAppend(",");
773    if(! Li.hasItem()) break;
774  }
775  for(i=0;i<pVariables;i++)
776  {
777    if(mark[i]==0)
778    {
779      StringAppend(currRing->names[i]);
780      cnt--;
781      if(cnt==0) break;
782      StringAppend(",");
783    }
784  }
785  return mstrdup(StringAppend(""));
786#else
787  return NULL;
788#endif
789}
790
791BOOLEAN singclap_isSqrFree(poly f)
792{
793  BOOLEAN b=FALSE;
794  Off(SW_RATIONAL);
795  //  Q / Fp
796  if (((nGetChar() == 0) || (nGetChar() > 1) )
797  &&(currRing->parameter==NULL))
798  {
799    setCharacteristic( nGetChar() );
800    CanonicalForm F( convSingPClapP( f ) );
801    if((nGetChar()>1)&&(!F.isUnivariate()))
802      goto err;
803    b=(BOOLEAN)isSqrFree(F);
804  }
805  // and over Q(a) / Fp(a)
806  else if (( nGetChar()==1 ) /* Q(a) */
807  || (nGetChar() <-1))       /* Fp(a) */
808  {
809    if (nGetChar()==1) setCharacteristic( 0 );
810    else               setCharacteristic( -nGetChar() );
811    //if (currRing->minpoly!=NULL)
812    //{
813    //  CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
814    //  Variable a=rootOf(mipo);
815    //  CanonicalForm F( convSingAPClapAP( f,a ) );
816    //  ...
817    //}
818    //else
819    {
820      CanonicalForm F( convSingTrPClapP( f ) );
821      b=(BOOLEAN)isSqrFree(F);
822    }
823    Off(SW_RATIONAL);
824  }
825  else
826  {
827err:
828    WerrorS( feNotImplemented );
829  }
830  return b;
831}
832
833poly singclap_det( const matrix m )
834{
835  int r=m->rows();
836  if (r!=m->cols())
837  {
838    Werror("det of %d x %d matrix",r,m->cols());
839    return NULL;
840  }
841  poly res=NULL;
842  if (( nGetChar() == 0 || nGetChar() > 1 )
843  && (currRing->parameter==NULL))
844  {
845    setCharacteristic( nGetChar() );
846    CFMatrix M(r,r);
847    int i,j;
848    for(i=r;i>0;i--)
849    {
850      for(j=r;j>0;j--)
851      {
852        M(i,j)=convSingPClapP(MATELEM(m,i,j));
853      }
854    }
855    res= convClapPSingP( determinant(M,r) ) ;
856  }
857  // and over Q(a) / Fp(a)
858  else if (( nGetChar()==1 ) /* Q(a) */
859  || (nGetChar() <-1))       /* Fp(a) */
860  {
861    if (nGetChar()==1) setCharacteristic( 0 );
862    else               setCharacteristic( -nGetChar() );
863    CFMatrix M(r,r);
864    poly res;
865    if (currRing->minpoly!=NULL)
866    {
867      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
868      Variable a=rootOf(mipo);
869      int i,j;
870      for(i=r;i>0;i--)
871      {
872        for(j=r;j>0;j--)
873        {
874          M(i,j)=convSingAPClapAP(MATELEM(m,i,j),a);
875        }
876      }
877      res= convClapAPSingAP( determinant(M,r) ) ;
878    }
879    else
880    {
881      int i,j;
882      for(i=r;i>0;i--)
883      {
884        for(j=r;j>0;j--)
885        {
886          M(i,j)=convSingTrPClapP(MATELEM(m,i,j));
887        }
888      }
889      res= convClapPSingTrP( determinant(M,r) );
890    }
891  }
892  else
893    WerrorS( feNotImplemented );
894  Off(SW_RATIONAL);
895  return res;
896}
897
898int singclap_det_i( intvec * m )
899{
900  setCharacteristic( 0 );
901  CFMatrix M(m->rows(),m->cols());
902  int i,j;
903  for(i=m->rows();i>0;i--)
904  {
905    for(j=m->cols();j>0;j--)
906    {
907      M(i,j)=IMATELEM(*m,i,j);
908    }
909  }
910  int res= convClapISingI( determinant(M,m->rows())) ;
911  Off(SW_RATIONAL);
912  return res;
913}
914/*==============================================================*/
915/* interpreter interface : */
916BOOLEAN jjGCD_P(leftv res, leftv u, leftv v)
917{
918  res->data=(void *)singclap_gcd((poly)(u->CopyD()),((poly)v->CopyD()));
919  return FALSE;
920}
921
922BOOLEAN jjFAC_P(leftv res, leftv u)
923{
924  intvec *v=NULL;
925  ideal f=singclap_factorize((poly)(u->Data()), &v, 0);
926  if (f==NULL) return TRUE;
927  lists l=(lists)Alloc(sizeof(slists));
928  l->Init(2);
929  l->m[0].rtyp=IDEAL_CMD;
930  l->m[0].data=(void *)f;
931  l->m[1].rtyp=INTVEC_CMD;
932  l->m[1].data=(void *)v;
933  res->data=(void *)l;
934  return FALSE;
935}
936
937BOOLEAN jjSQR_FREE_DEC(leftv res, leftv u,leftv dummy)
938{
939  intvec *v=NULL;
940  int sw=(int)dummy->Data();
941  ideal f=singclap_factorize((poly)(u->Data()), &v, sw);
942  switch(sw)
943  {
944    case 0:
945    case 2:
946    {
947      lists l=(lists)Alloc(sizeof(slists));
948      l->Init(2);
949      l->m[0].rtyp=IDEAL_CMD;
950      l->m[0].data=(void *)f;
951      l->m[1].rtyp=INTVEC_CMD;
952      l->m[1].data=(void *)v;
953      res->data=(void *)l;
954      res->rtyp=LIST_CMD;
955      return FALSE;
956    }
957    case 1:
958      res->data=(void *)f;
959      return f==NULL;
960  }
961  WerrorS("invalid switch");
962  return TRUE;
963}
964
965#if 0
966BOOLEAN jjIS_SQR_FREE(leftv res, leftv u)
967{
968  BOOLEAN b=singclap_factorize((poly)(u->Data()), &v, 0);
969  res->data=(void *)b;
970}
971#endif
972
973BOOLEAN jjEXTGCD_P(leftv res, leftv u, leftv v)
974{
975  res->data=singclap_extgcd((poly)u->Data(),(poly)v->Data());
976  return (res->data==NULL);
977}
978BOOLEAN jjRESULTANT(leftv res, leftv u, leftv v, leftv w)
979{
980  res->data=singclap_resultant((poly)u->Data(),(poly)v->Data(), (poly)w->Data());
981  return errorreported;
982}
983BOOLEAN jjCHARSERIES(leftv res, leftv u)
984{
985  res->data=singclap_irrCharSeries((ideal)u->Data());
986  return (res->data==NULL);
987}
988
989alg singclap_alglcm ( alg f, alg g )
990{
991 FACTORY_ALGOUT( "f = ", f );
992 FACTORY_ALGOUT( "g = ", g );
993
994 // over Q(a) / Fp(a)
995 if (nGetChar()==1) setCharacteristic( 0 );
996 else               setCharacteristic( -nGetChar() );
997 alg res;
998
999 if (currRing->minpoly!=NULL)
1000 {
1001   CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1002   Variable a=rootOf(mipo);
1003   CanonicalForm F( convSingAClapA( f,a ) ), G( convSingAClapA( g,a ) );
1004   CanonicalForm GCD;
1005
1006   TIMING_START( algLcmTimer );
1007   // calculate gcd
1008#ifdef FACTORY_GCD_TEST
1009   GCD = CFPrimitiveGcdUtil::gcd( F, G );
1010#else
1011   GCD = gcd( F, G );
1012#endif
1013   TIMING_END( algLcmTimer );
1014
1015   FACTORY_CFAOUT( "gcd = ", GCD );
1016
1017   // calculate lcm
1018   res= convClapASingA( (F/GCD)*G );
1019 }
1020 else
1021 {
1022   CanonicalForm F( convSingTrClapP( f ) ), G( convSingTrClapP( g ) );
1023   CanonicalForm GCD;
1024   TIMING_START( algLcmTimer );
1025   // calculate gcd
1026#ifdef FACTORY_GCD_TEST
1027   GCD = CFPrimitiveGcdUtil::gcd( F, G );
1028#else
1029   GCD = gcd( F, G );
1030#endif
1031   TIMING_END( algLcmTimer );
1032
1033   FACTORY_CFTROUT( "gcd = ", GCD );
1034
1035   // calculate lcm
1036   res= convClapPSingTr( (F/GCD)*G );
1037 }
1038
1039 Off(SW_RATIONAL);
1040 return res;
1041}
1042
1043void singclap_algdividecontent ( alg f, alg g, alg &ff, alg &gg )
1044{
1045 FACTORY_ALGOUT( "f = ", f );
1046 FACTORY_ALGOUT( "g = ", g );
1047
1048 // over Q(a) / Fp(a)
1049 if (nGetChar()==1) setCharacteristic( 0 );
1050 else               setCharacteristic( -nGetChar() );
1051 ff=gg=NULL;
1052
1053 if (currRing->minpoly!=NULL)
1054 {
1055   CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1056   Variable a=rootOf(mipo);
1057   CanonicalForm F( convSingAClapA( f,a ) ), G( convSingAClapA( g,a ) );
1058   CanonicalForm GCD;
1059
1060   TIMING_START( algContentTimer );
1061#ifdef FACTORY_GCD_TEST
1062   GCD=CFPrimitiveGcdUtil::gcd( F, G );
1063#else
1064   GCD=gcd( F, G );
1065#endif
1066   TIMING_END( algContentTimer );
1067
1068   FACTORY_CFAOUT( "gcd = ", GCD );
1069
1070   if (GCD!=1)
1071   {
1072     ff= convClapASingA( F/ GCD );
1073     gg= convClapASingA( G/ GCD );
1074   }
1075 }
1076 else
1077 {
1078   CanonicalForm F( convSingTrClapP( f ) ), G( convSingTrClapP( g ) );
1079   CanonicalForm GCD;
1080
1081   TIMING_START( algContentTimer );
1082#ifdef FACTORY_GCD_TEST
1083   GCD=CFPrimitiveGcdUtil::gcd( F, G );
1084#else
1085   GCD=gcd( F, G );
1086#endif
1087   TIMING_END( algContentTimer );
1088
1089   FACTORY_CFTROUT( "gcd = ", GCD );
1090
1091   if (GCD!=1)
1092   {
1093     ff= convClapPSingTr( F/ GCD );
1094     gg= convClapPSingTr( G/ GCD );
1095   }
1096 }
1097
1098 Off(SW_RATIONAL);
1099}
1100#endif
Note: See TracBrowser for help on using the repository browser.