/* -*-ePiX-*- */
#include "epix.h"
using namespace ePiX;

const double M(10);
const int N(40);
double f(double p, double q)
{
  return sqrt((M + p*p)/(M + q*q));
}

int main()
{
  picture(P(-1.5, -1), P(1.5, 1.5), "4.2 x 3.5in");

  begin();

  for (int i=1; i<N; ++i)
    for (int j = ceil(-1.5*i); j <= floor(1.5*i); ++j)
      if (gcd(i, j) == 1 && fabs(fabs(j*1.0/i) - f(j, i)) > 0.025)
	ddot(P(j*1.0/i, f(j, i)));

  for (int i=-1.5*N; i<=1.5*N; ++i)
    {
      double t(i*1.0/(N+1));
      ddot(P(t, fabs(t)));
    }

  bbold();
  line(P(ymin(), ymin()), P(xmax(), xmax()));

  bold();
  line(P(xmin(), 0), P(xmax(), 0));
  line(P(0, ymin()), P(0, ymax()));

  tikz_format();
  end();
}