source: git/kernel/clapsing.cc @ d1f231

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