source: git/kernel/clapsing.cc @ 320251

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