source: git/libpolys/polys/clapconv.cc @ fbdfd4

spielwiese
Last change on this file since fbdfd4 was fbdfd4, checked in by Hans Schoenemann <hannes@…>, 11 years ago
fix: enforce integer coeffs in rational functions
  • Property mode set to 100644
File size: 11.0 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2/****************************************
3*  Computer Algebra System SINGULAR     *
4****************************************/
5/*
6* ABSTRACT: convert data between Singular and factory
7*/
8
9
10#ifdef HAVE_CONFIG_H
11#include "libpolysconfig.h"
12#endif /* HAVE_CONFIG_H */
13#include <misc/auxiliary.h>
14
15#ifdef HAVE_FACTORY
16#define SI_DONT_HAVE_GLOBAL_VARS
17#include <factory/factory.h>
18
19#include <omalloc/omalloc.h>
20#include <coeffs/coeffs.h>
21#include <coeffs/longrat.h>
22#include <coeffs/modulop.h>
23#include <polys/monomials/p_polys.h>
24#include <polys/sbuckets.h>
25#include <polys/clapconv.h>
26
27#include "simpleideals.h"
28
29#define TRANSEXT_PRIVATES
30#include <polys/ext_fields/transext.h>
31
32void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
33
34static void conv_RecPP ( const CanonicalForm & f, int * exp, sBucket_pt result, ring r );
35
36static void convRecTrP ( const CanonicalForm & f, int * exp, poly & result, int offs, const ring r );
37
38static void convRecGFGF ( const CanonicalForm & f, int * exp, poly & result );
39
40static number convFactoryNSingAN( const CanonicalForm &f, const ring r);
41
42poly convFactoryPSingP ( const CanonicalForm & f, const ring r )
43{
44  int n = rVar(r)+1;
45  /* ASSERT( level( f ) <= pVariables, "illegal number of variables" ); */
46  int * exp = (int*)omAlloc0(n*sizeof(int));
47  sBucket_pt result_bucket=sBucketCreate(r);
48  conv_RecPP( f, exp, result_bucket, r );
49  poly result; int dummy;
50  sBucketDestroyMerge(result_bucket,&result,&dummy);
51  omFreeSize((ADDRESS)exp,n*sizeof(int));
52  return result;
53}
54
55static void conv_RecPP ( const CanonicalForm & f, int * exp, sBucket_pt result, ring r )
56{
57  if (f.isZero())
58    return;
59  if ( ! f.inCoeffDomain() )
60  {
61    int l = f.level();
62    for ( CFIterator i = f; i.hasTerms(); i++ )
63    {
64      exp[l] = i.exp();
65      conv_RecPP( i.coeff(), exp, result, r );
66    }
67    exp[l] = 0;
68  }
69  else
70  {
71    poly term = p_Init(r);
72    pNext( term ) = NULL;
73    for ( int i = 1; i <= r->N; i++ )
74      p_SetExp( term, i, exp[i], r);
75    pGetCoeff( term )=r->cf->convFactoryNSingN(f, r->cf);
76    p_Setm( term, r );
77    if ( n_IsZero(pGetCoeff(term), r->cf) )
78    {
79      p_Delete(&term,r);
80    }
81    else
82    {
83      sBucket_Merge_p(result,term,1);
84    }
85  }
86}
87
88
89CanonicalForm convSingPFactoryP( poly p, const ring r )
90{
91  CanonicalForm result = 0;
92  int e, n = rVar(r);
93  BOOLEAN setChar=TRUE;
94
95  while ( p!=NULL )
96  {
97    CanonicalForm term;
98    term=r->cf->convSingNFactoryN(pGetCoeff( p ),setChar, r->cf);
99    if (errorreported) break;
100    setChar=FALSE;
101    for ( int i = n; i >0; i-- )
102    {
103      if ( (e = p_GetExp( p, i, r)) != 0 )
104        term *= power( Variable( i ), e );
105    }
106    result += term;
107    pIter( p );
108  }
109 return result;
110}
111
112int convFactoryISingI( const CanonicalForm & f)
113{
114  if (!f.isImm()) WerrorS("int overflow in det");
115  return f.intval();
116}
117
118CanonicalForm convSingAPFactoryAP ( poly p , const Variable & a, const ring r)
119{
120  CanonicalForm result = 0;
121  int e, n = r-> N;
122  int off=rPar(r);
123
124  if (!rField_is_Zp_a(r))
125    On(SW_RATIONAL);
126  while ( p!=NULL)
127  {
128    CanonicalForm term=convSingAFactoryA(((poly)p_GetCoeff(p, r->cf->extRing)),a, r);
129    for ( int i = 1; i <= n; i++ )
130    {
131      if ( (e = p_GetExp( p, i, r )) != 0 )
132        term *= power( Variable( i + off), e );
133    }
134    result += term;
135    pIter( p );
136  }
137  return result;
138}
139
140static void
141convRecAP_R ( const CanonicalForm & f, int * exp, poly & result, int par_start, int var_start, const ring r) ;
142
143poly convFactoryAPSingAP_R ( const CanonicalForm & f, int par_start, int var_start, const ring r )
144{
145  int n = rVar(r)+rPar(r)+1;
146  int * exp = (int *)omAlloc0(n*sizeof(int));
147  poly result = NULL;
148  convRecAP_R( f, exp, result,par_start, var_start, r );
149  omFreeSize((ADDRESS)exp,n*sizeof(int));
150  return result;
151}
152
153poly convFactoryAPSingAP ( const CanonicalForm & f, const ring r )
154{
155  return convFactoryAPSingAP_R(f,0,rPar(r),r);
156}
157
158static void convRecAP_R ( const CanonicalForm & f, int * exp, poly & result, int par_start, int var_start, const ring r )
159{
160  if (f.isZero())
161    return;
162  if ( ! f.inCoeffDomain() )
163  {
164    int l = f.level();
165    for ( CFIterator i = f; i.hasTerms(); i++ )
166    {
167      exp[l] = i.exp();
168      convRecAP_R( i.coeff(), exp, result, par_start, var_start, r);
169    }
170    exp[l] = 0;
171  }
172  else
173  {
174    poly z=(poly)convFactoryASingA( f,r );
175    if (z!=NULL)
176    {
177      poly term = p_Init(r);
178      pNext( term ) = NULL;
179      int i;
180      for ( i = rVar(r); i>0 ; i-- )
181        p_SetExp( term, i , exp[i+var_start],r);
182      //if (rRing_has_Comp(currRing->extRing)) p_SetComp(term, 0, currRing->extRing); // done by pInit
183      if (par_start==0)
184      {
185        for ( i = 1; i <= var_start; i++ )
186        //z->e[i-1]+=exp[i];
187          p_AddExp(z,i,exp[i],r->cf->extRing);
188      }
189      else
190      {
191        for ( i = par_start+1; i <= var_start+rPar(r); i++ )
192        //z->e[i-1]+=exp[i];
193          p_AddExp(z,i,exp[i-par_start],r->cf->extRing);
194      }
195      pGetCoeff(term)=(number)ALLOC0_RNUMBER();
196      p_GetCoeff(term, r->cf->extRing)=(number) z;
197      p_Setm( term,r );
198      result = p_Add_q( result, term, r );
199    }
200  }
201}
202
203CanonicalForm convSingAFactoryA ( poly p , const Variable & a, const ring r )
204{
205  CanonicalForm result = 0;
206  int e;
207
208  while ( p!=NULL )
209  {
210    CanonicalForm term;
211    if ( rField_is_Zp_a(r) )
212    {
213      term = n_Int( p_GetCoeff( p, r->cf->extRing ), r->cf->extRing->cf );
214    }
215    else
216    {
217      if ( SR_HDL(p_GetCoeff( p, r->cf->extRing )) & SR_INT )
218        term = SR_TO_INT(p_GetCoeff( p, r->cf->extRing )) ;
219      else
220      {
221        if ( p_GetCoeff( p, r->cf->extRing )->s == 3 )
222        {
223          mpz_t dummy;
224          mpz_init_set( dummy, (p_GetCoeff( p,r->cf->extRing )->z) );
225          term = make_cf( dummy );
226        }
227        else
228        {
229          // assume s==0 or s==1
230          mpz_t num, den;
231          On(SW_RATIONAL);
232          mpz_init_set( num, (p_GetCoeff( p, r->cf->extRing )->z) );
233          mpz_init_set( den, (p_GetCoeff( p, r->cf->extRing )->n) );
234          term = make_cf( num, den, ( p_GetCoeff( p, r->cf->extRing )->s != 1 ));
235        }
236      }
237    }
238    if ( (e = p_GetExp( p, 1, r->cf->extRing )) != 0 )
239      term *= power( a , e );
240    result += term;
241    p = pNext( p );
242  }
243  return result;
244}
245
246static number convFactoryNSingAN( const CanonicalForm &f, const ring r)
247{
248  if ( f.isImm() )
249  {
250    long longf=f.intval();
251    int intf=(int) longf;
252    if((long)intf==longf)
253    {
254      assume (r->cf->extRing != NULL);
255      return n_Init(f.intval(),r->cf->extRing->cf);
256    }
257    else return nlRInit( longf );
258  }
259  else
260  {
261    number z=ALLOC_RNUMBER();
262#if defined(LDEBUG)
263    z->debug=123456;
264#endif
265    gmp_numerator( f, z->z );
266    if ( f.den().isOne() )
267    {
268      z->s = 3;
269    }
270    else
271    {
272      gmp_denominator( f, z->n );
273      z->s = 0;
274      nlNormalize(z,r->cf->extRing->cf);
275    }
276    /*#ifdef LDEBUG
277    nlTest(z,r->cf->extRing->cf);
278    #endif*/
279    return z;
280  }
281}
282
283poly convFactoryASingA ( const CanonicalForm & f, const ring r )
284{
285  poly a=NULL;
286  poly t;
287  for( CFIterator i=f; i.hasTerms(); i++)
288  {
289    t= p_Init (r->cf->extRing);
290    p_GetCoeff(t, r->cf->extRing)= convFactoryNSingAN( i.coeff(), r );
291    if (n_IsZero(p_GetCoeff(t,r->cf->extRing),r->cf->extRing->cf))
292    {
293      p_Delete(&t,r->cf->extRing);
294    }
295    else
296    {
297      p_SetExp(t,1,i.exp(),r->cf->extRing);
298      a=p_Add_q(a,t,r->cf->extRing);
299    }
300  }
301  if (a!=NULL)
302  {
303    if( r->cf->extRing != NULL )
304      if (r->cf->extRing->qideal->m[0]!=NULL)
305      {
306        poly l=r->cf->extRing->qideal->m[0];
307        if (p_GetExp(a,1,r->cf->extRing) >= p_GetExp(l,1,r->cf->extRing))
308          a = p_PolyDiv (a, l, FALSE, r->cf->extRing); // ???
309      }
310  }
311  return a;
312}
313
314CanonicalForm convSingTrPFactoryP ( poly p, const ring r )
315{
316  CanonicalForm result = 0;
317  int e, n = rVar(r);
318  int offs = rPar(r);
319
320  while ( p!=NULL )
321  {
322    n_Normalize(p_GetCoeff(p, r), r->cf);
323
324    // test if denominator is constant
325    if (!p_IsConstantPoly(DEN (p_GetCoeff (p,r)),r->cf->extRing) && !errorreported)
326      WerrorS("conversion error: denominator!= 1");
327
328    CanonicalForm term=convSingPFactoryP(NUM (p_GetCoeff(p, r)),r->cf->extRing);
329
330    // if denominator is not NULL it should be a constant at this point
331    if (DEN (p_GetCoeff(p,r)) != NULL)
332    {
333      CanonicalForm den= convSingPFactoryP(DEN (p_GetCoeff(p, r)),r->cf->extRing);
334      if (rChar (r) == 0)
335        On (SW_RATIONAL);
336      term /= den;
337    }
338
339    for ( int i = n; i > 0; i-- )
340    {
341      if ( (e = p_GetExp( p, i,r )) != 0 )
342        term = term * power( Variable( i + offs ), e );
343    }
344    result += term;
345    p = pNext( p );
346  }
347  return result;
348}
349
350poly convFactoryPSingTrP ( const CanonicalForm & f, const ring r )
351{
352  int n = rVar(r)+1;
353  int * exp = (int*)omAlloc0(n*sizeof(int));
354  poly result = NULL;
355  convRecTrP( f, exp, result , rPar(r), r );
356  omFreeSize((ADDRESS)exp,n*sizeof(int));
357  return result;
358}
359
360static void
361convRecTrP ( const CanonicalForm & f, int * exp, poly & result , int offs, const ring r)
362{
363  if (f.isZero())
364    return;
365  if ( f.level() > offs )
366  {
367    int l = f.level();
368    for ( CFIterator i = f; i.hasTerms(); i++ )
369    {
370      exp[l-offs] = i.exp();
371      convRecTrP( i.coeff(), exp, result, offs, r );
372    }
373    exp[l-offs] = 0;
374  }
375  else
376  {
377    poly term = p_Init(r);
378    pNext( term ) = NULL;
379    for ( int i = rVar(r); i>0; i-- )
380      p_SetExp( term, i ,exp[i], r);
381    //if (rRing_has_Comp(currRing)) p_SetComp(term, 0, currRing); // done by pInit
382    pGetCoeff(term)=ntInit(convFactoryPSingP( f, r->cf->extRing ), r->cf);
383    p_Setm( term,r );
384    result = p_Add_q( result, term,r );
385  }
386}
387
388#if 0
389CanonicalForm
390convSingGFFactoryGF( poly p )
391{
392  CanonicalForm result=CanonicalForm(0);
393  int e, n = pVariables;
394
395  while ( p != NULL )
396  {
397    CanonicalForm term;
398    term = make_cf_from_gf( (int)(long)pGetCoeff( p ) );
399    //int * A=(int *)&term;
400    //Print("term=%x, == 0 ?: %d\n",*A, term.isZero());
401    for ( int i = 1; i <= n; i++ )
402    {
403      if ( (e = pGetExp( p, i )) != 0 )
404         term *= power( Variable( i ), e );
405    }
406    result += term;
407    p = pNext( p );
408  }
409  return result;
410}
411
412poly
413convFactoryGFSingGF ( const CanonicalForm & f )
414{
415//    cerr << " f = " << f << endl;
416  int n = pVariables+1;
417  /* ASSERT( level( f ) <= pVariables, "illegal number of variables" ); */
418  int * exp = (int*)omAlloc0(n*sizeof(int));
419  poly result = NULL;
420  convRecGFGF( f, exp, result );
421  omFreeSize((ADDRESS)exp,n*sizeof(int));
422  return result;
423}
424
425static void
426convRecGFGF ( const CanonicalForm & f, int * exp, poly & result )
427{
428  if (f.isZero())
429    return;
430  if ( ! f.inCoeffDomain() )
431  {
432    int l = f.level();
433    for ( CFIterator i = f; i.hasTerms(); i++ )
434    {
435      exp[l] = i.exp();
436      convRecGFGF( i.coeff(), exp, result );
437    }
438    exp[l] = 0;
439  }
440  else
441  {
442    poly term = pInit();
443    pNext( term ) = NULL;
444    for ( int i = 1; i <= pVariables; i++ )
445      pSetExp( term, i, exp[i]);
446    //if (rRing_has_Comp(currRing)) p_SetComp(term, 0, currRing); // done by pInit
447    pGetCoeff( term ) = (number) gf_value (f);
448    pSetm( term );
449    result = pAdd( result, term );
450  }
451}
452
453#endif
454#endif /* HAVE_FACTORY */
Note: See TracBrowser for help on using the repository browser.