Advertisement
aaaaaa123456789

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. }
Advertisement
RAW Paste Data Copied
Advertisement