My Project
Loading...
Searching...
No Matches
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 );
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:4078
CanonicalForm cf
Definition: cfModGcd.cc:4083
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:575
#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:389
#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:104
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:380
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:713
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:500
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:506
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
#define IDELEMS(i)
Definition: simpleideals.h:23