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

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