source: git/kernel/clapsing.cc @ 7c05ef

spielwiese
Last change on this file since 7c05ef was 7c05ef, checked in by Hans Schönemann <hannes@…>, 18 years ago
*hannes: bug_44 git-svn-id: file:///usr/local/Singular/svn/trunk@9437 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 37.7 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2/****************************************
3*  Computer Algebra System SINGULAR     *
4****************************************/
5// $Id: clapsing.cc,v 1.14 2006-10-02 14:47:51 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=pGetCoeff(T_F_conv);
895      number n_f=pGetCoeff(f);
896      poly n_f_m=pMult_nn(pCopy(f),n_T);
897      T_F_conv=pMult_nn(T_F_conv,n_f);
898      T_F_conv=pSub(T_F_conv,n_f_m);
899      if (T_F_conv!=NULL)
900      {
901        if (singclap_factorize_retry<3)
902        {
903          singclap_factorize_retry++;
904          //if( si_factor_reminder) Print("problem with factorize, retrying\n");
905        //#define FEHLER_FACTORIZE
906        #ifdef FEHLER_FACTORIZE
907          Print("Problem....:");pWrite(f);
908          J=L;
909          for ( ; J.hasItem(); J++ )
910          {
911            if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
912              pWrite0( convClapPSingP( J.getItem().factor() ));
913            else if (rField_is_Extension())     /* Q(a), Fp(a) */
914            {
915              if (currRing->minpoly==NULL)
916                pWrite0( convClapPSingTrP( J.getItem().factor() ));
917              else
918                pWrite0( convClapAPSingAP( J.getItem().factor() ));
919            }
920            Print(" exp: %d\n", J.getItem().exp());
921          }
922          Print("mult:");
923          if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
924            pWrite( convClapPSingP( T ));
925          else if (rField_is_Extension())     /* Q(a), Fp(a) */
926          {
927            if (currRing->minpoly==NULL)
928              pWrite( convClapPSingTrP( T ));
929            else
930              pWrite( convClapAPSingAP( T ));
931          }
932          Print("diff: sing:"); pWrite(T_F_conv);
933          Print("diff: factory:");
934          if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
935            pWrite( convClapPSingP( T_F ));
936          else if (rField_is_Extension())     /* Q(a), Fp(a) */
937          {
938            if (currRing->minpoly==NULL)
939              pWrite( convClapPSingTrP( T_F ));
940            else
941              pWrite( convClapAPSingAP( T_F ));
942          }
943        #endif
944          ideal T_i=singclap_factorize ( f, v , with_exps);
945          if (N!=NULL) nDelete(&N);
946          pDelete(&T_F_conv);
947          return T_i;
948        }
949        else
950        {
951          singclap_factorize_retry=0;
952          WarnS("problem with factorize: irreducibility assumed");
953          ideal T_i=idInit(2,1);
954          T_i->m[0]=pOne();
955          T_i->m[1]=pCopy(f);
956          if (N!=NULL) nDelete(&N);
957          pDelete(&T_F_conv);
958          if (with_exps!=1)
959          {
960            (*v)=new intvec(2);
961            (**v)[0]=1;
962            (**v)[1]=1;
963          }
964          return T_i;
965        }
966      }
967    }
968    J=L;
969    int j=0;
970    if (with_exps!=1)
971    {
972      if ((with_exps==2)&&(n>1))
973      {
974        n--;
975        J++;
976      }
977      *v = new intvec( n );
978    }
979    res = idInit( n ,1);
980    for ( ; J.hasItem(); J++, j++ )
981    {
982      if (with_exps!=1) (**v)[j] = J.getItem().exp();
983      if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
984        res->m[j] = convClapPSingP( J.getItem().factor() );
985      #if 0
986      else if (rField_is_GF())
987        res->m[j] = convClapGFSingGF( J.getItem().factor() );
988      #endif
989      else if (rField_is_Extension())     /* Q(a), Fp(a) */
990      {
991        if (currRing->minpoly==NULL)
992          res->m[j] = convClapPSingTrP( J.getItem().factor() );
993        else
994          res->m[j] = convClapAPSingAP( J.getItem().factor() );
995      }
996    }
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]=0;
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          *v = new intvec( si_max(n-j,1) );
1043          for (i=0,j=0;i<w->length();i++)
1044          {
1045            if((*w)[i]!=0)
1046            {
1047              (**v)[j]=(*w)[i]; j++;
1048            }
1049          }
1050          delete w;
1051        }
1052      }
1053      if (res->m[0]==NULL)
1054      {
1055        res->m[0]=pOne();
1056      }
1057    }
1058  }
1059  if (rField_is_Q_a() && (currRing->minpoly!=NULL))
1060  {
1061    int i=IDELEMS(res)-1;
1062    int stop=1;
1063    if (with_exps!=0) stop=0;
1064    for(;i>=stop;i--)
1065    {
1066      pNorm(res->m[i]);
1067    }
1068    if (with_exps==0) pSetCoeff(res->m[0],old_lead_coeff);
1069    else nDelete(&old_lead_coeff);
1070  }
1071  else
1072    nDelete(&old_lead_coeff);
1073notImpl:
1074  if (res==NULL)
1075    WerrorS( feNotImplemented );
1076  if (NN!=NULL)
1077  {
1078    pMult_nn(f,NN);
1079    nDelete(&NN);
1080  }
1081  if (N!=NULL)
1082  {
1083    nDelete(&N);
1084  }
1085  //PrintS("......S\n");
1086  return res;
1087}
1088matrix singclap_irrCharSeries ( ideal I)
1089{
1090#ifdef HAVE_LIBFAC_P
1091  if (idIs0(I)) return mpNew(1,1);
1092
1093  // for now there is only the possibility to handle polynomials over
1094  // Q and Fp ...
1095  matrix res=NULL;
1096  int i;
1097  Off(SW_RATIONAL);
1098  On(SW_SYMMETRIC_FF);
1099  CFList L;
1100  ListCFList LL;
1101  if (((nGetChar() == 0) || (nGetChar() > 1) )
1102  && (currRing->parameter==NULL))
1103  {
1104    setCharacteristic( nGetChar() );
1105    for(i=0;i<IDELEMS(I);i++)
1106    {
1107      poly p=I->m[i];
1108      if (p!=NULL)
1109      {
1110        p=pCopy(p);
1111        pCleardenom(p);
1112        L.append(convSingPClapP(p));
1113      }
1114    }
1115  }
1116  // and over Q(a) / Fp(a)
1117  else if (( nGetChar()==1 ) /* Q(a) */
1118  || (nGetChar() <-1))       /* Fp(a) */
1119  {
1120    if (nGetChar()==1) setCharacteristic( 0 );
1121    else               setCharacteristic( -nGetChar() );
1122    for(i=0;i<IDELEMS(I);i++)
1123    {
1124      poly p=I->m[i];
1125      if (p!=NULL)
1126      {
1127        p=pCopy(p);
1128        pCleardenom(p);
1129        L.append(convSingTrPClapP(p));
1130      }
1131    }
1132  }
1133  else
1134  {
1135    WerrorS( feNotImplemented );
1136    return res;
1137  }
1138
1139  // a very bad work-around --- FIX IT in libfac
1140  // should be fixed as of 2001/6/27
1141  int tries=0;
1142  int m,n;
1143  ListIterator<CFList> LLi;
1144  loop
1145  {
1146    LL=IrrCharSeries(L);
1147    m= LL.length(); // Anzahl Zeilen
1148    n=0;
1149    for ( LLi = LL; LLi.hasItem(); LLi++ )
1150    {
1151      n = si_max(LLi.getItem().length(),n);
1152    }
1153    if ((m!=0) && (n!=0)) break;
1154    tries++;
1155    if (tries>=5) break;
1156  }
1157  if ((m==0) || (n==0))
1158  {
1159    Warn("char_series returns %d x %d matrix from %d input polys (%d)",
1160      m,n,IDELEMS(I)+1,LL.length());
1161    iiWriteMatrix((matrix)I,"I",2,0);
1162    m=si_max(m,1);
1163    n=si_max(n,1);
1164  }
1165  res=mpNew(m,n);
1166  CFListIterator Li;
1167  for ( m=1, LLi = LL; LLi.hasItem(); LLi++, m++ )
1168  {
1169    for (n=1, Li = LLi.getItem(); Li.hasItem(); Li++, n++)
1170    {
1171      if ( (nGetChar() == 0) || (nGetChar() > 1) )
1172        MATELEM(res,m,n)=convClapPSingP(Li.getItem());
1173      else
1174        MATELEM(res,m,n)=convClapPSingTrP(Li.getItem());
1175    }
1176  }
1177  Off(SW_RATIONAL);
1178  return res;
1179#else
1180  return NULL;
1181#endif
1182}
1183
1184char* singclap_neworder ( ideal I)
1185{
1186#ifdef HAVE_LIBFAC_P
1187  int i;
1188  Off(SW_RATIONAL);
1189  On(SW_SYMMETRIC_FF);
1190  CFList L;
1191  if (((nGetChar() == 0) || (nGetChar() > 1) )
1192  && (currRing->parameter==NULL))
1193  {
1194    setCharacteristic( nGetChar() );
1195    for(i=0;i<IDELEMS(I);i++)
1196    {
1197      L.append(convSingPClapP(I->m[i]));
1198    }
1199  }
1200  // and over Q(a) / Fp(a)
1201  else if (( nGetChar()==1 ) /* Q(a) */
1202  || (nGetChar() <-1))       /* Fp(a) */
1203  {
1204    if (nGetChar()==1) setCharacteristic( 0 );
1205    else               setCharacteristic( -nGetChar() );
1206    for(i=0;i<IDELEMS(I);i++)
1207    {
1208      L.append(convSingTrPClapP(I->m[i]));
1209    }
1210  }
1211  else
1212  {
1213    WerrorS( feNotImplemented );
1214    return NULL;
1215  }
1216
1217  List<int> IL=neworderint(L);
1218  ListIterator<int> Li;
1219  StringSetS("");
1220  Li = IL;
1221  int offs=rPar(currRing);
1222  int* mark=(int*)omAlloc0((pVariables+offs)*sizeof(int));
1223  int cnt=pVariables+offs;
1224  loop
1225  {
1226    if(! Li.hasItem()) break;
1227    BOOLEAN done=TRUE;
1228    i=Li.getItem()-1;
1229    mark[i]=1;
1230    if (i<offs)
1231    {
1232      done=FALSE;
1233      //StringAppendS(currRing->parameter[i]);
1234    }
1235    else
1236    {
1237      StringAppendS(currRing->names[i-offs]);
1238    }
1239    Li++;
1240    cnt--;
1241    if(cnt==0) break;
1242    if (done) StringAppendS(",");
1243  }
1244  for(i=0;i<pVariables+offs;i++)
1245  {
1246    BOOLEAN done=TRUE;
1247    if(mark[i]==0)
1248    {
1249      if (i<offs)
1250      {
1251        done=FALSE;
1252        //StringAppendS(currRing->parameter[i]);
1253      }
1254      else
1255      {
1256        StringAppendS(currRing->names[i-offs]);
1257      }
1258      cnt--;
1259      if(cnt==0) break;
1260      if (done) StringAppendS(",");
1261    }
1262  }
1263  char * s=omStrDup(StringAppendS(""));
1264  if (s[strlen(s)-1]==',') s[strlen(s)-1]='\0';
1265  return s;
1266#else
1267  return NULL;
1268#endif
1269}
1270
1271BOOLEAN singclap_isSqrFree(poly f)
1272{
1273  BOOLEAN b=FALSE;
1274  Off(SW_RATIONAL);
1275  //  Q / Fp
1276  if (((nGetChar() == 0) || (nGetChar() > 1) )
1277  &&(currRing->parameter==NULL))
1278  {
1279    setCharacteristic( nGetChar() );
1280    CanonicalForm F( convSingPClapP( f ) );
1281    if((nGetChar()>1)&&(!F.isUnivariate()))
1282      goto err;
1283    b=(BOOLEAN)isSqrFree(F);
1284  }
1285  // and over Q(a) / Fp(a)
1286  else if (( nGetChar()==1 ) /* Q(a) */
1287  || (nGetChar() <-1))       /* Fp(a) */
1288  {
1289    if (nGetChar()==1) setCharacteristic( 0 );
1290    else               setCharacteristic( -nGetChar() );
1291    //if (currRing->minpoly!=NULL)
1292    //{
1293    //  CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1294    //  Variable a=rootOf(mipo);
1295    //  CanonicalForm F( convSingAPClapAP( f,a ) );
1296    //  ...
1297    //}
1298    //else
1299    {
1300      CanonicalForm F( convSingTrPClapP( f ) );
1301      b=(BOOLEAN)isSqrFree(F);
1302    }
1303    Off(SW_RATIONAL);
1304  }
1305  else
1306  {
1307err:
1308    WerrorS( feNotImplemented );
1309  }
1310  return b;
1311}
1312
1313poly singclap_det( const matrix m )
1314{
1315  int r=m->rows();
1316  if (r!=m->cols())
1317  {
1318    Werror("det of %d x %d matrix",r,m->cols());
1319    return NULL;
1320  }
1321  poly res=NULL;
1322  if (( nGetChar() == 0 || nGetChar() > 1 )
1323  && (currRing->parameter==NULL))
1324  {
1325    setCharacteristic( nGetChar() );
1326    CFMatrix M(r,r);
1327    int i,j;
1328    for(i=r;i>0;i--)
1329    {
1330      for(j=r;j>0;j--)
1331      {
1332        M(i,j)=convSingPClapP(MATELEM(m,i,j));
1333      }
1334    }
1335    res= convClapPSingP( determinant(M,r) ) ;
1336  }
1337  // and over Q(a) / Fp(a)
1338  else if (( nGetChar()==1 ) /* Q(a) */
1339  || (nGetChar() <-1))       /* Fp(a) */
1340  {
1341    if (nGetChar()==1) setCharacteristic( 0 );
1342    else               setCharacteristic( -nGetChar() );
1343    CFMatrix M(r,r);
1344    poly res;
1345    if (currRing->minpoly!=NULL)
1346    {
1347      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1348      Variable a=rootOf(mipo);
1349      int i,j;
1350      for(i=r;i>0;i--)
1351      {
1352        for(j=r;j>0;j--)
1353        {
1354          M(i,j)=convSingAPClapAP(MATELEM(m,i,j),a);
1355        }
1356      }
1357      res= convClapAPSingAP( determinant(M,r) ) ;
1358    }
1359    else
1360    {
1361      int i,j;
1362      for(i=r;i>0;i--)
1363      {
1364        for(j=r;j>0;j--)
1365        {
1366          M(i,j)=convSingTrPClapP(MATELEM(m,i,j));
1367        }
1368      }
1369      res= convClapPSingTrP( determinant(M,r) );
1370    }
1371  }
1372  else
1373    WerrorS( feNotImplemented );
1374  Off(SW_RATIONAL);
1375  return res;
1376}
1377
1378int singclap_det_i( intvec * m )
1379{
1380  setCharacteristic( 0 );
1381  CFMatrix M(m->rows(),m->cols());
1382  int i,j;
1383  for(i=m->rows();i>0;i--)
1384  {
1385    for(j=m->cols();j>0;j--)
1386    {
1387      M(i,j)=IMATELEM(*m,i,j);
1388    }
1389  }
1390  int res= convClapISingI( determinant(M,m->rows())) ;
1391  Off(SW_RATIONAL);
1392  return res;
1393}
1394napoly singclap_alglcm ( napoly f, napoly g )
1395{
1396 FACTORY_ALGOUT( "f", f );
1397 FACTORY_ALGOUT( "g", g );
1398
1399 // over Q(a) / Fp(a)
1400 if (nGetChar()==1) setCharacteristic( 0 );
1401 else               setCharacteristic( -nGetChar() );
1402 napoly res;
1403
1404 if (currRing->minpoly!=NULL)
1405 {
1406   CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1407   Variable a=rootOf(mipo);
1408   CanonicalForm F( convSingAClapA( f,a ) ), G( convSingAClapA( g,a ) );
1409   CanonicalForm GCD;
1410
1411   TIMING_START( algLcmTimer );
1412   // calculate gcd
1413#ifdef FACTORY_GCD_TEST
1414   GCD = CFPrimitiveGcdUtil::gcd( F, G );
1415#else
1416   GCD = gcd( F, G );
1417#endif
1418   TIMING_END( algLcmTimer );
1419
1420   FACTORY_CFAOUT( "d", GCD );
1421   FACTORY_GCDSTAT( "alcm:", F, G, GCD );
1422
1423   // calculate lcm
1424   res= convClapASingA( (F/GCD)*G );
1425 }
1426 else
1427 {
1428   CanonicalForm F( convSingTrClapP( f ) ), G( convSingTrClapP( g ) );
1429   CanonicalForm GCD;
1430   TIMING_START( algLcmTimer );
1431   // calculate gcd
1432#ifdef FACTORY_GCD_TEST
1433   GCD = CFPrimitiveGcdUtil::gcd( F, G );
1434#else
1435   GCD = gcd( F, G );
1436#endif
1437   TIMING_END( algLcmTimer );
1438
1439   FACTORY_CFTROUT( "d", GCD );
1440   FACTORY_GCDSTAT( "alcm:", F, G, GCD );
1441
1442   // calculate lcm
1443   res= convClapPSingTr( (F/GCD)*G );
1444 }
1445
1446 Off(SW_RATIONAL);
1447 return res;
1448}
1449
1450void singclap_algdividecontent ( napoly f, napoly g, napoly &ff, napoly &gg )
1451{
1452 FACTORY_ALGOUT( "f", f );
1453 FACTORY_ALGOUT( "g", g );
1454
1455 // over Q(a) / Fp(a)
1456 if (nGetChar()==1) setCharacteristic( 0 );
1457 else               setCharacteristic( -nGetChar() );
1458 ff=gg=NULL;
1459 On(SW_RATIONAL);
1460
1461 if (currRing->minpoly!=NULL)
1462 {
1463   CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1464   Variable a=rootOf(mipo);
1465   CanonicalForm F( convSingAClapA( f,a ) ), G( convSingAClapA( g,a ) );
1466   CanonicalForm GCD;
1467
1468   TIMING_START( algContentTimer );
1469#ifdef FACTORY_GCD_TEST
1470   GCD=CFPrimitiveGcdUtil::gcd( F, G );
1471#else
1472   GCD=gcd( F, G );
1473#endif
1474   TIMING_END( algContentTimer );
1475
1476   FACTORY_CFAOUT( "d", GCD );
1477   FACTORY_GCDSTAT( "acnt:", F, G, GCD );
1478
1479   if ((GCD!=1) && (GCD!=0))
1480   {
1481     ff= convClapASingA( F/ GCD );
1482     gg= convClapASingA( G/ GCD );
1483   }
1484 }
1485 else
1486 {
1487   CanonicalForm F( convSingTrClapP( f ) ), G( convSingTrClapP( g ) );
1488   CanonicalForm GCD;
1489
1490   TIMING_START( algContentTimer );
1491#ifdef FACTORY_GCD_TEST
1492   GCD=CFPrimitiveGcdUtil::gcd( F, G );
1493#else
1494   GCD=gcd( F, G );
1495#endif
1496   TIMING_END( algContentTimer );
1497
1498   FACTORY_CFTROUT( "d", GCD );
1499   FACTORY_GCDSTAT( "acnt:", F, G, GCD );
1500
1501   if ((GCD!=1) && (GCD!=0))
1502   {
1503     ff= convClapPSingTr( F/ GCD );
1504     gg= convClapPSingTr( G/ GCD );
1505   }
1506 }
1507
1508 Off(SW_RATIONAL);
1509}
1510
1511#if 0
1512lists singclap_chineseRemainder(lists x, lists q)
1513{
1514  //assume(x->nr == q->nr);
1515  //assume(x->nr >= 0);
1516  int n=x->nr+1;
1517  if ((x->nr<0) || (x->nr!=q->nr))
1518  {
1519    WerrorS("list are empty or not of equal length");
1520    return NULL;
1521  }
1522  lists res=(lists)omAlloc0Bin(slists_bin);
1523  CFArray X(1,n), Q(1,n);
1524  int i;
1525  for(i=0; i<n; i++)
1526  {
1527    if (x->m[i-1].Typ()==INT_CMD)
1528    {
1529      X[i]=(int)x->m[i-1].Data();
1530    }
1531    else if (x->m[i-1].Typ()==NUMBER_CMD)
1532    {
1533      number N=(number)x->m[i-1].Data();
1534      X[i]=convSingNClapN(N);
1535    }
1536    else
1537    {
1538      WerrorS("illegal type in chineseRemainder");
1539      omFreeBin(res,slists_bin);
1540      return NULL;
1541    }
1542    if (q->m[i-1].Typ()==INT_CMD)
1543    {
1544      Q[i]=(int)q->m[i-1].Data();
1545    }
1546    else if (q->m[i-1].Typ()==NUMBER_CMD)
1547    {
1548      number N=(number)x->m[i-1].Data();
1549      Q[i]=convSingNClapN(N);
1550    }
1551    else
1552    {
1553      WerrorS("illegal type in chineseRemainder");
1554      omFreeBin(res,slists_bin);
1555      return NULL;
1556    }
1557  }
1558  CanonicalForm r, prod;
1559  chineseRemainder( X, Q, r, prod );
1560  res->Init(2);
1561  res->m[0].rtyp=NUMBER_CMD;
1562  res->m[1].rtyp=NUMBER_CMD;
1563  res->m[0].data=(char *)convClapNSingN( r );
1564  res->m[1].data=(char *)convClapNSingN( prod );
1565  return res;
1566}
1567#endif
1568#endif
Note: See TracBrowser for help on using the repository browser.