# Square root approximation code

Jan 2nd, 2020
154
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. /* This code is hereby released to the public domain.
2.  * ~aaaaaa123456789, 2020-01-01                       */
3.
4. #include <stdio.h>
5. #include <stdlib.h>
6. #include <stdint.h>
7. #include <math.h>
8. #include <time.h>
9.
10. #define FRACTIONAL_BITS 20
11.
12. int32_t integer_inverse (unsigned value) {
13.   if (!value) return INT32_MAX;
14.   return ((2 << FRACTIONAL_BITS) / value + 1) >> 1;
15. }
16.
17. int32_t multiply (int32_t first, int32_t second) {
18.   int64_t result = ((int64_t) first * (int64_t) second + ((int64_t) 1 << (FRACTIONAL_BITS - 1))) >> FRACTIONAL_BITS;
19.   if (result > INT32_MAX) return INT32_MAX;
20.   if (result < INT32_MIN) return INT32_MIN;
21.   return result;
22. }
23.
24. int32_t approximate_square_root (int32_t operand, unsigned iterations) {
25.   unsigned iteration;
26.   int32_t result = (operand + (1 << FRACTIONAL_BITS)) >> 1;
27.   int32_t factor = ((1 << FRACTIONAL_BITS) - operand) >> 1;
28.   int32_t power = factor;
29.   int32_t coefficient = 1 << FRACTIONAL_BITS;
30.   for (iteration = 2; iteration <= iterations; iteration ++) {
31.     coefficient = multiply(coefficient, (2 * iteration - 3) << FRACTIONAL_BITS);
32.     coefficient = multiply(coefficient, integer_inverse(iteration));
33.     power = multiply(power, factor);
34.     result -= multiply(power, coefficient);
35.   }
36.   return result;
37. }
38.
39. int32_t true_square_root (int32_t operand) {
40.   #if FRACTIONAL_BITS & 1
41.     #error This function requires an even number of fractional bits!
42.   #endif
43.   if (operand < 0) return 0;
44.   return ldexp(sqrt(operand), FRACTIONAL_BITS / 2) + 0.5;
45. }
46.
47. unsigned calculate_error (int32_t approximate, int32_t real) {
48.   double error = (double) approximate / (double) real - 1.0;
49.   error = fabs(error);
50.   return error * 10000 + 0.5;
51. }
52.
53. int main (void) {
54.   unsigned tests = 100, limit = 6, current;
55.   srand(clock() ^ time(NULL));
56.   fputs("  value    sqrt        1 term    ", stdout);
57.   for (current = 2; current <= limit; current ++) printf("     %2u terms   ", current);
58.   fputs("\n-------- -------- ---------------", stdout);
59.   for (current = 2; current <= limit; current ++) fputs(" ---------------", stdout);
60.   putchar('\n');
61.   while (tests --) {
62.     int32_t value = ((rand() ^ (rand() << 12)) >> 4) & ((1 << (FRACTIONAL_BITS + 1)) - 1); // random between 0.0 and 2.0
63.     int32_t true_root = true_square_root(value);
64.     printf("%.6f %.6f", ldexp(value, -FRACTIONAL_BITS), ldexp(true_root, -FRACTIONAL_BITS));
65.     for (current = 1; current <= limit; current ++) {
66.       int32_t root = approximate_square_root(value, current);
67.       printf(" %.6f", ldexp(root, -FRACTIONAL_BITS));
68.       unsigned error = calculate_error(root, true_root);
69.       if (error >= 10000)
70.         printf(" %5u%%", (error + 50) / 100);
71.       else
72.         printf(" %2u.%02u%%", error / 100, error % 100);
73.     }
74.     putchar('\n');
75.   }
76.   return 0;
77. }