source: git/kernel/clapsing.cc @ 4e35a89

spielwiese
Last change on this file since 4e35a89 was a41623, checked in by Hans Schoenemann <hannes@…>, 13 years ago
removed dead code, unused variables etc git-svn-id: file:///usr/local/Singular/svn/trunk@14249 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 41.9 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
618BOOLEAN count_Factors(ideal I, intvec *v,int j, poly &f, poly fac)
619{
620  pTest(f);
621  pTest(fac);
622  int e=0;
623  if (!pIsConstantPoly(fac))
624  {
625#ifdef FACTORIZE2_DEBUG
626    printf("start count_Factors(%d), Fdeg=%d, factor deg=%d\n",j,pTotaldegree(f),pTotaldegree(fac));
627    p_wrp(fac,currRing);PrintLn();
628#endif
629    On(SW_RATIONAL);
630    CanonicalForm F, FAC,Q,R;
631    Variable a;
632    if (rField_is_Zp() || rField_is_Q())
633    {
634      F=convSingPFactoryP( f );
635      FAC=convSingPFactoryP( fac );
636    }
637    else if (rField_is_Extension())
638    {
639      if (currRing->minpoly!=NULL)
640      {
641        CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
642                                    currRing->algring);
643        a=rootOf(mipo);
644        F=convSingAPFactoryAP( f,a,currRing );
645        FAC=convSingAPFactoryAP( fac,a,currRing );
646      }
647      else
648      {
649        F=convSingTrPFactoryP( f );
650        FAC=convSingTrPFactoryP( fac );
651      }
652    }
653    else
654      WerrorS( feNotImplemented );
655
656    poly q;
657    loop
658    {
659      Q=F;
660      Q/=FAC;
661      R=Q;
662      R*=FAC;
663      R-=F;
664      if (R.isZero())
665      {
666        if (rField_is_Zp() || rField_is_Q())
667        {
668          q = convFactoryPSingP( Q );
669        }
670        else if (rField_is_Extension())
671        {
672          if (currRing->minpoly!=NULL)
673          {
674            q= convFactoryAPSingAP(  Q,currRing  );
675          }
676          else
677          {
678            q= convFactoryPSingTrP(  Q  );
679          }
680        }
681        e++; pDelete(&f); f=q; q=NULL; F=Q;
682      }
683      else
684      {
685        break;
686      }
687    }
688    if (e==0)
689    {
690      Off(SW_RATIONAL);
691      return FALSE;
692    }
693  }
694  else e=1;
695  I->m[j]=fac;
696  if (v!=NULL) (*v)[j]=e;
697  Off(SW_RATIONAL);
698  return TRUE;
699}
700
701int singclap_factorize_retry;
702extern int libfac_interruptflag;
703
704ideal singclap_factorize ( poly f, intvec ** v , int with_exps)
705/* destroys f, sets *v */
706{
707  pTest(f);
708#ifdef FACTORIZE2_DEBUG
709  printf("singclap_factorize, degree %ld\n",pTotaldegree(f));
710#endif
711  // with_exps: 3,1 return only true factors, no exponents
712  //            2 return true factors and exponents
713  //            0 return coeff, factors and exponents
714  BOOLEAN save_errorreported=errorreported;
715
716  ideal res=NULL;
717
718  // handle factorize(0) =========================================
719  if (f==NULL)
720  {
721    res=idInit(1,1);
722    if (with_exps!=1)
723    {
724      (*v)=new intvec(1);
725      (**v)[0]=1;
726    }
727    return res;
728  }
729  // handle factorize(mon) =========================================
730  if (pNext(f)==NULL)
731  {
732    int i=0;
733    int n=0;
734    int e;
735    for(i=pVariables;i>0;i--) if(pGetExp(f,i)!=0) n++;
736    if (with_exps==0) n++; // with coeff
737    res=idInit(si_max(n,1),1);
738    switch(with_exps)
739    {
740      case 0: // with coef & exp.
741        res->m[0]=pNSet(nCopy(pGetCoeff(f)));
742        // no break
743      case 2: // with exp.
744        (*v)=new intvec(si_max(1,n));
745        (**v)[0]=1;
746        // no break
747      case 1: ;
748      #ifdef TEST
749      default: ;
750      #endif
751    }
752    if (n==0)
753    {
754      res->m[0]=pOne();
755      // (**v)[0]=1; is already done
756    }
757    else
758    {
759      for(i=pVariables;i>0;i--)
760      {
761        e=pGetExp(f,i);
762        if(e!=0)
763        {
764          n--;
765          poly p=pOne();
766          pSetExp(p,i,1);
767          pSetm(p);
768          res->m[n]=p;
769          if (with_exps!=1) (**v)[n]=e;
770        }
771      }
772    }
773    pDelete(&f);
774    return res;
775  }
776  //PrintS("S:");pWrite(f);PrintLn();
777  // use factory/libfac in general ==============================
778  Off(SW_RATIONAL);
779  On(SW_SYMMETRIC_FF);
780  #ifdef HAVE_NTL
781  extern int prime_number;
782  if(rField_is_Q()) prime_number=0;
783  #endif
784  CFFList L;
785  number N=NULL;
786  number NN=NULL;
787  number old_lead_coeff=nCopy(pGetCoeff(f));
788
789  if (!rField_is_Zp() && !rField_is_Zp_a()) /* Q, Q(a) */
790  {
791    //if (f!=NULL) // already tested at start of routine
792    {
793      number n0=nCopy(pGetCoeff(f));
794      if (with_exps==0)
795        N=nCopy(n0);
796      p_Cleardenom(f, currRing);
797      NN=nDiv(n0,pGetCoeff(f));
798      nDelete(&n0);
799      if (with_exps==0)
800      {
801        nDelete(&N);
802        N=nCopy(NN);
803      }
804    }
805  }
806  else if (rField_is_Zp_a())
807  {
808    //if (f!=NULL) // already tested at start of routine
809    if (singclap_factorize_retry==0)
810    {
811      number n0=nCopy(pGetCoeff(f));
812      if (with_exps==0)
813        N=nCopy(n0);
814      pNorm(f);
815      p_Cleardenom(f, currRing);
816      NN=nDiv(n0,pGetCoeff(f));
817      nDelete(&n0);
818      if (with_exps==0)
819      {
820        nDelete(&N);
821        N=nCopy(NN);
822      }
823    }
824  }
825  if (rField_is_Q() || rField_is_Zp())
826  {
827    setCharacteristic( nGetChar() );
828    CanonicalForm F( convSingPFactoryP( f ) );
829    L = factorize( F );
830  }
831  #if 0
832  else if (rField_is_GF())
833  {
834    int c=rChar(currRing);
835    setCharacteristic( c, primepower(c) );
836    CanonicalForm F( convSingGFFactoryGF( f ) );
837    if (F.isUnivariate())
838    {
839      L = factorize( F );
840    }
841    else
842    {
843      goto notImpl;
844    }
845  }
846  #endif
847  // and over Q(a) / Fp(a)
848  else if (rField_is_Extension())
849  {
850    if (rField_is_Q_a()) setCharacteristic( 0 );
851    else                 setCharacteristic( -nGetChar() );
852    if (currRing->minpoly!=NULL)
853    {
854      CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
855                   currRing->algring);
856      Variable a=rootOf(mipo);
857      CanonicalForm F( convSingAPFactoryAP( f,a,currRing ) );
858      if (rField_is_Zp_a())
859      {
860        L = factorize( F, a );
861      }
862      else
863      {
864        //  over Q(a)
865        if (F.isUnivariate())
866        {
867          L= factorize (F, a);
868        }
869        else
870        {
871          CanonicalForm G( convSingTrPFactoryP( f ) );
872          do
873          {
874            libfac_interruptflag=0;
875            L=Factorize2(G, mipo);
876          }
877          while ((libfac_interruptflag!=0) ||(L.isEmpty()));
878          #ifdef FACTORIZE2_DEBUG
879          printf("while okay\n");
880          #endif
881          libfac_interruptflag=0;
882        }
883      }
884    }
885    else
886    {
887      CanonicalForm F( convSingTrPFactoryP( f ) );
888      L = factorize( F );
889    }
890  }
891  else
892  {
893    goto notImpl;
894  }
895  {
896    poly ff=pCopy(f); // a copy for the retry stuff
897    // the first factor should be a constant
898    if ( ! L.getFirst().factor().inCoeffDomain() )
899      L.insert(CFFactor(1,1));
900    // convert into ideal
901    int n = L.length();
902    if (n==0) n=1;
903    CFFListIterator J=L;
904    int j=0;
905    if (with_exps!=1)
906    {
907      if ((with_exps==2)&&(n>1))
908      {
909        n--;
910        J++;
911      }
912      *v = new intvec( n );
913    }
914    res = idInit( n ,1);
915    for ( ; J.hasItem(); J++, j++ )
916    {
917      if (with_exps!=1) (**v)[j] = J.getItem().exp();
918      if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
919      {
920        //count_Factors(res,*v,f, j, convFactoryPSingP( J.getItem().factor() );
921        res->m[j] = convFactoryPSingP( J.getItem().factor() );
922      }
923      #if 0
924      else if (rField_is_GF())
925        res->m[j] = convFactoryGFSingGF( J.getItem().factor() );
926      #endif
927      else if (rField_is_Extension())     /* Q(a), Fp(a) */
928      {
929        intvec *w=NULL;
930        if (v!=NULL) w=*v;
931        if (currRing->minpoly==NULL)
932        {
933          if(!count_Factors(res,w,j,ff,convFactoryPSingTrP( J.getItem().factor() )))
934          {
935            if (w!=NULL)
936              (*w)[j]=1;
937            res->m[j]=pOne();
938          }
939        }
940        else
941        {
942          if (!count_Factors(res,w,j,ff,convFactoryAPSingAP( J.getItem().factor(),currRing )))
943          {
944            if (w!=NULL)
945              (*w)[j]=1;
946            res->m[j]=pOne();
947          }
948        }
949      }
950    }
951    if (rField_is_Extension() && (!pIsConstantPoly(ff)))
952    {
953      singclap_factorize_retry++;
954      if (singclap_factorize_retry<3)
955      {
956        int jj;
957        #ifdef FACTORIZE2_DEBUG
958        printf("factorize_retry\n");
959        #endif
960        intvec *ww=NULL;
961        idTest(res);
962        ideal h=singclap_factorize ( ff, &ww , with_exps);
963        idTest(h);
964        int l=(*v)->length();
965        (*v)->resize(l+ww->length());
966        for(jj=0;jj<ww->length();jj++)
967          (**v)[jj+l]=(*ww)[jj];
968        delete ww;
969        ideal hh=idInit(IDELEMS(res)+IDELEMS(h),1);
970        for(jj=IDELEMS(res)-1;jj>=0;jj--)
971        {
972          hh->m[jj]=res->m[jj];
973          res->m[jj]=NULL;
974        }
975        for(jj=IDELEMS(h)-1;jj>=0;jj--)
976        {
977          hh->m[jj+IDELEMS(res)]=h->m[jj];
978          h->m[jj]=NULL;
979        }
980        idDelete(&res);
981        idDelete(&h);
982        res=hh;
983        idTest(res);
984        ff=NULL;
985      }
986      else
987      {
988        WarnS("problem with factorize");
989        #if 0
990        pWrite(ff);
991        idShow(res);
992        #endif
993        idDelete(&res);
994        res=idInit(2,1);
995        res->m[0]=pOne();
996        res->m[1]=ff; ff=NULL;
997      }
998    }
999    pDelete(&ff);
1000    if (N!=NULL)
1001    {
1002      pMult_nn(res->m[0],N);
1003      nDelete(&N);
1004      N=NULL;
1005    }
1006    // delete constants
1007    if (res!=NULL)
1008    {
1009      int i=IDELEMS(res)-1;
1010      int j=0;
1011      for(;i>=0;i--)
1012      {
1013        if ((res->m[i]!=NULL)
1014        && (pNext(res->m[i])==NULL)
1015        && (pIsConstant(res->m[i])))
1016        {
1017          if (with_exps!=0)
1018          {
1019            pDelete(&(res->m[i]));
1020            if ((v!=NULL) && ((*v)!=NULL))
1021              (**v)[i]=0;
1022            j++;
1023          }
1024          else if (i!=0)
1025          {
1026            while ((v!=NULL) && ((*v)!=NULL) && ((**v)[i]>1))
1027            {
1028              res->m[0]=pMult(res->m[0],pCopy(res->m[i]));
1029              (**v)[i]--;
1030            }
1031            res->m[0]=pMult(res->m[0],res->m[i]);
1032            res->m[i]=NULL;
1033            if ((v!=NULL) && ((*v)!=NULL))
1034              (**v)[i]=1;
1035            j++;
1036          }
1037        }
1038      }
1039      if (j>0)
1040      {
1041        idSkipZeroes(res);
1042        if ((v!=NULL) && ((*v)!=NULL))
1043        {
1044          intvec *w=*v;
1045          int len=IDELEMS(res);
1046          *v = new intvec( len );
1047          for (i=0,j=0;i<si_min(w->length(),len);i++)
1048          {
1049            if((*w)[i]!=0)
1050            {
1051              (**v)[j]=(*w)[i]; j++;
1052            }
1053          }
1054          delete w;
1055        }
1056      }
1057      if (res->m[0]==NULL)
1058      {
1059        res->m[0]=pOne();
1060      }
1061    }
1062  }
1063  if (rField_is_Q_a() && (currRing->minpoly!=NULL))
1064  {
1065    int i=IDELEMS(res)-1;
1066    int stop=1;
1067    if (with_exps!=0) stop=0;
1068    for(;i>=stop;i--)
1069    {
1070      pNorm(res->m[i]);
1071    }
1072    if (with_exps==0) pSetCoeff(res->m[0],old_lead_coeff);
1073    else nDelete(&old_lead_coeff);
1074  }
1075  else
1076    nDelete(&old_lead_coeff);
1077  errorreported=save_errorreported;
1078notImpl:
1079  if (res==NULL)
1080    WerrorS( feNotImplemented );
1081  if (NN!=NULL)
1082  {
1083    nDelete(&NN);
1084  }
1085  if (N!=NULL)
1086  {
1087    nDelete(&N);
1088  }
1089  if (f!=NULL) pDelete(&f);
1090  //PrintS("......S\n");
1091  return res;
1092}
1093ideal singclap_sqrfree ( poly f)
1094{
1095  pTest(f);
1096#ifdef FACTORIZE2_DEBUG
1097  printf("singclap_sqrfree, degree %d\n",pTotaldegree(f));
1098#endif
1099  // with_exps: 3,1 return only true factors, no exponents
1100  //            2 return true factors and exponents
1101  //            0 return coeff, factors and exponents
1102  BOOLEAN save_errorreported=errorreported;
1103
1104  ideal res=NULL;
1105
1106  // handle factorize(0) =========================================
1107  if (f==NULL)
1108  {
1109    res=idInit(1,1);
1110    return res;
1111  }
1112  // handle factorize(mon) =========================================
1113  if (pNext(f)==NULL)
1114  {
1115    int i=0;
1116    int n=0;
1117    int e;
1118    for(i=pVariables;i>0;i--) if(pGetExp(f,i)!=0) n++;
1119    n++; // with coeff
1120    res=idInit(si_max(n,1),1);
1121    res->m[0]=pNSet(nCopy(pGetCoeff(f)));
1122    if (n==0)
1123    {
1124      res->m[0]=pOne();
1125      // (**v)[0]=1; is already done
1126      return res;
1127    }
1128    for(i=pVariables;i>0;i--)
1129    {
1130      e=pGetExp(f,i);
1131      if(e!=0)
1132      {
1133        n--;
1134        poly p=pOne();
1135        pSetExp(p,i,1);
1136        pSetm(p);
1137        res->m[n]=p;
1138      }
1139    }
1140    return res;
1141  }
1142  //PrintS("S:");pWrite(f);PrintLn();
1143  // use factory/libfac in general ==============================
1144  Off(SW_RATIONAL);
1145  On(SW_SYMMETRIC_FF);
1146  #ifdef HAVE_NTL
1147  extern int prime_number;
1148  if(rField_is_Q()) prime_number=0;
1149  #endif
1150  CFFList L;
1151
1152  if (!rField_is_Zp() && !rField_is_Zp_a()) /* Q, Q(a) */
1153  {
1154    //if (f!=NULL) // already tested at start of routine
1155    {
1156      p_Cleardenom(f, currRing);
1157    }
1158  }
1159  else if (rField_is_Zp_a())
1160  {
1161    //if (f!=NULL) // already tested at start of routine
1162    if (singclap_factorize_retry==0)
1163    {
1164      pNorm(f);
1165      p_Cleardenom(f, currRing);
1166    }
1167  }
1168  if (rField_is_Q() || rField_is_Zp())
1169  {
1170    setCharacteristic( nGetChar() );
1171    CanonicalForm F( convSingPFactoryP( f ) );
1172    L = sqrFree( F );
1173  }
1174  #if 0
1175  else if (rField_is_GF())
1176  {
1177    int c=rChar(currRing);
1178    setCharacteristic( c, primepower(c) );
1179    CanonicalForm F( convSingGFFactoryGF( f ) );
1180    if (F.isUnivariate())
1181    {
1182      L = factorize( F );
1183    }
1184    else
1185    {
1186      goto notImpl;
1187    }
1188  }
1189  #endif
1190  // and over Q(a) / Fp(a)
1191  else if (rField_is_Extension())
1192  {
1193    if (rField_is_Q_a()) setCharacteristic( 0 );
1194    else                 setCharacteristic( -nGetChar() );
1195    if (currRing->minpoly!=NULL)
1196    {
1197      CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
1198                    currRing->algring);
1199      Variable a=rootOf(mipo);
1200      CanonicalForm F( convSingAPFactoryAP( f,a,currRing ) );
1201      CFFList SqrFreeMV( const CanonicalForm & f , const CanonicalForm & mipo=0) ;
1202
1203      L = SqrFreeMV( F,mipo );
1204      //WarnS("L = sqrFree( F,mipo );");
1205      //L = sqrFree( F );
1206    }
1207    else
1208    {
1209      CanonicalForm F( convSingTrPFactoryP( f ) );
1210      L = sqrFree( F );
1211    }
1212  }
1213  else
1214  {
1215    goto notImpl;
1216  }
1217  {
1218    // convert into ideal
1219    int n = L.length();
1220    if (n==0) n=1;
1221    CFFListIterator J=L;
1222    int j=0;
1223    res = idInit( n ,1);
1224    for ( ; J.hasItem(); J++, j++ )
1225    {
1226      if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
1227        //count_Factors(res,*v,f, j, convFactoryPSingP( J.getItem().factor() );
1228        res->m[j] = convFactoryPSingP( J.getItem().factor() );
1229      #if 0
1230      else if (rField_is_GF())
1231        res->m[j] = convFactoryGFSingGF( J.getItem().factor() );
1232      #endif
1233      else if (rField_is_Extension())     /* Q(a), Fp(a) */
1234      {
1235        if (currRing->minpoly==NULL)
1236          res->m[j]=convFactoryPSingTrP( J.getItem().factor() );
1237        else
1238          res->m[j]=convFactoryAPSingAP( J.getItem().factor(),currRing );
1239      }
1240    }
1241    if (res->m[0]==NULL)
1242    {
1243      res->m[0]=pOne();
1244    }
1245  }
1246  pDelete(&f);
1247  errorreported=save_errorreported;
1248notImpl:
1249  if (res==NULL)
1250    WerrorS( feNotImplemented );
1251  return res;
1252}
1253matrix singclap_irrCharSeries ( ideal I)
1254{
1255  if (idIs0(I)) return mpNew(1,1);
1256
1257  // for now there is only the possibility to handle polynomials over
1258  // Q and Fp ...
1259  matrix res=NULL;
1260  int i;
1261  Off(SW_RATIONAL);
1262  On(SW_SYMMETRIC_FF);
1263  CFList L;
1264  ListCFList LL;
1265  if (((nGetChar() == 0) || (nGetChar() > 1) )
1266  && (currRing->parameter==NULL))
1267  {
1268    setCharacteristic( nGetChar() );
1269    for(i=0;i<IDELEMS(I);i++)
1270    {
1271      poly p=I->m[i];
1272      if (p!=NULL)
1273      {
1274        p=pCopy(p);
1275        p_Cleardenom(p, currRing);
1276        L.append(convSingPFactoryP(p));
1277      }
1278    }
1279  }
1280  // and over Q(a) / Fp(a)
1281  else if (( nGetChar()==1 ) /* Q(a) */
1282  || (nGetChar() <-1))       /* Fp(a) */
1283  {
1284    if (nGetChar()==1) setCharacteristic( 0 );
1285    else               setCharacteristic( -nGetChar() );
1286    for(i=0;i<IDELEMS(I);i++)
1287    {
1288      poly p=I->m[i];
1289      if (p!=NULL)
1290      {
1291        p=pCopy(p);
1292        p_Cleardenom(p, currRing);
1293        L.append(convSingTrPFactoryP(p));
1294      }
1295    }
1296  }
1297  else
1298  {
1299    WerrorS( feNotImplemented );
1300    return res;
1301  }
1302
1303  // a very bad work-around --- FIX IT in libfac
1304  // should be fixed as of 2001/6/27
1305  int tries=0;
1306  int m,n;
1307  ListIterator<CFList> LLi;
1308  loop
1309  {
1310    LL=IrrCharSeries(L);
1311    m= LL.length(); // Anzahl Zeilen
1312    n=0;
1313    for ( LLi = LL; LLi.hasItem(); LLi++ )
1314    {
1315      n = si_max(LLi.getItem().length(),n);
1316    }
1317    if ((m!=0) && (n!=0)) break;
1318    tries++;
1319    if (tries>=5) break;
1320  }
1321  if ((m==0) || (n==0))
1322  {
1323    Warn("char_series returns %d x %d matrix from %d input polys (%d)",
1324      m,n,IDELEMS(I)+1,LL.length());
1325    iiWriteMatrix((matrix)I,"I",2,0);
1326    m=si_max(m,1);
1327    n=si_max(n,1);
1328  }
1329  res=mpNew(m,n);
1330  CFListIterator Li;
1331  for ( m=1, LLi = LL; LLi.hasItem(); LLi++, m++ )
1332  {
1333    for (n=1, Li = LLi.getItem(); Li.hasItem(); Li++, n++)
1334    {
1335      if ( (nGetChar() == 0) || (nGetChar() > 1) )
1336        MATELEM(res,m,n)=convFactoryPSingP(Li.getItem());
1337      else
1338        MATELEM(res,m,n)=convFactoryPSingTrP(Li.getItem());
1339    }
1340  }
1341  Off(SW_RATIONAL);
1342  return res;
1343}
1344
1345char* singclap_neworder ( ideal I)
1346{
1347  int i;
1348  Off(SW_RATIONAL);
1349  On(SW_SYMMETRIC_FF);
1350  CFList L;
1351  if (((nGetChar() == 0) || (nGetChar() > 1) )
1352  && (currRing->parameter==NULL))
1353  {
1354    setCharacteristic( nGetChar() );
1355    for(i=0;i<IDELEMS(I);i++)
1356    {
1357      L.append(convSingPFactoryP(I->m[i]));
1358    }
1359  }
1360  // and over Q(a) / Fp(a)
1361  else if (( nGetChar()==1 ) /* Q(a) */
1362  || (nGetChar() <-1))       /* Fp(a) */
1363  {
1364    if (nGetChar()==1) setCharacteristic( 0 );
1365    else               setCharacteristic( -nGetChar() );
1366    for(i=0;i<IDELEMS(I);i++)
1367    {
1368      L.append(convSingTrPFactoryP(I->m[i]));
1369    }
1370  }
1371  else
1372  {
1373    WerrorS( feNotImplemented );
1374    return NULL;
1375  }
1376
1377  List<int> IL=neworderint(L);
1378  ListIterator<int> Li;
1379  StringSetS("");
1380  Li = IL;
1381  int offs=rPar(currRing);
1382  int* mark=(int*)omAlloc0((pVariables+offs)*sizeof(int));
1383  int cnt=pVariables+offs;
1384  loop
1385  {
1386    if(! Li.hasItem()) break;
1387    BOOLEAN done=TRUE;
1388    i=Li.getItem()-1;
1389    mark[i]=1;
1390    if (i<offs)
1391    {
1392      done=FALSE;
1393      //StringAppendS(currRing->parameter[i]);
1394    }
1395    else
1396    {
1397      StringAppendS(currRing->names[i-offs]);
1398    }
1399    Li++;
1400    cnt--;
1401    if(cnt==0) break;
1402    if (done) StringAppendS(",");
1403  }
1404  for(i=0;i<pVariables+offs;i++)
1405  {
1406    BOOLEAN done=TRUE;
1407    if(mark[i]==0)
1408    {
1409      if (i<offs)
1410      {
1411        done=FALSE;
1412        //StringAppendS(currRing->parameter[i]);
1413      }
1414      else
1415      {
1416        StringAppendS(currRing->names[i-offs]);
1417      }
1418      cnt--;
1419      if(cnt==0) break;
1420      if (done) StringAppendS(",");
1421    }
1422  }
1423  char * s=omStrDup(StringAppendS(""));
1424  if (s[strlen(s)-1]==',') s[strlen(s)-1]='\0';
1425  return s;
1426}
1427
1428BOOLEAN singclap_isSqrFree(poly f)
1429{
1430  BOOLEAN b=FALSE;
1431  Off(SW_RATIONAL);
1432  //  Q / Fp
1433  if (((nGetChar() == 0) || (nGetChar() > 1) )
1434  &&(currRing->parameter==NULL))
1435  {
1436    setCharacteristic( nGetChar() );
1437    CanonicalForm F( convSingPFactoryP( f ) );
1438    if((nGetChar()>1)&&(!F.isUnivariate()))
1439      goto err;
1440    b=(BOOLEAN)isSqrFree(F);
1441  }
1442  // and over Q(a) / Fp(a)
1443  else if (( nGetChar()==1 ) /* Q(a) */
1444  || (nGetChar() <-1))       /* Fp(a) */
1445  {
1446    if (nGetChar()==1) setCharacteristic( 0 );
1447    else               setCharacteristic( -nGetChar() );
1448    //if (currRing->minpoly!=NULL)
1449    //{
1450    //  CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
1451    //                                             currRing->algring);
1452    //  Variable a=rootOf(mipo);
1453    //  CanonicalForm F( convSingAPFactoryAP( f,a ) );
1454    //  ...
1455    //}
1456    //else
1457    {
1458      CanonicalForm F( convSingTrPFactoryP( f ) );
1459      b=(BOOLEAN)isSqrFree(F);
1460    }
1461    Off(SW_RATIONAL);
1462  }
1463  else
1464  {
1465err:
1466    WerrorS( feNotImplemented );
1467  }
1468  return b;
1469}
1470
1471poly singclap_det( const matrix m )
1472{
1473  int r=m->rows();
1474  if (r!=m->cols())
1475  {
1476    Werror("det of %d x %d matrix",r,m->cols());
1477    return NULL;
1478  }
1479  poly res=NULL;
1480  if (( nGetChar() == 0 || nGetChar() > 1 )
1481  && (currRing->parameter==NULL))
1482  {
1483    setCharacteristic( nGetChar() );
1484    CFMatrix M(r,r);
1485    int i,j;
1486    for(i=r;i>0;i--)
1487    {
1488      for(j=r;j>0;j--)
1489      {
1490        M(i,j)=convSingPFactoryP(MATELEM(m,i,j));
1491      }
1492    }
1493    res= convFactoryPSingP( determinant(M,r) ) ;
1494  }
1495  // and over Q(a) / Fp(a)
1496  else if (( nGetChar()==1 ) /* Q(a) */
1497  || (nGetChar() <-1))       /* Fp(a) */
1498  {
1499    if (nGetChar()==1) setCharacteristic( 0 );
1500    else               setCharacteristic( -nGetChar() );
1501    CFMatrix M(r,r);
1502    poly res;
1503    if (currRing->minpoly!=NULL)
1504    {
1505      CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
1506                                            currRing->algring);
1507      Variable a=rootOf(mipo);
1508      int i,j;
1509      for(i=r;i>0;i--)
1510      {
1511        for(j=r;j>0;j--)
1512        {
1513          M(i,j)=convSingAPFactoryAP(MATELEM(m,i,j),a,currRing);
1514        }
1515      }
1516      res= convFactoryAPSingAP( determinant(M,r),currRing ) ;
1517    }
1518    else
1519    {
1520      int i,j;
1521      for(i=r;i>0;i--)
1522      {
1523        for(j=r;j>0;j--)
1524        {
1525          M(i,j)=convSingTrPFactoryP(MATELEM(m,i,j));
1526        }
1527      }
1528      res= convFactoryPSingTrP( determinant(M,r) );
1529    }
1530  }
1531  else
1532    WerrorS( feNotImplemented );
1533  Off(SW_RATIONAL);
1534  return res;
1535}
1536
1537int singclap_det_i( intvec * m )
1538{
1539  setCharacteristic( 0 );
1540  CFMatrix M(m->rows(),m->cols());
1541  int i,j;
1542  for(i=m->rows();i>0;i--)
1543  {
1544    for(j=m->cols();j>0;j--)
1545    {
1546      M(i,j)=IMATELEM(*m,i,j);
1547    }
1548  }
1549  int res= convFactoryISingI( determinant(M,m->rows())) ;
1550  Off(SW_RATIONAL);
1551  return res;
1552}
1553#ifdef HAVE_NTL
1554matrix singntl_HNF(matrix  m )
1555{
1556  int r=m->rows();
1557  if (r!=m->cols())
1558  {
1559    Werror("HNF of %d x %d matrix",r,m->cols());
1560    return NULL;
1561  }
1562  matrix res=mpNew(r,r);
1563  if (rField_is_Q(currRing))
1564  {
1565
1566    CFMatrix M(r,r);
1567    int i,j;
1568    for(i=r;i>0;i--)
1569    {
1570      for(j=r;j>0;j--)
1571      {
1572        M(i,j)=convSingPFactoryP(MATELEM(m,i,j));
1573      }
1574    }
1575    CFMatrix *MM=cf_HNF(M);
1576    for(i=r;i>0;i--)
1577    {
1578      for(j=r;j>0;j--)
1579      {
1580        MATELEM(res,i,j)=convFactoryPSingP((*MM)(i,j));
1581      }
1582    }
1583    delete MM;
1584  }
1585  return res;
1586}
1587intvec* singntl_HNF(intvec*  m )
1588{
1589  int r=m->rows();
1590  if (r!=m->cols())
1591  {
1592    Werror("HNF of %d x %d matrix",r,m->cols());
1593    return NULL;
1594  }
1595  setCharacteristic( 0 );
1596  CFMatrix M(r,r);
1597  int i,j;
1598  for(i=r;i>0;i--)
1599  {
1600    for(j=r;j>0;j--)
1601    {
1602      M(i,j)=IMATELEM(*m,i,j);
1603    }
1604  }
1605  CFMatrix *MM=cf_HNF(M);
1606  intvec *mm=ivCopy(m);
1607  for(i=r;i>0;i--)
1608  {
1609    for(j=r;j>0;j--)
1610    {
1611      IMATELEM(*mm,i,j)=convFactoryISingI((*MM)(i,j));
1612    }
1613  }
1614  delete MM;
1615  return mm;
1616}
1617matrix singntl_LLL(matrix  m )
1618{
1619  int r=m->rows();
1620  int c=m->cols();
1621  matrix res=mpNew(r,c);
1622  if (rField_is_Q(currRing))
1623  {
1624    CFMatrix M(r,c);
1625    int i,j;
1626    for(i=r;i>0;i--)
1627    {
1628      for(j=c;j>0;j--)
1629      {
1630        M(i,j)=convSingPFactoryP(MATELEM(m,i,j));
1631      }
1632    }
1633    CFMatrix *MM=cf_LLL(M);
1634    for(i=r;i>0;i--)
1635    {
1636      for(j=c;j>0;j--)
1637      {
1638        MATELEM(res,i,j)=convFactoryPSingP((*MM)(i,j));
1639      }
1640    }
1641    delete MM;
1642  }
1643  return res;
1644}
1645intvec* singntl_LLL(intvec*  m )
1646{
1647  int r=m->rows();
1648  int c=m->cols();
1649  setCharacteristic( 0 );
1650  CFMatrix M(r,c);
1651  int i,j;
1652  for(i=r;i>0;i--)
1653  {
1654    for(j=r;j>0;j--)
1655    {
1656      M(i,j)=IMATELEM(*m,i,j);
1657    }
1658  }
1659  CFMatrix *MM=cf_LLL(M);
1660  intvec *mm=ivCopy(m);
1661  for(i=r;i>0;i--)
1662  {
1663    for(j=c;j>0;j--)
1664    {
1665      IMATELEM(*mm,i,j)=convFactoryISingI((*MM)(i,j));
1666    }
1667  }
1668  delete MM;
1669  return mm;
1670}
1671#endif
1672
1673napoly singclap_alglcm ( napoly f, napoly g )
1674{
1675
1676 // over Q(a) / Fp(a)
1677 if (nGetChar()==1) setCharacteristic( 0 );
1678 else               setCharacteristic( -nGetChar() );
1679 napoly res;
1680
1681 if (currRing->minpoly!=NULL)
1682 {
1683   CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
1684                                         currRing->algring);
1685   Variable a=rootOf(mipo);
1686   CanonicalForm F( convSingAFactoryA( f,a, currRing ) ),
1687                 G( convSingAFactoryA( g,a, currRing ) );
1688   CanonicalForm GCD;
1689
1690   // calculate gcd
1691   GCD = gcd( F, G );
1692
1693   // calculate lcm
1694   res= convFactoryASingA( (F/GCD)*G,currRing );
1695 }
1696 else
1697 {
1698   CanonicalForm F( convSingPFactoryP( f,currRing->algring ) ),
1699                 G( convSingPFactoryP( g,currRing->algring ) );
1700   CanonicalForm GCD;
1701   // calculate gcd
1702   GCD = gcd( F, G );
1703
1704   // calculate lcm
1705   res= convFactoryPSingP( (F/GCD)*G, currRing->algring );
1706 }
1707
1708 Off(SW_RATIONAL);
1709 return res;
1710}
1711
1712void singclap_algdividecontent ( napoly f, napoly g, napoly &ff, napoly &gg )
1713{
1714 // over Q(a) / Fp(a)
1715 if (nGetChar()==1) setCharacteristic( 0 );
1716 else               setCharacteristic( -nGetChar() );
1717 ff=gg=NULL;
1718 On(SW_RATIONAL);
1719
1720 if (currRing->minpoly!=NULL)
1721 {
1722   CanonicalForm mipo=convSingPFactoryP(((lnumber)currRing->minpoly)->z,
1723                         currRing->algring);
1724   Variable a=rootOf(mipo);
1725   CanonicalForm F( convSingAFactoryA( f,a, currRing ) ),
1726                 G( convSingAFactoryA( g,a, currRing ) );
1727   CanonicalForm GCD;
1728
1729   GCD=gcd( F, G );
1730
1731   if ((GCD!=1) && (GCD!=0))
1732   {
1733     ff= convFactoryASingA( F/ GCD, currRing );
1734     gg= convFactoryASingA( G/ GCD, currRing );
1735   }
1736 }
1737 else
1738 {
1739   CanonicalForm F( convSingPFactoryP( f,currRing->algring ) ),
1740                 G( convSingPFactoryP( g,currRing->algring ) );
1741   CanonicalForm GCD;
1742
1743   GCD=gcd( F, G );
1744
1745   if ((GCD!=1) && (GCD!=0))
1746   {
1747     ff= convFactoryPSingP( F/ GCD, currRing->algring );
1748     gg= convFactoryPSingP( G/ GCD, currRing->algring );
1749   }
1750 }
1751
1752 Off(SW_RATIONAL);
1753}
1754
1755#if 0
1756lists singclap_chineseRemainder(lists x, lists q)
1757{
1758  //assume(x->nr == q->nr);
1759  //assume(x->nr >= 0);
1760  int n=x->nr+1;
1761  if ((x->nr<0) || (x->nr!=q->nr))
1762  {
1763    WerrorS("list are empty or not of equal length");
1764    return NULL;
1765  }
1766  lists res=(lists)omAlloc0Bin(slists_bin);
1767  CFArray X(1,n), Q(1,n);
1768  int i;
1769  for(i=0; i<n; i++)
1770  {
1771    if (x->m[i-1].Typ()==INT_CMD)
1772    {
1773      X[i]=(int)x->m[i-1].Data();
1774    }
1775    else if (x->m[i-1].Typ()==NUMBER_CMD)
1776    {
1777      number N=(number)x->m[i-1].Data();
1778      X[i]=convSingNFactoryN(N);
1779    }
1780    else
1781    {
1782      WerrorS("illegal type in chineseRemainder");
1783      omFreeBin(res,slists_bin);
1784      return NULL;
1785    }
1786    if (q->m[i-1].Typ()==INT_CMD)
1787    {
1788      Q[i]=(int)q->m[i-1].Data();
1789    }
1790    else if (q->m[i-1].Typ()==NUMBER_CMD)
1791    {
1792      number N=(number)x->m[i-1].Data();
1793      Q[i]=convSingNFactoryN(N);
1794    }
1795    else
1796    {
1797      WerrorS("illegal type in chineseRemainder");
1798      omFreeBin(res,slists_bin);
1799      return NULL;
1800    }
1801  }
1802  CanonicalForm r, prod;
1803  chineseRemainder( X, Q, r, prod );
1804  res->Init(2);
1805  res->m[0].rtyp=NUMBER_CMD;
1806  res->m[1].rtyp=NUMBER_CMD;
1807  res->m[0].data=(char *)convFactoryNSingN( r );
1808  res->m[1].data=(char *)convFactoryNSingN( prod );
1809  return res;
1810}
1811#endif
1812#endif
Note: See TracBrowser for help on using the repository browser.