Arbitrary precision formula experiment (#103)

* examples: allow custom formula filename in cmake

* [#7] example experiment: formula with arbitrary precision support

* MP example: add swap idiom and other improvements

* add debug information to cpp examples compilation

examples/README.md

@@ -152,4 +152,13 @@ Execute:
 ```
 examples/cpp/custom_formula.sh
 ```
+Then you should see a new file under `examples/output`.
+
+### Creating a simple mandelbrot using a custom formula with MPFR
+In this example we did an experiment using [mpfr library](https://www.mpfr.org/). We replaced the double types in the formula from the previos example with a custom type with operator overloading using a multiple precision support from mpfr.
+Of course this is not changing much the final result as for that we'd need to change the fract4dc engine (sources under /model) to support multiple precision. It's just an experiment about one possible approach to tackle the arbitrary precision support which would also need to change the compiler, among many other parts of the program.
+Execute:
+```
+examples/cpp/mp_mandelbrot.sh
+```
 Then you should see a new file under `examples/output`.

examples/cpp/CMakeLists.txt

@@ -3,6 +3,12 @@ project(simple_mandelbrot)
 
 set(CMAKE_CXX_STANDARD 17)
 
+# add some debug information: elapsed time and compiler information
+set_property(GLOBAL PROPERTY RULE_LAUNCH_COMPILE "${CMAKE_COMMAND} -E time")
+set_property(GLOBAL PROPERTY RULE_LAUNCH_LINK "${CMAKE_COMMAND} -E time")
+set(CMAKE_VERBOSE_MAKEFILE on)
+message("Compiler: ${CMAKE_CXX_COMPILER_ID} ${CMAKE_CXX_COMPILER_VERSION}")
+
 include_directories(${PROJECT_SOURCE_DIR})
 include_directories(${PROJECT_SOURCE_DIR}/model)
 
@@ -42,7 +48,30 @@ set_target_properties(fract_stdlib PROPERTIES SUFFIX ".so")
 
 # FORMULA LIB
 
+get_filename_component(FORMULA_FILE $ENV{formula_source} NAME)
+
 add_library(formula SHARED
-    ${PROJECT_SOURCE_DIR}/formula.c)
+        ${PROJECT_SOURCE_DIR}/${FORMULA_FILE})
+
+if("$ENV{gmp_support}" STREQUAL "1")
+    message("gmp support enabled")
+
+    # https://github.com/symengine/symengine/blob/master/cmake/LibFindMacros.cmake
+    find_library(MPFR_LIBRARIES NAMES mpfr)
+    find_library(GMP_LIBRARIES NAMES gmpxx)
+
+    find_path(MPFR_INCLUDES NAMES mpfr.h)
+    find_path(GMP_INCLUDE_DIR NAMES gmpxx.h)
+
+    include_directories(${GMP_INCLUDE_DIR})
+    target_link_libraries(formula ${GMP_LIBRARIES})
+    include_directories(${MPFR_INCLUDES})
+    target_link_libraries(formula ${MPFR_LIBRARIES})
+
+else()
+    message("gmp support disabled")
+
+endif()
+
 set_target_properties(formula PROPERTIES SUFFIX ".so")
 set_target_properties(formula PROPERTIES PREFIX "")

examples/cpp/Dockerfile

@@ -5,17 +5,23 @@ ENV DEBIAN_FRONTEND=noninteractive
 
 RUN apt-get update && \
 	apt-get install -y build-essential git cmake autoconf libtool pkg-config libpng-dev libjpeg-dev
+RUN apt-get install -y libmpfr-dev libgmp-dev
 
 WORKDIR /src
 
 ARG main_source
 ARG formula_source
+ARG gmp_support
 
 COPY examples/cpp/${main_source} ./main.cpp
 COPY examples/cpp/CMakeLists.txt ./
 COPY fract4d/c/fract_stdlib.cpp fract4d/c/fract_stdlib.h fract4d/c/pf.h  ./
 COPY fract4d/c/model/ ./model
-COPY ${formula_source} ./formula.c
+COPY ${formula_source} ./
+
+ENV formula_source=$formula_source
+ENV gmp_support=$gmp_support
+COPY examples/cpp/mp_double.h ./mp_double.h
 
 RUN cmake . && make
 

examples/cpp/mp_double.h

@@ -0,0 +1,155 @@
+#include <utility>
+#include <cstdio>
+
+#include <mpfr.h>
+
+class MpDouble final {
+    mpfr_t value;
+public:
+
+    // copy constructor
+    MpDouble(const MpDouble &original) noexcept {
+        mpfr_init(value);
+        mpfr_set(value, original.value, MPFR_RNDN);
+    }
+    // move constructor
+    MpDouble(MpDouble &&original) noexcept {
+        mpfr_init(value); // this is needed to keep original is a safe state for destruction
+        swap(original);
+    }
+    // implicit constructor
+    MpDouble(const double literal) noexcept {
+        mpfr_init(value);
+        mpfr_set_d(value, literal, MPFR_RNDN);
+    }
+    // explicit constructors
+    explicit MpDouble(const char *literal) noexcept {
+        mpfr_init(value);
+        mpfr_set_str(value, literal, 10, MPFR_RNDN);
+    }
+    explicit MpDouble(mpfr_t val) noexcept {
+        mpfr_swap(value, val);
+    }
+    // swap
+    inline void swap(MpDouble &other) noexcept {
+        mpfr_swap(value, other.value);
+    }
+    // destructor
+    ~MpDouble() {
+        mpfr_clear(value);
+    }
+
+    static void setPreccisionInBits(int bits) {
+        mpfr_set_default_prec(bits);
+    }
+
+    static void cleanUp() {
+        mpfr_free_cache();
+    }
+
+    // explicit static cast to double
+    explicit operator double() const {
+        return mpfr_get_d(value, MPFR_RNDN);
+    }
+
+    // assignment operation
+    MpDouble &operator=(const MpDouble &other) noexcept {
+        if(this == &other) return *this;
+        mpfr_set(value, other.value, MPFR_RNDN);
+        return *this;
+    }
+    // move assignment operation
+    MpDouble &operator=(MpDouble &&other) noexcept {
+        if(this == &other) return *this;
+        swap(other);
+        return *this;
+    }
+
+    // arithmetic operations
+    MpDouble &operator+=(const MpDouble &other) {
+        mpfr_add(value, value, other.value, MPFR_RNDN);
+        return *this;
+    }
+    MpDouble &operator-=(const MpDouble &other) {
+        mpfr_sub(value, value, other.value, MPFR_RNDN);
+        return *this;
+    }
+    MpDouble &operator/=(const MpDouble &other) {
+        mpfr_div(value, value, other.value, MPFR_RNDN);
+        return *this;
+    }
+    MpDouble &operator*=(const MpDouble &other) {
+        mpfr_mul(value, value, other.value, MPFR_RNDN);
+        return *this;
+    }
+
+    friend MpDouble operator+(const MpDouble &a, const MpDouble &b);
+    friend MpDouble operator-(const MpDouble &a, const MpDouble &b);
+    friend MpDouble operator/(const MpDouble &a, const MpDouble &b);
+    friend MpDouble operator*(const MpDouble &a, const MpDouble &b);
+
+    friend bool operator>(const MpDouble &a, const MpDouble &b);
+    friend bool operator<(const MpDouble &a, const MpDouble &b);
+    friend bool operator>=(const MpDouble &a, const MpDouble &b);
+    friend bool operator<=(const MpDouble &a, const MpDouble &b);
+    friend bool operator!=(const MpDouble &a, const MpDouble &b);
+    friend bool operator==(const MpDouble &a, const MpDouble &b);
+
+    friend void swap(MpDouble &a, MpDouble &b);
+};
+
+void swap(MpDouble &a, MpDouble &b) {
+    a.swap(b);
+}
+
+MpDouble operator+(const MpDouble &a, const MpDouble &b) {
+    mpfr_t result;
+    mpfr_init(result);
+    mpfr_add(result, a.value, b.value, MPFR_RNDN);
+    return MpDouble {result};
+}
+
+MpDouble operator-(const MpDouble &a, const MpDouble &b) {
+    mpfr_t result;
+    mpfr_init(result);
+    mpfr_sub(result, a.value, b.value, MPFR_RNDN);
+    return MpDouble {result};
+}
+
+MpDouble operator/(const MpDouble &a, const MpDouble &b) {
+    mpfr_t result;
+    mpfr_init(result);
+    mpfr_div(result, a.value, b.value, MPFR_RNDN);
+    return MpDouble {result};
+}
+
+MpDouble operator*(const MpDouble &a, const MpDouble &b) {
+    mpfr_t result;
+    mpfr_init(result);
+    mpfr_mul(result, a.value, b.value, MPFR_RNDN);
+    return MpDouble {result};
+}
+
+bool operator>(const MpDouble &a, const MpDouble &b) {
+    return mpfr_greater_p(a.value, b.value);
+}
+
+bool operator<(const MpDouble &a, const MpDouble &b) {
+    return mpfr_less_p(a.value, b.value);
+}
+
+bool operator>=(const MpDouble &a, const MpDouble &b) {
+    return mpfr_greaterequal_p(a.value, b.value);
+}
+
+bool operator<=(const MpDouble &a, const MpDouble &b) {
+    return mpfr_lessequal_p(a.value, b.value);
+}
+
+bool operator==(const MpDouble &a, const MpDouble &b) {
+    return mpfr_equal_p(a.value, b.value);
+}
+
+bool operator!=(const MpDouble &a, const MpDouble &b) {
+    return mpfr_lessgreater_p(a.value, b.value);
+}