source: git/Singular/clapsing.cc @ 54035b

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