# How to mark rational points on a sphere

I found this picture on mathoverflow, which I find very intriguing and so I like to know how to draw such an image with a simple computer program.

To calculate the rational point, I can draw a line from P_0(0,0,1) and P_1(u,v,0) and calculate the intersection with the sphere as follows:

$$x=\frac{2u}{u^2+v^2+1};y=\frac{2v}{u^2+v^2+1};z=\frac{u^2+v^2-1}{u^2+v^2+1}$$

So the intersection coordinates are rational numbers where

$$a=2u;b=2v;c=u^2+v^2-1;d=u^2+v^2+1$$

with

$$a^2+b^2+c^2=d^2$$

If I understand correctly the guy who posted the picture he would mark the points depending on the value of d. So a value of d below a certain threshold would override the pixel color to white.

Now inside a program we would loop through all pixels along the x and y axis using the following algorithm:

for(int y = 0; y < height; ++y)
{
for(int x = 0; x < width; ++x)
{
// does pixel ray intersect with sphere?
// using orthgraphic projection ( all rays are parallel ) of
// sphere with radius of image height
if(intersect(sphere,x,y))
{
// calculate the distance of the intersection of
// pixel ray and sphere
double z = calc_distance(sphere,x,y);

// the further away the darker the color
rgb_color color = make_color(z);

// calculate d
// the gcd of a,b,c,d will be used to make d as small as possible
int d = calculate_denominator(x,y)

if(d < 100)
{
color = white;
}

set_pixel(img,x,y,color);
}
}
}


Could someone help in correcting my algorithm?

-

I'd iterate over $u$ and $v$ instead, since the point which actually has easy rational coordinates is not exactly at a pixel position.

double d;
for (int u = 0; u <= 100; ++u) {
for (int v = 0; (d = u*u+v*v+1) <= 100; ++v) {
markBlack(2*u/d, 2*v/d); // arguments are double in [0,1]
}
}


or, if you want to stay in integer arithmetic:

int d;
for (int u = 0; u <= 100; ++u) {
for (int v = 0; (d = u*u+v*v+1) <= 100; ++v) {
int dh = d/2; // rounded down
int x = (size*2*u + dh)/d;
int y = (size*2*v + dh)/d;
set_pixel(img,x,y,black);
}
}


In both cases, I'm ommitting the third coordinate since we are projecting onto a coordinate plane.

The picture still doesn't look like the one you indicated, though. That's likely because you'd have to consider fractional $u,v$ as well. You might be better off using Pythagorean quadruples instead. Then you can use nested loops there.

int size = 2047, limit = 100;
int dm, dn, dp, d;
for (int m = 0; (dm = m*m) <= limit; ++m)
for (int n = 0; (dn = dm + n*n) <= limit; ++n)
for (int p = 0; (dp = dn + p*p) <= limit; ++p)
for (int q = 0; (d = (dp + q*q) <= limit; ++q) {
if (d == 0) continue;
int a = abs(m*m + n*n - p*p - q*q);
int b =     2*(m*q + n*p) ;
int c = abs(2*(n*q - m*p));
int dh = d/2;
int x = (size*a + dh)/d;
int y = (size*b + dh)/d;
int z = (size*c + dh)/d;
set_pixel(img,x,y,black);
set_pixel(img,x,z,black);
set_pixel(img,y,x,black);
set_pixel(img,y,z,black);
set_pixel(img,z,x,black);
set_pixel(img,z,y,black);
}


However, instead of marking pixels either black or white, you might also consider marking them with a suitable gray level, with the gray value depending on the number of rational points within a given pixel. If you do that, then you want to make sure that you mark every rational point exactly once. I did so in the code used for the animation below, building on this paper.

-
Thanks for your solution. Once I understand it I'll mark your answer. Could you post the whole code so I can test it out? I tried the the second loop but don't know what size means. – chhenning Aug 31 '14 at 16:24
@chhenning: size is the size of the image, in pixels, Actually one less, since coordinates range from 0 to size inclusive. I never ran this C code, but I ran its Python analog in sage, and posted that here. – MvG Aug 31 '14 at 19:23