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

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