#include <iostream>
#include <Singular/libsingular.h>
#include <omp.h>

int main() {
   //initialize singular
   siInit((char *)"/usr/lib/libSingular.so");

   #pragma omp parallel
   {
      // construct the ring Z/32003[x,y,z]
      // the variable names
      char **n=(char**)omalloc(3*sizeof(char*));
      n[0]=omStrDup("x");
      n[1]=omStrDup("y");
      n[2]=omStrDup("z");

      ring R=rDefault(32003,3,n);
      // make R the default ring: not needed if p_<...> functions are used
      //rChangeCurrRing(R);

      // create the polynomial 1
      poly p1=p_ISet(1,R);

      // create the polynomial 2*x^3*z^2
      poly p2=p_ISet(2,R);
      p_SetExp(p2,1,3,R);
      p_SetExp(p2,3,2,R);
      p_Setm(p2,R);

      // print p1 + p2
      #pragma omp critical
      {
         p_Write(p1,R); printf(" + \n"); p_Write(p2,R); printf("\n");
      }
   }

   return 0;
}

// ***dError: Bug reported: p_Procs is NULL
occurred at ./templates/p_Procs_Set.h,170
// ** Singular will work properly, but much slower
// ** If you chose a coef ring, it may not work at all