source: git/kernel/clapconv.cc @ 599326

spielwiese
Last change on this file since 599326 was 599326, checked in by Kai Krüger <krueger@…>, 14 years ago
Anne, Kai, Frank: - changes to #include "..." statements to allow cleaner build structure - affected directories: omalloc, kernel, Singular - not yet done: IntergerProgramming git-svn-id: file:///usr/local/Singular/svn/trunk@13032 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 13.2 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2/****************************************
3*  Computer Algebra System SINGULAR     *
4****************************************/
5// $Id$
6/*
7* ABSTRACT: convert data between Singular and factory
8*/
9
10
11#include <kernel/mod2.h>
12#ifdef HAVE_FACTORY
13#include <omalloc.h>
14#include <kernel/structs.h>
15#define SI_DONT_HAVE_GLOBAL_VARS
16#include <kernel/clapconv.h>
17#include <kernel/numbers.h>
18#include <kernel/modulop.h>
19#include <kernel/longalg.h>
20#include <kernel/polys.h>
21#include <kernel/febase.h>
22#include <kernel/ring.h>
23#include <kernel/sbuckets.h>
24#include <kernel/ffields.h>
25
26typedef __mpz_struct lint;
27
28void out_cf(char *s1,const CanonicalForm &f,char *s2);
29
30static void convRec( const CanonicalForm & f, int * exp, poly & result );
31
32static void convRecAlg( const CanonicalForm & f, int * exp, napoly & result );
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
42CanonicalForm convSingNFactoryN( number n, const ring r )
43{
44  CanonicalForm term;
45  /* does only work for Zp, Q */
46  if ( getCharacteristic() != 0 )
47  {
48    term = npInt( n,r );
49  }
50  else
51  {
52    if ( SR_HDL(n) & SR_INT )
53    {
54      term = SR_TO_INT(n);
55    }
56    else
57    {
58      if ( n->s == 3 )
59      {
60        lint dummy;
61        mpz_init_set( &dummy,n->z );
62        term = make_cf( dummy );
63      }
64      else
65      {
66        // assume s==0 or s==1
67        lint num, den;
68        On(SW_RATIONAL);
69        mpz_init_set( &num, n->z );
70        mpz_init_set( &den, n->n );
71        term = make_cf( num, den, ( n->s != 1 ));
72      }
73    }
74  }
75  return term;
76}
77
78number convFactoryNSingN( const CanonicalForm & n)
79{
80  if (n.isImm())
81    return nInit(n.intval());
82  else
83  {
84    number z=(number)omAllocBin(rnumber_bin);
85#if defined(LDEBUG)
86    z->debug=123456;
87#endif
88    gmp_numerator( n, z->z );
89    if ( n.den().isOne() )
90      z->s = 3;
91    else
92    {
93      gmp_denominator( n, z->n );
94      z->s = 0;
95    }
96    return z;
97  }
98}
99
100poly convFactoryPSingP ( const CanonicalForm & f, const ring r )
101{
102//    cerr << " f = " << f << endl;
103  int n = rVar(r)+1;
104  /* ASSERT( level( f ) <= pVariables, "illegal number of variables" ); */
105  int * exp = (int*)omAlloc0(n*sizeof(int));
106  sBucket_pt result_bucket=sBucketCreate(r);
107  conv_RecPP( f, exp, result_bucket, r );
108  poly result; int dummy;
109  sBucketDestroyMerge(result_bucket,&result,&dummy);
110  omFreeSize((ADDRESS)exp,n*sizeof(int));
111  return result;
112}
113
114static void conv_RecPP ( const CanonicalForm & f, int * exp, sBucket_pt result, ring r )
115{
116  if (f.isZero())
117    return;
118  if ( ! f.inCoeffDomain() )
119  {
120    int l = f.level();
121    for ( CFIterator i = f; i.hasTerms(); i++ )
122    {
123      exp[l] = i.exp();
124      conv_RecPP( i.coeff(), exp, result, r );
125    }
126    exp[l] = 0;
127  }
128  else
129  {
130    poly term = p_Init(r);
131    pNext( term ) = NULL;
132    for ( int i = 1; i <= r->N; i++ )
133      p_SetExp( term, i, exp[i], r);
134    //if (rRing_has_Comp(r)) p_SetComp(term, 0, r); // done by p_Init
135    if ( f.isImm() )
136      pGetCoeff( term ) = n_Init( f.intval(), r );
137    else
138    {
139      number z=(number)omAllocBin(rnumber_bin);
140#if defined(LDEBUG)
141      z->debug=123456;
142#endif
143      gmp_numerator( f, z->z );
144      if ( f.den().isOne() )
145        z->s = 3;
146      else
147      {
148        gmp_denominator( f, z->n );
149        z->s = 0;
150        nlNormalize(z);
151      }
152      pGetCoeff( term ) = z;
153    }
154    p_Setm( term, r );
155    if ( n_IsZero(pGetCoeff(term), r) )
156    {
157      p_Delete(&term,r);
158    }
159    else
160    {
161      sBucket_Merge_p(result,term,1);
162    }
163  }
164}
165
166
167CanonicalForm convSingPFactoryP( poly p, const ring r )
168{
169  CanonicalForm result = 0;
170  int e, n = rVar(r);
171
172  while ( p!=NULL)
173  {
174    CanonicalForm term;
175    /* does only work for finite fields */
176    if ( rField_is_Zp(r) )
177    {
178      term = npInt( pGetCoeff( p ),r );
179    }
180    else if (rField_is_Q(r))
181    {
182      if ( SR_HDL(pGetCoeff( p )) & SR_INT )
183        term = SR_TO_INT( pGetCoeff( p ) );
184      else
185      {
186        if ( pGetCoeff( p )->s == 3 )
187        {
188          lint dummy;
189          mpz_init_set( &dummy, (pGetCoeff( p )->z) );
190          term = make_cf( dummy );
191        }
192        else
193        {
194          // assume s==0 or s==1
195          lint num, den;
196          On(SW_RATIONAL);
197          mpz_init_set( &num, (pGetCoeff( p )->z) );
198          mpz_init_set( &den, (pGetCoeff( p )->n) );
199          term = make_cf( num, den, ( pGetCoeff( p )->s != 1 ));
200        }
201      }
202    }
203    else
204    {
205      WerrorS("conversion error");
206      return result;
207    }
208    for ( int i = n; i >0; i-- )
209    {
210      if ( (e = p_GetExp( p, i, r)) != 0 )
211        term *= power( Variable( i ), e );
212    }
213    result += term;
214    pIter( p );
215  }
216  return result;
217}
218
219CanonicalForm convSingAPFactoryAP ( poly p , const Variable & a, const ring r)
220{
221  CanonicalForm result = 0;
222  int e, n = pVariables;
223  int off=rPar(currRing);
224
225  On(SW_RATIONAL);
226  while ( p!=NULL)
227  {
228    CanonicalForm term=convSingAFactoryA(((lnumber)pGetCoeff(p))->z,a, r);
229    for ( int i = 1; i <= n; i++ )
230    {
231      if ( (e = pGetExp( p, i )) != 0 )
232        term *= power( Variable( i + off), e );
233    }
234    result += term;
235    pIter( p );
236  }
237  return result;
238}
239
240static void
241convRecAP_R ( const CanonicalForm & f, int * exp, poly & result, int par_start, int var_start, const ring r) ;
242
243poly convFactoryAPSingAP_R ( const CanonicalForm & f, int par_start, int var_start, const ring r )
244{
245  int n = rVar(r)+rPar(r)+1;
246  int * exp = (int *)omAlloc0(n*sizeof(int));
247  poly result = NULL;
248  convRecAP_R( f, exp, result,par_start, var_start, r );
249  omFreeSize((ADDRESS)exp,n*sizeof(int));
250  return result;
251}
252
253poly convFactoryAPSingAP ( const CanonicalForm & f, const ring r )
254{
255  return convFactoryAPSingAP_R(f,0,rPar(r),r);
256}
257
258static void convRecAP_R ( const CanonicalForm & f, int * exp, poly & result, int par_start, int var_start, const ring r )
259{
260  if (f.isZero())
261    return;
262  if ( ! f.inCoeffDomain() )
263  {
264    int l = f.level();
265    for ( CFIterator i = f; i.hasTerms(); i++ )
266    {
267      exp[l] = i.exp();
268      convRecAP_R( i.coeff(), exp, result, par_start, var_start, r);
269    }
270    exp[l] = 0;
271  }
272  else
273  {
274    napoly z=(napoly)convFactoryASingA( f,r );
275    if (z!=NULL)
276    {
277      poly term = p_Init(r);
278      pNext( term ) = NULL;
279      int i;
280      for ( i = rVar(r); i>0 ; i-- )
281        p_SetExp( term, i , exp[i+var_start],r);
282      //if (rRing_has_Comp(currRing->algring)) p_SetComp(term, 0, currRing->algring); // done by pInit
283      if (par_start==0)
284      {
285        for ( i = 1; i <= var_start; i++ )
286        //z->e[i-1]+=exp[i];
287          p_AddExp(z,i,exp[i],r->algring);
288      }
289      else
290      {
291        for ( i = par_start+1; i <= var_start+rPar(currRing); i++ )
292        //z->e[i-1]+=exp[i];
293          p_AddExp(z,i,exp[i-par_start],r->algring);
294      }
295      pGetCoeff(term)=(number)omAlloc0Bin(rnumber_bin);
296      ((lnumber)pGetCoeff(term))->z=z;
297      p_Setm( term,r );
298      result = p_Add_q( result, term, r );
299    }
300  }
301}
302
303CanonicalForm convSingAFactoryA ( napoly p , const Variable & a, const ring r )
304{
305  CanonicalForm result = 0;
306  int e;
307
308  while ( p!=NULL )
309  {
310    CanonicalForm term;
311    /* does only work for finite fields:Z/p */
312    if ( rField_is_Zp_a() )
313    {
314      term = npInt( napGetCoeff( p ),r );
315    }
316    else
317    {
318      if ( SR_HDL(napGetCoeff( p )) & SR_INT )
319        term = SR_TO_INT(napGetCoeff( p )) ;
320      else
321      {
322        if ( napGetCoeff( p )->s == 3 )
323        {
324          lint dummy;
325          mpz_init_set( &dummy, (napGetCoeff( p )->z) );
326          term = make_cf( dummy );
327        }
328        else
329        {
330          // assume s==0 or s==1
331          lint num, den;
332          On(SW_RATIONAL);
333          mpz_init_set( &num, (napGetCoeff( p )->z) );
334          mpz_init_set( &den, (napGetCoeff( p )->n) );
335          term = make_cf( num, den, ( napGetCoeff( p )->s != 1 ));
336        }
337      }
338    }
339    if ( (e = p_GetExp( p, 1, r->algring )) != 0 )
340      term *= power( a , e );
341    result += term;
342    p = pNext( p );
343  }
344  return result;
345}
346
347static number convFactoryNSingAN( const CanonicalForm &f, const ring r)
348{
349  if ( f.isImm() )
350    return n_Init( f.intval(), r->algring );
351  else
352  {
353    number z=(number)omAllocBin(rnumber_bin);
354#if defined(LDEBUG)
355    z->debug=123456;
356#endif
357    gmp_numerator( f, z->z );
358    if ( f.den().isOne() )
359    {
360      z->s = 3;
361    }
362    else
363    {
364      gmp_denominator( f, z->n );
365      z->s = 0;
366      nlNormalize(z);
367    }
368    #ifdef LDEBUG
369    nlDBTest(z,__FILE__,__LINE__);
370    #endif
371    return z;
372  }
373}
374
375napoly convFactoryASingA ( const CanonicalForm & f, const ring r )
376{
377  poly a=NULL;
378  poly t;
379  for( CFIterator i=f; i.hasTerms(); i++)
380  {
381    t=napNew();
382    // pNext( t ) = NULL; //already done by napNew
383    pGetCoeff(t)=convFactoryNSingAN( i.coeff(), r );
384    if (n_IsZero(napGetCoeff(t),r->algring))
385    {
386      p_Delete(&t,r->algring);
387    }
388    else
389    {
390      p_SetExp(t,1,i.exp(),r->algring);
391      a=p_Add_q(a,t,r->algring);
392    }
393  }
394  if (a!=NULL)
395  {
396    if (r->minpoly!=NULL)
397    {
398      lnumber l=(lnumber)r->minpoly;
399      if (p_GetExp(a,1,r->algring) >= p_GetExp(l->z,1,r->algring))
400        a = napRemainder( a, l->z);
401    }
402  }
403  return a;
404}
405
406CanonicalForm convSingTrPFactoryP ( poly p, const ring r )
407{
408  CanonicalForm result = 0;
409  int e, n = rVar(r);
410  int offs = rPar(r);
411
412  while ( p!=NULL )
413  {
414    n_Normalize(pGetCoeff(p),r);
415    CanonicalForm term=convSingPFactoryP(((lnumber)pGetCoeff(p))->z,r->algring);
416
417    if ((((lnumber)pGetCoeff(p))->n!=NULL)
418    && (!errorreported))
419    {
420      WerrorS("conversion error: denominator!= 1");
421    }
422
423    for ( int i = n; i > 0; i-- )
424    {
425      if ( (e = p_GetExp( p, i,r )) != 0 )
426        term = term * power( Variable( i + offs ), e );
427    }
428    result += term;
429    p = pNext( p );
430  }
431  return result;
432}
433
434poly convFactoryPSingTrP ( const CanonicalForm & f, const ring r )
435{
436  int n = rVar(r)+1;
437  int * exp = (int*)omAlloc0(n*sizeof(int));
438  poly result = NULL;
439  convRecTrP( f, exp, result , rPar(r), r );
440  omFreeSize((ADDRESS)exp,n*sizeof(int));
441  return result;
442}
443
444static void
445convRecTrP ( const CanonicalForm & f, int * exp, poly & result , int offs, const ring r)
446{
447  if (f.isZero())
448    return;
449  if ( f.level() > offs )
450  {
451    int l = f.level();
452    for ( CFIterator i = f; i.hasTerms(); i++ )
453    {
454      exp[l-offs] = i.exp();
455      convRecTrP( i.coeff(), exp, result, offs, r );
456    }
457    exp[l-offs] = 0;
458  }
459  else
460  {
461    poly term = p_Init(r);
462    pNext( term ) = NULL;
463    for ( int i = rVar(r); i>0; i-- )
464      p_SetExp( term, i ,exp[i], r);
465    //if (rRing_has_Comp(currRing)) p_SetComp(term, 0, currRing); // done by pInit
466    pGetCoeff(term)=(number)omAlloc0Bin(rnumber_bin);
467    ((lnumber)pGetCoeff(term))->z=convFactoryPSingP( f, r->algring );
468    p_Setm( term,r );
469    result = p_Add_q( result, term,r );
470  }
471}
472
473number   nlChineseRemainder(number *x, number *q,int rl)
474// elemenst in the array are x[0..(rl-1)], q[0..(rl-1)]
475{
476#ifdef HAVE_FACTORY
477  setCharacteristic( 0 ); // only in char 0
478  CFArray X(rl), Q(rl);
479  int i;
480  for(i=rl-1;i>=0;i--)
481  {
482    X[i]=convSingNFactoryN(x[i],currRing); // may be larger MAX_INT
483    Q[i]=convSingNFactoryN(q[i],currRing); // may be larger MAX_INT
484  }
485  CanonicalForm xnew,qnew;
486  chineseRemainder(X,Q,xnew,qnew);
487  number n=convFactoryNSingN(xnew);
488  number p=convFactoryNSingN(qnew);
489  number p2=nlIntDiv(p,nlInit(2, currRing));
490  if (nlGreater(n,p2))
491  {
492     number n2=nlSub(n,p);
493     nlDelete(&n,currRing);
494     n=n2;
495  }
496  nlDelete(&p,currRing);
497  nlDelete(&p2,currRing);
498  return n;
499#else
500  WerrorS("not implemented");
501  return nlInit(0);
502#endif
503}
504
505CanonicalForm
506convSingGFFactoryGF( poly p )
507{
508  CanonicalForm result=CanonicalForm(0);
509  int e, n = pVariables;
510
511  while ( p != NULL )
512  {
513    CanonicalForm term;
514    term = make_cf_from_gf( (int)(long)pGetCoeff( p ) );
515    //int * A=(int *)&term;
516    //Print("term=%x, == 0 ?: %d\n",*A, term.isZero());
517    for ( int i = 1; i <= n; i++ )
518    {
519      if ( (e = pGetExp( p, i )) != 0 )
520         term *= power( Variable( i ), e );
521    }
522    result += term;
523    p = pNext( p );
524  }
525  return result;
526}
527
528poly
529convFactoryGFSingGF ( const CanonicalForm & f )
530{
531//    cerr << " f = " << f << endl;
532  int n = pVariables+1;
533  /* ASSERT( level( f ) <= pVariables, "illegal number of variables" ); */
534  int * exp = (int*)omAlloc0(n*sizeof(int));
535  poly result = NULL;
536  convRecGFGF( f, exp, result );
537  omFreeSize((ADDRESS)exp,n*sizeof(int));
538  return result;
539}
540
541static void
542convRecGFGF ( const CanonicalForm & f, int * exp, poly & result )
543{
544  if (f.isZero())
545    return;
546  if ( ! f.inCoeffDomain() )
547  {
548    int l = f.level();
549    for ( CFIterator i = f; i.hasTerms(); i++ )
550    {
551      exp[l] = i.exp();
552      convRecGFGF( i.coeff(), exp, result );
553    }
554    exp[l] = 0;
555  }
556  else
557  {
558    poly term = pInit();
559    pNext( term ) = NULL;
560    for ( int i = 1; i <= pVariables; i++ )
561      pSetExp( term, i, exp[i]);
562    //if (rRing_has_Comp(currRing)) p_SetComp(term, 0, currRing); // done by pInit
563    pGetCoeff( term ) = (number) gf_value (f);
564    pSetm( term );
565    result = pAdd( result, term );
566  }
567}
568
569int convFactoryISingI( const CanonicalForm & f)
570{
571  if (!f.isImm()) WerrorS("int overflow in det");
572  return f.intval();
573}
574#endif
Note: See TracBrowser for help on using the repository browser.