Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- /* This code is hereby released to the public domain.
- * ~aaaaaa123456789, 2020-01-01 */
- #include <stdio.h>
- #include <stdlib.h>
- #include <stdint.h>
- #include <math.h>
- #include <time.h>
- #define FRACTIONAL_BITS 20
- int32_t integer_inverse (unsigned value) {
- if (!value) return INT32_MAX;
- return ((2 << FRACTIONAL_BITS) / value + 1) >> 1;
- }
- int32_t multiply (int32_t first, int32_t second) {
- int64_t result = ((int64_t) first * (int64_t) second + ((int64_t) 1 << (FRACTIONAL_BITS - 1))) >> FRACTIONAL_BITS;
- if (result > INT32_MAX) return INT32_MAX;
- if (result < INT32_MIN) return INT32_MIN;
- return result;
- }
- int32_t approximate_square_root (int32_t operand, unsigned iterations) {
- unsigned iteration;
- int32_t result = (operand + (1 << FRACTIONAL_BITS)) >> 1;
- int32_t factor = ((1 << FRACTIONAL_BITS) - operand) >> 1;
- int32_t power = factor;
- int32_t coefficient = 1 << FRACTIONAL_BITS;
- for (iteration = 2; iteration <= iterations; iteration ++) {
- coefficient = multiply(coefficient, (2 * iteration - 3) << FRACTIONAL_BITS);
- coefficient = multiply(coefficient, integer_inverse(iteration));
- power = multiply(power, factor);
- result -= multiply(power, coefficient);
- }
- return result;
- }
- int32_t true_square_root (int32_t operand) {
- #if FRACTIONAL_BITS & 1
- #error This function requires an even number of fractional bits!
- #endif
- if (operand < 0) return 0;
- return ldexp(sqrt(operand), FRACTIONAL_BITS / 2) + 0.5;
- }
- unsigned calculate_error (int32_t approximate, int32_t real) {
- double error = (double) approximate / (double) real - 1.0;
- error = fabs(error);
- return error * 10000 + 0.5;
- }
- int main (void) {
- unsigned tests = 100, limit = 6, current;
- srand(clock() ^ time(NULL));
- fputs(" value sqrt 1 term ", stdout);
- for (current = 2; current <= limit; current ++) printf(" %2u terms ", current);
- fputs("\n-------- -------- ---------------", stdout);
- for (current = 2; current <= limit; current ++) fputs(" ---------------", stdout);
- putchar('\n');
- while (tests --) {
- int32_t value = ((rand() ^ (rand() << 12)) >> 4) & ((1 << (FRACTIONAL_BITS + 1)) - 1); // random between 0.0 and 2.0
- int32_t true_root = true_square_root(value);
- printf("%.6f %.6f", ldexp(value, -FRACTIONAL_BITS), ldexp(true_root, -FRACTIONAL_BITS));
- for (current = 1; current <= limit; current ++) {
- int32_t root = approximate_square_root(value, current);
- printf(" %.6f", ldexp(root, -FRACTIONAL_BITS));
- unsigned error = calculate_error(root, true_root);
- if (error >= 10000)
- printf(" %5u%%", (error + 50) / 100);
- else
- printf(" %2u.%02u%%", error / 100, error % 100);
- }
- putchar('\n');
- }
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement