source: git/kernel/clapsing.cc @ f2b6f0b

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