Re: Incorporating Pari into a C program

Le sam. 17 mai 2025 à 09:55, Gordon Royle <[email protected]> a
écrit :

> Hi Pari Users
>
>
>
> I do a lot of things with graph polynomials and their roots and for that I
> use GP a lot, but only in a naïve way – essentially just calling
> polrootsreal on billions of polynomials (using a Perl script to actually
> format the relevant GP commands).
>
>
>
> Now I need to do something different and although I have tried to read the
> manual, I am lost.
>
>
>
> What I need is a C function that will operate as a “plugin” to an existing
> C program (not written by me).
>
>
>
> The C program traverses a huge search tree, and at every node it calls the
> plugin function (whose name is given to the outer program at compile time)
> and whose arguments are fixed. So I write a routine, something like
>
>
>
> int decision(int *a, int n)
>
>
>
> which is called by main program at every node (and there will be billions
> of nodes).
>
>
>
> The response given by my plugin determines whether that branch of the
> search tree should be extended or pruned.
>
>
>
> What does my function do?
>
>
>
> It uses the input array to construct an n x n real symmetric matrix, and
> tests whether its largest eigenvalue is greater than 2, and then returns 1
> or 0 accordingly.
>
>
>
> The matrices are of moderate-ish size (20x20 – 30x30) and the entries are
> small rational numbers.
>
>
>
> At the moment, (just in the testing phase) I create the matrix using
> floating point arithmetic and use a simple eigenvalue program that had made
> available on some academic’s home page.
>
>
>
> Now I want to do it rigorously – that is, I need to make sure that when my
> routine says “the largest eigenvalue is greater than 2” it is *definitely
> *correct.  I am willing to have a routine that makes an occasional
> mistake in the other direction – it can incorrectly report that the
> eigenvalue is less than 2, but only very occasionally.
>
>
>
> (I also want to do it faster (of course) and the obvious place to gain
> speed is to only calculate the largest eigenvalue.)
>
>
>
> There are various ways that I could try using Pari, but I cannot
> understand the documentation about using Pari from a C program – in
> particular, the instructions about initializing the Pari stack and
> allocating and freeing memory are beyond me.
>
>
>
> Can anyone provide a simple example of how I could use a Pari routine from
> within my plugin function?   And how I should manage the initialization
> step… the plugin mechanism does allow me define a second function that is
> called just once at the beginning of the search, and presumably I would put
> the pari-initialize stuff in there?
>

      You'll have to put

#include <pari/pari.h>

    at the beginning. In your decision(int a*, int n) you'll have to do a
variant of

    ***Since I do not know how the elements of the matrix are encoded, you
said
  "small rational numbers" but provide a pointer a an array of integers,
I'll make
  the (wrong) assumption that a is actually int ***a, where a[i-1][j-1][0]
is the
  numerator of the (i,j)-entry and a[i-1][j-1][1] is the denominator.
You'll have to
  translate that into the linear versione (maybe
a[2*n*n*(i-1)+2*n*(j-1)+k], where
  k is 0 or 1 above).***

    GEN M, col;
    long i, j;

    pari_init(2*1024*1024,500*1000);
    paristack_setsize(2*1024*1024,1024uL*1024uL*1024uL);
   M = cgetg(n+1, t_MAT);
   for (j = 1; j < n; j++) {
        GEN col = cgetg(n+1, t_COL);
        for (i = 1; i < n; i++)
            gel(col, i) = gdiv(ltoi(a[i-1][j-1][0]),ltoi(a[i-1][j-1][1]));
        gel(M, j) = col;
   }
   /* Compute with M, store result in r */
   pari_close();
   return r;



> Unfortunately, my C is sufficiently rusty that I am not sure how the
> plugin can access things created by the “initialize plugin” function.
>
>
>
> I have also considered that perhaps Pari is overkill for what I need, and
> that it would be better to use, say GMP, to create my rational matrix, then
> simply run a few iterations of the power method starting with some vector
> v0, so v1=Bv0, v2=Bv1, v3=Bv2, … and then check if the Rayleigh quotient
> v3^TBv^3/v3^Tv is greater than 2.
>

     This tells you whether the largest *modulus* of (complex) eigenvalues
is greater
  than 2. (Except that you might miss it!) Do you need the *real*
eigenvalue 2?

     You'll probably be better off multiplying the matrix to get integer
entries (and obvisously
  adjust 2).

         Hope this helps,

                Loïc