xref: /btstack/test/maths/sqrt_test.cpp (revision a64cbea79fa5fb9daa58e135a19c6cce10f3e642)
1 #include <math.h>
2 #include <sys/time.h>
3 #include <stdio.h>
4 #include <stdlib.h>
5 
6 #include "CppUTest/TestHarness.h"
7 #include "CppUTest/CommandLineTestRunner.h"
8 
9 static struct timeval init_tv;
10 
11 union {
12     int i;
13     float x;
14 } u;
15 
16 // taken from http://www.codeproject.com/Articles/69941/Best-Square-Root-Method-Algorithm-Function-Precisi
17 // Algorithm: Babylonian Method + some manipulations on IEEE 32 bit floating point representation
18 float sqrt1(const float x){
19     u.x = x;
20     u.i = (1<<29) + (u.i >> 1) - (1<<22);
21 
22     // Two Babylonian Steps (simplified from:)
23     // u.x = 0.5f * (u.x + x/u.x);
24     // u.x = 0.5f * (u.x + x/u.x);
25     u.x =       u.x + x/u.x;
26     u.x = 0.25f*u.x + x/u.x;
27 
28     return u.x;
29 }
30 
31 // Algorithm: Log base 2 approximation and Newton's Method
32 float sqrt3(const float x){
33     u.x = x;
34     u.i = (1<<29) + (u.i >> 1) - (1<<22);
35     return u.x;
36 }
37 
38 float sqrt2(const float n) {
39     /*using n itself as initial approximation => improve */
40     float x = n;
41     float y = 1;
42     float e = 0.001; /* e decides the accuracy level*/
43     while(x - y > e){
44         x = (x + y)/2;
45         y = n/x;
46     }
47     return x;
48 }
49 
50 static uint32_t get_time_ms(void){
51     struct timeval tv;
52     gettimeofday(&tv, NULL);
53     uint32_t time_ms = (uint32_t)((tv.tv_sec  - init_tv.tv_sec) * 1000) + (tv.tv_usec / 1000);
54     return time_ms;
55 }
56 
57 static int values_len = 100000;
58 
59 TEST_GROUP(SqrtTest){
60     void setup(void){
61     }
62 
63     void test_method(float (*my_sqrt)(const float x)){
64         int i, j;
65         float precision = 0;
66         int ta = 0;
67         int te = 0;
68 
69         for (j=0; j<100; j++){
70             for (i=0; i<values_len; i++){
71                 int t1 = get_time_ms();
72                 float expected = sqrt(i);
73                 int t2 = get_time_ms();
74                 float actual = my_sqrt(i);
75                 int t3 = get_time_ms();
76                 te += t2 - t1;
77                 ta += t3 - t2;
78                 precision += fabs(expected - actual);
79             }
80         }
81         printf("Precision: %f, Time: (%d, %d)ms\n", precision/values_len, te, ta);
82     }
83 };
84 
85 TEST(SqrtTest, Sqrt1){
86     printf("\nsqrt1: ");
87     test_method(sqrt1);
88 }
89 
90 TEST(SqrtTest, Sqrt2){
91     printf("\nsqrt2: ");
92     test_method(sqrt2);
93 }
94 
95 TEST(SqrtTest, Sqrt3){
96     printf("\nsqrt3: ");
97     test_method(sqrt3);
98 }
99 
100 int main (int argc, const char * argv[]){
101     return CommandLineTestRunner::RunAllTests(argc, argv);
102 }
103 
104 // TODO: check http://www.embedded.com/electronics-blogs/programmer-s-toolbox/4219659/Integer-Square-Roots