The maximum radius you can have, using the premises you have described is 100000
, which leads to a squared modulus of 10 000 000 000
, which requires (as all numbers are positive) 35
bit integer (of unsigned
square radius) to be represented.
Based on these premises, and the fact that you have no easy way to get 64bit integers, and having quite low extra bits, we can scale the results four bits in norm (two bits in the source coordinates) to achieve full capacity to handle upto 100000
coordinates in a 32bit unsigned integer.
In my first edition of this answer, I assumed using only one shift in the coordinates was enough to handle the full set of values (two bits in the calculated norm), and losing 1 bit of precision was considered, but I was wrong and one extra bit was needed. It is needed to shift the results at least three bits to hold the full set of inputs, so I decided to scale the coordinates two bits, and so the results would be scaled by four. As I decided also to always scale, and return the fraction of a square unit as an integer ranging from 0
to 15
(in sixteenths of a square unit). So you will achieve exact results by comparing first the integer parts of the two points and use the fractional parts, in case the integer parts match. This makes the computation and the meaning of the results returned back more coherent than earlier, and gives you complete exactitude with integer coordinates.
You requested a working implementation, so I have posted one for you below:
#include <stdio.h>
#include <stdint.h>
/* calculate the square of a divided by four number and
* accumulate the fraction (in sixteenths of a square unit)
* into the reference pointed by frac_p. */
uint32_t
square_of_div16(uint32_t x, int *frac_p)
{
/* we use (IP + FP)^2 = IP^2 + 2*IP*FP + FP^2 */
uint32_t int_part = x >> 2; /* divide by four */
uint32_t frac_part = x & 0x3; /* mod 4 */
uint32_t int_result = int_part * int_part; /* square of IP */
int frac_result = frac_part * frac_part; /* square of FP */
uint32_t mixed_prod = int_part * frac_part; /* IP*FP */
int_result += mixed_prod >> 1;
frac_result += (mixed_prod & 1) << 3;
if (frac_result >= 0x10) { /* carry process */
int_result += frac_result >> 4;
frac_result &= 0x0f;
}
if (frac_p) *frac_p += frac_result; /* accumulate */
return int_result;
}
/* this calculates the squared norm scaled to one sixteenth
* of the original coordinates (scaled by one fourth).
* The ref_fraction pointer is a reference of a variable to
* accumulate the fraction sixteenths of a square unit. If
* you are not interested in the fraction value, you can just
* pass NULL as parameter. */
uint32_t
norm_scaled(uint32_t x, uint32_t y, int *ref_fraction)
{
int fraction = 0;
uint32_t result = 0;
result += square_of_div16(x, &fraction);
result += square_of_div16(y, &fraction);
if (ref_fraction)
*ref_fraction += fraction; /* the excess */
return result;
}
/* TEST MAIN PROGRAM. Just input pairs of coordinates in the
* same line (separated by spaces) and calculate the squared
* norm of the vector, scaled by 1/16 (accumulating the
* fraction of the value in 1/16s of a square unit in the
* location referenced. This is done using double floating
* point numbers and uint32_t integers. */
int main()
{
char line[256];
while (fgets(line, sizeof line, stdin) != NULL) {
int x = 0, y = 0, fraction = 0;
sscanf(line, "%u%u", &x, &y);
uint32_t norm_16th = norm_scaled(x, y, &fraction);
printf("Trying (%u, %u) => %u (fraction = %d/16)\n",
x, y, norm_16th, fraction);
double norm_sq_16th
= (double) x/4.0 * (double)x/4.0
+ (double) y/4.0 * (double)y/4.0;
printf("squared norm scaled: %.8f\n", norm_sq_16th);
}
printf("Program finished\n");
}
The function square_of_div16
calculates a scaled modulus divided by 16 of a number, so we can reuse it to calculate the squares of x
and y
coordinates. The function takes a pointer frac_p to an integer variable to store the fraction part (in sixteenths of a square unit)
The function norm_scaled
then calculates the norm, by using the square_of_div16
function and adding both results. The fractional part is accumulated for both calls and the result is accumulated to the referred variable by pointer ref_fraction
. A carry processing is done in here, to give correct results.
Finally a main()
routine is in charge of querying the user to input pairs of coordinates and calculate the scaled norm of the resulting vector by calling the function and using the squares pithagorean formula applied to double
values. The results should be the same in all cases.
x^2 + y^2 < r^2
is true, you can test if eitherx
ory
is greater thanr
fairly quickly, and whether bothx
andy
are less thanr/sqrt(2)
. Depending on how yourx
,y
andr
are distributed, this could be a useful fast test.x^2 + y^2 <= r^2
, because the addition could overflow even whenx<r
andy<r
. Instead computex^2 <= r^2 - y^2
. As long asy <= r
the subtraction is guaranteed to work.