source: git/kernel/clapsing.cc @ d1eca41

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