source: git/Singular/clapsing.cc @ 425e8b8

fieker-DuValspielwiese
Last change on this file since 425e8b8 was 425e8b8, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes/siebert: fixed bug_4 + examples git-svn-id: file:///usr/local/Singular/svn/trunk@2555 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 26.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.42 1998-10-09 12:29:37 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!=NULL) pCleardenom(f);
101  if (g!=NULL) 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  number NN=NULL;
500
501  if (( (nGetChar() == 0) || (nGetChar() > 1) )
502  && (currRing->parameter==NULL))
503  {
504    setCharacteristic( nGetChar() );
505    if (nGetChar()==0) /* Q */
506    {
507      if (f!=NULL)
508      {
509        number n0=nCopy(pGetCoeff(f));
510        if (with_exps==0)
511          N=nCopy(n0);
512        pCleardenom(f);
513        NN=nDiv(n0,pGetCoeff(f));
514        nDelete(&n0);
515        if (with_exps==0)
516        {
517          nDelete(&N);
518          N=nCopy(NN);
519        }
520      }
521    }
522    CanonicalForm F( convSingPClapP( f ) );
523    if (nGetChar()==0) /* Q */
524    {
525      L = factorize( F );
526    }
527    else /* Fp */
528    {
529#ifdef HAVE_LIBFAC_P
530      L = Factorize( F );
531#else
532      goto notImpl;
533#endif
534    }
535  }
536  // and over Q(a) / Fp(a)
537  else if (( nGetChar()==1 ) /* Q(a) */
538  || (nGetChar() <-1))       /* Fp(a) */
539  {
540    if (nGetChar()==1) setCharacteristic( 0 );
541    else               setCharacteristic( -nGetChar() );
542    if ((currRing->minpoly!=NULL)
543    && (nGetChar()<(-1)))
544    {
545      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
546      Variable a=rootOf(mipo);
547      CanonicalForm F( convSingAPClapAP( f,a ) );
548      if (F.isUnivariate())
549      {
550        L = factorize( F, a );
551      }
552      else
553      {
554        CanonicalForm G( convSingTrPClapP( f ) );
555        L = factorize( G );
556      }
557    }
558    else
559    {
560      CanonicalForm F( convSingTrPClapP( f ) );
561      if (nGetChar()==1) /* Q(a) */
562      {
563        L = factorize( F );
564      }
565      else /* Fp(a) */
566      {
567#ifdef HAVE_LIBFAC_P
568        L = Factorize( F );
569#else
570        goto notImpl;
571#endif
572      }
573    }
574  }
575  else
576  {
577    goto notImpl;
578  }
579  {
580    // the first factor should be a constant
581    if ( ! L.getFirst().factor().inCoeffDomain() )
582      L.insert(CFFactor(1,1));
583    // convert into ideal
584    int n = L.length();
585    CFFListIterator J=L;
586    int j=0;
587    if (with_exps!=1)
588    {
589      if ((with_exps==2)&&(n>1))
590      {
591        n--;
592        J++;
593      }
594      *v = new intvec( n );
595    }
596    res = idInit( n ,1);
597    for ( ; J.hasItem(); J++, j++ )
598    {
599      if (with_exps!=1) (**v)[j] = J.getItem().exp();
600      if ((nGetChar()==0)||(nGetChar()>1))           /* Q, Fp */
601        res->m[j] = convClapPSingP( J.getItem().factor() );
602      else if ((nGetChar()==1)||(nGetChar()<-1))     /* Q(a), Fp(a) */
603      {
604        if (currRing->minpoly==NULL)
605          res->m[j] = convClapPSingTrP( J.getItem().factor() );
606        else
607          res->m[j] = convClapAPSingAP( J.getItem().factor() );
608      }
609    }
610    if (N!=NULL)
611    {
612      pMultN(res->m[0],N);
613      nDelete(&N);
614      N=NULL;
615    }
616    // delete constants
617    if ((with_exps!=0) && (res!=NULL))
618    {
619      int i=IDELEMS(res)-1;
620      int j=0;
621      for(;i>=0;i--)
622      {
623        if (pIsConstant(res->m[i]))
624        {
625          pDelete(&(res->m[i]));
626          if ((v!=NULL) && ((*v)!=NULL))
627            (**v)[i]=0;
628          j++;
629        }
630      }
631      if (j>0)
632      {
633        idSkipZeroes(res);
634        if ((v!=NULL) && ((*v)!=NULL))
635        {
636          intvec *w=*v;
637          *v = new intvec( max(n-j,1) );
638          for (i=0,j=0;i<w->length();i++)
639          {
640            if((*w)[i]!=0)
641            {
642              (**v)[j]=(*w)[i]; j++;
643            }
644          }
645          delete w;
646        }
647      }
648      if (res->m[0]==NULL)
649      {
650        res->m[0]=pOne();
651      }
652    }
653  }
654notImpl:
655  if (res==NULL)
656    WerrorS( feNotImplemented );
657  if (NN!=NULL)
658  {
659    pMultN(f,NN);
660    nDelete(&NN);
661  }
662  if (N!=NULL)
663  {
664    nDelete(&N);
665  }
666  return res;
667}
668
669matrix singclap_irrCharSeries ( ideal I)
670{
671#ifdef HAVE_LIBFAC_P
672  // for now there is only the possibility to handle polynomials over
673  // Q and Fp ...
674  matrix res=NULL;
675  int i;
676  Off(SW_RATIONAL);
677  On(SW_SYMMETRIC_FF);
678  CFList L;
679  ListCFList LL;
680  if (((nGetChar() == 0) || (nGetChar() > 1) )
681  && (currRing->parameter==NULL))
682  {
683    setCharacteristic( nGetChar() );
684    for(i=0;i<IDELEMS(I);i++)
685    {
686      L.append(convSingPClapP(I->m[i]));
687    }
688  }
689  // and over Q(a) / Fp(a)
690  else if (( nGetChar()==1 ) /* Q(a) */
691  || (nGetChar() <-1))       /* Fp(a) */
692  {
693    if (nGetChar()==1) setCharacteristic( 0 );
694    else               setCharacteristic( -nGetChar() );
695    for(i=0;i<IDELEMS(I);i++)
696    {
697      L.append(convSingTrPClapP(I->m[i]));
698    }
699  }
700  else
701  {
702    WerrorS( feNotImplemented );
703    return res;
704  }
705
706  LL=IrrCharSeries(L);
707  int m= LL.length(); // Anzahl Zeilen
708  int n=0;
709  ListIterator<CFList> LLi;
710  CFListIterator Li;
711  for ( LLi = LL; LLi.hasItem(); LLi++ )
712  {
713    n = max(LLi.getItem().length(),n);
714  }
715  res=mpNew(m,n);
716  if ((m==0) || (n==0))
717  {
718    Warn("char_series returns %d x %d matrix from %d input polys (%d)\n",m,n,IDELEMS(I)+1,LL.length());
719    iiWriteMatrix((matrix)I,"I",2,0);
720  }
721  for ( m=1, LLi = LL; LLi.hasItem(); LLi++, m++ )
722  {
723    for (n=1, Li = LLi.getItem(); Li.hasItem(); Li++, n++)
724    {
725      if ( (nGetChar() == 0) || (nGetChar() > 1) )
726        MATELEM(res,m,n)=convClapPSingP(Li.getItem());
727      else
728        MATELEM(res,m,n)=convClapPSingTrP(Li.getItem());
729    }
730  }
731  Off(SW_RATIONAL);
732  return res;
733#else
734  return NULL;
735#endif
736}
737
738char* singclap_neworder ( ideal I)
739{
740#ifdef HAVE_LIBFAC_P
741  int i;
742  Off(SW_RATIONAL);
743  On(SW_SYMMETRIC_FF);
744  CFList L;
745  if (((nGetChar() == 0) || (nGetChar() > 1) )
746  && (currRing->parameter==NULL))
747  {
748    setCharacteristic( nGetChar() );
749    for(i=0;i<IDELEMS(I);i++)
750    {
751      L.append(convSingPClapP(I->m[i]));
752    }
753  }
754  // and over Q(a) / Fp(a)
755  else if (( nGetChar()==1 ) /* Q(a) */
756  || (nGetChar() <-1))       /* Fp(a) */
757  {
758    if (nGetChar()==1) setCharacteristic( 0 );
759    else               setCharacteristic( -nGetChar() );
760    for(i=0;i<IDELEMS(I);i++)
761    {
762      L.append(convSingTrPClapP(I->m[i]));
763    }
764  }
765  else
766  {
767    WerrorS( feNotImplemented );
768    return NULL;
769  }
770
771  List<int> IL=neworderint(L);
772  ListIterator<int> Li;
773  StringSet("");
774  Li = IL;
775  int offs=rPar(currRing);
776  int* mark=(int*)Alloc0((pVariables+offs)*sizeof(int));
777  int cnt=pVariables+offs;
778  loop
779  {
780    i=Li.getItem()-1;
781    mark[i]=1;
782    if (i<offs)
783    {
784      StringAppend(currRing->parameter[i]);
785    }
786    else
787    {
788      StringAppend(currRing->names[i-offs]);
789    }
790    Li++;
791    cnt--;
792    if(cnt==0) break;
793    StringAppend(",");
794    if(! Li.hasItem()) break;
795  }
796  for(i=0;i<pVariables+offs;i++)
797  {
798    if(mark[i]==0)
799    {
800      if (i<offs)
801      {
802        StringAppend(currRing->parameter[i]);
803      }
804      else
805      {
806        StringAppend(currRing->names[i-offs]);
807      }
808      cnt--;
809      if(cnt==0) break;
810      StringAppend(",");
811    }
812  }
813  return mstrdup(StringAppend(""));
814#else
815  return NULL;
816#endif
817}
818
819BOOLEAN singclap_isSqrFree(poly f)
820{
821  BOOLEAN b=FALSE;
822  Off(SW_RATIONAL);
823  //  Q / Fp
824  if (((nGetChar() == 0) || (nGetChar() > 1) )
825  &&(currRing->parameter==NULL))
826  {
827    setCharacteristic( nGetChar() );
828    CanonicalForm F( convSingPClapP( f ) );
829    if((nGetChar()>1)&&(!F.isUnivariate()))
830      goto err;
831    b=(BOOLEAN)isSqrFree(F);
832  }
833  // and over Q(a) / Fp(a)
834  else if (( nGetChar()==1 ) /* Q(a) */
835  || (nGetChar() <-1))       /* Fp(a) */
836  {
837    if (nGetChar()==1) setCharacteristic( 0 );
838    else               setCharacteristic( -nGetChar() );
839    //if (currRing->minpoly!=NULL)
840    //{
841    //  CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
842    //  Variable a=rootOf(mipo);
843    //  CanonicalForm F( convSingAPClapAP( f,a ) );
844    //  ...
845    //}
846    //else
847    {
848      CanonicalForm F( convSingTrPClapP( f ) );
849      b=(BOOLEAN)isSqrFree(F);
850    }
851    Off(SW_RATIONAL);
852  }
853  else
854  {
855err:
856    WerrorS( feNotImplemented );
857  }
858  return b;
859}
860
861poly singclap_det( const matrix m )
862{
863  int r=m->rows();
864  if (r!=m->cols())
865  {
866    Werror("det of %d x %d matrix",r,m->cols());
867    return NULL;
868  }
869  poly res=NULL;
870  if (( nGetChar() == 0 || nGetChar() > 1 )
871  && (currRing->parameter==NULL))
872  {
873    setCharacteristic( nGetChar() );
874    CFMatrix M(r,r);
875    int i,j;
876    for(i=r;i>0;i--)
877    {
878      for(j=r;j>0;j--)
879      {
880        M(i,j)=convSingPClapP(MATELEM(m,i,j));
881      }
882    }
883    res= convClapPSingP( determinant(M,r) ) ;
884  }
885  // and over Q(a) / Fp(a)
886  else if (( nGetChar()==1 ) /* Q(a) */
887  || (nGetChar() <-1))       /* Fp(a) */
888  {
889    if (nGetChar()==1) setCharacteristic( 0 );
890    else               setCharacteristic( -nGetChar() );
891    CFMatrix M(r,r);
892    poly res;
893    if (currRing->minpoly!=NULL)
894    {
895      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
896      Variable a=rootOf(mipo);
897      int i,j;
898      for(i=r;i>0;i--)
899      {
900        for(j=r;j>0;j--)
901        {
902          M(i,j)=convSingAPClapAP(MATELEM(m,i,j),a);
903        }
904      }
905      res= convClapAPSingAP( determinant(M,r) ) ;
906    }
907    else
908    {
909      int i,j;
910      for(i=r;i>0;i--)
911      {
912        for(j=r;j>0;j--)
913        {
914          M(i,j)=convSingTrPClapP(MATELEM(m,i,j));
915        }
916      }
917      res= convClapPSingTrP( determinant(M,r) );
918    }
919  }
920  else
921    WerrorS( feNotImplemented );
922  Off(SW_RATIONAL);
923  return res;
924}
925
926int singclap_det_i( intvec * m )
927{
928  setCharacteristic( 0 );
929  CFMatrix M(m->rows(),m->cols());
930  int i,j;
931  for(i=m->rows();i>0;i--)
932  {
933    for(j=m->cols();j>0;j--)
934    {
935      M(i,j)=IMATELEM(*m,i,j);
936    }
937  }
938  int res= convClapISingI( determinant(M,m->rows())) ;
939  Off(SW_RATIONAL);
940  return res;
941}
942/*==============================================================*/
943/* interpreter interface : */
944BOOLEAN jjGCD_P(leftv res, leftv u, leftv v)
945{
946  res->data=(void *)singclap_gcd((poly)(u->CopyD()),((poly)v->CopyD()));
947  return FALSE;
948}
949
950BOOLEAN jjFAC_P(leftv res, leftv u)
951{
952  intvec *v=NULL;
953  ideal f=singclap_factorize((poly)(u->Data()), &v, 0);
954  if (f==NULL) return TRUE;
955  lists l=(lists)Alloc(sizeof(slists));
956  l->Init(2);
957  l->m[0].rtyp=IDEAL_CMD;
958  l->m[0].data=(void *)f;
959  l->m[1].rtyp=INTVEC_CMD;
960  l->m[1].data=(void *)v;
961  res->data=(void *)l;
962  return FALSE;
963}
964
965BOOLEAN jjSQR_FREE_DEC(leftv res, leftv u,leftv dummy)
966{
967  intvec *v=NULL;
968  int sw=(int)dummy->Data();
969  ideal f=singclap_factorize((poly)(u->Data()), &v, sw);
970  if (f==NULL)
971    return TRUE;
972  switch(sw)
973  {
974    case 0:
975    case 2:
976    {
977      lists l=(lists)Alloc(sizeof(slists));
978      l->Init(2);
979      l->m[0].rtyp=IDEAL_CMD;
980      l->m[0].data=(void *)f;
981      l->m[1].rtyp=INTVEC_CMD;
982      l->m[1].data=(void *)v;
983      res->data=(void *)l;
984      res->rtyp=LIST_CMD;
985      return FALSE;
986    }
987    case 1:
988      res->data=(void *)f;
989      return FALSE;
990    case 3:
991      {
992        poly p=f->m[0];
993        int i=IDELEMS(f);
994        f->m[0]=NULL;
995        while(i>1)
996        {
997          i--;
998          p=pMult(p,f->m[i]);
999          f->m[i]=NULL;
1000        }
1001        res->data=(void *)p;
1002        res->rtyp=POLY_CMD;
1003      }
1004      return FALSE;
1005  }
1006  WerrorS("invalid switch");
1007  return TRUE;
1008}
1009
1010#if 0
1011BOOLEAN jjIS_SQR_FREE(leftv res, leftv u)
1012{
1013  BOOLEAN b=singclap_factorize((poly)(u->Data()), &v, 0);
1014  res->data=(void *)b;
1015}
1016#endif
1017
1018BOOLEAN jjEXTGCD_P(leftv res, leftv u, leftv v)
1019{
1020  res->data=singclap_extgcd((poly)u->Data(),(poly)v->Data());
1021  return (res->data==NULL);
1022}
1023BOOLEAN jjRESULTANT(leftv res, leftv u, leftv v, leftv w)
1024{
1025  res->data=singclap_resultant((poly)u->Data(),(poly)v->Data(), (poly)w->Data());
1026  return errorreported;
1027}
1028BOOLEAN jjCHARSERIES(leftv res, leftv u)
1029{
1030  res->data=singclap_irrCharSeries((ideal)u->Data());
1031  return (res->data==NULL);
1032}
1033
1034alg singclap_alglcm ( alg f, alg g )
1035{
1036 FACTORY_ALGOUT( "f = ", f );
1037 FACTORY_ALGOUT( "g = ", g );
1038
1039 // over Q(a) / Fp(a)
1040 if (nGetChar()==1) setCharacteristic( 0 );
1041 else               setCharacteristic( -nGetChar() );
1042 alg res;
1043
1044 if (currRing->minpoly!=NULL)
1045 {
1046   CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1047   Variable a=rootOf(mipo);
1048   CanonicalForm F( convSingAClapA( f,a ) ), G( convSingAClapA( g,a ) );
1049   CanonicalForm GCD;
1050
1051   TIMING_START( algLcmTimer );
1052   // calculate gcd
1053#ifdef FACTORY_GCD_TEST
1054   GCD = CFPrimitiveGcdUtil::gcd( F, G );
1055#else
1056   GCD = gcd( F, G );
1057#endif
1058   TIMING_END( algLcmTimer );
1059
1060   FACTORY_CFAOUT( "gcd = ", GCD );
1061
1062   // calculate lcm
1063   res= convClapASingA( (F/GCD)*G );
1064 }
1065 else
1066 {
1067   CanonicalForm F( convSingTrClapP( f ) ), G( convSingTrClapP( g ) );
1068   CanonicalForm GCD;
1069   TIMING_START( algLcmTimer );
1070   // calculate gcd
1071#ifdef FACTORY_GCD_TEST
1072   GCD = CFPrimitiveGcdUtil::gcd( F, G );
1073#else
1074   GCD = gcd( F, G );
1075#endif
1076   TIMING_END( algLcmTimer );
1077
1078   FACTORY_CFTROUT( "gcd = ", GCD );
1079
1080   // calculate lcm
1081   res= convClapPSingTr( (F/GCD)*G );
1082 }
1083
1084 Off(SW_RATIONAL);
1085 return res;
1086}
1087
1088void singclap_algdividecontent ( alg f, alg g, alg &ff, alg &gg )
1089{
1090 FACTORY_ALGOUT( "f = ", f );
1091 FACTORY_ALGOUT( "g = ", g );
1092
1093 // over Q(a) / Fp(a)
1094 if (nGetChar()==1) setCharacteristic( 0 );
1095 else               setCharacteristic( -nGetChar() );
1096 ff=gg=NULL;
1097
1098 if (currRing->minpoly!=NULL)
1099 {
1100   CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1101   Variable a=rootOf(mipo);
1102   CanonicalForm F( convSingAClapA( f,a ) ), G( convSingAClapA( g,a ) );
1103   CanonicalForm GCD;
1104
1105   TIMING_START( algContentTimer );
1106#ifdef FACTORY_GCD_TEST
1107   GCD=CFPrimitiveGcdUtil::gcd( F, G );
1108#else
1109   GCD=gcd( F, G );
1110#endif
1111   TIMING_END( algContentTimer );
1112
1113   FACTORY_CFAOUT( "gcd = ", GCD );
1114
1115   if (GCD!=1)
1116   {
1117     ff= convClapASingA( F/ GCD );
1118     gg= convClapASingA( G/ GCD );
1119   }
1120 }
1121 else
1122 {
1123   CanonicalForm F( convSingTrClapP( f ) ), G( convSingTrClapP( g ) );
1124   CanonicalForm GCD;
1125
1126   TIMING_START( algContentTimer );
1127#ifdef FACTORY_GCD_TEST
1128   GCD=CFPrimitiveGcdUtil::gcd( F, G );
1129#else
1130   GCD=gcd( F, G );
1131#endif
1132   TIMING_END( algContentTimer );
1133
1134   FACTORY_CFTROUT( "gcd = ", GCD );
1135
1136   if (GCD!=1)
1137   {
1138     ff= convClapPSingTr( F/ GCD );
1139     gg= convClapPSingTr( G/ GCD );
1140   }
1141 }
1142
1143 Off(SW_RATIONAL);
1144}
1145#endif
Note: See TracBrowser for help on using the repository browser.