interpolation.h File Reference

#include "polys/monomials/ring.h"
#include "misc/intvec.h"
#include <vector>

Go to the source code of this file.

Functions
ideal	interpolation (const std::vector< ideal > &L, intvec *v)

Function Documentation

interpolation()

ideal interpolation	(	const std::vector< ideal > &	L,
		intvec *	v
	)

Definition at line 1484 of file interpolation.cc.

{
  protocol=TEST_OPT_PROT;  // should be set if option(prot) is enabled
 
  bool data_ok=true;
 
  // reading the ring data ***************************************************
  if ((currRing==NULL) || ((!rField_is_Zp (currRing))&&(!rField_is_Q (currRing))))
  {
     WerrorS("coefficient field should be Zp or Q!");
     return NULL;
  }
  if ((currRing->qideal)!=NULL)
  {
     WerrorS("quotient ring not supported!");
     return NULL;
  }
  if ((currRing->OrdSgn)!=1)
  {
     WerrorS("ordering must be global!");
     return NULL;
  }
  n_points=v->length ();
  if (n_points!=L.size())
  {
     WerrorS("list and intvec must have the same length!");
     return NULL;
  }
  assume( n_points > 0 );
  variables=currRing->N;
  only_modp=rField_is_Zp(currRing);
  if (only_modp) myp=rChar(currRing);
  // ring data read **********************************************************
 
 
  multiplicity=(int*)omAlloc(sizeof(int)*n_points);
  int i;
  for (i=0;i<n_points;i++) multiplicity[i]=(*v)[i];
 
  final_base_dim = CalcBaseDim ();
 
#ifdef writemsg
  Print("number of variables: %d\n", variables);
  Print("number of points: %d\n", n_points);
  PrintS("multiplicities: ");
  for (i=0;i<n_points;i++) Print("%d ", multiplicity[i]);
  PrintLn();
  Print("general initialization for dimension %d ...\n", final_base_dim);
#endif
 
  GeneralInit ();
 
// reading coordinates of points from ideals **********************************
  mpq_t divisor;
  if (!only_modp) mpq_init(divisor);
  int j;
  for(i=0; i<L.size();i++)
  {
    ideal I = L[i];
    for(j=0;j<IDELEMS(I);j++)
    {
      poly p=I->m[j];
      if (p!=NULL)
      {
        poly ph=pHead(p);
        int pcvar=pVar(ph);
        if (pcvar!=0)
        {
          pcvar--;
          if (coord_exist[i][pcvar])
          {
             Print("coordinate %d for point %d initialized twice!\n",pcvar+1,i+1);
             data_ok=false;
          }
          number m;
          m=pGetCoeff(p); // possible coefficient standing by a leading monomial
          if (!only_modp) ResolveCoeff (divisor,m);
          number n;
          if (pNext(p)!=NULL) n=pGetCoeff(pNext(p));
          else n=nInit(0);
          if (only_modp)
          {
            n=nInpNeg(n);
            n=nDiv(n,m);
            modp_points[i][pcvar]=(int)((long)n);
          }
          else
          {
            ResolveCoeff (q_points[i][pcvar],n);
            mpq_neg(q_points[i][pcvar],q_points[i][pcvar]);
            mpq_div(q_points[i][pcvar],q_points[i][pcvar],divisor);
          }
          coord_exist[i][pcvar]=true;
        }
        else
        {
          PrintS("not a variable? ");
          wrp(p);
          PrintLn();
          data_ok=false;
        }
        pDelete(&ph);
      }
    }
  }
  if (!only_modp) mpq_clear(divisor);
  // data from ideal read *******************************************************
 
  // ckecking if all coordinates are initialized
  for (i=0;i<n_points;i++)
  {
      for (j=0;j<variables;j++)
      {
          if (!coord_exist[i][j])
          {
             Print("coordinate %d for point %d not known!\n",j+1,i+1);
             data_ok=false;
          }
      }
  }
 
  if (!data_ok)
  {
     GeneralDone();
     WerrorS("data structure is invalid");
     return NULL;
  }
 
  if (!only_modp) IntegerPoints ();
  MakeConditions ();
#ifdef writemsg
  PrintS("done.\n");
#else
  if (protocol) Print("[vdim %d]",final_base_dim);
#endif
 
 
// main procedure *********************************************************************
  int modp_cycles=10;
  bool correct_gen=false;
  if (only_modp) modp_cycles=1;
  myp_index=cf_getNumSmallPrimes ();
 
  while ((!correct_gen)&&(myp_index>1))
  {
#ifdef writemsg
        Print("trying %d cycles mod p...\n",modp_cycles);
#else
        if (protocol) Print("(%d)",modp_cycles);
#endif
        while ((n_results<modp_cycles)&&(myp_index>1))  // some computations mod p
        {
            if (!only_modp) myp=TakePrime (myp);
            NewResultEntry ();
            InitProcData ();
            if (only_modp) modp_PrepareProducts (); else int_PrepareProducts ();
 
            modp_Main ();
 
            if (!only_modp)
            {
               MultGenerators ();
               CheckColumnSequence ();
            }
            else
            {
               modp_SetColumnNames ();
            }
            FreeProcData ();
        }
 
        if (!only_modp)
        {
           PrepareChinese (modp_cycles);
           correct_gen=true;
           for (i=0;i<generic_n_generators;i++)
           {
               ReconstructGenerator (i,modp_cycles);
               correct_gen=CheckGenerator ();
               if (!correct_gen)
               {
#ifdef writemsg
                  PrintS("wrong generator!\n");
#else
//                  if (protocol) PrintS("!g");
#endif
                  ClearGenList ();
                  break;
               }
               else
               {
                  UpdateGenList ();
               }
           }
#ifdef checksize
           Print("maximal size of output: %d, precision bound: %d.\n",maximalsize,mpz_sizeinbase(bigcongr,10));
#endif
           CloseChinese ();
           modp_cycles=modp_cycles+10;
        }
        else
        {
           correct_gen=true;
        }
  }
// end of main procedure ************************************************************************************
 
#ifdef writemsg
  PrintS("computations finished.\n");
#else
  if (protocol) PrintLn();
#endif
 
  if (!correct_gen)
  {
     GeneralDone ();
     ClearGenList ();
     WerrorS("internal error - coefficient too big!");
     return NULL;
  }
 
// passing data to ideal *************************************************************************************
  ideal ret;
 
  if (only_modp)
  {
    mono_type mon;
    ret=idInit(modp_result->n_generators,1);
    generator_entry *cur_gen=modp_result->generator;
    for(i=0;i<IDELEMS(ret);i++)
    {
      poly p,sum;
      sum=NULL;
      int a;
      int cf;
      for (a=final_base_dim;a>=0;a--)
      {
        if (a==final_base_dim) cf=cur_gen->ltcoef; else cf=cur_gen->coef[a];
        if (cf!=0)
        {
            p=pISet(cf);
            if (a==final_base_dim) mon=cur_gen->lt; else mon=generic_column_name[a];
            for (j=0;j<variables;j++) pSetExp(p,j+1,mon[j]);
            pSetm(p);
            sum=pAdd(sum,p);
        }
      }
      ret->m[i]=sum;
      cur_gen=cur_gen->next;
    }
  }
  else
  {
    ret=idInit(generic_n_generators,1);
    gen_list_entry *temp=gen_list;
    for(i=0;i<IDELEMS(ret);i++)
    {
      poly p,sum;
      sum=NULL;
      int a;
      for (a=final_base_dim;a>=0;a--) // build one polynomial
      {
          if (mpz_sgn(temp->polycoef[a])!=0)
          {
             number n=ALLOC_RNUMBER();
#ifdef LDEBUG
             n->debug=123456;
#endif
             mpz_init_set(n->z,temp->polycoef[a]);
             n->s=3;
             n_Normalize(n, currRing->cf);
             p=pNSet(n); //a monomial
             for (j=0;j<variables;j++) pSetExp(p,j+1,temp->polyexp[a][j]);
             pSetm(p); // after all pSetExp
             sum=pAdd(sum,p);
          }
      }
      ret->m[i]=sum;
      temp=temp->next;
    }
  }
// data transferred ****************************************************************************
 
 
  GeneralDone ();
  ClearGenList ();
  return ret;
}