/*  Program in C  */

/* begin source file */
#include <math.h>
#include <assert.h>
#include <stdlib.h>
#include <malloc.h>
/*
	Randy Thomson
	Math 365
	2-29-98
*/

int gcd(int x, int y)
{
/*
	Simple Euclidean style gcd algorithm.
*/

  int g;
  if (x<0)
    x=-x;
  if (y<0)
    y=-y;
  assert((x+y)!=0);
  g=y;
  while (x>0) {
    g=x;
    x=y%x;
    y=g;
  };
  return g;
};



int jacobi(int a, int b)
{
/*	
	This program recursively calculates the Jacobi symbol (a,b).
	It reduces the pair recursively according to a series of
	Jacobi identities, similar to the Euclidean gcd routine.
*/

  int g;
  assert((b%2)==1);
  if (a>=b) a%=b;
  if (a==0) return 0;
  if (a==1) return 1;

  if (a<0)
    {
      if (((b-1)/2 % 2)==0)
	{
	  return jacobi(-a,b);
	}
      else
	{
	  return -jacobi(-a,b);
	};
    };
  if ((a%2)==0)
    {
      if(((b*b-1)/8)%2 ==0)
	{
	  return jacobi(a/2,b);
	}
      else
	{
	  return -jacobi(a/2,b);
	};
    };
  g=gcd(a,b);
  assert((a%2)==1);
  
  if (g==a)
    return 0;
  else if (g!=1)
    return jacobi(g,b)*jacobi(a/g,b);
  else if (((a-1)*(b-1)/4) % 2 == 0)
    return jacobi(b,a);
  else
    return -jacobi(b,a);
};


int fill_with_primes(int *array,int lo,int hi)
{
  /*	I use this very slow, but memory efficient function, to
	initialize an array to all the prime values between lo
	and hi.   Memory was a consideration, but not cpu time,
	at least for this step.

	Basically, it uses simple trial division to test for
	primality.  It's a simple, deterministic method.  I
	considered Pollard-rho and Miller-Rabin, but they
	were more trouble than it was worth, and I wanted
	a guaranteed deterministic algorithm.
  */
  int i,j;
  int pindex;
  int sn;
  int isprime;
  pindex=0;
  if (lo==2)
    {
      array[pindex]=2;
      ++pindex;
    };
  if ((lo % 2)==0)
    ++lo;
  for (i=lo;i<=hi;i+=2)
    {
      sn=sqrt(i);
      isprime=1;
      for (j=3;j<=sn;j+=2)
      {
	if ((i%j)==0)
	  {
	    isprime=0;
	    break;
	  };
      };
      if (isprime)
	{
	  array[pindex]=i;
	  ++pindex;
	};
    };
  return pindex;
};

int main(int argc,char **argv)
{
  int lo,hi;
  int residues[50];
  int *primes;
  int n;
  int num_primes,num_residues;
  int i,j;
  int p;
  int num_satisfied;
  int v;
  int sp;
  int step;
  lo=atoi(argv[1]);
  hi=atoi(argv[2]);

  residues[0]=-1;
  /*Save time, and only do the computation with primes between 2 and 50,
    plus -1.*/
  num_residues=fill_with_primes(residues+1,2,50)+1;
  n=(hi-lo+1);
  for (i=lo;i<hi;++i)
    {
	/*Use congruences to only consider the numbers that are 1 mod 24.
	This speeds things up considerably.*/
      if (!(i%24==1))
	{
	  if ((i%24)==0)
	    continue;
	  else
	    i+=25-(i%24);
	};
      p=i;
      num_satisfied=0;
	/*Check the congruence for each of the appropriate v.  Break
	out of the loop as soon as one fails.  This saves a lot of time.
	*/
      for (j=0;j<num_residues;++j)
	{
	  
	  v=(residues[j]+p) % p;
	  if (jacobi(v,p)==1)
	    ++num_satisfied;
	  else
	    break;
	};
      if (num_satisfied==num_residues)
      {
	sp=sqrt(p);
	/*Print only non-square p*/
	if (sp*sp!=p)
	  printf("%d\n",p);
      }
    };
  return 0;


};
*/