source: git/kernel/clapsing.cc @ 599326

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