My Project  debian-1:4.1.2-p1+ds-2
Functions
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 1483 of file interpolation.cc.

1484 {
1485  protocol=TEST_OPT_PROT; // should be set if option(prot) is enabled
1486 
1487  bool data_ok=true;
1488 
1489  // reading the ring data ***************************************************
1490  if ((currRing==NULL) || ((!rField_is_Zp (currRing))&&(!rField_is_Q (currRing))))
1491  {
1492  WerrorS("coefficient field should be Zp or Q!");
1493  return NULL;
1494  }
1495  if ((currRing->qideal)!=NULL)
1496  {
1497  WerrorS("quotient ring not supported!");
1498  return NULL;
1499  }
1500  if ((currRing->OrdSgn)!=1)
1501  {
1502  WerrorS("ordering must be global!");
1503  return NULL;
1504  }
1505  n_points=v->length ();
1506  if (n_points!=L.size())
1507  {
1508  WerrorS("list and intvec must have the same length!");
1509  return NULL;
1510  }
1511  assume( n_points > 0 );
1512  variables=currRing->N;
1514  if (only_modp) myp=rChar(currRing);
1515  // ring data read **********************************************************
1516 
1517 
1518  multiplicity=(int*)malloc(sizeof(int)*n_points); // TODO: use omalloc!
1519  int i;
1520  for (i=0;i<n_points;i++) multiplicity[i]=(*v)[i];
1521 
1523 
1524 #ifdef writemsg
1525  Print("number of variables: %d\n", variables);
1526  Print("number of points: %d\n", n_points);
1527  PrintS("multiplicities: ");
1528  for (i=0;i<n_points;i++) Print("%d ", multiplicity[i]);
1529  PrintLn();
1530  Print("general initialization for dimension %d ...\n", final_base_dim);
1531 #endif
1532 
1533  GeneralInit ();
1534 
1535 // reading coordinates of points from ideals **********************************
1536  mpq_t divisor;
1537  if (!only_modp) mpq_init(divisor);
1538  int j;
1539  for(i=0; i<L.size();i++)
1540  {
1541  ideal I = L[i];
1542  for(j=0;j<IDELEMS(I);j++)
1543  {
1544  poly p=I->m[j];
1545  if (p!=NULL)
1546  {
1547  poly ph=pHead(p);
1548  int pcvar=pVar(ph);
1549  if (pcvar!=0)
1550  {
1551  pcvar--;
1552  if (coord_exist[i][pcvar])
1553  {
1554  Print("coordinate %d for point %d initialized twice!\n",pcvar+1,i+1);
1555  data_ok=false;
1556  }
1557  number m;
1558  m=pGetCoeff(p); // possible coefficient standing by a leading monomial
1559  if (!only_modp) ResolveCoeff (divisor,m);
1560  number n;
1561  if (pNext(p)!=NULL) n=pGetCoeff(pNext(p));
1562  else n=nInit(0);
1563  if (only_modp)
1564  {
1565  n=nInpNeg(n);
1566  n=nDiv(n,m);
1567  modp_points[i][pcvar]=(int)((long)n);
1568  }
1569  else
1570  {
1571  ResolveCoeff (q_points[i][pcvar],n);
1572  mpq_neg(q_points[i][pcvar],q_points[i][pcvar]);
1573  mpq_div(q_points[i][pcvar],q_points[i][pcvar],divisor);
1574  }
1575  coord_exist[i][pcvar]=true;
1576  }
1577  else
1578  {
1579  PrintS("not a variable? ");
1580  wrp(p);
1581  PrintLn();
1582  data_ok=false;
1583  }
1584  pDelete(&ph);
1585  }
1586  }
1587  }
1588  if (!only_modp) mpq_clear(divisor);
1589  // data from ideal read *******************************************************
1590 
1591  // ckecking if all coordinates are initialized
1592  for (i=0;i<n_points;i++)
1593  {
1594  for (j=0;j<variables;j++)
1595  {
1596  if (!coord_exist[i][j])
1597  {
1598  Print("coordinate %d for point %d not known!\n",j+1,i+1);
1599  data_ok=false;
1600  }
1601  }
1602  }
1603 
1604  if (!data_ok)
1605  {
1606  GeneralDone();
1607  WerrorS("data structure is invalid");
1608  return NULL;
1609  }
1610 
1611  if (!only_modp) IntegerPoints ();
1612  MakeConditions ();
1613 #ifdef writemsg
1614  PrintS("done.\n");
1615 #else
1616  if (protocol) Print("[vdim %d]",final_base_dim);
1617 #endif
1618 
1619 
1620 // main procedure *********************************************************************
1621  int modp_cycles=10;
1622  bool correct_gen=false;
1623  if (only_modp) modp_cycles=1;
1625 
1626  while ((!correct_gen)&&(myp_index>1))
1627  {
1628 #ifdef writemsg
1629  Print("trying %d cycles mod p...\n",modp_cycles);
1630 #else
1631  if (protocol) Print("(%d)",modp_cycles);
1632 #endif
1633  while ((n_results<modp_cycles)&&(myp_index>1)) // some computations mod p
1634  {
1635  if (!only_modp) myp=TakePrime (myp);
1636  NewResultEntry ();
1637  InitProcData ();
1639 
1640  modp_Main ();
1641 
1642  if (!only_modp)
1643  {
1644  MultGenerators ();
1646  }
1647  else
1648  {
1650  }
1651  FreeProcData ();
1652  }
1653 
1654  if (!only_modp)
1655  {
1656  PrepareChinese (modp_cycles);
1657  correct_gen=true;
1658  for (i=0;i<generic_n_generators;i++)
1659  {
1660  ReconstructGenerator (i,modp_cycles);
1661  correct_gen=CheckGenerator ();
1662  if (!correct_gen)
1663  {
1664 #ifdef writemsg
1665  PrintS("wrong generator!\n");
1666 #else
1667 // if (protocol) PrintS("!g");
1668 #endif
1669  ClearGenList ();
1670  break;
1671  }
1672  else
1673  {
1674  UpdateGenList ();
1675  }
1676  }
1677 #ifdef checksize
1678  Print("maximal size of output: %d, precision bound: %d.\n",maximalsize,mpz_sizeinbase(bigcongr,10));
1679 #endif
1680  CloseChinese ();
1681  modp_cycles=modp_cycles+10;
1682  }
1683  else
1684  {
1685  correct_gen=true;
1686  }
1687  }
1688 // end of main procedure ************************************************************************************
1689 
1690 #ifdef writemsg
1691  PrintS("computations finished.\n");
1692 #else
1693  if (protocol) PrintLn();
1694 #endif
1695 
1696  if (!correct_gen)
1697  {
1698  GeneralDone ();
1699  ClearGenList ();
1700  WerrorS("internal error - coefficient too big!");
1701  return NULL;
1702  }
1703 
1704 // passing data to ideal *************************************************************************************
1705  ideal ret;
1706 
1707  if (only_modp)
1708  {
1709  mono_type mon;
1710  ret=idInit(modp_result->n_generators,1);
1711  generator_entry *cur_gen=modp_result->generator;
1712  for(i=0;i<IDELEMS(ret);i++)
1713  {
1714  poly p,sum;
1715  sum=NULL;
1716  int a;
1717  int cf;
1718  for (a=final_base_dim;a>=0;a--)
1719  {
1720  if (a==final_base_dim) cf=cur_gen->ltcoef; else cf=cur_gen->coef[a];
1721  if (cf!=0)
1722  {
1723  p=pISet(cf);
1724  if (a==final_base_dim) mon=cur_gen->lt; else mon=generic_column_name[a];
1725  for (j=0;j<variables;j++) pSetExp(p,j+1,mon[j]);
1726  pSetm(p);
1727  sum=pAdd(sum,p);
1728  }
1729  }
1730  ret->m[i]=sum;
1731  cur_gen=cur_gen->next;
1732  }
1733  }
1734  else
1735  {
1737  gen_list_entry *temp=gen_list;
1738  for(i=0;i<IDELEMS(ret);i++)
1739  {
1740  poly p,sum;
1741  sum=NULL;
1742  int a;
1743  for (a=final_base_dim;a>=0;a--) // build one polynomial
1744  {
1745  if (mpz_sgn(temp->polycoef[a])!=0)
1746  {
1747  number n=ALLOC_RNUMBER();
1748 #ifdef LDEBUG
1749  n->debug=123456;
1750 #endif
1751  mpz_init_set(n->z,temp->polycoef[a]);
1752  n->s=3;
1753  n_Normalize(n, currRing->cf);
1754  p=pNSet(n); //a monomial
1755  for (j=0;j<variables;j++) pSetExp(p,j+1,temp->polyexp[a][j]);
1756  pSetm(p); // after all pSetExp
1757  sum=pAdd(sum,p);
1758  }
1759  }
1760  ret->m[i]=sum;
1761  temp=temp->next;
1762  }
1763  }
1764 // data transferred ****************************************************************************
1765 
1766 
1767  GeneralDone ();
1768  ClearGenList ();
1769  return ret;
1770 }
int m
Definition: cfEzgcd.cc:121
int i
Definition: cfEzgcd.cc:125
int p
Definition: cfModGcd.cc:4019
CanonicalForm cf
Definition: cfModGcd.cc:4024
int cf_getNumSmallPrimes()
Definition: cf_primes.cc:34
#define ALLOC_RNUMBER()
Definition: coeffs.h:87
static FORCE_INLINE void n_Normalize(number &n, const coeffs r)
inplace-normalization of n; produces some canonical representation of n;
Definition: coeffs.h:578
#define Print
Definition: emacs.cc:80
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
int j
Definition: facHensel.cc:105
void WerrorS(const char *s)
Definition: feFopen.cc:24
static modp_number TakePrime(modp_number)
STATIC_VAR int * multiplicity
static void ClearGenList()
static int CalcBaseDim()
STATIC_VAR int myp_index
static void UpdateGenList()
static void FreeProcData()
STATIC_VAR mpz_t bigcongr
static void ResolveCoeff(mpq_t c, number m)
static void IntegerPoints()
static void MakeConditions()
static void CloseChinese()
static void NewResultEntry()
static void PrepareChinese(int n)
STATIC_VAR modp_number myp
STATIC_VAR q_coordinates * q_points
static void ReconstructGenerator(int ngen, int n)
STATIC_VAR modp_result_entry * modp_result
static void modp_PrepareProducts()
STATIC_VAR int n_points
STATIC_VAR modp_coordinates * modp_points
STATIC_VAR bool only_modp
static bool CheckGenerator()
static void InitProcData()
exponent * mono_type
STATIC_VAR int final_base_dim
static void MultGenerators()
STATIC_VAR mono_type * generic_column_name
static void GeneralInit()
STATIC_VAR gen_list_entry * gen_list
STATIC_VAR int variables
static void modp_Main()
static void int_PrepareProducts()
static void GeneralDone()
static void CheckColumnSequence()
STATIC_VAR int generic_n_generators
STATIC_VAR bool protocol
STATIC_VAR int n_results
STATIC_VAR coord_exist_table * coord_exist
static void modp_SetColumnNames()
#define assume(x)
Definition: mod2.h:390
#define pNext(p)
Definition: monomials.h:36
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
#define nDiv(a, b)
Definition: numbers.h:32
#define nInpNeg(n)
Definition: numbers.h:21
#define nInit(i)
Definition: numbers.h:24
#define NULL
Definition: omList.c:12
void * malloc(size_t size)
Definition: omalloc.c:92
#define TEST_OPT_PROT
Definition: options.h:101
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
#define pAdd(p, q)
Definition: polys.h:199
#define pDelete(p_ptr)
Definition: polys.h:182
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL
Definition: polys.h:67
#define pSetm(p)
Definition: polys.h:267
#define pNSet(n)
Definition: polys.h:309
#define pVar(m)
Definition: polys.h:377
void wrp(poly p)
Definition: polys.h:306
#define pSetExp(p, i, v)
Definition: polys.h:42
#define pISet(i)
Definition: polys.h:308
void PrintS(const char *s)
Definition: reporter.cc:284
void PrintLn()
Definition: reporter.cc:310
int rChar(ring r)
Definition: ring.cc:713
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:495
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:501
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
#define IDELEMS(i)
Definition: simpleideals.h:23