source: git/libpolys/polys/clapsing.cc @ a934fc8

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