// check modular polynomial solutions for N = 3
//
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define p 5     // the prime modulus
#define M 9     // coefficients a[0]...a[M]
#define N 3     // variables x[0]...x[N-1]

void init( int a[], int n)
{
  for( int i = 0; i < n; ++i)
  {
    a[i] = rand() % p; // general case
    // a[i] = 1 + rand() % (p-1); // all non-zero coefficients
  }
}

void print( const char *msg, int a[], int n)
{
  printf( "%s = (", msg);
  for( int i = 0; i < n; ++i) printf( "%d%s", a[i], i < n-1 ? "," : "");
  printf( ")\n");
}

int f( int a[], int x[])
{
  return (a[0] + a[1]*x[0] + a[2]*x[1] + a[3]*x[2]
    + a[4]*x[0]*x[0] + a[5]*x[1]*x[1] + a[6]*x[2]*x[2]
    + a[7]*x[0]*x[1] + a[8]*x[0]*x[2] + a[9]*x[1]*x[2]) % p;
}

int main( void)
{
  int a[M+1], x[N], s;

  srand(time(0));

  init( a, M+1);

  print( "a", a, M+1);

  s = 0;

  for( x[0] = 0; x[0] < p; ++x[0])
   for( x[1] = 0; x[1] < p; ++x[1])
    for( x[2] = 0; x[2] < p; ++x[2])
    {
        if( f(a,x) == 0)
        {
          ++s; print( "x", x, N);
        }
    }

  printf( "p = %d, %d solutions, %d%%%d = %d\n", p, s, s, p, s%p);

  return 0;
}