#include <stdio.h>
#include <math.h>
float sqrt_(float x) {
static const float T[7] = { 0.545253866f, 0.594603558f, 0.648419777f, 0.707106781f, 0.771105413f, 0.840896415f, 0.917004043f };
static const float C0[2][8] = {
{2.075945786f, 1.987932496f, 1.903650680f, 1.822942136f, 1.745655370f, 1.671645309f, 1.600773032f, 1.532905508f },
{1.467915343f, 1.405680548f, 1.346084305f, 1.289014746f, 1.234364750f, 1.182031734f, 1.131917466f, 1.083927879f },
};
static const float C1[2][8] = {
{-1.324459632f, -1.163042545f, -1.021297991f, -0.896828402f, -0.787528409f, -0.691549235f, -0.607267419f, -0.533257358f },
{-0.936534387f, -0.822395271f, -0.722166735f, -0.634153445f, -0.556866678f, -0.488999154f, -0.429402910f, -0.377069894f },
};
if (x < 0.f) {
return NAN;
}
int e;
float f = frexpf(x, &e);
int e_div_2 = e >> 1;
int e_mod_2 = e & 1;
int n = 3;
int m = 4;
for (int i = 0; i < 3; i++) {
if (f < T[n]) {
m >>= 1;
n -= m;
}
else {
m = (m + 1 >> 1);
n += m;
}
}
x = C0[e_mod_2][n] + C1[e_mod_2][n] * f;
f = ldexpf(f, e_mod_2);
x += .5f * x * (1.f - f * x * x);
return ldexpf(f * x, e_div_2);
}
int main() {
float TABLE1[] = { 0.958202441f, 0.99999994f, 1.f };
float TABLE2[] = { .5f, 0.545253866f, 0.594603558f, 0.648419777f, 0.707106781f, 0.771105413f, 0.840896415f, 0.917004043f };
float TABLE3[] = { 0.125f, 0.25f, 0.5f, 1.f, 2.f, 4.f, 8.f, 16.f, 32.f };
for (float x3 : TABLE3) {
for (float x2 : TABLE2) {
for (float x1 : TABLE1) {
float x = x1 * x2 * x3;
float y1 = sqrt_(x);
float y0 = sqrt(x);
printf("%e\n", (double)y1 / (double)y0 - 1.);
}
}
}
return 0;
}
I2luY2x1ZGUgPHN0ZGlvLmg+CiNpbmNsdWRlIDxtYXRoLmg+CgpmbG9hdCBzcXJ0XyhmbG9hdCB4KSB7CglzdGF0aWMgY29uc3QgZmxvYXQgVFs3XSA9IHsgMC41NDUyNTM4NjZmLCAwLjU5NDYwMzU1OGYsIDAuNjQ4NDE5Nzc3ZiwgMC43MDcxMDY3ODFmLCAwLjc3MTEwNTQxM2YsIDAuODQwODk2NDE1ZiwgMC45MTcwMDQwNDNmIH07CglzdGF0aWMgY29uc3QgZmxvYXQgQzBbMl1bOF0gPSB7CgkJezIuMDc1OTQ1Nzg2ZiwgMS45ODc5MzI0OTZmLCAxLjkwMzY1MDY4MGYsIDEuODIyOTQyMTM2ZiwgMS43NDU2NTUzNzBmLCAxLjY3MTY0NTMwOWYsIDEuNjAwNzczMDMyZiwgMS41MzI5MDU1MDhmIH0sCgkJezEuNDY3OTE1MzQzZiwgMS40MDU2ODA1NDhmLCAxLjM0NjA4NDMwNWYsIDEuMjg5MDE0NzQ2ZiwgMS4yMzQzNjQ3NTBmLCAxLjE4MjAzMTczNGYsIDEuMTMxOTE3NDY2ZiwgMS4wODM5Mjc4NzlmIH0sCgl9OwoJc3RhdGljIGNvbnN0IGZsb2F0IEMxWzJdWzhdID0gewoJCXstMS4zMjQ0NTk2MzJmLCAtMS4xNjMwNDI1NDVmLCAtMS4wMjEyOTc5OTFmLCAtMC44OTY4Mjg0MDJmLCAtMC43ODc1Mjg0MDlmLCAtMC42OTE1NDkyMzVmLCAtMC42MDcyNjc0MTlmLCAtMC41MzMyNTczNThmIH0sCgkJey0wLjkzNjUzNDM4N2YsIC0wLjgyMjM5NTI3MWYsIC0wLjcyMjE2NjczNWYsIC0wLjYzNDE1MzQ0NWYsIC0wLjU1Njg2NjY3OGYsIC0wLjQ4ODk5OTE1NGYsIC0wLjQyOTQwMjkxMGYsIC0wLjM3NzA2OTg5NGYgfSwKCX07CgoJaWYgKHggPCAwLmYpIHsKCQlyZXR1cm4gTkFOOwoJfQoKCWludCBlOwoJZmxvYXQgZiA9IGZyZXhwZih4LCAmZSk7CglpbnQgZV9kaXZfMiA9IGUgPj4gMTsKCWludCBlX21vZF8yID0gZSAmIDE7CglpbnQgbiA9IDM7CglpbnQgbSA9IDQ7Cglmb3IgKGludCBpID0gMDsgaSA8IDM7IGkrKykgewoJCWlmIChmIDwgVFtuXSkgewoJCQltID4+PSAxOwoJCQluIC09IG07CgkJfQoJCWVsc2UgewoJCQltID0gKG0gKyAxID4+IDEpOwoJCQluICs9IG07CgkJfQoJfQoJeCA9IEMwW2VfbW9kXzJdW25dICsgQzFbZV9tb2RfMl1bbl0gKiBmOwoJZiA9IGxkZXhwZihmLCBlX21vZF8yKTsKCXggKz0gLjVmICogeCAqICgxLmYgLSBmICogeCAqIHgpOwoJcmV0dXJuIGxkZXhwZihmICogeCwgZV9kaXZfMik7Cn0KCmludCBtYWluKCkgewoJZmxvYXQgVEFCTEUxW10gPSB7IDAuOTU4MjAyNDQxZiwgMC45OTk5OTk5NGYsIDEuZiB9OwoJZmxvYXQgVEFCTEUyW10gPSB7IC41ZiwgMC41NDUyNTM4NjZmLCAwLjU5NDYwMzU1OGYsIDAuNjQ4NDE5Nzc3ZiwgMC43MDcxMDY3ODFmLCAwLjc3MTEwNTQxM2YsIDAuODQwODk2NDE1ZiwgMC45MTcwMDQwNDNmIH07CglmbG9hdCBUQUJMRTNbXSA9IHsgMC4xMjVmLCAwLjI1ZiwgMC41ZiwgMS5mLCAyLmYsIDQuZiwgOC5mLCAxNi5mLCAzMi5mIH07Cglmb3IgKGZsb2F0IHgzIDogVEFCTEUzKSB7CgkJZm9yIChmbG9hdCB4MiA6IFRBQkxFMikgewoJCQlmb3IgKGZsb2F0IHgxIDogVEFCTEUxKSB7CgkJCQlmbG9hdCB4ID0geDEgKiB4MiAqIHgzOwoJCQkJZmxvYXQgeTEgPSBzcXJ0Xyh4KTsKCQkJCWZsb2F0IHkwID0gc3FydCh4KTsKCQkJCXByaW50ZigiJWVcbiIsIChkb3VibGUpeTEgLyAoZG91YmxlKXkwIC0gMS4pOwoJCQl9CgkJfQoJfQoJcmV0dXJuIDA7Cn0K