Loading...
Searching...
No Matches
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 *******************************************************
// checking 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;
}
m
int m
Definition cfEzgcd.cc:128
i
int i
Definition cfEzgcd.cc:132
p
int p
Definition cfModGcd.cc:4086
cf
CanonicalForm cf
Definition cfModGcd.cc:4091
cf_getNumSmallPrimes
int cf_getNumSmallPrimes()
Definition cf_primes.cc:34
ALLOC_RNUMBER
#define ALLOC_RNUMBER()
Definition coeffs.h:94
n_Normalize
static FORCE_INLINE void n_Normalize(number &n, const coeffs r)
inplace-normalization of n; produces some canonical representation of n;
Definition coeffs.h:581
Print
#define Print
Definition emacs.cc:80
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39
j
int j
Definition facHensel.cc:110
WerrorS
void WerrorS(const char *s)
Definition feFopen.cc:24
TakePrime
static modp_number TakePrime(modp_number)
Definition interpolation.cc:1038
multiplicity
STATIC_VAR int * multiplicity
Definition interpolation.cc:86
ClearGenList
static void ClearGenList()
Definition interpolation.cc:1325
CalcBaseDim
static int CalcBaseDim()
Definition interpolation.cc:250
myp_index
STATIC_VAR int myp_index
Definition interpolation.cc:34
UpdateGenList
static void UpdateGenList()
Definition interpolation.cc:1344
FreeProcData
static void FreeProcData()
Definition interpolation.cc:501
bigcongr
STATIC_VAR mpz_t bigcongr
Definition interpolation.cc:111
ResolveCoeff
static void ResolveCoeff(mpq_t c, number m)
Definition interpolation.cc:1462
IntegerPoints
static void IntegerPoints()
Definition interpolation.cc:599
MakeConditions
static void MakeConditions()
Definition interpolation.cc:665
CloseChinese
static void CloseChinese()
Definition interpolation.cc:1071
NewResultEntry
static void NewResultEntry()
Definition interpolation.cc:937
PrepareChinese
static void PrepareChinese(int n)
Definition interpolation.cc:1044
myp
STATIC_VAR modp_number myp
Definition interpolation.cc:33
q_points
STATIC_VAR q_coordinates * q_points
Definition interpolation.cc:92
ReconstructGenerator
static void ReconstructGenerator(int ngen, int n)
Definition interpolation.cc:1101
modp_result
STATIC_VAR modp_result_entry * modp_result
Definition interpolation.cc:107
modp_PrepareProducts
static void modp_PrepareProducts()
Definition interpolation.cc:651
n_points
STATIC_VAR int n_points
Definition interpolation.cc:85
modp_points
STATIC_VAR modp_coordinates * modp_points
Definition interpolation.cc:91
only_modp
STATIC_VAR bool only_modp
Definition interpolation.cc:89
CheckGenerator
static bool CheckGenerator()
Definition interpolation.cc:1298
InitProcData
static void InitProcData()
Definition interpolation.cc:399
mono_type
exponent * mono_type
Definition interpolation.cc:47
final_base_dim
STATIC_VAR int final_base_dim
Definition interpolation.cc:83
MultGenerators
static void MultGenerators()
Definition interpolation.cc:1000
generic_column_name
STATIC_VAR mono_type * generic_column_name
Definition interpolation.cc:124
GeneralInit
static void GeneralInit()
Definition interpolation.cc:338
gen_list
STATIC_VAR gen_list_entry * gen_list
Definition interpolation.cc:121
variables
STATIC_VAR int variables
Definition interpolation.cc:87
modp_Main
static void modp_Main()
Definition interpolation.cc:1405
int_PrepareProducts
static void int_PrepareProducts()
Definition interpolation.cc:627
GeneralDone
static void GeneralDone()
Definition interpolation.cc:449
CheckColumnSequence
static void CheckColumnSequence()
Definition interpolation.cc:1222
generic_n_generators
STATIC_VAR int generic_n_generators
Definition interpolation.cc:123
protocol
STATIC_VAR bool protocol
Definition interpolation.cc:136
n_results
STATIC_VAR int n_results
Definition interpolation.cc:105
coord_exist
STATIC_VAR coord_exist_table * coord_exist
Definition interpolation.cc:94
modp_SetColumnNames
static void modp_SetColumnNames()
Definition interpolation.cc:1216
assume
#define assume(x)
Definition mod2.h:389
pNext
#define pNext(p)
Definition monomials.h:36
pGetCoeff
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
nDiv
#define nDiv(a, b)
Definition numbers.h:32
nInpNeg
#define nInpNeg(n)
Definition numbers.h:21
nInit
#define nInit(i)
Definition numbers.h:24
omAlloc
#define omAlloc(size)
Definition omAllocDecl.h:210
NULL
#define NULL
Definition omList.c:12
TEST_OPT_PROT
#define TEST_OPT_PROT
Definition options.h:105
currRing
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition polys.cc:13
pAdd
#define pAdd(p, q)
Definition polys.h:204
pDelete
#define pDelete(p_ptr)
Definition polys.h:187
pHead
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL
Definition polys.h:68
pSetm
#define pSetm(p)
Definition polys.h:272
pNSet
#define pNSet(n)
Definition polys.h:314
pVar
#define pVar(m)
Definition polys.h:381
wrp
void wrp(poly p)
Definition polys.h:311
pSetExp
#define pSetExp(p, i, v)
Definition polys.h:43
pISet
#define pISet(i)
Definition polys.h:313
PrintS
void PrintS(const char *s)
Definition reporter.cc:288
PrintLn
void PrintLn()
Definition reporter.cc:314
rChar
int rChar(ring r)
Definition ring.cc:719
rField_is_Zp
static BOOLEAN rField_is_Zp(const ring r)
Definition ring.h:506
rField_is_Q
static BOOLEAN rField_is_Q(const ring r)
Definition ring.h:512
idInit
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition simpleideals.cc:35
IDELEMS
#define IDELEMS(i)
Definition simpleideals.h:23