source: git/kernel/clapsing.cc @ 6f6c1b

spielwiese
Last change on this file since 6f6c1b was 6f6c1b, checked in by Hans Schönemann <hannes@…>, 18 years ago
*hannes: fixed case: not LIBFAC git-svn-id: file:///usr/local/Singular/svn/trunk@9031 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 37.6 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2/****************************************
3*  Computer Algebra System SINGULAR     *
4****************************************/
5// $Id: clapsing.cc,v 1.12 2006-03-21 15:27:20 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  if (( nGetChar() == 0 || nGetChar() > 1 )
212  && (currRing->parameter==NULL))
213  {
214    setCharacteristic( nGetChar() );
215    CanonicalForm F( convSingPClapP( f ) ), G( convSingPClapP( g ) );
216    res=convClapPSingP( gcd( F, G ) );
217    Off(SW_RATIONAL);
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    Off(SW_RATIONAL);
257  }
258  #if 0
259  else if (( nGetChar()>1 )&&(currRing->parameter!=NULL)) /* GF(q) */
260  {
261    int p=rChar(currRing);
262    int n=2;
263    int t=p*p;
264    while (t!=nChar) { t*=p;n++; }
265    setCharacteristic(p,n,'a');
266    CanonicalForm F( convSingGFClapGF( f ) ), G( convSingGFClapGF( g ) );
267    res= convClapGFSingGF( gcd( F, G ) );
268  }
269  #endif
270  else
271    WerrorS( feNotImplemented );
272
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    if (currRing->minpoly!=NULL)
323    {
324      Variable X(i);
325      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
326      Variable a=rootOf(mipo);
327      CanonicalForm F( convSingAPClapAP( f,a ) ), G( convSingAPClapAP( g,a ) );
328      res= convClapAPSingAP( resultant( F, G, X ) );
329    }
330    else
331    {
332      Variable X(i+rPar(currRing));
333      number nf,ng;
334      pCleardenom_n(f,nf);pCleardenom_n(g,ng);
335      int ef,eg;
336      ef=pGetExp_Var(f,i);
337      eg=pGetExp_Var(g,i);
338      CanonicalForm F( convSingTrPClapP( f ) ), G( convSingTrPClapP( g ) );
339      res= convClapPSingTrP( resultant( F, G, X ) );
340      if ((nf!=NULL)&&(!nIsOne(nf))&&(!nIsZero(nf)))
341      {
342        number n=nInvers(nf);
343        while(eg>0)
344        {
345          res=pMult_nn(res,n);
346          eg--;
347        }
348        nDelete(&n);
349      }
350      nDelete(&nf);
351      if ((ng!=NULL)&&(!nIsOne(ng))&&(!nIsZero(ng)))
352      {
353        number n=nInvers(ng);
354        while(ef>0)
355        {
356          res=pMult_nn(res,n);
357          ef--;
358        }
359        nDelete(&n);
360      }
361      nDelete(&ng);
362    }
363    Off(SW_RATIONAL);
364    return res;
365  }
366  else
367    WerrorS( feNotImplemented );
368  return NULL;
369}
370//poly singclap_resultant ( poly f, poly g , poly x)
371//{
372//  int i=pVar(x);
373//  if (i==0)
374//  {
375//    WerrorS("ringvar expected");
376//    return NULL;
377//  }
378//  ideal I=idInit(1,1);
379//
380//  // get the coeffs von f wrt. x:
381//  I->m[0]=pCopy(f);
382//  matrix ffi=mpCoeffs(I,i);
383//  ffi->rank=1;
384//  ffi->ncols=ffi->nrows;
385//  ffi->nrows=1;
386//  ideal fi=(ideal)ffi;
387//
388//  // get the coeffs von g wrt. x:
389//  I->m[0]=pCopy(g);
390//  matrix ggi=mpCoeffs(I,i);
391//  ggi->rank=1;
392//  ggi->ncols=ggi->nrows;
393//  ggi->nrows=1;
394//  ideal gi=(ideal)ggi;
395//
396//  // contruct the matrix:
397//  int fn=IDELEMS(fi); //= deg(f,x)+1
398//  int gn=IDELEMS(gi); //= deg(g,x)+1
399//  matrix m=mpNew(fn+gn-2,fn+gn-2);
400//  if(m==NULL)
401//  {
402//    return NULL;
403//  }
404//
405//  // enter the coeffs into m:
406//  int j;
407//  for(i=0;i<gn-1;i++)
408//  {
409//    for(j=0;j<fn;j++)
410//    {
411//      MATELEM(m,i+1,fn-j+i)=pCopy(fi->m[j]);
412//    }
413//  }
414//  for(i=0;i<fn-1;i++)
415//  {
416//    for(j=0;j<gn;j++)
417//    {
418//      MATELEM(m,gn+i,gn-j+i)=pCopy(gi->m[j]);
419//    }
420//  }
421//
422//  poly r=mpDet(m);
423//
424//  idDelete(&fi);
425//  idDelete(&gi);
426//  idDelete((ideal *)&m);
427//  return r;
428//}
429
430BOOLEAN singclap_extgcd ( poly f, poly g, poly &res, poly &pa, poly &pb )
431{
432  // for now there is only the possibility to handle univariate
433  // polynomials over
434  // Q and Fp ...
435  res=NULL;pa=NULL;pb=NULL;
436  On(SW_SYMMETRIC_FF);
437  if (( nGetChar() == 0 || nGetChar() > 1 )
438  && (currRing->parameter==NULL))
439  {
440    setCharacteristic( nGetChar() );
441    CanonicalForm F( convSingPClapP( f ) ), G( convSingPClapP( g ) );
442    CanonicalForm FpG=F+G;
443    if (!(FpG.isUnivariate()|| FpG.inCoeffDomain()))
444    //if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
445    {
446      Off(SW_RATIONAL);
447      WerrorS("not univariate");
448      return TRUE;
449    }
450    CanonicalForm Fa,Gb;
451    On(SW_RATIONAL);
452    res=convClapPSingP( extgcd( F, G, Fa, Gb ) );
453    pa=convClapPSingP(Fa);
454    pb=convClapPSingP(Gb);
455    Off(SW_RATIONAL);
456  }
457  // and over Q(a) / Fp(a)
458  else if (( nGetChar()==1 ) /* Q(a) */
459  || (nGetChar() <-1))       /* Fp(a) */
460  {
461    if (nGetChar()==1) setCharacteristic( 0 );
462    else               setCharacteristic( -nGetChar() );
463    CanonicalForm Fa,Gb;
464    if (currRing->minpoly!=NULL)
465    {
466      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
467      Variable a=rootOf(mipo);
468      CanonicalForm F( convSingAPClapAP( f,a ) ), G( convSingAPClapAP( g,a ) );
469      CanonicalForm FpG=F+G;
470      if (!(FpG.isUnivariate()|| FpG.inCoeffDomain()))
471      //if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
472      {
473        WerrorS("not univariate");
474        return TRUE;
475      }
476      res= convClapAPSingAP( extgcd( F, G, Fa, Gb ) );
477      pa=convClapAPSingAP(Fa);
478      pb=convClapAPSingAP(Gb);
479    }
480    else
481    {
482      CanonicalForm F( convSingTrPClapP( f ) ), G( convSingTrPClapP( g ) );
483      CanonicalForm FpG=F+G;
484      if (!(FpG.isUnivariate()|| FpG.inCoeffDomain()))
485      //if (!F.isUnivariate() || !G.isUnivariate() || F.mvar()!=G.mvar())
486      {
487        Off(SW_RATIONAL);
488        WerrorS("not univariate");
489        return TRUE;
490      }
491      res= convClapPSingTrP( extgcd( F, G, Fa, Gb ) );
492      pa=convClapPSingTrP(Fa);
493      pb=convClapPSingTrP(Gb);
494    }
495    Off(SW_RATIONAL);
496  }
497  else
498  {
499    WerrorS( feNotImplemented );
500    return TRUE;
501  }
502  return FALSE;
503}
504
505poly singclap_pdivide ( poly f, poly g )
506{
507  // for now there is only the possibility to handle polynomials over
508  // Q and Fp ...
509  poly res=NULL;
510  On(SW_RATIONAL);
511  if (( nGetChar() == 0 || nGetChar() > 1 )
512  && (currRing->parameter==NULL))
513  {
514    setCharacteristic( nGetChar() );
515    CanonicalForm F( convSingPClapP( f ) ), G( convSingPClapP( g ) );
516    res = convClapPSingP( F / G );
517  }
518  // and over Q(a) / Fp(a)
519  else if (( nGetChar()==1 ) /* Q(a) */
520  || (nGetChar() <-1))       /* Fp(a) */
521  {
522    if (nGetChar()==1) setCharacteristic( 0 );
523    else               setCharacteristic( -nGetChar() );
524    if (currRing->minpoly!=NULL)
525    {
526      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
527      Variable a=rootOf(mipo);
528      CanonicalForm F( convSingAPClapAP( f,a ) ), G( convSingAPClapAP( g,a ) );
529      res= convClapAPSingAP(  F / G  );
530    }
531    else
532    {
533      CanonicalForm F( convSingTrPClapP( f ) ), G( convSingTrPClapP( g ) );
534      res= convClapPSingTrP(  F / G  );
535    }
536  }
537  else
538    WerrorS( feNotImplemented );
539  Off(SW_RATIONAL);
540  return res;
541}
542
543void singclap_divide_content ( poly f )
544{
545  if ( f==NULL )
546  {
547    return;
548  }
549  else  if ( pNext( f ) == NULL )
550  {
551    pSetCoeff( f, nInit( 1 ) );
552    return;
553  }
554  else
555  {
556    if ( nGetChar() == 1 )
557      setCharacteristic( 0 );
558    else  if ( nGetChar() == -1 )
559      return; /* not implemented for R */
560    else  if ( nGetChar() < 0 )
561      setCharacteristic( -nGetChar() );
562    else if (currRing->parameter==NULL) /* not GF(q) */
563      setCharacteristic( nGetChar() );
564    else
565      return; /* not implemented*/
566
567    CFList L;
568    CanonicalForm g, h;
569    poly p = pNext(f);
570
571    // first attemp: find 2 smallest g:
572
573    number g1=pGetCoeff(f);
574    number g2=pGetCoeff(p); // p==pNext(f);
575    pIter(p);
576    int sz1=nSize(g1);
577    int sz2=nSize(g2);
578    if (sz1>sz2)
579    {
580      number gg=g1;
581      g1=g2; g2=gg;
582      int sz=sz1;
583      sz1=sz2; sz2=sz;
584    }
585    while (p!=NULL)
586    {
587      int n_sz=nSize(pGetCoeff(p));
588      if (n_sz<sz1)
589      {
590        sz2=sz1;
591        g2=g1;
592        g1=pGetCoeff(p);
593        sz1=n_sz;
594        if (sz1<=3) break;
595      }
596      else if(n_sz<sz2)
597      {
598        sz2=n_sz;
599        g2=pGetCoeff(p);
600        sz2=n_sz;
601      }
602      pIter(p);
603    }
604    FACTORY_ALGOUT( "G", ((lnumber)g1)->z );
605    g = convSingTrClapP( ((lnumber)g1)->z );
606    g = gcd( g, convSingTrClapP( ((lnumber)g2)->z ));
607
608    // second run: gcd's
609
610    p = f;
611    TIMING_START( contentTimer );
612    while ( (p != NULL) && (g != 1)  && ( g != 0))
613    {
614      FACTORY_ALGOUT( "h", (((lnumber)pGetCoeff(p))->z) );
615      h = convSingTrClapP( ((lnumber)pGetCoeff(p))->z );
616      pIter( p );
617#ifdef FACTORY_GCD_STAT
618      // save g
619      CanonicalForm gOld = g;
620#endif
621
622#ifdef FACTORY_GCD_TEST
623      g = CFPrimitiveGcdUtil::gcd( g, h );
624#else
625      g = gcd( g, h );
626#endif
627
628      FACTORY_GCDSTAT( "gcnt:", gOld, h, g );
629      FACTORY_CFTROUT( "g", g );
630      L.append( h );
631    }
632    TIMING_END( contentTimer );
633    FACTORY_CONTSTAT( "cont:", g );
634    if (( g == 1 ) || (g == 0))
635    {
636      // pTest(f);
637      return;
638    }
639    else
640    {
641      CFListIterator i;
642      for ( i = L, p = f; i.hasItem(); i++, p=pNext(p) )
643      {
644        lnumber c=(lnumber)pGetCoeff(p);
645        napDelete(&c->z);
646        c->z=convClapPSingTr( i.getItem() / g );
647        //nTest((number)c);
648        //#ifdef LDEBUG
649        //number cn=(number)c;
650        //StringSetS(""); nWrite(nt); StringAppend(" ==> ");
651        //nWrite(cn);PrintS(StringAppend("\n"));
652        //#endif
653      }
654    }
655    // pTest(f);
656  }
657}
658
659static int primepower(int c)
660{
661  int p=1;
662  int cc=c;
663  while(cc!= rInternalChar(currRing)) { cc*=c; p++; }
664  return p;
665}
666
667int singclap_factorize_retry;
668//extern int si_factor_reminder;
669
670ideal singclap_factorize ( poly f, intvec ** v , int with_exps)
671{
672  // with_exps: 3,1 return only true factors, no exponents
673  //            2 return true factors and exponents
674  //            0 return coeff, factors and exponents
675
676
677  ideal res=NULL;
678
679  // handle factorize(0) =========================================
680  if (f==NULL)
681  {
682    res=idInit(1,1);
683    if (with_exps!=1)
684    {
685      (*v)=new intvec(1);
686      (**v)[0]=1;
687    }
688    return res;
689  }
690  // handle factorize(mon) =========================================
691  if (pNext(f)==NULL)
692  {
693    int i=0;
694    int n=0;
695    int e;
696    for(i=pVariables;i>0;i--) if(pGetExp(f,i)!=0) n++;
697    if (with_exps==0) n++; // with coeff
698    res=idInit(si_max(n,1),1);
699    switch(with_exps)
700    {
701      case 0: // with coef & exp.
702        res->m[0]=pOne();
703        pSetCoeff(res->m[0],nCopy(pGetCoeff(f)));
704        // no break
705      case 2: // with exp.
706        (*v)=new intvec(si_max(1,n));
707        (**v)[0]=1;
708        // no break
709      case 1: ;
710      #ifdef TEST
711      default: ;
712      #endif
713    }
714    if (n==0)
715    {
716      res->m[0]=pOne();
717      // (**v)[0]=1; is already done
718      return res;
719    }
720    for(i=pVariables;i>0;i--)
721    {
722      e=pGetExp(f,i);
723      if(e!=0)
724      {
725        n--;
726        poly p=pOne();
727        pSetExp(p,i,1);
728        pSetm(p);
729        res->m[n]=p;
730        if (with_exps!=1) (**v)[n]=e;
731      }
732    }
733    return res;
734  }
735  //PrintS("S:");pWrite(f);PrintLn();
736  // use factory/libfac in general ==============================
737  Off(SW_RATIONAL);
738  On(SW_SYMMETRIC_FF);
739  #ifdef HAVE_NTL
740  extern int prime_number;
741  if(rField_is_Q()) prime_number=0;
742  #endif
743  CFFList L;
744  number N=NULL;
745  number NN=NULL;
746  CanonicalForm T_F(0);
747  number old_lead_coeff=nCopy(pGetCoeff(f));
748
749  if (!rField_is_Zp()) /* Q, Q(a), Zp(a) */
750  {
751    //if (f!=NULL) // already tested at start of routine
752    {
753      number n0=nCopy(pGetCoeff(f));
754      if (with_exps==0)
755        N=nCopy(n0);
756      pCleardenom(f);
757      NN=nDiv(n0,pGetCoeff(f));
758      nDelete(&n0);
759      if (with_exps==0)
760      {
761        nDelete(&N);
762        N=nCopy(NN);
763      }
764    }
765  }
766  if (rField_is_Q() || rField_is_Zp())
767  {
768    setCharacteristic( nGetChar() );
769    CanonicalForm F( convSingPClapP( f ) );
770    T_F=F;
771    if (nGetChar()==0) /* Q */
772    {
773      L = factorize( F );
774    }
775    else /* Fp */
776    {
777#ifdef HAVE_LIBFAC_P
778      L = Factorize( F );
779#else
780      goto notImpl;
781#endif
782    }
783  }
784  #if 0
785  else if (rField_is_GF())
786  {
787    int c=rChar(currRing);
788    setCharacteristic( c, primepower(c) );
789    CanonicalForm F( convSingGFClapGF( f ) );
790    T_F=F;
791    if (F.isUnivariate())
792    {
793      L = factorize( F );
794    }
795    else
796    {
797      goto notImpl;
798    }
799  }
800  #endif
801  // and over Q(a) / Fp(a)
802  else if (rField_is_Extension())
803  {
804    if (rField_is_Q_a()) setCharacteristic( 0 );
805    else                 setCharacteristic( -nGetChar() );
806    if (currRing->minpoly!=NULL)
807    {
808      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
809      Variable a=rootOf(mipo);
810      CanonicalForm F( convSingAPClapAP( f,a ) );
811      T_F=F;
812      L.insert(F);
813      if (rField_is_Zp_a() && F.isUnivariate())
814      {
815        L = factorize( F, a );
816      }
817      else
818      {
819        CanonicalForm G( convSingTrPClapP( f ) );
820        T_F=G;
821#ifdef HAVE_LIBFAC_P
822        //  over Q(a) / Fp(a)
823        L=Factorize(G, mipo);
824#else
825        WarnS("complete factorization only for univariate polynomials");
826        if (rField_is_Q_a() ||(!F.isUnivariate())) /* Q(a) */
827        {
828          L = factorize( G );
829        }
830        else
831        {
832          L = factorize( G, a );
833        }
834#endif
835      }
836    }
837    else
838    {
839      CanonicalForm F( convSingTrPClapP( f ) );
840      T_F=F;
841      if (rField_is_Q_a())
842      {
843        L = factorize( F );
844      }
845      else /* Fp(a) */
846      {
847#ifdef HAVE_LIBFAC_P
848        L = Factorize( F );
849#else
850        goto notImpl;
851#endif
852      }
853    }
854  }
855  else
856  {
857    goto notImpl;
858  }
859  {
860    // the first factor should be a constant
861    if ( ! L.getFirst().factor().inCoeffDomain() )
862      L.insert(CFFactor(1,1));
863    // convert into ideal
864    int n = L.length();
865    CFFListIterator J=L;
866    CanonicalForm T=1;
867    for ( ; J.hasItem(); J++ )
868    {
869      int T_e = J.getItem().exp();
870      while(T_e>0)  { T *= J.getItem().factor(); T_e--; }
871    }
872    T_F-=T;
873    if (!T_F.isZero())
874    {
875      poly T_F_conv=pOne();
876      J=L;
877      for ( ; J.hasItem(); J++ )
878      {
879        poly p;
880        int T_e = J.getItem().exp();
881        if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
882          p=( convClapPSingP( J.getItem().factor() ));
883        else if (rField_is_Extension())     /* Q(a), Fp(a) */
884        {
885          if (currRing->minpoly==NULL)
886            p=( convClapPSingTrP( J.getItem().factor() ));
887          else
888            p=( convClapAPSingAP( J.getItem().factor() ));
889        }
890        while(T_e>0)  { T_F_conv=pMult(T_F_conv,pCopy(p)); T_e--; }
891        pDelete(&p);
892      }
893      number n_T=pGetCoeff(T_F_conv);
894      number n_f=pGetCoeff(f);
895      poly n_f_m=pMult_nn(pCopy(f),n_T);
896      T_F_conv=pMult_nn(T_F_conv,n_f);
897      T_F_conv=pSub(T_F_conv,n_f_m);
898      if (T_F_conv!=NULL)
899      {
900        if (singclap_factorize_retry<3)
901        {
902          singclap_factorize_retry++;
903          //if( si_factor_reminder) Print("problem with factorize, retrying\n");
904        #ifdef FEHLER_FACTORIZE
905          Print("Problem....:");pWrite(f);
906          J=L;
907          for ( ; J.hasItem(); J++ )
908          {
909            if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
910              pWrite0( convClapPSingP( J.getItem().factor() ));
911            else if (rField_is_Extension())     /* Q(a), Fp(a) */
912            {
913              if (currRing->minpoly==NULL)
914                pWrite0( convClapPSingTrP( J.getItem().factor() ));
915              else
916                pWrite0( convClapAPSingAP( J.getItem().factor() ));
917            }
918            Print(" exp: %d\n", J.getItem().exp());
919          }
920          Print("mult:");
921          if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
922            pWrite( convClapPSingP( T ));
923          else if (rField_is_Extension())     /* Q(a), Fp(a) */
924          {
925            if (currRing->minpoly==NULL)
926              pWrite( convClapPSingTrP( T ));
927            else
928              pWrite( convClapAPSingAP( T ));
929          }
930          Print("diff: sing:"); pWrite(T_F_conv);
931          Print("diff: factory:");
932          if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
933            pWrite( convClapPSingP( T_F ));
934          else if (rField_is_Extension())     /* Q(a), Fp(a) */
935          {
936            if (currRing->minpoly==NULL)
937              pWrite( convClapPSingTrP( T_F ));
938            else
939              pWrite( convClapAPSingAP( T_F ));
940          }
941        #endif
942          ideal T_i=singclap_factorize ( f, v , with_exps);
943          if (N!=NULL) nDelete(&N);
944          pDelete(&T_F_conv);
945          return T_i;
946        }
947        else
948        {
949          singclap_factorize_retry=0;
950          WarnS("problem with factorize: irreducibility assumed");
951          ideal T_i=idInit(2,1);
952          T_i->m[0]=pOne();
953          T_i->m[1]=pCopy(f);
954          if (N!=NULL) nDelete(&N);
955          pDelete(&T_F_conv);
956          if (with_exps!=1)
957          {
958            (*v)=new intvec(2);
959            (**v)[0]=1;
960            (**v)[1]=1;
961          }
962          return T_i;
963        }
964      }
965    }
966    J=L;
967    int j=0;
968    if (with_exps!=1)
969    {
970      if ((with_exps==2)&&(n>1))
971      {
972        n--;
973        J++;
974      }
975      *v = new intvec( n );
976    }
977    res = idInit( n ,1);
978    for ( ; J.hasItem(); J++, j++ )
979    {
980      if (with_exps!=1) (**v)[j] = J.getItem().exp();
981      if (rField_is_Zp() || rField_is_Q())           /* Q, Fp */
982        res->m[j] = convClapPSingP( J.getItem().factor() );
983      #if 0
984      else if (rField_is_GF())
985        res->m[j] = convClapGFSingGF( J.getItem().factor() );
986      #endif
987      else if (rField_is_Extension())     /* Q(a), Fp(a) */
988      {
989        if (currRing->minpoly==NULL)
990          res->m[j] = convClapPSingTrP( J.getItem().factor() );
991        else
992          res->m[j] = convClapAPSingAP( J.getItem().factor() );
993      }
994    }
995    if (N!=NULL)
996    {
997      pMult_nn(res->m[0],N);
998      nDelete(&N);
999      N=NULL;
1000    }
1001    // delete constants
1002    if (res!=NULL)
1003    {
1004      int i=IDELEMS(res)-1;
1005      int j=0;
1006      for(;i>=0;i--)
1007      {
1008        if ((res->m[i]!=NULL)
1009        && (pNext(res->m[i])==NULL)
1010        && (pIsConstant(res->m[i])))
1011        {
1012          if (with_exps!=0)
1013          {
1014            pDelete(&(res->m[i]));
1015            if ((v!=NULL) && ((*v)!=NULL))
1016              (**v)[i]=0;
1017            j++;
1018          }
1019          else if (i!=0)
1020          {
1021            while ((v!=NULL) && ((*v)!=NULL) && ((**v)[i]>1))
1022            {
1023              res->m[0]=pMult(res->m[0],pCopy(res->m[i]));
1024              (**v)[i]--;
1025            }
1026            res->m[0]=pMult(res->m[0],res->m[i]);
1027            res->m[i]=NULL;
1028            if ((v!=NULL) && ((*v)!=NULL))
1029              (**v)[i]=0;
1030            j++;
1031          }
1032        }
1033      }
1034      if (j>0)
1035      {
1036        idSkipZeroes(res);
1037        if ((v!=NULL) && ((*v)!=NULL))
1038        {
1039          intvec *w=*v;
1040          *v = new intvec( si_max(n-j,1) );
1041          for (i=0,j=0;i<w->length();i++)
1042          {
1043            if((*w)[i]!=0)
1044            {
1045              (**v)[j]=(*w)[i]; j++;
1046            }
1047          }
1048          delete w;
1049        }
1050      }
1051      if (res->m[0]==NULL)
1052      {
1053        res->m[0]=pOne();
1054      }
1055    }
1056  }
1057  if (rField_is_Q_a() && (currRing->minpoly!=NULL))
1058  {
1059    int i=IDELEMS(res)-1;
1060    int stop=1;
1061    if (with_exps!=0) stop=0;
1062    for(;i>=stop;i--)
1063    {
1064      pNorm(res->m[i]);
1065    }
1066    if (with_exps==0) pSetCoeff(res->m[0],old_lead_coeff);
1067    else nDelete(&old_lead_coeff);
1068  }
1069  else
1070    nDelete(&old_lead_coeff);
1071notImpl:
1072  if (res==NULL)
1073    WerrorS( feNotImplemented );
1074  if (NN!=NULL)
1075  {
1076    pMult_nn(f,NN);
1077    nDelete(&NN);
1078  }
1079  if (N!=NULL)
1080  {
1081    nDelete(&N);
1082  }
1083  //PrintS("......S\n");
1084  return res;
1085}
1086matrix singclap_irrCharSeries ( ideal I)
1087{
1088#ifdef HAVE_LIBFAC_P
1089  if (idIs0(I)) return mpNew(1,1);
1090
1091  // for now there is only the possibility to handle polynomials over
1092  // Q and Fp ...
1093  matrix res=NULL;
1094  int i;
1095  Off(SW_RATIONAL);
1096  On(SW_SYMMETRIC_FF);
1097  CFList L;
1098  ListCFList LL;
1099  if (((nGetChar() == 0) || (nGetChar() > 1) )
1100  && (currRing->parameter==NULL))
1101  {
1102    setCharacteristic( nGetChar() );
1103    for(i=0;i<IDELEMS(I);i++)
1104    {
1105      poly p=I->m[i];
1106      if (p!=NULL)
1107      {
1108        p=pCopy(p);
1109        pCleardenom(p);
1110        L.append(convSingPClapP(p));
1111      }
1112    }
1113  }
1114  // and over Q(a) / Fp(a)
1115  else if (( nGetChar()==1 ) /* Q(a) */
1116  || (nGetChar() <-1))       /* Fp(a) */
1117  {
1118    if (nGetChar()==1) setCharacteristic( 0 );
1119    else               setCharacteristic( -nGetChar() );
1120    for(i=0;i<IDELEMS(I);i++)
1121    {
1122      poly p=I->m[i];
1123      if (p!=NULL)
1124      {
1125        p=pCopy(p);
1126        pCleardenom(p);
1127        L.append(convSingTrPClapP(p));
1128      }
1129    }
1130  }
1131  else
1132  {
1133    WerrorS( feNotImplemented );
1134    return res;
1135  }
1136
1137  // a very bad work-around --- FIX IT in libfac
1138  // should be fixed as of 2001/6/27
1139  int tries=0;
1140  int m,n;
1141  ListIterator<CFList> LLi;
1142  loop
1143  {
1144    LL=IrrCharSeries(L);
1145    m= LL.length(); // Anzahl Zeilen
1146    n=0;
1147    for ( LLi = LL; LLi.hasItem(); LLi++ )
1148    {
1149      n = si_max(LLi.getItem().length(),n);
1150    }
1151    if ((m!=0) && (n!=0)) break;
1152    tries++;
1153    if (tries>=5) break;
1154  }
1155  if ((m==0) || (n==0))
1156  {
1157    Warn("char_series returns %d x %d matrix from %d input polys (%d)",
1158      m,n,IDELEMS(I)+1,LL.length());
1159    iiWriteMatrix((matrix)I,"I",2,0);
1160    m=si_max(m,1);
1161    n=si_max(n,1);
1162  }
1163  res=mpNew(m,n);
1164  CFListIterator Li;
1165  for ( m=1, LLi = LL; LLi.hasItem(); LLi++, m++ )
1166  {
1167    for (n=1, Li = LLi.getItem(); Li.hasItem(); Li++, n++)
1168    {
1169      if ( (nGetChar() == 0) || (nGetChar() > 1) )
1170        MATELEM(res,m,n)=convClapPSingP(Li.getItem());
1171      else
1172        MATELEM(res,m,n)=convClapPSingTrP(Li.getItem());
1173    }
1174  }
1175  Off(SW_RATIONAL);
1176  return res;
1177#else
1178  return NULL;
1179#endif
1180}
1181
1182char* singclap_neworder ( ideal I)
1183{
1184#ifdef HAVE_LIBFAC_P
1185  int i;
1186  Off(SW_RATIONAL);
1187  On(SW_SYMMETRIC_FF);
1188  CFList L;
1189  if (((nGetChar() == 0) || (nGetChar() > 1) )
1190  && (currRing->parameter==NULL))
1191  {
1192    setCharacteristic( nGetChar() );
1193    for(i=0;i<IDELEMS(I);i++)
1194    {
1195      L.append(convSingPClapP(I->m[i]));
1196    }
1197  }
1198  // and over Q(a) / Fp(a)
1199  else if (( nGetChar()==1 ) /* Q(a) */
1200  || (nGetChar() <-1))       /* Fp(a) */
1201  {
1202    if (nGetChar()==1) setCharacteristic( 0 );
1203    else               setCharacteristic( -nGetChar() );
1204    for(i=0;i<IDELEMS(I);i++)
1205    {
1206      L.append(convSingTrPClapP(I->m[i]));
1207    }
1208  }
1209  else
1210  {
1211    WerrorS( feNotImplemented );
1212    return NULL;
1213  }
1214
1215  List<int> IL=neworderint(L);
1216  ListIterator<int> Li;
1217  StringSetS("");
1218  Li = IL;
1219  int offs=rPar(currRing);
1220  int* mark=(int*)omAlloc0((pVariables+offs)*sizeof(int));
1221  int cnt=pVariables+offs;
1222  loop
1223  {
1224    if(! Li.hasItem()) break;
1225    BOOLEAN done=TRUE;
1226    i=Li.getItem()-1;
1227    mark[i]=1;
1228    if (i<offs)
1229    {
1230      done=FALSE;
1231      //StringAppendS(currRing->parameter[i]);
1232    }
1233    else
1234    {
1235      StringAppendS(currRing->names[i-offs]);
1236    }
1237    Li++;
1238    cnt--;
1239    if(cnt==0) break;
1240    if (done) StringAppendS(",");
1241  }
1242  for(i=0;i<pVariables+offs;i++)
1243  {
1244    BOOLEAN done=TRUE;
1245    if(mark[i]==0)
1246    {
1247      if (i<offs)
1248      {
1249        done=FALSE;
1250        //StringAppendS(currRing->parameter[i]);
1251      }
1252      else
1253      {
1254        StringAppendS(currRing->names[i-offs]);
1255      }
1256      cnt--;
1257      if(cnt==0) break;
1258      if (done) StringAppendS(",");
1259    }
1260  }
1261  char * s=omStrDup(StringAppendS(""));
1262  if (s[strlen(s)-1]==',') s[strlen(s)-1]='\0';
1263  return s;
1264#else
1265  return NULL;
1266#endif
1267}
1268
1269BOOLEAN singclap_isSqrFree(poly f)
1270{
1271  BOOLEAN b=FALSE;
1272  Off(SW_RATIONAL);
1273  //  Q / Fp
1274  if (((nGetChar() == 0) || (nGetChar() > 1) )
1275  &&(currRing->parameter==NULL))
1276  {
1277    setCharacteristic( nGetChar() );
1278    CanonicalForm F( convSingPClapP( f ) );
1279    if((nGetChar()>1)&&(!F.isUnivariate()))
1280      goto err;
1281    b=(BOOLEAN)isSqrFree(F);
1282  }
1283  // and over Q(a) / Fp(a)
1284  else if (( nGetChar()==1 ) /* Q(a) */
1285  || (nGetChar() <-1))       /* Fp(a) */
1286  {
1287    if (nGetChar()==1) setCharacteristic( 0 );
1288    else               setCharacteristic( -nGetChar() );
1289    //if (currRing->minpoly!=NULL)
1290    //{
1291    //  CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1292    //  Variable a=rootOf(mipo);
1293    //  CanonicalForm F( convSingAPClapAP( f,a ) );
1294    //  ...
1295    //}
1296    //else
1297    {
1298      CanonicalForm F( convSingTrPClapP( f ) );
1299      b=(BOOLEAN)isSqrFree(F);
1300    }
1301    Off(SW_RATIONAL);
1302  }
1303  else
1304  {
1305err:
1306    WerrorS( feNotImplemented );
1307  }
1308  return b;
1309}
1310
1311poly singclap_det( const matrix m )
1312{
1313  int r=m->rows();
1314  if (r!=m->cols())
1315  {
1316    Werror("det of %d x %d matrix",r,m->cols());
1317    return NULL;
1318  }
1319  poly res=NULL;
1320  if (( nGetChar() == 0 || nGetChar() > 1 )
1321  && (currRing->parameter==NULL))
1322  {
1323    setCharacteristic( nGetChar() );
1324    CFMatrix M(r,r);
1325    int i,j;
1326    for(i=r;i>0;i--)
1327    {
1328      for(j=r;j>0;j--)
1329      {
1330        M(i,j)=convSingPClapP(MATELEM(m,i,j));
1331      }
1332    }
1333    res= convClapPSingP( determinant(M,r) ) ;
1334  }
1335  // and over Q(a) / Fp(a)
1336  else if (( nGetChar()==1 ) /* Q(a) */
1337  || (nGetChar() <-1))       /* Fp(a) */
1338  {
1339    if (nGetChar()==1) setCharacteristic( 0 );
1340    else               setCharacteristic( -nGetChar() );
1341    CFMatrix M(r,r);
1342    poly res;
1343    if (currRing->minpoly!=NULL)
1344    {
1345      CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1346      Variable a=rootOf(mipo);
1347      int i,j;
1348      for(i=r;i>0;i--)
1349      {
1350        for(j=r;j>0;j--)
1351        {
1352          M(i,j)=convSingAPClapAP(MATELEM(m,i,j),a);
1353        }
1354      }
1355      res= convClapAPSingAP( determinant(M,r) ) ;
1356    }
1357    else
1358    {
1359      int i,j;
1360      for(i=r;i>0;i--)
1361      {
1362        for(j=r;j>0;j--)
1363        {
1364          M(i,j)=convSingTrPClapP(MATELEM(m,i,j));
1365        }
1366      }
1367      res= convClapPSingTrP( determinant(M,r) );
1368    }
1369  }
1370  else
1371    WerrorS( feNotImplemented );
1372  Off(SW_RATIONAL);
1373  return res;
1374}
1375
1376int singclap_det_i( intvec * m )
1377{
1378  setCharacteristic( 0 );
1379  CFMatrix M(m->rows(),m->cols());
1380  int i,j;
1381  for(i=m->rows();i>0;i--)
1382  {
1383    for(j=m->cols();j>0;j--)
1384    {
1385      M(i,j)=IMATELEM(*m,i,j);
1386    }
1387  }
1388  int res= convClapISingI( determinant(M,m->rows())) ;
1389  Off(SW_RATIONAL);
1390  return res;
1391}
1392napoly singclap_alglcm ( napoly f, napoly g )
1393{
1394 FACTORY_ALGOUT( "f", f );
1395 FACTORY_ALGOUT( "g", g );
1396
1397 // over Q(a) / Fp(a)
1398 if (nGetChar()==1) setCharacteristic( 0 );
1399 else               setCharacteristic( -nGetChar() );
1400 napoly res;
1401
1402 if (currRing->minpoly!=NULL)
1403 {
1404   CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1405   Variable a=rootOf(mipo);
1406   CanonicalForm F( convSingAClapA( f,a ) ), G( convSingAClapA( g,a ) );
1407   CanonicalForm GCD;
1408
1409   TIMING_START( algLcmTimer );
1410   // calculate gcd
1411#ifdef FACTORY_GCD_TEST
1412   GCD = CFPrimitiveGcdUtil::gcd( F, G );
1413#else
1414   GCD = gcd( F, G );
1415#endif
1416   TIMING_END( algLcmTimer );
1417
1418   FACTORY_CFAOUT( "d", GCD );
1419   FACTORY_GCDSTAT( "alcm:", F, G, GCD );
1420
1421   // calculate lcm
1422   res= convClapASingA( (F/GCD)*G );
1423 }
1424 else
1425 {
1426   CanonicalForm F( convSingTrClapP( f ) ), G( convSingTrClapP( g ) );
1427   CanonicalForm GCD;
1428   TIMING_START( algLcmTimer );
1429   // calculate gcd
1430#ifdef FACTORY_GCD_TEST
1431   GCD = CFPrimitiveGcdUtil::gcd( F, G );
1432#else
1433   GCD = gcd( F, G );
1434#endif
1435   TIMING_END( algLcmTimer );
1436
1437   FACTORY_CFTROUT( "d", GCD );
1438   FACTORY_GCDSTAT( "alcm:", F, G, GCD );
1439
1440   // calculate lcm
1441   res= convClapPSingTr( (F/GCD)*G );
1442 }
1443
1444 Off(SW_RATIONAL);
1445 return res;
1446}
1447
1448void singclap_algdividecontent ( napoly f, napoly g, napoly &ff, napoly &gg )
1449{
1450 FACTORY_ALGOUT( "f", f );
1451 FACTORY_ALGOUT( "g", g );
1452
1453 // over Q(a) / Fp(a)
1454 if (nGetChar()==1) setCharacteristic( 0 );
1455 else               setCharacteristic( -nGetChar() );
1456 ff=gg=NULL;
1457 On(SW_RATIONAL);
1458
1459 if (currRing->minpoly!=NULL)
1460 {
1461   CanonicalForm mipo=convSingTrClapP(((lnumber)currRing->minpoly)->z);
1462   Variable a=rootOf(mipo);
1463   CanonicalForm F( convSingAClapA( f,a ) ), G( convSingAClapA( g,a ) );
1464   CanonicalForm GCD;
1465
1466   TIMING_START( algContentTimer );
1467#ifdef FACTORY_GCD_TEST
1468   GCD=CFPrimitiveGcdUtil::gcd( F, G );
1469#else
1470   GCD=gcd( F, G );
1471#endif
1472   TIMING_END( algContentTimer );
1473
1474   FACTORY_CFAOUT( "d", GCD );
1475   FACTORY_GCDSTAT( "acnt:", F, G, GCD );
1476
1477   if ((GCD!=1) && (GCD!=0))
1478   {
1479     ff= convClapASingA( F/ GCD );
1480     gg= convClapASingA( G/ GCD );
1481   }
1482 }
1483 else
1484 {
1485   CanonicalForm F( convSingTrClapP( f ) ), G( convSingTrClapP( g ) );
1486   CanonicalForm GCD;
1487
1488   TIMING_START( algContentTimer );
1489#ifdef FACTORY_GCD_TEST
1490   GCD=CFPrimitiveGcdUtil::gcd( F, G );
1491#else
1492   GCD=gcd( F, G );
1493#endif
1494   TIMING_END( algContentTimer );
1495
1496   FACTORY_CFTROUT( "d", GCD );
1497   FACTORY_GCDSTAT( "acnt:", F, G, GCD );
1498
1499   if ((GCD!=1) && (GCD!=0))
1500   {
1501     ff= convClapPSingTr( F/ GCD );
1502     gg= convClapPSingTr( G/ GCD );
1503   }
1504 }
1505
1506 Off(SW_RATIONAL);
1507}
1508
1509#if 0
1510lists singclap_chineseRemainder(lists x, lists q)
1511{
1512  //assume(x->nr == q->nr);
1513  //assume(x->nr >= 0);
1514  int n=x->nr+1;
1515  if ((x->nr<0) || (x->nr!=q->nr))
1516  {
1517    WerrorS("list are empty or not of equal length");
1518    return NULL;
1519  }
1520  lists res=(lists)omAlloc0Bin(slists_bin);
1521  CFArray X(1,n), Q(1,n);
1522  int i;
1523  for(i=0; i<n; i++)
1524  {
1525    if (x->m[i-1].Typ()==INT_CMD)
1526    {
1527      X[i]=(int)x->m[i-1].Data();
1528    }
1529    else if (x->m[i-1].Typ()==NUMBER_CMD)
1530    {
1531      number N=(number)x->m[i-1].Data();
1532      X[i]=convSingNClapN(N);
1533    }
1534    else
1535    {
1536      WerrorS("illegal type in chineseRemainder");
1537      omFreeBin(res,slists_bin);
1538      return NULL;
1539    }
1540    if (q->m[i-1].Typ()==INT_CMD)
1541    {
1542      Q[i]=(int)q->m[i-1].Data();
1543    }
1544    else if (q->m[i-1].Typ()==NUMBER_CMD)
1545    {
1546      number N=(number)x->m[i-1].Data();
1547      Q[i]=convSingNClapN(N);
1548    }
1549    else
1550    {
1551      WerrorS("illegal type in chineseRemainder");
1552      omFreeBin(res,slists_bin);
1553      return NULL;
1554    }
1555  }
1556  CanonicalForm r, prod;
1557  chineseRemainder( X, Q, r, prod );
1558  res->Init(2);
1559  res->m[0].rtyp=NUMBER_CMD;
1560  res->m[1].rtyp=NUMBER_CMD;
1561  res->m[0].data=(char *)convClapNSingN( r );
1562  res->m[1].data=(char *)convClapNSingN( prod );
1563  return res;
1564}
1565#endif
1566#endif
Note: See TracBrowser for help on using the repository browser.