source: git/kernel/clapsing.cc @ b3af93

spielwiese
Last change on this file since b3af93 was 1629ab, checked in by Hans Schönemann <hannes@…>, 15 years ago
*hannes: cleanup longalg git-svn-id: file:///usr/local/Singular/svn/trunk@12071 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 39.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.44 2009-08-14 17:14:30 Singular Exp $
6/*
7* ABSTRACT: interface between Singular and factory
8*/
9
10//#define FACTORIZE2_DEBUG
11#include "mod2.h"
12#include "omalloc.h"
13#ifdef HAVE_FACTORY
14#define SI_DONT_HAVE_GLOBAL_VARS
15#include "structs.h"
16#include "clapsing.h"
17#include "numbers.h"
18#include "ring.h"
19#include "ffields.h"
20#include <factory.h>
21#include "clapconv.h"
22#include <factor.h>
23//CanonicalForm algcd(const CanonicalForm & F, const CanonicalForm & g, const CFList & as, const Varlist & order);
24CanonicalForm alg_gcd(const CanonicalForm &, const CanonicalForm &, const CFList &);
25#include "ring.h"
26
27void out_cf(char *s1,const CanonicalForm &f,char *s2);
28
29poly singclap_gcd_r ( poly f, poly g, const ring r )
30{
31  // assume pCleardenom is done
32  // assume f!=0, g!=0
33  poly res=NULL;
34
35  if (p_IsConstantPoly(f,r) || p_IsConstantPoly(g,r))
36  {
37    return p_One(r);
38  }
39
40  // for now there is only the possibility to handle polynomials over
41  // Q and Fp ...
42  Off(SW_RATIONAL);
43  if (rField_is_Q(r) || (rField_is_Zp(r)))
44  {
45    CanonicalForm newGCD(const CanonicalForm & A, const CanonicalForm & B);
46    setCharacteristic( n_GetChar(r) );
47    CanonicalForm F( convSingPFactoryP( f,r ) ), G( convSingPFactoryP( g, r ) );
48    //if (nGetChar() > 1 )
49    //{
50    //  res=convFactoryPSingP( newGCD( F,G ));
51    //  if (!nGreaterZero(pGetCoeff(res))) res=pNeg(res);
52    //}
53    //else
54      res=convFactoryPSingP( gcd( F, G ) , r);
55  }
56  // and over Q(a) / Fp(a)
57  else if (( n_GetChar(r)==1 ) /* Q(a) */
58  || (n_GetChar(r) <-1))       /* Fp(a) */
59  {
60    if (n_GetChar(r)==1) setCharacteristic( 0 );
61    else                 setCharacteristic( -n_GetChar(r) );
62    if (r->minpoly!=NULL)
63    {
64    #if 0
65      if (( n_GetChar(r)==1 ) /* Q(a) */ && (!isOn(SW_USE_QGCD)))
66      {
67      //  WerrorS( feNotImplemented );
68        CanonicalForm mipo=convSingPFactoryP(((lnumber)r->minpoly)->z,r->algring);
69        //Varlist ord;
70        //ord.append(mipo.mvar());
71        CFList as(mipo);
72        Variable a=rootOf(mipo);
73        //CanonicalForm F( convSingAPFactoryAP( f,a,r ) ), G( convSingAPFactoryAP( g,a,r ) );
74        CanonicalForm F( convSingTrPFactoryP(f,r) ), G( convSingTrPFactoryP(g,r) );
75        //res= convFactoryAPSingAP( algcd( F, G, as, ord) );
76        //res= convFactoryAPSingAP( alg_gcd( F, G, as) );
77        res= convFactoryAPSingAP( alg_gcd( F, G, as),r );
78      }
79      else
80      #endif
81      {
82        bool b=isOn(SW_USE_QGCD);
83        if ( rField_is_Q_a() ) On(SW_USE_QGCD);
84        CanonicalForm mipo=convSingPFactoryP(((lnumber)r->minpoly)->z,
85                                           r->algring);
86        Variable a=rootOf(mipo);
87        CanonicalForm F( convSingAPFactoryAP( f,a ) ),
88                      G( convSingAPFactoryAP( g,a ) );
89        res= convFactoryAPSingAP( gcd( F, G ) );
90        if (!b) Off(SW_USE_QGCD);
91      }
92    }
93    else
94    {
95      CanonicalForm F( convSingTrPFactoryP( f,r ) ), G( convSingTrPFactoryP( g,r ) );
96      res= convFactoryPSingTrP( gcd( F, G ),r );
97    }
98  }
99  #if 0
100  else if (( n_GetChar(r)>1 )&&(r->parameter!=NULL)) /* GF(q) */
101  {
102    int p=rChar(r);
103    int n=2;
104    int t=p*p;
105    while (t!=n_Char(r)) { t*=p;n++; }
106    setCharacteristic(p,n,'a');
107    CanonicalForm F( convSingGFFactoryGF( f,r ) ), G( convSingGFFactoryGF( g,r ) );
108    res= convFactoryGFSingGF( gcd( F, G ),r );
109  }
110  #endif
111  else
112    WerrorS( feNotImplemented );
113
114  Off(SW_RATIONAL);
115  return res;
116}
117
118poly singclap_gcd ( poly f, poly g )
119{
120  poly res=NULL;
121
122  if (f!=NULL) pCleardenom(f);
123  if (g!=NULL) pCleardenom(g);
124  else         return f; // g==0 => gcd=f (but do a pCleardenom)
125  if (f==NULL) return g; // f==0 => gcd=g (but do a pCleardenom)
126
127  res=singclap_gcd_r(f,g,currRing);
128  pDelete(&f);
129  pDelete(&g);
130  return res;
131}
132
133/*2 find the maximal exponent of var(i) in poly p*/
134int pGetExp_Var(poly p, int i)
135{
136  int m=0;
137  int mm;
138  while (p!=NULL)
139  {
140    mm=pGetExp(p,i);
141    if (mm>m) m=mm;
142    pIter(p);
143  }
144  return m;
145}
146
147poly singclap_resultant ( poly f, poly g , poly x)
148{
149  int i=pIsPurePower(x);
150  if (i==0)
151  {
152    WerrorS("3rd argument must be a ring variable");
153    return NULL;
154  }
155  if ((f==NULL) || (g==NULL))
156    return NULL;
157  // for now there is only the possibility to handle polynomials over
158  // Q and Fp ...
159  if (rField_is_Zp() || rField_is_Q())
160  {
161    Variable X(i);
162    setCharacteristic( nGetChar() );
163    CanonicalForm F( convSingPFactoryP( f ) ), G( convSingPFactoryP( g ) );
164    poly res=convFactoryPSingP( resultant( F, G, X ) );
165    Off(SW_RATIONAL);
166    return res;
167  }
168  // and over Q(a) / Fp(a)
169  else if (rField_is_Extension())
170  {
171    if (rField_is_Q_a()) setCharacteristic( 0 );
172    else               setCharacteristic( -nGetChar() );
173    poly res;
174    Variable X(i+rPar(currRing));
175    if (currRing->minpoly!=NULL)
176    {
177      //Variable X(i);
178      CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
179                                                 currRing->algring);
180      Variable a=rootOf(mipo);
181      CanonicalForm F( convSingAPFactoryAP( f,a ) ), G( convSingAPFactoryAP( g,a ) );
182      res= convFactoryAPSingAP( resultant( F, G, X ) );
183    }
184    else
185    {
186      //Variable X(i+rPar(currRing));
187      number nf,ng;
188      pCleardenom_n(f,nf);pCleardenom_n(g,ng);
189      int ef,eg;
190      ef=pGetExp_Var(f,i);
191      eg=pGetExp_Var(g,i);
192      CanonicalForm F( convSingTrPFactoryP( f ) ), G( convSingTrPFactoryP( g ) );
193      res= convFactoryPSingTrP( resultant( F, G, X ) );
194      if ((nf!=NULL)&&(!nIsOne(nf))&&(!nIsZero(nf)))
195      {
196        number n=nInvers(nf);
197        while(eg>0)
198        {
199          res=pMult_nn(res,n);
200          eg--;
201        }
202        nDelete(&n);
203      }
204      nDelete(&nf);
205      if ((ng!=NULL)&&(!nIsOne(ng))&&(!nIsZero(ng)))
206      {
207        number n=nInvers(ng);
208        while(ef>0)
209        {
210          res=pMult_nn(res,n);
211          ef--;
212        }
213        nDelete(&n);
214      }
215      nDelete(&ng);
216    }
217    Off(SW_RATIONAL);
218    return res;
219  }
220  else
221    WerrorS( feNotImplemented );
222  return NULL;
223}
224//poly singclap_resultant ( poly f, poly g , poly x)
225//{
226//  int i=pVar(x);
227//  if (i==0)
228//  {
229//    WerrorS("ringvar expected");
230//    return NULL;
231//  }
232//  ideal I=idInit(1,1);
233//
234//  // get the coeffs von f wrt. x:
235//  I->m[0]=pCopy(f);
236//  matrix ffi=mpCoeffs(I,i);
237//  ffi->rank=1;
238//  ffi->ncols=ffi->nrows;
239//  ffi->nrows=1;
240//  ideal fi=(ideal)ffi;
241//
242//  // get the coeffs von g wrt. x:
243//  I->m[0]=pCopy(g);
244//  matrix ggi=mpCoeffs(I,i);
245//  ggi->rank=1;
246//  ggi->ncols=ggi->nrows;
247//  ggi->nrows=1;
248//  ideal gi=(ideal)ggi;
249//
250//  // contruct the matrix:
251//  int fn=IDELEMS(fi); //= deg(f,x)+1
252//  int gn=IDELEMS(gi); //= deg(g,x)+1
253//  matrix m=mpNew(fn+gn-2,fn+gn-2);
254//  if(m==NULL)
255//  {
256//    return NULL;
257//  }
258//
259//  // enter the coeffs into m:
260//  int j;
261//  for(i=0;i<gn-1;i++)
262//  {
263//    for(j=0;j<fn;j++)
264//    {
265//      MATELEM(m,i+1,fn-j+i)=pCopy(fi->m[j]);
266//    }
267//  }
268//  for(i=0;i<fn-1;i++)
269//  {
270//    for(j=0;j<gn;j++)
271//    {
272//      MATELEM(m,gn+i,gn-j+i)=pCopy(gi->m[j]);
273//    }
274//  }
275//
276//  poly r=mpDet(m);
277//
278//  idDelete(&fi);
279//  idDelete(&gi);
280//  idDelete((ideal *)&m);
281//  return r;
282//}
283
284BOOLEAN singclap_extgcd ( poly f, poly g, poly &res, poly &pa, poly &pb )
285{
286  // for now there is only the possibility to handle univariate
287  // polynomials over
288  // Q and Fp ...
289  res=NULL;pa=NULL;pb=NULL;
290  On(SW_SYMMETRIC_FF);
291  if (rField_is_Zp() || rField_is_Q())
292  {
293    setCharacteristic( nGetChar() );
294    CanonicalForm F( convSingPFactoryP( f ) ), G( convSingPFactoryP( g ) );
295    CanonicalForm FpG=F+G;
296    if (!(FpG.isUnivariate()|| FpG.inCoeffDomain()))
297    //if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
298    {
299      Off(SW_RATIONAL);
300      WerrorS("not univariate");
301      return TRUE;
302    }
303    CanonicalForm Fa,Gb;
304    On(SW_RATIONAL);
305    res=convFactoryPSingP( extgcd( F, G, Fa, Gb ) );
306    pa=convFactoryPSingP(Fa);
307    pb=convFactoryPSingP(Gb);
308    Off(SW_RATIONAL);
309  }
310  // and over Q(a) / Fp(a)
311  else if (rField_is_Extension())
312  {
313    if (rField_is_Q_a()) setCharacteristic( 0 );
314    else               setCharacteristic( -nGetChar() );
315    CanonicalForm Fa,Gb;
316    if (currRing->minpoly!=NULL)
317    {
318      CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
319                                                 currRing->algring);
320      Variable a=rootOf(mipo);
321      CanonicalForm F( convSingAPFactoryAP( f,a ) ), G( convSingAPFactoryAP( g,a ) );
322      CanonicalForm FpG=F+G;
323      if (!(FpG.isUnivariate()|| FpG.inCoeffDomain()))
324      //if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
325      {
326        WerrorS("not univariate");
327        return TRUE;
328      }
329      res= convFactoryAPSingAP( extgcd( F, G, Fa, Gb ) );
330      pa=convFactoryAPSingAP(Fa);
331      pb=convFactoryAPSingAP(Gb);
332    }
333    else
334    {
335      CanonicalForm F( convSingTrPFactoryP( f ) ), G( convSingTrPFactoryP( g ) );
336      CanonicalForm FpG=F+G;
337      if (!(FpG.isUnivariate()|| FpG.inCoeffDomain()))
338      //if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
339      {
340        Off(SW_RATIONAL);
341        WerrorS("not univariate");
342        return TRUE;
343      }
344      res= convFactoryPSingTrP( extgcd( F, G, Fa, Gb ) );
345      pa=convFactoryPSingTrP(Fa);
346      pb=convFactoryPSingTrP(Gb);
347    }
348    Off(SW_RATIONAL);
349  }
350  else
351  {
352    WerrorS( feNotImplemented );
353    return TRUE;
354  }
355  return FALSE;
356}
357
358BOOLEAN singclap_extgcd_r ( poly f, poly g, poly &res, poly &pa, poly &pb, const ring r )
359{
360  // for now there is only the possibility to handle univariate
361  // polynomials over
362  // Q and Fp ...
363  res=NULL;pa=NULL;pb=NULL;
364  On(SW_SYMMETRIC_FF);
365  if (( n_GetChar(r) == 0 || n_GetChar(r) > 1 )
366  && (r->parameter==NULL))
367  {
368    setCharacteristic( n_GetChar(r) );
369    CanonicalForm F( convSingPFactoryP( f,r ) ), G( convSingPFactoryP( g,r) );
370    CanonicalForm FpG=F+G;
371    if (!(FpG.isUnivariate()|| FpG.inCoeffDomain()))
372    //if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
373    {
374      Off(SW_RATIONAL);
375      WerrorS("not univariate");
376      return TRUE;
377    }
378    CanonicalForm Fa,Gb;
379    On(SW_RATIONAL);
380    res=convFactoryPSingP( extgcd( F, G, Fa, Gb ),r );
381    pa=convFactoryPSingP(Fa,r);
382    pb=convFactoryPSingP(Gb,r);
383    Off(SW_RATIONAL);
384  }
385  // and over Q(a) / Fp(a)
386  else if (( n_GetChar(r)==1 ) /* Q(a) */
387  || (n_GetChar(r) <-1))       /* Fp(a) */
388  {
389    if (n_GetChar(r)==1) setCharacteristic( 0 );
390    else                 setCharacteristic( -n_GetChar(r) );
391    CanonicalForm Fa,Gb;
392    if (r->minpoly!=NULL)
393    {
394      CanonicalForm mipo=convSingPFactoryP(((lnumber)r->minpoly)->z,
395                                            r->algring);
396      Variable a=rootOf(mipo);
397      CanonicalForm F( convSingAPFactoryAP( f,a ) ),
398                    G( convSingAPFactoryAP( g,a ) );
399      CanonicalForm FpG=F+G;
400      if (!(FpG.isUnivariate()|| FpG.inCoeffDomain()))
401      //if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
402      {
403        WerrorS("not univariate");
404        return TRUE;
405      }
406      res= convFactoryAPSingAP( extgcd( F, G, Fa, Gb ) );
407      pa=convFactoryAPSingAP(Fa);
408      pb=convFactoryAPSingAP(Gb);
409    }
410    else
411    {
412      CanonicalForm F( convSingTrPFactoryP( f, r ) ), G( convSingTrPFactoryP( g, r ) );
413      CanonicalForm FpG=F+G;
414      if (!(FpG.isUnivariate()|| FpG.inCoeffDomain()))
415      //if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
416      {
417        Off(SW_RATIONAL);
418        WerrorS("not univariate");
419        return TRUE;
420      }
421      res= convFactoryPSingTrP( extgcd( F, G, Fa, Gb ), r );
422      pa=convFactoryPSingTrP(Fa, r);
423      pb=convFactoryPSingTrP(Gb, r);
424    }
425    Off(SW_RATIONAL);
426  }
427  else
428  {
429    WerrorS( feNotImplemented );
430    return TRUE;
431  }
432  return FALSE;
433}
434
435poly singclap_pdivide ( poly f, poly g )
436{
437  poly res=NULL;
438  On(SW_RATIONAL);
439  if (rField_is_Zp() || rField_is_Q())
440  {
441    setCharacteristic( nGetChar() );
442    CanonicalForm F( convSingPFactoryP( f ) ), G( convSingPFactoryP( g ) );
443    res = convFactoryPSingP( F / G );
444  }
445  else if (rField_is_Extension())
446  {
447    if (rField_is_Q_a()) setCharacteristic( 0 );
448    else               setCharacteristic( -nGetChar() );
449    if (currRing->minpoly!=NULL)
450    {
451      CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
452                                                 currRing->algring);
453      Variable a=rootOf(mipo);
454      CanonicalForm F( convSingAPFactoryAP( f,a ) ), G( convSingAPFactoryAP( g,a ) );
455      res= convFactoryAPSingAP(  F / G  );
456    }
457    else
458    {
459      CanonicalForm F( convSingTrPFactoryP( f ) ), G( convSingTrPFactoryP( g ) );
460      res= convFactoryPSingTrP(  F / G  );
461    }
462  }
463  #if 0 // not yet working
464  else if (rField_is_GF())
465  {
466    //Print("GF(%d^%d)\n",nfCharP,nfMinPoly[0]);
467    setCharacteristic( nfCharP,nfMinPoly[0], currRing->parameter[0][0] );
468    CanonicalForm F( convSingGFFactoryGF( f ) ), G( convSingGFFactoryGF( g ) );
469    res = convFactoryGFSingGF( F / G );
470  }
471  #endif
472  else
473    WerrorS( feNotImplemented );
474  Off(SW_RATIONAL);
475  return res;
476}
477
478void singclap_divide_content ( poly f )
479{
480  if ( f==NULL )
481  {
482    return;
483  }
484  else  if ( pNext( f ) == NULL )
485  {
486    pSetCoeff( f, nInit( 1 ) );
487    return;
488  }
489  else
490  {
491    if ( rField_is_Q_a() )
492      setCharacteristic( 0 );
493    else  if ( rField_is_Zp_a() )
494      setCharacteristic( -nGetChar() );
495    else
496      return; /* not implemented*/
497
498    CFList L;
499    CanonicalForm g, h;
500    poly p = pNext(f);
501
502    // first attemp: find 2 smallest g:
503
504    number g1=pGetCoeff(f);
505    number g2=pGetCoeff(p); // p==pNext(f);
506    pIter(p);
507    int sz1=nSize(g1);
508    int sz2=nSize(g2);
509    if (sz1>sz2)
510    {
511      number gg=g1;
512      g1=g2; g2=gg;
513      int sz=sz1;
514      sz1=sz2; sz2=sz;
515    }
516    while (p!=NULL)
517    {
518      int n_sz=nSize(pGetCoeff(p));
519      if (n_sz<sz1)
520      {
521        sz2=sz1;
522        g2=g1;
523        g1=pGetCoeff(p);
524        sz1=n_sz;
525        if (sz1<=3) break;
526      }
527      else if(n_sz<sz2)
528      {
529        sz2=n_sz;
530        g2=pGetCoeff(p);
531        sz2=n_sz;
532      }
533      pIter(p);
534    }
535    g = convSingPFactoryP( ((lnumber)g1)->z, currRing->algring );
536    g = gcd( g, convSingPFactoryP( ((lnumber)g2)->z , currRing->algring));
537
538    // second run: gcd's
539
540    p = f;
541    while ( (p != NULL) && (g != 1)  && ( g != 0))
542    {
543      h = convSingPFactoryP( ((lnumber)pGetCoeff(p))->z, currRing->algring );
544      pIter( p );
545
546      g = gcd( g, h );
547
548      L.append( h );
549    }
550    if (( g == 1 ) || (g == 0))
551    {
552      // pTest(f);
553      return;
554    }
555    else
556    {
557      CFListIterator i;
558      for ( i = L, p = f; i.hasItem(); i++, p=pNext(p) )
559      {
560        lnumber c=(lnumber)pGetCoeff(p);
561        p_Delete(&c->z,nacRing);
562        c->z=convFactoryPSingP( i.getItem() / g, currRing->algring );
563        //nTest((number)c);
564        //#ifdef LDEBUG
565        //number cn=(number)c;
566        //StringSetS(""); nWrite(nt); StringAppend(" ==> ");
567        //nWrite(cn);PrintS(StringAppend("\n"));
568        //#endif
569      }
570    }
571    // pTest(f);
572  }
573}
574
575static int primepower(int c)
576{
577  int p=1;
578  int cc=c;
579  while(cc!= rInternalChar(currRing)) { cc*=c; p++; }
580  return p;
581}
582
583BOOLEAN count_Factors(ideal I, intvec *v,int j, poly &f, poly fac)
584{
585  pTest(f);
586  pTest(fac);
587  int e=0;
588  if (!pIsConstantPoly(fac))
589  {
590#ifdef FACTORIZE2_DEBUG
591    printf("start count_Factors(%d), Fdeg=%d, factor deg=%d\n",j,pTotaldegree(f),pTotaldegree(fac));
592    p_wrp(fac,currRing);PrintLn();
593#endif
594    On(SW_RATIONAL);
595    CanonicalForm F, FAC,Q,R;
596    Variable a;
597    if (rField_is_Zp() || rField_is_Q())
598    {
599      F=convSingPFactoryP( f );
600      FAC=convSingPFactoryP( fac );
601    }
602    else if (rField_is_Extension())
603    {
604      if (currRing->minpoly!=NULL)
605      {
606        CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
607                                    currRing->algring);
608        a=rootOf(mipo);
609        F=convSingAPFactoryAP( f,a );
610        FAC=convSingAPFactoryAP( fac,a );
611      }
612      else
613      {
614        F=convSingTrPFactoryP( f );
615        FAC=convSingTrPFactoryP( fac );
616      }
617    }
618    else
619      WerrorS( feNotImplemented );
620
621    poly q;
622    loop
623    {
624      Q=F;
625      Q/=FAC;
626      R=Q;
627      R*=FAC;
628      R-=F;
629      if (R.isZero())
630      {
631        if (rField_is_Zp() || rField_is_Q())
632        {
633          q = convFactoryPSingP( Q );
634        }
635        else if (rField_is_Extension())
636        {
637          if (currRing->minpoly!=NULL)
638          {
639            q= convFactoryAPSingAP(  Q  );
640          }
641          else
642          {
643            q= convFactoryPSingTrP(  Q  );
644          }
645        }
646        e++; pDelete(&f); f=q; q=NULL; F=Q;
647      }
648      else
649      {
650        break;
651      }
652    }
653    if (e==0)
654    {
655      Off(SW_RATIONAL);
656      return FALSE;
657    }
658  }
659  else e=1;
660  I->m[j]=fac;
661  if (v!=NULL) (*v)[j]=e;
662  Off(SW_RATIONAL);
663  return TRUE;
664}
665
666int singclap_factorize_retry;
667extern int libfac_interruptflag;
668
669ideal singclap_factorize ( poly f, intvec ** v , int with_exps)
670{
671  pTest(f);
672#ifdef FACTORIZE2_DEBUG
673  printf("singclap_factorize, degree %d\n",pTotaldegree(f));
674#endif
675  // with_exps: 3,1 return only true factors, no exponents
676  //            2 return true factors and exponents
677  //            0 return coeff, factors and exponents
678  BOOLEAN save_errorreported=errorreported;
679
680  ideal res=NULL;
681
682  // handle factorize(0) =========================================
683  if (f==NULL)
684  {
685    res=idInit(1,1);
686    if (with_exps!=1)
687    {
688      (*v)=new intvec(1);
689      (**v)[0]=1;
690    }
691    return res;
692  }
693  // handle factorize(mon) =========================================
694  if (pNext(f)==NULL)
695  {
696    int i=0;
697    int n=0;
698    int e;
699    for(i=pVariables;i>0;i--) if(pGetExp(f,i)!=0) n++;
700    if (with_exps==0) n++; // with coeff
701    res=idInit(si_max(n,1),1);
702    switch(with_exps)
703    {
704      case 0: // with coef & exp.
705        res->m[0]=pOne();
706        pSetCoeff(res->m[0],nCopy(pGetCoeff(f)));
707        // no break
708      case 2: // with exp.
709        (*v)=new intvec(si_max(1,n));
710        (**v)[0]=1;
711        // no break
712      case 1: ;
713      #ifdef TEST
714      default: ;
715      #endif
716    }
717    if (n==0)
718    {
719      res->m[0]=pOne();
720      // (**v)[0]=1; is already done
721      return res;
722    }
723    for(i=pVariables;i>0;i--)
724    {
725      e=pGetExp(f,i);
726      if(e!=0)
727      {
728        n--;
729        poly p=pOne();
730        pSetExp(p,i,1);
731        pSetm(p);
732        res->m[n]=p;
733        if (with_exps!=1) (**v)[n]=e;
734      }
735    }
736    return res;
737  }
738  //PrintS("S:");pWrite(f);PrintLn();
739  // use factory/libfac in general ==============================
740  Off(SW_RATIONAL);
741  On(SW_SYMMETRIC_FF);
742  #ifdef HAVE_NTL
743  extern int prime_number;
744  if(rField_is_Q()) prime_number=0;
745  #endif
746  CFFList L;
747  number N=NULL;
748  number NN=NULL;
749  number old_lead_coeff=nCopy(pGetCoeff(f));
750
751  if (!rField_is_Zp() && !rField_is_Zp_a()) /* Q, Q(a) */
752  {
753    //if (f!=NULL) // already tested at start of routine
754    {
755      number n0=nCopy(pGetCoeff(f));
756      if (with_exps==0)
757        N=nCopy(n0);
758      pCleardenom(f);
759      NN=nDiv(n0,pGetCoeff(f));
760      nDelete(&n0);
761      if (with_exps==0)
762      {
763        nDelete(&N);
764        N=nCopy(NN);
765      }
766    }
767  }
768  else if (rField_is_Zp_a())
769  {
770    //if (f!=NULL) // already tested at start of routine
771    if (singclap_factorize_retry==0)
772    {
773      number n0=nCopy(pGetCoeff(f));
774      if (with_exps==0)
775        N=nCopy(n0);
776      pNorm(f);
777      pCleardenom(f);
778      NN=nDiv(n0,pGetCoeff(f));
779      nDelete(&n0);
780      if (with_exps==0)
781      {
782        nDelete(&N);
783        N=nCopy(NN);
784      }
785    }
786  }
787  if (rField_is_Q() || rField_is_Zp())
788  {
789    setCharacteristic( nGetChar() );
790    CanonicalForm F( convSingPFactoryP( f ) );
791    if (rField_is_Q())
792    {
793      L = factorize( F );
794    }
795    else /* Fp */
796    {
797      do
798      {
799        libfac_interruptflag=0;
800        L = Factorize( F );
801      }
802      while ((libfac_interruptflag!=0) ||(L.isEmpty()));
803    }
804  }
805  #if 0
806  else if (rField_is_GF())
807  {
808    int c=rChar(currRing);
809    setCharacteristic( c, primepower(c) );
810    CanonicalForm F( convSingGFFactoryGF( f ) );
811    if (F.isUnivariate())
812    {
813      L = factorize( F );
814    }
815    else
816    {
817      goto notImpl;
818    }
819  }
820  #endif
821  // and over Q(a) / Fp(a)
822  else if (rField_is_Extension())
823  {
824    if (rField_is_Q_a()) setCharacteristic( 0 );
825    else                 setCharacteristic( -nGetChar() );
826    if (currRing->minpoly!=NULL)
827    {
828      CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
829                   currRing->algring);
830      Variable a=rootOf(mipo);
831      CanonicalForm F( convSingAPFactoryAP( f,a ) );
832      if (rField_is_Zp_a() && F.isUnivariate())
833      {
834        L = factorize( F, a );
835      }
836      else
837      {
838        CanonicalForm G( convSingTrPFactoryP( f ) );
839        //  over Q(a) / multivariate over Fp(a)
840        do
841        {
842          libfac_interruptflag=0;
843          L=Factorize2(G, mipo);
844        }
845        while ((libfac_interruptflag!=0) ||(L.isEmpty()));
846        #ifdef FACTORIZE2_DEBUG
847        printf("while okay\n");
848        #endif
849        libfac_interruptflag=0;
850      }
851    }
852    else
853    {
854      CanonicalForm F( convSingTrPFactoryP( f ) );
855      if (rField_is_Q_a())
856      {
857        L = factorize( F );
858      }
859      else /* Fp(a) */
860      {
861        L = Factorize( F );
862      }
863    }
864  }
865  else
866  {
867    goto notImpl;
868  }
869  {
870    poly ff=pCopy(f); // a copy for the retry stuff
871    // the first factor should be a constant
872    if ( ! L.getFirst().factor().inCoeffDomain() )
873      L.insert(CFFactor(1,1));
874    // convert into ideal
875    int n = L.length();
876    if (n==0) n=1;
877    CFFListIterator J=L;
878    int j=0;
879    if (with_exps!=1)
880    {
881      if ((with_exps==2)&&(n>1))
882      {
883        n--;
884        J++;
885      }
886      *v = new intvec( n );
887    }
888    res = idInit( n ,1);
889    for ( ; J.hasItem(); J++, j++ )
890    {
891      poly p;
892      if (with_exps!=1) (**v)[j] = J.getItem().exp();
893      if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
894      {
895        //count_Factors(res,*v,f, j, convFactoryPSingP( J.getItem().factor() );
896        res->m[j] = convFactoryPSingP( J.getItem().factor() );
897      }
898      #if 0
899      else if (rField_is_GF())
900        res->m[j] = convFactoryGFSingGF( J.getItem().factor() );
901      #endif
902      else if (rField_is_Extension())     /* Q(a), Fp(a) */
903      {
904        intvec *w=NULL;
905        if (v!=NULL) w=*v;
906        if (currRing->minpoly==NULL)
907        {
908          if(!count_Factors(res,w,j,ff,convFactoryPSingTrP( J.getItem().factor() )))
909          {
910            if (w!=NULL)
911              (*w)[j]=1;
912            res->m[j]=pOne();
913          }
914        }
915        else
916        {
917          if (!count_Factors(res,w,j,ff,convFactoryAPSingAP( J.getItem().factor() )))
918          {
919            if (w!=NULL)
920              (*w)[j]=1;
921            res->m[j]=pOne();
922          }
923        }
924      }
925    }
926    if (rField_is_Extension() && (!pIsConstantPoly(ff)))
927    {
928      singclap_factorize_retry++;
929      if (singclap_factorize_retry<3)
930      {
931        int jj;
932        #ifdef FACTORIZE2_DEBUG
933        printf("factorize_retry\n");
934        #endif
935        intvec *ww=NULL;
936        idTest(res);
937        ideal h=singclap_factorize ( ff, &ww , with_exps);
938        idTest(h);
939        int l=(*v)->length();
940        (*v)->resize(l+ww->length());
941        for(jj=0;jj<ww->length();jj++)
942          (**v)[jj+l]=(*ww)[jj];
943        delete ww;
944        ideal hh=idInit(IDELEMS(res)+IDELEMS(h),1);
945        for(jj=IDELEMS(res)-1;jj>=0;jj--)
946        {
947          hh->m[jj]=res->m[jj];
948          res->m[jj]=NULL;
949        }
950        for(jj=IDELEMS(h)-1;jj>=0;jj--)
951        {
952          hh->m[jj+IDELEMS(res)]=h->m[jj];
953          h->m[jj]=NULL;
954        }
955        idDelete(&res);
956        idDelete(&h);
957        res=hh;
958        idTest(res);
959        ff=NULL;
960      }
961      else
962      {
963        WarnS("problem with factorize");
964        #if 0
965        pWrite(ff);
966        idShow(res);
967        #endif
968        idDelete(&res);
969        res=idInit(2,1);
970        res->m[0]=pOne();
971        res->m[1]=ff; ff=NULL;
972      }
973    }
974    pDelete(&ff);
975    if (N!=NULL)
976    {
977      pMult_nn(res->m[0],N);
978      nDelete(&N);
979      N=NULL;
980    }
981    // delete constants
982    if (res!=NULL)
983    {
984      int i=IDELEMS(res)-1;
985      int j=0;
986      for(;i>=0;i--)
987      {
988        if ((res->m[i]!=NULL)
989        && (pNext(res->m[i])==NULL)
990        && (pIsConstant(res->m[i])))
991        {
992          if (with_exps!=0)
993          {
994            pDelete(&(res->m[i]));
995            if ((v!=NULL) && ((*v)!=NULL))
996              (**v)[i]=0;
997            j++;
998          }
999          else if (i!=0)
1000          {
1001            while ((v!=NULL) && ((*v)!=NULL) && ((**v)[i]>1))
1002            {
1003              res->m[0]=pMult(res->m[0],pCopy(res->m[i]));
1004              (**v)[i]--;
1005            }
1006            res->m[0]=pMult(res->m[0],res->m[i]);
1007            res->m[i]=NULL;
1008            if ((v!=NULL) && ((*v)!=NULL))
1009              (**v)[i]=1;
1010            j++;
1011          }
1012        }
1013      }
1014      if (j>0)
1015      {
1016        idSkipZeroes(res);
1017        if ((v!=NULL) && ((*v)!=NULL))
1018        {
1019          intvec *w=*v;
1020          int len=IDELEMS(res);
1021          *v = new intvec( len );
1022          for (i=0,j=0;i<si_min(w->length(),len);i++)
1023          {
1024            if((*w)[i]!=0)
1025            {
1026              (**v)[j]=(*w)[i]; j++;
1027            }
1028          }
1029          delete w;
1030        }
1031      }
1032      if (res->m[0]==NULL)
1033      {
1034        res->m[0]=pOne();
1035      }
1036    }
1037  }
1038  if (rField_is_Q_a() && (currRing->minpoly!=NULL))
1039  {
1040    int i=IDELEMS(res)-1;
1041    int stop=1;
1042    if (with_exps!=0) stop=0;
1043    for(;i>=stop;i--)
1044    {
1045      pNorm(res->m[i]);
1046    }
1047    if (with_exps==0) pSetCoeff(res->m[0],old_lead_coeff);
1048    else nDelete(&old_lead_coeff);
1049  }
1050  else
1051    nDelete(&old_lead_coeff);
1052  errorreported=save_errorreported;
1053notImpl:
1054  if (res==NULL)
1055    WerrorS( feNotImplemented );
1056  if (NN!=NULL)
1057  {
1058    nDelete(&NN);
1059  }
1060  if (N!=NULL)
1061  {
1062    nDelete(&N);
1063  }
1064  //if (f!=NULL) pDelete(&f);
1065  //PrintS("......S\n");
1066  return res;
1067}
1068ideal singclap_sqrfree ( poly f)
1069{
1070  pTest(f);
1071#ifdef FACTORIZE2_DEBUG
1072  printf("singclap_sqrfree, degree %d\n",pTotaldegree(f));
1073#endif
1074  // with_exps: 3,1 return only true factors, no exponents
1075  //            2 return true factors and exponents
1076  //            0 return coeff, factors and exponents
1077  BOOLEAN save_errorreported=errorreported;
1078
1079  ideal res=NULL;
1080
1081  // handle factorize(0) =========================================
1082  if (f==NULL)
1083  {
1084    res=idInit(1,1);
1085    return res;
1086  }
1087  // handle factorize(mon) =========================================
1088  if (pNext(f)==NULL)
1089  {
1090    int i=0;
1091    int n=0;
1092    int e;
1093    for(i=pVariables;i>0;i--) if(pGetExp(f,i)!=0) n++;
1094    n++; // with coeff
1095    res=idInit(si_max(n,1),1);
1096    res->m[0]=pOne();
1097    pSetCoeff(res->m[0],nCopy(pGetCoeff(f)));
1098    if (n==0)
1099    {
1100      res->m[0]=pOne();
1101      // (**v)[0]=1; is already done
1102      return res;
1103    }
1104    for(i=pVariables;i>0;i--)
1105    {
1106      e=pGetExp(f,i);
1107      if(e!=0)
1108      {
1109        n--;
1110        poly p=pOne();
1111        pSetExp(p,i,1);
1112        pSetm(p);
1113        res->m[n]=p;
1114      }
1115    }
1116    return res;
1117  }
1118  //PrintS("S:");pWrite(f);PrintLn();
1119  // use factory/libfac in general ==============================
1120  Off(SW_RATIONAL);
1121  On(SW_SYMMETRIC_FF);
1122  #ifdef HAVE_NTL
1123  extern int prime_number;
1124  if(rField_is_Q()) prime_number=0;
1125  #endif
1126  CFFList L;
1127
1128  if (!rField_is_Zp() && !rField_is_Zp_a()) /* Q, Q(a) */
1129  {
1130    //if (f!=NULL) // already tested at start of routine
1131    {
1132      pCleardenom(f);
1133    }
1134  }
1135  else if (rField_is_Zp_a())
1136  {
1137    //if (f!=NULL) // already tested at start of routine
1138    if (singclap_factorize_retry==0)
1139    {
1140      pNorm(f);
1141      pCleardenom(f);
1142    }
1143  }
1144  if (rField_is_Q() || rField_is_Zp())
1145  {
1146    setCharacteristic( nGetChar() );
1147    CanonicalForm F( convSingPFactoryP( f ) );
1148    L = sqrFree( F );
1149  }
1150  #if 0
1151  else if (rField_is_GF())
1152  {
1153    int c=rChar(currRing);
1154    setCharacteristic( c, primepower(c) );
1155    CanonicalForm F( convSingGFFactoryGF( f ) );
1156    if (F.isUnivariate())
1157    {
1158      L = factorize( F );
1159    }
1160    else
1161    {
1162      goto notImpl;
1163    }
1164  }
1165  #endif
1166  // and over Q(a) / Fp(a)
1167  else if (rField_is_Extension())
1168  {
1169    if (rField_is_Q_a()) setCharacteristic( 0 );
1170    else                 setCharacteristic( -nGetChar() );
1171    if (currRing->minpoly!=NULL)
1172    {
1173      CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
1174                    currRing->algring);
1175      Variable a=rootOf(mipo);
1176      CanonicalForm F( convSingAPFactoryAP( f,a ) );
1177      CFFList SqrFreeMV( const CanonicalForm & f , const CanonicalForm & mipo=0) ;
1178
1179      L = SqrFreeMV( F,mipo );
1180      //WarnS("L = sqrFree( F,mipo );");
1181      //L = sqrFree( F );
1182    }
1183    else
1184    {
1185      CanonicalForm F( convSingTrPFactoryP( f ) );
1186      L = sqrFree( F );
1187    }
1188  }
1189  else
1190  {
1191    goto notImpl;
1192  }
1193  {
1194    poly ff=pCopy(f); // a copy for the retry stuff
1195    // convert into ideal
1196    int n = L.length();
1197    if (n==0) n=1;
1198    CFFListIterator J=L;
1199    int j=0;
1200    res = idInit( n ,1);
1201    for ( ; J.hasItem(); J++, j++ )
1202    {
1203      poly p;
1204      if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
1205        //count_Factors(res,*v,f, j, convFactoryPSingP( J.getItem().factor() );
1206        res->m[j] = convFactoryPSingP( J.getItem().factor() );
1207      #if 0
1208      else if (rField_is_GF())
1209        res->m[j] = convFactoryGFSingGF( J.getItem().factor() );
1210      #endif
1211      else if (rField_is_Extension())     /* Q(a), Fp(a) */
1212      {
1213        if (currRing->minpoly==NULL)
1214          res->m[j]=convFactoryPSingTrP( J.getItem().factor() );
1215        else
1216          res->m[j]=convFactoryAPSingAP( J.getItem().factor() );
1217      }
1218    }
1219    if (res->m[0]==NULL)
1220    {
1221      res->m[0]=pOne();
1222    }
1223  }
1224  errorreported=save_errorreported;
1225notImpl:
1226  if (res==NULL)
1227    WerrorS( feNotImplemented );
1228  return res;
1229}
1230matrix singclap_irrCharSeries ( ideal I)
1231{
1232  if (idIs0(I)) return mpNew(1,1);
1233
1234  // for now there is only the possibility to handle polynomials over
1235  // Q and Fp ...
1236  matrix res=NULL;
1237  int i;
1238  Off(SW_RATIONAL);
1239  On(SW_SYMMETRIC_FF);
1240  CFList L;
1241  ListCFList LL;
1242  if (((nGetChar() == 0) || (nGetChar() > 1) )
1243  && (currRing->parameter==NULL))
1244  {
1245    setCharacteristic( nGetChar() );
1246    for(i=0;i<IDELEMS(I);i++)
1247    {
1248      poly p=I->m[i];
1249      if (p!=NULL)
1250      {
1251        p=pCopy(p);
1252        pCleardenom(p);
1253        L.append(convSingPFactoryP(p));
1254      }
1255    }
1256  }
1257  // and over Q(a) / Fp(a)
1258  else if (( nGetChar()==1 ) /* Q(a) */
1259  || (nGetChar() <-1))       /* Fp(a) */
1260  {
1261    if (nGetChar()==1) setCharacteristic( 0 );
1262    else               setCharacteristic( -nGetChar() );
1263    for(i=0;i<IDELEMS(I);i++)
1264    {
1265      poly p=I->m[i];
1266      if (p!=NULL)
1267      {
1268        p=pCopy(p);
1269        pCleardenom(p);
1270        L.append(convSingTrPFactoryP(p));
1271      }
1272    }
1273  }
1274  else
1275  {
1276    WerrorS( feNotImplemented );
1277    return res;
1278  }
1279
1280  // a very bad work-around --- FIX IT in libfac
1281  // should be fixed as of 2001/6/27
1282  int tries=0;
1283  int m,n;
1284  ListIterator<CFList> LLi;
1285  loop
1286  {
1287    LL=IrrCharSeries(L);
1288    m= LL.length(); // Anzahl Zeilen
1289    n=0;
1290    for ( LLi = LL; LLi.hasItem(); LLi++ )
1291    {
1292      n = si_max(LLi.getItem().length(),n);
1293    }
1294    if ((m!=0) && (n!=0)) break;
1295    tries++;
1296    if (tries>=5) break;
1297  }
1298  if ((m==0) || (n==0))
1299  {
1300    Warn("char_series returns %d x %d matrix from %d input polys (%d)",
1301      m,n,IDELEMS(I)+1,LL.length());
1302    iiWriteMatrix((matrix)I,"I",2,0);
1303    m=si_max(m,1);
1304    n=si_max(n,1);
1305  }
1306  res=mpNew(m,n);
1307  CFListIterator Li;
1308  for ( m=1, LLi = LL; LLi.hasItem(); LLi++, m++ )
1309  {
1310    for (n=1, Li = LLi.getItem(); Li.hasItem(); Li++, n++)
1311    {
1312      if ( (nGetChar() == 0) || (nGetChar() > 1) )
1313        MATELEM(res,m,n)=convFactoryPSingP(Li.getItem());
1314      else
1315        MATELEM(res,m,n)=convFactoryPSingTrP(Li.getItem());
1316    }
1317  }
1318  Off(SW_RATIONAL);
1319  return res;
1320}
1321
1322char* singclap_neworder ( ideal I)
1323{
1324  int i;
1325  Off(SW_RATIONAL);
1326  On(SW_SYMMETRIC_FF);
1327  CFList L;
1328  if (((nGetChar() == 0) || (nGetChar() > 1) )
1329  && (currRing->parameter==NULL))
1330  {
1331    setCharacteristic( nGetChar() );
1332    for(i=0;i<IDELEMS(I);i++)
1333    {
1334      L.append(convSingPFactoryP(I->m[i]));
1335    }
1336  }
1337  // and over Q(a) / Fp(a)
1338  else if (( nGetChar()==1 ) /* Q(a) */
1339  || (nGetChar() <-1))       /* Fp(a) */
1340  {
1341    if (nGetChar()==1) setCharacteristic( 0 );
1342    else               setCharacteristic( -nGetChar() );
1343    for(i=0;i<IDELEMS(I);i++)
1344    {
1345      L.append(convSingTrPFactoryP(I->m[i]));
1346    }
1347  }
1348  else
1349  {
1350    WerrorS( feNotImplemented );
1351    return NULL;
1352  }
1353
1354  List<int> IL=neworderint(L);
1355  ListIterator<int> Li;
1356  StringSetS("");
1357  Li = IL;
1358  int offs=rPar(currRing);
1359  int* mark=(int*)omAlloc0((pVariables+offs)*sizeof(int));
1360  int cnt=pVariables+offs;
1361  loop
1362  {
1363    if(! Li.hasItem()) break;
1364    BOOLEAN done=TRUE;
1365    i=Li.getItem()-1;
1366    mark[i]=1;
1367    if (i<offs)
1368    {
1369      done=FALSE;
1370      //StringAppendS(currRing->parameter[i]);
1371    }
1372    else
1373    {
1374      StringAppendS(currRing->names[i-offs]);
1375    }
1376    Li++;
1377    cnt--;
1378    if(cnt==0) break;
1379    if (done) StringAppendS(",");
1380  }
1381  for(i=0;i<pVariables+offs;i++)
1382  {
1383    BOOLEAN done=TRUE;
1384    if(mark[i]==0)
1385    {
1386      if (i<offs)
1387      {
1388        done=FALSE;
1389        //StringAppendS(currRing->parameter[i]);
1390      }
1391      else
1392      {
1393        StringAppendS(currRing->names[i-offs]);
1394      }
1395      cnt--;
1396      if(cnt==0) break;
1397      if (done) StringAppendS(",");
1398    }
1399  }
1400  char * s=omStrDup(StringAppendS(""));
1401  if (s[strlen(s)-1]==',') s[strlen(s)-1]='\0';
1402  return s;
1403}
1404
1405BOOLEAN singclap_isSqrFree(poly f)
1406{
1407  BOOLEAN b=FALSE;
1408  Off(SW_RATIONAL);
1409  //  Q / Fp
1410  if (((nGetChar() == 0) || (nGetChar() > 1) )
1411  &&(currRing->parameter==NULL))
1412  {
1413    setCharacteristic( nGetChar() );
1414    CanonicalForm F( convSingPFactoryP( f ) );
1415    if((nGetChar()>1)&&(!F.isUnivariate()))
1416      goto err;
1417    b=(BOOLEAN)isSqrFree(F);
1418  }
1419  // and over Q(a) / Fp(a)
1420  else if (( nGetChar()==1 ) /* Q(a) */
1421  || (nGetChar() <-1))       /* Fp(a) */
1422  {
1423    if (nGetChar()==1) setCharacteristic( 0 );
1424    else               setCharacteristic( -nGetChar() );
1425    //if (currRing->minpoly!=NULL)
1426    //{
1427    //  CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
1428    //                                             currRing->algring);
1429    //  Variable a=rootOf(mipo);
1430    //  CanonicalForm F( convSingAPFactoryAP( f,a ) );
1431    //  ...
1432    //}
1433    //else
1434    {
1435      CanonicalForm F( convSingTrPFactoryP( f ) );
1436      b=(BOOLEAN)isSqrFree(F);
1437    }
1438    Off(SW_RATIONAL);
1439  }
1440  else
1441  {
1442err:
1443    WerrorS( feNotImplemented );
1444  }
1445  return b;
1446}
1447
1448poly singclap_det( const matrix m )
1449{
1450  int r=m->rows();
1451  if (r!=m->cols())
1452  {
1453    Werror("det of %d x %d matrix",r,m->cols());
1454    return NULL;
1455  }
1456  poly res=NULL;
1457  if (( nGetChar() == 0 || nGetChar() > 1 )
1458  && (currRing->parameter==NULL))
1459  {
1460    setCharacteristic( nGetChar() );
1461    CFMatrix M(r,r);
1462    int i,j;
1463    for(i=r;i>0;i--)
1464    {
1465      for(j=r;j>0;j--)
1466      {
1467        M(i,j)=convSingPFactoryP(MATELEM(m,i,j));
1468      }
1469    }
1470    res= convFactoryPSingP( determinant(M,r) ) ;
1471  }
1472  // and over Q(a) / Fp(a)
1473  else if (( nGetChar()==1 ) /* Q(a) */
1474  || (nGetChar() <-1))       /* Fp(a) */
1475  {
1476    if (nGetChar()==1) setCharacteristic( 0 );
1477    else               setCharacteristic( -nGetChar() );
1478    CFMatrix M(r,r);
1479    poly res;
1480    if (currRing->minpoly!=NULL)
1481    {
1482      CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
1483                                            currRing->algring);
1484      Variable a=rootOf(mipo);
1485      int i,j;
1486      for(i=r;i>0;i--)
1487      {
1488        for(j=r;j>0;j--)
1489        {
1490          M(i,j)=convSingAPFactoryAP(MATELEM(m,i,j),a);
1491        }
1492      }
1493      res= convFactoryAPSingAP( determinant(M,r) ) ;
1494    }
1495    else
1496    {
1497      int i,j;
1498      for(i=r;i>0;i--)
1499      {
1500        for(j=r;j>0;j--)
1501        {
1502          M(i,j)=convSingTrPFactoryP(MATELEM(m,i,j));
1503        }
1504      }
1505      res= convFactoryPSingTrP( determinant(M,r) );
1506    }
1507  }
1508  else
1509    WerrorS( feNotImplemented );
1510  Off(SW_RATIONAL);
1511  return res;
1512}
1513
1514int singclap_det_i( intvec * m )
1515{
1516  setCharacteristic( 0 );
1517  CFMatrix M(m->rows(),m->cols());
1518  int i,j;
1519  for(i=m->rows();i>0;i--)
1520  {
1521    for(j=m->cols();j>0;j--)
1522    {
1523      M(i,j)=IMATELEM(*m,i,j);
1524    }
1525  }
1526  int res= convFactoryISingI( determinant(M,m->rows())) ;
1527  Off(SW_RATIONAL);
1528  return res;
1529}
1530napoly singclap_alglcm ( napoly f, napoly g )
1531{
1532
1533 // over Q(a) / Fp(a)
1534 if (nGetChar()==1) setCharacteristic( 0 );
1535 else               setCharacteristic( -nGetChar() );
1536 napoly res;
1537
1538 if (currRing->minpoly!=NULL)
1539 {
1540   CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
1541                                         currRing->algring);
1542   Variable a=rootOf(mipo);
1543   CanonicalForm F( convSingAFactoryA( f,a ) ), G( convSingAFactoryA( g,a ) );
1544   CanonicalForm GCD;
1545
1546   // calculate gcd
1547   GCD = gcd( F, G );
1548
1549   // calculate lcm
1550   res= convFactoryASingA( (F/GCD)*G );
1551 }
1552 else
1553 {
1554   CanonicalForm F( convSingPFactoryP( f,currRing->algring ) ),
1555                 G( convSingPFactoryP( g,currRing->algring ) );
1556   CanonicalForm GCD;
1557   // calculate gcd
1558   GCD = gcd( F, G );
1559
1560   // calculate lcm
1561   res= convFactoryPSingP( (F/GCD)*G, currRing->algring );
1562 }
1563
1564 Off(SW_RATIONAL);
1565 return res;
1566}
1567
1568void singclap_algdividecontent ( napoly f, napoly g, napoly &ff, napoly &gg )
1569{
1570 // over Q(a) / Fp(a)
1571 if (nGetChar()==1) setCharacteristic( 0 );
1572 else               setCharacteristic( -nGetChar() );
1573 ff=gg=NULL;
1574 On(SW_RATIONAL);
1575
1576 if (currRing->minpoly!=NULL)
1577 {
1578   CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
1579                         currRing->algring);
1580   Variable a=rootOf(mipo);
1581   CanonicalForm F( convSingAFactoryA( f,a ) ), G( convSingAFactoryA( g,a ) );
1582   CanonicalForm GCD;
1583
1584   GCD=gcd( F, G );
1585
1586   if ((GCD!=1) && (GCD!=0))
1587   {
1588     ff= convFactoryASingA( F/ GCD );
1589     gg= convFactoryASingA( G/ GCD );
1590   }
1591 }
1592 else
1593 {
1594   CanonicalForm F( convSingPFactoryP( f,currRing->algring ) ),
1595                 G( convSingPFactoryP( g,currRing->algring ) );
1596   CanonicalForm GCD;
1597
1598   GCD=gcd( F, G );
1599
1600   if ((GCD!=1) && (GCD!=0))
1601   {
1602     ff= convFactoryPSingP( F/ GCD, currRing->algring );
1603     gg= convFactoryPSingP( G/ GCD, currRing->algring );
1604   }
1605 }
1606
1607 Off(SW_RATIONAL);
1608}
1609
1610#if 0
1611lists singclap_chineseRemainder(lists x, lists q)
1612{
1613  //assume(x->nr == q->nr);
1614  //assume(x->nr >= 0);
1615  int n=x->nr+1;
1616  if ((x->nr<0) || (x->nr!=q->nr))
1617  {
1618    WerrorS("list are empty or not of equal length");
1619    return NULL;
1620  }
1621  lists res=(lists)omAlloc0Bin(slists_bin);
1622  CFArray X(1,n), Q(1,n);
1623  int i;
1624  for(i=0; i<n; i++)
1625  {
1626    if (x->m[i-1].Typ()==INT_CMD)
1627    {
1628      X[i]=(int)x->m[i-1].Data();
1629    }
1630    else if (x->m[i-1].Typ()==NUMBER_CMD)
1631    {
1632      number N=(number)x->m[i-1].Data();
1633      X[i]=convSingNFactoryN(N);
1634    }
1635    else
1636    {
1637      WerrorS("illegal type in chineseRemainder");
1638      omFreeBin(res,slists_bin);
1639      return NULL;
1640    }
1641    if (q->m[i-1].Typ()==INT_CMD)
1642    {
1643      Q[i]=(int)q->m[i-1].Data();
1644    }
1645    else if (q->m[i-1].Typ()==NUMBER_CMD)
1646    {
1647      number N=(number)x->m[i-1].Data();
1648      Q[i]=convSingNFactoryN(N);
1649    }
1650    else
1651    {
1652      WerrorS("illegal type in chineseRemainder");
1653      omFreeBin(res,slists_bin);
1654      return NULL;
1655    }
1656  }
1657  CanonicalForm r, prod;
1658  chineseRemainder( X, Q, r, prod );
1659  res->Init(2);
1660  res->m[0].rtyp=NUMBER_CMD;
1661  res->m[1].rtyp=NUMBER_CMD;
1662  res->m[0].data=(char *)convFactoryNSingN( r );
1663  res->m[1].data=(char *)convFactoryNSingN( prod );
1664  return res;
1665}
1666#endif
1667#endif
Note: See TracBrowser for help on using the repository browser.