Re: Arbitrary precision modulo operation

/*
** Hopefully, I am not annoying anyone with all of this.
** This is a brief demonstration of how Newton's division algorithm
works:
*/
#include <stdio.h>
double divide(double x, double y)
{
unsigned __int64 est;
double y1 = 1.0 / y;
double y2;
puts("\nForming reciprocal:");
est = y1 * 10000;
y1 = est / 10000.0; /* make an estimate good to 4 decimal
places... */
printf("x=%.7f, y=%.7f, y1=%.7f\n", x, y, y1);
y1 *= (2 - y * y1); /* should be 8 places */
printf("x=%.14f, y=%.14f, y1=%.14f\n", x, y, y1);
y1 *= (2 - y * y1); /* In theory ~16 digits */
printf("x=%.20f, y=%.20f, y1=%.20f\n", x, y, y1);
y1 *= (2 - y * y1); /* One more is really needed here... */
printf("x=%.20f, y=%.20f, y1=%.20f\n", x, y, y1);
puts("\ndividing:");
printf("result of division is:%f\n", x * y1);
return 0;
}

int main(void)
{
divide(1.0, 2.0);
divide(1.0, 17.0);
divide(3.0, 19.0);
divide(19.0, 3.0);
return 0;
}