source: git/kernel/clapsing.cc @ ffdd694

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