Sample use of gbasis in C++

Librairie C++ de calcul formel/ C++ symbolic computation library

Modérateur : xcasadmin

Répondre
jocaps
Messages : 118
Inscription : lun. avr. 17, 2017 4:32 pm

Sample use of gbasis in C++

Message par jocaps » dim. janv. 28, 2018 3:41 pm

Hi,

Is there anywhere I could see a minimal example of how I can compute groebner basis using giac in C++?
For instance suppose I have a vecpoly consisting of generator of an ideal (say some katsura polynomial system) and now I want to use gbasis. I see that in solve.cc there is a method called giac_gbasis (and another method simply called gbasis which looks like it uses the option if I have CoCoa). This method is given by

Code : Tout sélectionner

static bool giac_gbasis(vectpoly & res,const gen & order_,environment * env,int modularcheck,int & rur,GIAC_CONTEXT,bool eliminate_flag)
I do not know very well what the following parameters are:
  • environment
  • modularcheck (I just want to work over rationals or complex numbers) ?
  • rur
  • eliminate_flag (I can guess what this is but an explanation would not hurt)
Jose
Dernière modification par jocaps le dim. janv. 28, 2018 10:44 pm, modifié 1 fois.

jocaps
Messages : 118
Inscription : lun. avr. 17, 2017 4:32 pm

Re: Sample use of gbasis in C++

Message par jocaps » dim. janv. 28, 2018 9:26 pm

I have written a sample code testing gbasis on katsura9 (or katsura10 depending on definition) polynomials system.

Code : Tout sélectionner

#include <config.h>
#include <giac.h>

using namespace std;
using namespace giac;


int main()
{
  giac::context ct;
  vecteur x;
  vectpoly ideal;
  stringstream sstr;
  int n=9;

  //variable vector
  for (int i=0; i<n+1;++i)
  {
    stringstream sstr("x"); sstr << i;
    x.push_back(identificateur(sstr.str()));
  }

  // x_0 + 2\sum_{i=1}^n x_i = 1
  gen f = x[0];
  for (int i=1; i<n+1;++i) f += 2*x[i];
  ideal.push_back(polynome(*f._POLYptr));

  // \sum_{j=-n}^{n-i} x_{|j|}x_{|j+i|} = x_{i} \qquad i=0,\dots,n-1
  for (int i=0; i<n;++i)
  {
    f=-x[i];
    for (int j=-n; j< n-i+1;++j) f += x[abs(j)]*x[abs(j+i)];
    f = simplifier(f, &ct);
    ideal.push_back(polynome(*f._POLYptr));
  }

  environment env;
  env.modulo = 0;
  int rur=0;

  gbasis(ideal, x, false, false, rur, &env, &ct, false);
}
This does not compile in Visual Studio because of unresolved external context0 and gbasis (though, this could work if I link statically with respect to the giac library). I will report once I tested this with static links. But in general, can anyone confirm that this compiles in Linux?

Edit: I figured that _gbasis does compile and link properly. Here I do not need to give environment parameter, you basically need the polynomial system and the variables in _gbasis (all enclosed in a vecteur I suppose). But with katsura9 I do not get a proper answer and I think I even get a crash. I am sure I am doing something wrong. Here is the modified code for main():

Code : Tout sélectionner

int main()
{
  giac::context ct;
  vecteur x;
  vecteur ideal;
  stringstream sstr;
  int n=9;
  
  //variable vector
  for (int i=0; i<n+1;++i)
  {
    stringstream sstr("x"); sstr << i;
    x.push_back(identificateur(sstr.str()));
  }
  
  // x_0 + 2\sum_{i=1}^n x_i = 1
  gen f = x[0];
  for (int i=1; i<n+1;++i) f = f+2*x[i];
  ideal.push_back(f);

  
  // \sum_{j=-n}^{n-i} x_{|j|}x_{|j+i|} = x_{i} \qquad i=0,\dots,n-1
  for (int i=0; i<n;++i)
  {
    f=-x[i];
    for (int j=-n; j< n-i+1;++j) f = f+x[abs(j)]*x[abs(j+i)];
    f = simplifier(f, &ct);
    ideal.push_back(f);
  }
  
  _gbasis(makesequence(ideal,x), &ct);  
}
Jose

parisse
Messages : 5731
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: Sample use of gbasis in C++

Message par parisse » lun. janv. 29, 2018 8:27 am

Considering the fact that almost all the time will be spent in computing the Groebner basis, I would recommend that you do not try to use internal functions (static in C++ means that the function is not exported, therefore you will not be able to link with it outside of it's source file). Moreover polynome() does not convert to an internal giac polynome.
It's much easier to run user commands into your C++ program, user commands are accessed with the same name + a _ prefix. Arguments should be grouped in a sequence with makesequence. They can be parsed from a commandline string, like

Code : Tout sélectionner

for (int i=0; i<n+1;++i)
  {
    x.push_back(gen("x"+print_INT_(i),&ct));
  }
The call to gbasis would be _gbasis(makesequence(ideal,x),&ct);

jocaps
Messages : 118
Inscription : lun. avr. 17, 2017 4:32 pm

Re: Sample use of gbasis in C++

Message par jocaps » lun. janv. 29, 2018 10:23 am

Thanks for the tip about polynome. Nevertheless,_gbasis still gives me "not found". I am not sure if this is because of my visual studio build. I will investigate this more as I do not know exactly where the problem lies. The same code in giacpy shows no problem:

Code : Tout sélectionner

from giacpy import giac, gbasis
from time import clock

#katsura n, change n as pleased..
n=9

x=[]
ideal=[]
for i in xrange(n+1):
  s= "x" + str(i)
  x.append(giac(s))

# x_0 + 2\sum_{i=1}^n x_i = 1
f=x[0]
for i in xrange(1,n+1):
  f += 2*x[i]
ideal.append(f-1)

# \sum_{j=-n}^{n-i} x_{|j|}x_{|j+i|} = x_{i} \qquad i=0,\dots,n-1
for i in xrange(n):
  f=-x[i]
  for j in xrange(-n,n-i+1):
    f += x[abs(j)]*x[abs(j+i)]
  ideal.append(f.regroup())

print ideal 

"""
For katsura 9, ideal = 
[x0+2*x1+2*x2+2*x3+2*x4+2*x5+2*x6+2*x7+2*x8+2*x9-1, x0**2+2*x1**2+2*x2**2+2*x3**2+2*x4**2+2*x5**2+2*x6**2+2*x7**2+2*x8**2+2*x9**2-x0, 2*x0*x1+2*x1*x2+2*x2*x3+2*x3*x4+2*x4*x5+2*x5*x6+2*x6*x7+2*x7*x8+2*x8*x9-x1, x1**2+2*x0*x2+2*x1*x3+2*x2*x4+2*x3*x5+2*x4*x6+2*x5*x7+2*x6*x8+2*x7*x9-x2, 2*x0*x3+2*x1*x2+2*x1*x4+2*x2*x5+2*x3*x6+2*x4*x7+2*x5*x8+2*x6*x9-x3, x2**2+2*x0*x4+2*x1*x3+2*x1*x5+2*x2*x6+2*x3*x7+2*x4*x8+2*x5*x9-x4, 2*x0*x5+2*x1*x4+2*x1*x6+2*x2*x3+2*x2*x7+2*x3*x8+2*x4*x9-x5, x3**2+2*x0*x6+2*x1*x5+2*x1*x7+2*x2*x4+2*x2*x8+2*x3*x9-x6, 2*x0*x7+2*x1*x6+2*x1*x8+2*x2*x5+2*x2*x9+2*x3*x4-x7, x4**2+2*x0*x8+2*x1*x7+2*x1*x9+2*x2*x6+2*x3*x5-x8]
"""

start = clock()
gbasis(ideal,x)
print clock()-start, "seconds used"
Edit: I confirm this is a problem from visual studio build. With mingw32 the groebner basis works. I will have to go through debugging to find the problem and report.

jocaps
Messages : 118
Inscription : lun. avr. 17, 2017 4:32 pm

Re: Sample use of gbasis in C++

Message par jocaps » sam. févr. 03, 2018 9:44 am

Is the heart of the groebner basis computation in solve.cc or cocoa.cc? For some reasone the Visual Studio build goes to through the cocoa.cc code and attempts to call convert (line ca. 12295). The reason why I know this is because I get the text "not found" in the console which comes from this line of code.

Jose

parisse
Messages : 5731
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: Sample use of gbasis in C++

Message par parisse » sam. févr. 03, 2018 10:23 am

jocaps a écrit :Is the heart of the groebner basis computation in solve.cc or cocoa.cc?
in cocoa.cc
For some reasone the Visual Studio build goes to through the cocoa.cc code and attempts to call convert (line ca. 12295). The reason why I know this is because I get the text "not found" in the console which comes from this line of code.

Jose
Can you replicate this as a gbasis command?

jocaps
Messages : 118
Inscription : lun. avr. 17, 2017 4:32 pm

Re: Sample use of gbasis in C++

Message par jocaps » dim. févr. 04, 2018 1:37 pm

parisse a écrit : Can you replicate this as a gbasis command?
You mean in Xcas or giacpy? I can only replicated this in visual studio with _gbasis (gbasis command for C++). I feel that the vector is not recognized here and so the conversion is not taking place. It will be a very long debug before I can figure out what the exact problem is. Are there some inbuilt giac configurations that help with debugging (e.g. some form of logging when performing grobner basis computation)?

Jose

parisse
Messages : 5731
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: Sample use of gbasis in C++

Message par parisse » dim. févr. 04, 2018 3:56 pm

I meant in Xcas.
You can enable debugging by setting giac::debug_infolevel to 2 (a little verbose) or 3 (much more).

jocaps
Messages : 118
Inscription : lun. avr. 17, 2017 4:32 pm

Re:

Message par jocaps » lun. avr. 30, 2018 1:08 pm

Hi Bernard,

I am not getting any wiser with debugging. But maybe there is a clue. I have noticed that in visual studio the preprocessor definition IMMEDIATE_VECTOR is used, is this relevant?

Jose

PS: Would it be possible for the admin of this forum to allow writing inline code. Sometimes writing a code in displaymode is not desirable. This can be activated in phpbb: https://www.phpbb.com/community/viewtop ... &t=2113419

parisse
Messages : 5731
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: Re:

Message par parisse » lun. avr. 30, 2018 1:55 pm

jocaps a écrit :Hi Bernard,

I am not getting any wiser with debugging. But maybe there is a clue. I have noticed that in visual studio the preprocessor definition IMMEDIATE_VECTOR is used, is this relevant?
I don't think so. IMMEDIATE_VECTOR is a trick to avoid calling a malloc for small vectors. It was especially useful on calculators.
PS: Would it be possible for the admin of this forum to allow writing inline code. Sometimes writing a code in displaymode is not desirable. This can be activated in phpbb: https://www.phpbb.com/community/viewtop ... &t=2113419
The forum will move to a new hardware/latest version of phpbb in a few days (or maybe weeks), it's not worth making changes here now...

jocaps
Messages : 118
Inscription : lun. avr. 17, 2017 4:32 pm

Message par jocaps » lun. avr. 30, 2018 8:21 pm

parisse a écrit : The forum will move to a new hardware/latest version of phpbb in a few days (or maybe weeks), it's not worth making changes here now...
Ah ok. Wonderful, thanks for informing.

I tried to continue with my debugging. I had to remove GIAC_VECTOR in my preprocessor because I think I could get to the cause of the problem (and the problem with GIAC_VECTOR seems to be deeper when using gbasis in Visual Studio). Without this predef I narrowed the problem to Add_gen in gausspol.cc. It seems that Add_gen is (line ca. 210 in gausspol.cc) is comparing iterators from different containers (I think the constant iterators a,b,a_end,b_end are not for the container new_coord and that lines compare it with new_coord.begin()). Visual Studio vector does not like this (see for instance https://stackoverflow.com/questions/842 ... compatible ). Am I correct in assuming that you are comparing iterators from different containers? Is there a simpler way to implement Add_gen here? I can try to modify and see if I can resolve this and I will report as soon as I find a solution.

Jose

parisse
Messages : 5731
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re:

Message par parisse » mar. mai 01, 2018 5:27 am

I don't understand why there would be a problem, a, a_end etc. are std::vector< monomial<gen> >::const_iterator, new_coord is a std::vector< monomial<gen> >, therefore comparing a to new_coord.begin() makes sense.

jocaps
Messages : 118
Inscription : lun. avr. 17, 2017 4:32 pm

Re:

Message par jocaps » mar. mai 01, 2018 8:01 am

parisse a écrit :I don't understand why there would be a problem, a, a_end etc. are std::vector< monomial<gen> >::const_iterator, new_coord is a std::vector< monomial<gen> >, therefore comparing a to new_coord.begin() makes sense.
I think this is compiler/linker dependent (unfortunately this can be noticed as late as when linking). Comparing them as pointers would be fine but I think some compilers do complain if they are iterators of different containers (regardless of the type they point at). Maybe there is a point I am missing here.

May I ask the purpose of the first -if- statement in the beggining of Add_gen? What are the instances for which the if statement is true? I am no expert, so I am not sure how other compilers/linkers (e.g. gcc) behave when you inquire new_coord.begin()==a if a is not the iterator of the container new_coord (I would assume this to be false). If what I assume is correct (i.e. new_coord.begin()==a iff a is an iterator of new_coord and a is pointing at new_coord.begin()) then this would never happen if new_coord is a different container of elements than the one iterated/pointed by a. The only caller of Add_gen is addpoly from gen.cc and Add_gen itself (there is a recursive call to Add_gen) both seem to use a new container for the new_coord parameter (not the one iterated by a or b). Am I missing something here? Is it false to assume that the first -if- statement of Add_gen never happens?

Sorry for such a complex question. Probably I am misunderstanding something.

Jose

parisse
Messages : 5731
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re:

Message par parisse » mar. mai 01, 2018 9:38 am

Looking at the callers of Add_gen, it seems that the test will never be true, but I don't want to change the code. There is probably a define for visualc++ to avoid this additional (stupid) check.

jocaps
Messages : 118
Inscription : lun. avr. 17, 2017 4:32 pm

Re:

Message par jocaps » mar. mai 01, 2018 10:54 am

parisse a écrit :Looking at the callers of Add_gen, it seems that the test will never be true, but I don't want to change the code. There is probably a define for visualc++ to avoid this additional (stupid) check.
No problem, I understand that this is very sensitive so we do not want to break the code. I was able to get pass this problem for now by setting the following predef in visual studio :

Code : Tout sélectionner

_HAS_ITERATOR_DEBUGGING=0
_SECURE_SCL=0
Debugging is easier now, but problems still remains (unsolvable gbasis as compared to mingwc builds). I will report as soon as I have something new. Thanks for the support.

Jose

Répondre