source: git/kernel/clapconv.cc @ 07625cb

spielwiese
Last change on this file since 07625cb was bba835, checked in by Hans Schönemann <hannes@…>, 17 years ago
*hannes: chinese remainder git-svn-id: file:///usr/local/Singular/svn/trunk@10011 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 18.0 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2/****************************************
3*  Computer Algebra System SINGULAR     *
4****************************************/
5// $Id: clapconv.cc,v 1.9 2007-05-03 13:27:45 Singular Exp $
6/*
7* ABSTRACT: convert data between Singular and factory
8*/
9
10
11#include "mod2.h"
12#ifdef HAVE_FACTORY
13#include "omalloc.h"
14#include "structs.h"
15#define SI_DONT_HAVE_GLOBAL_VARS
16#include "clapconv.h"
17#include "numbers.h"
18#include "modulop.h"
19#include "longalg.h"
20#include "polys.h"
21#include "febase.h"
22#include "ring.h"
23
24static void convRec( const CanonicalForm & f, int * exp, poly & result );
25
26static void convRecAlg( const CanonicalForm & f, int * exp, napoly & result );
27
28static void convRecPP ( const CanonicalForm & f, int * exp, poly & result );
29static void conv_RecPP ( const CanonicalForm & f, int * exp, poly & result, ring r );
30
31static void convRecPTr ( const CanonicalForm & f, int * exp, napoly & result );
32
33static void convRecTrP ( const CanonicalForm & f, int * exp, poly & result, int offs );
34
35static void convRecGFGF ( const CanonicalForm & f, int * exp, poly & result );
36
37static number convClapNSingAN( const CanonicalForm &f);
38
39CanonicalForm
40convSingNClapN( number n )
41{
42  CanonicalForm term;
43  /* does only work for Zp, Q */
44  if ( getCharacteristic() != 0 )
45  {
46    term = npInt( n );
47  }
48  else
49  {
50    if ( ((long)n) & 1L )
51    {
52      term = ((long)n) >>2;
53    }
54    else
55    {
56      if ( n->s == 3 )
57      {
58        MP_INT dummy;
59        mpz_init_set( &dummy, &(n->z) );
60        term = make_cf( dummy );
61      }
62      else  if ( n->s == 1 )
63      {
64        MP_INT num, den;
65        On(SW_RATIONAL);
66        mpz_init_set( &num, &(n->z) );
67        mpz_init_set( &den, &(n->n) );
68        term = make_cf( num, den, false );
69      }
70      else
71      { // assume s == 0
72        MP_INT num, den;
73        mpz_init_set( &num, &(n->z) );
74        mpz_init_set( &den, &(n->n) );
75        term = make_cf( num, den, true );
76      }
77    }
78  }
79  return term;
80}
81
82number
83convClapNSingN( const CanonicalForm & n)
84{
85  if (n.isImm())
86    return nInit(n.intval());
87  else
88  {
89    number z=(number)omAllocBin(rnumber_bin);
90#if defined(LDEBUG)
91    z->debug=123456;
92#endif
93    z->z = gmp_numerator( n );
94    if ( n.den() == 1 )
95      z->s = 3;
96    else
97    {
98      z->n = gmp_denominator( n );
99      z->s = 0;
100    }
101    return z;
102  }
103}
104
105CanonicalForm conv_SingPClapP( poly p, ring r )
106{
107  CanonicalForm result = 0;
108  int e, n = r->N;
109  assume( rPar(r)==0);
110
111  while ( p != NULL )
112  {
113    CanonicalForm term;
114    /* does only work for finite fields */
115    if ( getCharacteristic() != 0 )
116    {
117      term = npInt( pGetCoeff( p ) );
118    }
119    else
120    {
121      if ( (long)(pGetCoeff( p )) & 1 )
122      {
123        term = ((long)( pGetCoeff( p ) )>>2);
124      }
125      else
126      {
127        if ( pGetCoeff( p )->s == 3 )
128        {
129          MP_INT dummy;
130          mpz_init_set( &dummy, &(pGetCoeff( p )->z) );
131          term = make_cf( dummy );
132        }
133        else  if ( pGetCoeff( p )->s == 1 )
134        {
135          MP_INT num, den;
136          On(SW_RATIONAL);
137          mpz_init_set( &num, &(pGetCoeff( p )->z) );
138          mpz_init_set( &den, &(pGetCoeff( p )->n) );
139          term = make_cf( num, den, false );
140        }
141        else
142        { // assume s == 0
143          MP_INT num, den;
144          mpz_init_set( &num, &(pGetCoeff( p )->z) );
145          mpz_init_set( &den, &(pGetCoeff( p )->n) );
146          term = make_cf( num, den, true );
147        }
148      }
149    }
150    for ( int i = 1; i <= n; i++ )
151    {
152      if ( (e = p_GetExp( p, i, r )) != 0 )
153         term *= power( Variable( i ), e );
154    }
155    result += term;
156    p = pNext( p );
157  }
158  return result;
159}
160
161poly convClapPSingP ( const CanonicalForm & f )
162{
163//    cerr << " f = " << f << endl;
164  int n = pVariables+1;
165  /* ASSERT( level( f ) <= pVariables, "illegal number of variables" ); */
166  int * exp = new int[n];
167  //for ( int i = 0; i < n; i++ ) exp[i] = 0;
168  memset(exp,0,n*sizeof(int));
169  poly result = NULL;
170  convRecPP( f, exp, result );
171  delete [] exp;
172  return result;
173}
174
175static void convRecPP ( const CanonicalForm & f, int * exp, poly & result )
176{
177  if (f == 0)
178    return;
179  if ( ! f.inCoeffDomain() )
180  {
181    int l = f.level();
182    for ( CFIterator i = f; i.hasTerms(); i++ )
183    {
184      exp[l] = i.exp();
185      convRecPP( i.coeff(), exp, result );
186    }
187    exp[l] = 0;
188  }
189  else
190  {
191    poly term = pInit();
192    pNext( term ) = NULL;
193    for ( int i = 1; i <= pVariables; i++ )
194      pSetExp( term, i, exp[i]);
195    pSetComp(term, 0);
196    if ( getCharacteristic() != 0 )
197    {
198      pGetCoeff( term ) = nInit( f.intval() );
199    }
200    else
201    {
202      if ( f.isImm() )
203        pGetCoeff( term ) = nInit( f.intval() );
204      else
205      {
206        number z=(number)omAllocBin(rnumber_bin);
207#if defined(LDEBUG)
208        z->debug=123456;
209#endif
210        z->z = gmp_numerator( f );
211        if ( f.den() == 1 )
212          z->s = 3;
213        else
214        {
215          z->n = gmp_denominator( f );
216          z->s = 0;
217        }
218        nlNormalize(z);
219        pGetCoeff( term ) = z;
220      }
221    }
222    pSetm( term );
223    if ( nIsZero(pGetCoeff(term)) )
224    {
225      pDelete(&term);
226    }
227    else
228    {
229      result = pAdd( result, term );
230    }
231  }
232}
233
234poly conv_ClapPSingP ( const CanonicalForm & f, ring r )
235{
236//    cerr << " f = " << f << endl;
237  int n = r->N+1;
238  /* ASSERT( level( f ) <= pVariables, "illegal number of variables" ); */
239  int * exp = new int[n];
240  //for ( int i = 0; i < n; i++ ) exp[i] = 0;
241  memset(exp,0,n*sizeof(int));
242  poly result = NULL;
243  conv_RecPP( f, exp, result, r );
244  delete [] exp;
245  return result;
246}
247
248static void
249conv_RecPP ( const CanonicalForm & f, int * exp, poly & result, ring r )
250{
251  if (f == 0)
252    return;
253  if ( ! f.inCoeffDomain() )
254  {
255    int l = f.level();
256    for ( CFIterator i = f; i.hasTerms(); i++ )
257    {
258      exp[l] = i.exp();
259      conv_RecPP( i.coeff(), exp, result, r );
260    }
261    exp[l] = 0;
262  }
263  else
264  {
265    poly term = p_Init(r);
266    pNext( term ) = NULL;
267    for ( int i = 1; i <= r->N; i++ )
268      p_SetExp( term, i, exp[i], r);
269    if (rRing_has_Comp(r)) p_SetComp(term, 0, r);
270    if ( getCharacteristic() != 0 )
271    {
272      pGetCoeff( term ) = n_Init( f.intval(), r );
273    }
274    else
275    {
276      if ( f.isImm() )
277        pGetCoeff( term ) = n_Init( f.intval(), r );
278      else
279      {
280        number z=(number)omAllocBin(rnumber_bin);
281#if defined(LDEBUG)
282        z->debug=123456;
283#endif
284        z->z = gmp_numerator( f );
285        if ( f.den() == 1 )
286          z->s = 3;
287        else
288        {
289          z->n = gmp_denominator( f );
290          z->s = 0;
291        }
292        pGetCoeff( term ) = z;
293        nlNormalize(z);
294      }
295    }
296    p_Setm( term, r );
297    if ( n_IsZero(pGetCoeff(term), r) )
298    {
299      p_Delete(&term,r);
300    }
301    else
302    {
303      result = p_Add_q( result, term, r );
304    }
305  }
306}
307
308
309CanonicalForm
310convSingTrClapP( napoly p )
311{
312  CanonicalForm result = 0;
313  int e, n = napVariables;
314
315  while ( p!=NULL)
316  {
317    CanonicalForm term;
318    /* does only work for finite fields */
319    if ( getCharacteristic() != 0 )
320    {
321      term = npInt( napGetCoeff( p ) );
322    }
323    else
324    {
325      //if ( (!(int)(napGetCoeff( p )) & 1 )
326      //&&  ( napGetCoeff( p )->s == 0))
327      //  naNormalize( naGetCoeff( p ) );
328      if ( (long)(napGetCoeff( p )) & 1L )
329        term = nlInt( napGetCoeff( p ) );
330      else
331      {
332        if ( napGetCoeff( p )->s == 3 )
333        {
334          MP_INT dummy;
335          mpz_init_set( &dummy, &(napGetCoeff( p )->z) );
336          term = make_cf( dummy );
337        }
338        else  if ( napGetCoeff( p )->s == 1 )
339        {
340          MP_INT num, den;
341          On(SW_RATIONAL);
342          mpz_init_set( &num, &(napGetCoeff( p )->z) );
343          mpz_init_set( &den, &(napGetCoeff( p )->n) );
344          term = make_cf( num, den, false );
345        }
346        else
347        { // assume s == 0
348          MP_INT num, den;
349          On(SW_RATIONAL);
350          mpz_init_set( &num, &(napGetCoeff( p )->z) );
351          mpz_init_set( &den, &(napGetCoeff( p )->n) );
352          term = make_cf( num, den, true );
353        }
354      }
355    }
356    for ( int i = 1; i <= n; i++ )
357    {
358      if ( (e = napGetExp( p, i )) != 0 )
359        term *= power( Variable( i ), e );
360    }
361    result += term;
362    p = napNext( p );
363  }
364  return result;
365}
366
367napoly
368convClapPSingTr ( const CanonicalForm & f )
369{
370//    cerr << " f = " << f << endl;
371  int n = napVariables+1;
372  /* ASSERT( level( f ) <= napVariables, "illegal number of variables" ); */
373  int * exp = new int[n];
374  // for ( int i = 0; i < n; i++ ) exp[i] = 0;
375  memset(exp,0,n*sizeof(int));
376  napoly result = NULL;
377  convRecPTr( f, exp, result );
378  delete [] exp;
379  return result;
380}
381
382static void
383convRecPTr ( const CanonicalForm & f, int * exp, napoly & result )
384{
385  if (f == 0)
386    return;
387  if ( ! f.inCoeffDomain() )
388  {
389    int l = f.level();
390    for ( CFIterator i = f; i.hasTerms(); i++ )
391    {
392        exp[l] = i.exp();
393        convRecPTr( i.coeff(), exp, result );
394    }
395    exp[l] = 0;
396  }
397  else
398  {
399    napoly term = napNew();
400    // napNext( term ) = NULL; //already done by napNew
401    for ( int i = 1; i <= napVariables; i++ )
402      napSetExp( term, i , exp[i]);
403    if ( getCharacteristic() != 0 )
404    {
405      napGetCoeff( term ) = npInit( f.intval() );
406    }
407    else
408    {
409      if ( f.isImm() )
410        napGetCoeff( term ) = nlInit( f.intval() );
411      else
412      {
413        number z=(number)omAllocBin(rnumber_bin);
414#if defined(LDEBUG)
415        z->debug=123456;
416#endif
417        z->z = gmp_numerator( f );
418        if ( f.den() == 1 )
419        {
420          z->s = 3;
421        }
422        else
423        {
424          z->n = gmp_denominator( f );
425          if (mpz_cmp_si(&z->n,(long)1)==0)
426          {
427            mpz_clear(&z->n);
428            z->s=3;
429          }
430          else
431          {
432            z->s = 0;
433          }
434          nacNormalize(z);
435        }
436        napGetCoeff( term ) = z;
437      }
438    }
439    if (nacIsZero(pGetCoeff(term)))
440    {
441      napDelete(&term);
442    }
443    else
444    {
445      result = napAdd( result, term );
446    }
447  }
448}
449
450CanonicalForm convSingAPClapAP ( poly p , const Variable & a)
451{
452  CanonicalForm result = 0;
453  int e, n = pVariables;
454  int off=rPar(currRing);
455
456  On(SW_RATIONAL);
457  while ( p!=NULL)
458  {
459    CanonicalForm term=convSingAClapA(((lnumber)pGetCoeff(p))->z,a);
460    for ( int i = 1; i <= n; i++ )
461    {
462      if ( (e = pGetExp( p, i )) != 0 )
463        term *= power( Variable( i + off), e );
464    }
465    result += term;
466    p = pNext( p );
467  }
468  return result;
469}
470
471static void
472convRecAP_R ( const CanonicalForm & f, int * exp, poly & result, int par_start, int var_start) ;
473
474poly convClapAPSingAP_R ( const CanonicalForm & f, int par_start, int var_start )
475{
476  int n = pVariables+1+rPar(currRing);
477  int * exp = new int[n];
478  // for ( int i = 0; i < n; i++ ) exp[i] = 0;
479  memset(exp,0,n*sizeof(int));
480  poly result = NULL;
481  convRecAP_R( f, exp, result,par_start, var_start );
482  delete [] exp;
483  return result;
484}
485poly convClapAPSingAP ( const CanonicalForm & f )
486{
487  return convClapAPSingAP_R(f,0,rPar(currRing));
488}
489
490static void
491convRecAP_R ( const CanonicalForm & f, int * exp, poly & result, int par_start, int var_start )
492{
493  if (f == 0)
494    return;
495  if ( ! f.inCoeffDomain() )
496  {
497    int l = f.level();
498    for ( CFIterator i = f; i.hasTerms(); i++ )
499    {
500      exp[l] = i.exp();
501      convRecAP_R( i.coeff(), exp, result, par_start, var_start );
502    }
503    exp[l] = 0;
504  }
505  else
506  {
507    napoly z=(napoly)convClapASingA( f );
508    if (z!=NULL)
509    {
510      poly term = pInit();
511      pNext( term ) = NULL;
512      int i;
513      for ( i = 1; i <= pVariables; i++ )
514        pSetExp( term, i , exp[i+var_start]);
515      pSetComp(term, 0);
516      if (par_start==0)
517      {
518        for ( i = 1; i <= var_start; i++ )
519        //z->e[i-1]+=exp[i];
520          napAddExp(z,i,exp[i]);
521      }
522      else
523      {
524        for ( i = par_start+1; i <= var_start+rPar(currRing); i++ )
525        //z->e[i-1]+=exp[i];
526          napAddExp(z,i,exp[i-par_start]);
527      }
528      pGetCoeff(term)=(number)omAlloc0Bin(rnumber_bin);
529      ((lnumber)pGetCoeff(term))->z=z;
530      pSetm( term );
531      result = pAdd( result, term );
532    }
533  }
534}
535
536CanonicalForm convSingAClapA ( napoly p , const Variable & a )
537{
538  CanonicalForm result = 0;
539  int e;
540
541  while ( p!=NULL )
542  {
543    CanonicalForm term;
544    /* does only work for finite fields:Z/p */
545    if ( rField_is_Zp_a() )
546    {
547      term = npInt( napGetCoeff( p ) );
548    }
549    else
550    {
551      //if ( (!(int)(napGetCoeff( p )) & 1 )
552      //&&  ( napGetCoeff( p )->s == 0))
553      //  naNormalize( naGetCoeff( p ) );
554      if ( (long)(napGetCoeff( p )) & SR_INT )
555        term = nlInt( napGetCoeff( p ) );
556        //term = SR_TO_INT(napGetCoeff( p )) ;
557      else
558      {
559        if ( napGetCoeff( p )->s == 3 )
560        {
561          MP_INT dummy;
562          mpz_init_set( &dummy, &(napGetCoeff( p )->z) );
563          term = make_cf( dummy );
564        }
565        else  if ( napGetCoeff( p )->s == 1 )
566        {
567          MP_INT num, den;
568          On(SW_RATIONAL);
569          mpz_init_set( &num, &(napGetCoeff( p )->z) );
570          mpz_init_set( &den, &(napGetCoeff( p )->n) );
571          term = make_cf( num, den, false );
572        }
573        else
574        { // assume s == 0
575          MP_INT num, den;
576          On(SW_RATIONAL);
577          mpz_init_set( &num, &(napGetCoeff( p )->z) );
578          mpz_init_set( &den, &(napGetCoeff( p )->n) );
579          term = make_cf( num, den, true );
580        }
581      }
582    }
583    if ( (e = napGetExp( p, 1 )) != 0 )
584      term *= power( a , e );
585    result += term;
586    p = napNext( p );
587  }
588  return result;
589}
590
591static number convClapNSingAN( const CanonicalForm &f)
592{
593  if ((getCharacteristic() != 0) || ( f.isImm() ))
594    return nacInit( f.intval() );
595  else
596  {
597    number z=(number)omAllocBin(rnumber_bin);
598#if defined(LDEBUG)
599    z->debug=123456;
600#endif
601    z->z = gmp_numerator( f );
602    if ( f.den() == 1 )
603    {
604      z->s = 3;
605    }
606    else
607    {
608      z->n = gmp_denominator( f );
609      z->s = 0;
610    }
611    nlNormalize(z);
612    #ifdef LDEBUG
613    nlDBTest(z,__FILE__,__LINE__);
614    #endif
615    return z;
616  }
617}
618
619napoly convClapASingA ( const CanonicalForm & f )
620{
621  napoly a=NULL;
622  napoly t;
623  for( CFIterator i=f; i.hasTerms(); i++)
624  {
625    t=napNew();
626    // napNext( t ) = NULL; //already done by napNew
627    napGetCoeff(t)=convClapNSingAN( i.coeff() );
628    if (nacIsZero(napGetCoeff(t)))
629    {
630      napDelete(&t);
631    }
632    else
633    {
634      napSetExp(t,1,i.exp());
635      a=napAdd(a,t);
636    }
637  }
638  if (a!=NULL)
639  {
640    if (naMinimalPoly!=NULL)
641    {
642      if (napGetExp(a,1) >= napGetExp(naMinimalPoly,1))
643        a = napRemainder( a, naMinimalPoly);
644    }
645  }
646  return a;
647}
648
649CanonicalForm convSingTrPClapP ( poly p )
650{
651  CanonicalForm result = 0;
652  int e, n = pVariables;
653  int offs = rPar(currRing);
654
655  while ( p!=NULL )
656  {
657    nNormalize(pGetCoeff(p));
658    CanonicalForm term=convSingTrClapP(((lnumber)pGetCoeff(p))->z);
659
660    if ((((lnumber)pGetCoeff(p))->n!=NULL)
661    && (!errorreported))
662    {
663      WerrorS("conversion error: denominator!= 1");
664    }
665
666    for ( int i = 1; i <= n; i++ )
667    {
668      if ( (e = pGetExp( p, i )) != 0 )
669        term = term * power( Variable( i + offs ), e );
670    }
671    result += term;
672    p = pNext( p );
673  }
674  return result;
675}
676
677poly convClapPSingTrP ( const CanonicalForm & f )
678{
679  int n = pVariables+1;
680  int * exp = new int[n];
681  // for ( int i = 0; i < n; i++ ) exp[i] = 0;
682  memset(exp,0,n*sizeof(int));
683  poly result = NULL;
684  convRecTrP( f, exp, result , rPar(currRing) );
685  delete [] exp;
686  return result;
687}
688
689static void
690convRecTrP ( const CanonicalForm & f, int * exp, poly & result , int offs)
691{
692  if (f == 0)
693    return;
694  if ( f.level() > offs )
695  {
696    int l = f.level();
697    for ( CFIterator i = f; i.hasTerms(); i++ )
698    {
699      exp[l-offs] = i.exp();
700      convRecTrP( i.coeff(), exp, result, offs );
701    }
702    exp[l-offs] = 0;
703  }
704  else
705  {
706    poly term = pInit();
707    pNext( term ) = NULL;
708    for ( int i = 1; i <= pVariables; i++ )
709      pSetExp( term, i ,exp[i]);
710    pSetComp(term, 0);
711    pGetCoeff(term)=(number)omAlloc0Bin(rnumber_bin);
712    ((lnumber)pGetCoeff(term))->z=convClapPSingTr( f );
713    pSetm( term );
714    result = pAdd( result, term );
715  }
716}
717
718number   nlChineseRemainder(number *x, number *q,int rl)
719// elemenst in the array are x[0..(rl-1)], q[0..(rl-1)]
720{
721#ifdef HAVE_FACTORY
722  CFArray X(rl), Q(rl);
723  int i;
724  for(i=rl-1;i>=0;i--)
725  {
726    X[i]=CanonicalForm(nlInt(x[i]));
727    Q[i]=CanonicalForm(nlInt(q[i]));
728  }
729  CanonicalForm xnew,qnew;
730  chineseRemainder(X,Q,xnew,qnew);
731  number n=convClapNSingN(xnew);
732  number p=convClapNSingN(qnew);
733  number p2=nlIntDiv(p,nlInit(2));
734  if (nlGreater(n,p2))
735  {
736     number n2=nlSub(n,p);
737     nlDelete(&n,currRing);
738     n=n2;
739  }
740  nlDelete(&p,currRing);
741  nlDelete(&p2,currRing);
742  return n;
743#else
744  WerrorS("not implemented");
745  return nlInit(0);
746#endif
747}
748
749#if 0
750CanonicalForm
751convSingGFClapGF( poly p )
752{
753  CanonicalForm result = 0;
754  int e, n = pVariables;
755
756  while ( p != NULL )
757  {
758    CanonicalForm term;
759    term = make_cf_from_gf( pGetCoeff( p ) );
760    for ( int i = 1; i <= n; i++ )
761    {
762      if ( (e = pGetExp( p, i )) != 0 )
763         term *= power( Variable( i ), e );
764    }
765    result += term;
766    p = pNext( p );
767  }
768  return result;
769}
770
771poly
772convClapGFSingGF ( const CanonicalForm & f )
773{
774//    cerr << " f = " << f << endl;
775  int n = pVariables+1;
776  /* ASSERT( level( f ) <= pVariables, "illegal number of variables" ); */
777  int * exp = new int[n];
778  //for ( int i = 0; i < n; i++ ) exp[i] = 0;
779  memset(exp,0,n*sizeof(int));
780  poly result = NULL;
781  convRecGFGF( f, exp, result );
782  delete [] exp;
783  return result;
784}
785
786static void
787convRecGFGF ( const CanonicalForm & f, int * exp, poly & result )
788{
789  if (f == 0)
790    return;
791  if ( ! f.inCoeffDomain() )
792  {
793    int l = f.level();
794    for ( CFIterator i = f; i.hasTerms(); i++ )
795    {
796      exp[l] = i.exp();
797      convRecGFGF( i.coeff(), exp, result );
798    }
799    exp[l] = 0;
800  }
801  else
802  {
803    poly term = pInit();
804    pNext( term ) = NULL;
805    for ( int i = 1; i <= pVariables; i++ )
806      pSetExp( term, i, exp[i]);
807    pSetComp(term, 0);
808    pGetCoeff( term ) = (number) gf_value (f);
809    pSetm( term );
810    result = pAdd( result, term );
811  }
812}
813#endif
814
815int convClapISingI( const CanonicalForm & f)
816{
817  if (!f.isImm()) WerrorS("int overflow in det");
818  return f.intval();
819}
820#endif
Note: See TracBrowser for help on using the repository browser.