| 26 | |
| 27 | template <typename T> |
| 28 | void newton_sqrt() |
| 29 | { |
| 30 | typedef exprtk::symbol_table<T> symbol_table_t; |
| 31 | typedef exprtk::expression<T> expression_t; |
| 32 | typedef exprtk::parser<T> parser_t; |
| 33 | typedef exprtk::function_compositor<T> compositor_t; |
| 34 | typedef typename compositor_t::function function_t; |
| 35 | |
| 36 | T x = T(0); |
| 37 | |
| 38 | symbol_table_t symbol_table; |
| 39 | |
| 40 | symbol_table.add_constants(); |
| 41 | symbol_table.add_variable("x",x); |
| 42 | |
| 43 | compositor_t compositor(symbol_table); |
| 44 | |
| 45 | compositor.add( |
| 46 | function_t("newton_sqrt") |
| 47 | .var("x") |
| 48 | .expression |
| 49 | ( |
| 50 | " switch " |
| 51 | " { " |
| 52 | " case x < 0 : null; " |
| 53 | " case x == 0 : 0; " |
| 54 | " case x == 1 : 1; " |
| 55 | " default: " |
| 56 | " { " |
| 57 | " var z := 100; " |
| 58 | " var sqrt_x := x / 2; " |
| 59 | " repeat " |
| 60 | " if (equal(sqrt_x^2, x)) " |
| 61 | " break[sqrt_x]; " |
| 62 | " else " |
| 63 | " sqrt_x := (1 / 2) * (sqrt_x + (x / sqrt_x)); " |
| 64 | " until ((z -= 1) <= 0); " |
| 65 | " }; " |
| 66 | " } " |
| 67 | )); |
| 68 | |
| 69 | const std::string expression_str = "newton_sqrt(x)"; |
| 70 | |
| 71 | expression_t expression; |
| 72 | expression.register_symbol_table(symbol_table); |
| 73 | |
| 74 | parser_t parser; |
| 75 | parser.compile(expression_str,expression); |
| 76 | |
| 77 | for (std::size_t i = 0; i < 1000; ++i) |
| 78 | { |
| 79 | x = static_cast<T>(i); |
| 80 | |
| 81 | const T result = expression.value(); |
| 82 | const T real = std::sqrt(x); |
| 83 | const T error = std::abs(result - real); |
| 84 | |
| 85 | printf("sqrt(%03d) - Result: %15.13f\tReal: %15.13f\tError: %18.16f\n", |