My Project
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 1484 of file interpolation.cc.

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