Added arbitrary-precision arithmetic to native module

3 files changed, 152 insertions, 94 deletions

diff --git a/native/include/edf/gel_pl.h b/native/include/edf/gel_pl.h
index 4093d95..1127d50 100644
--- a/native/include/edf/gel_pl.h
+++ b/native/include/edf/gel_pl.h
@@ -13,17 +13,21 @@ class GELPl
     int no_cpus;
     const TaskSet& tasks;
     int rounds;
-    std::vector<double> S_i;
-    std::vector<double> G_i;
+    std::vector<fractional_t> S_i;
+    std::vector<fractional_t> G_i;
 
     // For faster lookups, to avoid too many conversions.
-    std::vector<double> utilizations;
+    std::vector<fractional_t> utilizations;
 
-    double compute_exact_s(double S, const std::vector<double>& Y_ints);
-    double compute_binsearch_s(double S, const std::vector<double>& Y_ints);
+    void compute_exact_s(const fractional_t& S,
+                         const std::vector<fractional_t>& Y_ints,
+                         fractional_t& s);
+    void compute_binsearch_s(const fractional_t& S,
+                             const std::vector<fractional_t>& Y_ints,
+                             fractional_t& s);
 
-    inline double compute_M(double s, double S,
-                                   const std::vector<double>& Y_ints);
+    inline bool M_lt_0(const fractional_t& s, const fractional_t& S,
+                       const std::vector<fractional_t>& Y_ints);
 
     // These are basically just structs that override operator< to allow
     // sort algorithms to work.
@@ -31,8 +35,8 @@ class GELPl
      public:
         unsigned int old_task;
         unsigned int new_task;
-        double location;
-        double old_task_utilization;
+        fractional_t location;
+        fractional_t old_task_utilization;
 
         bool operator<(const ReplacementType& other) const {
             return (location < other.location)
@@ -44,7 +48,7 @@ class GELPl
     class TaggedValue {
      public:
         unsigned int task;
-        double value;
+        fractional_t value;
 
         //Order is reversed - we are going to want the largest, rather than the
         //smallest, values.
@@ -68,12 +72,14 @@ class GELPl
         return bounds[index];
    }
 
+   // Converted to double for the sake of Python
    double get_Si(unsigned int index) {
-        return S_i[index];
+        return S_i[index].get_d();
    }
 
+   // Converted to double for the sake of Python
    double get_Gi(unsigned int index) {
-        return G_i[index];
+        return G_i[index].get_d();
    }
 };
 
diff --git a/native/src/edf/gel_pl.cpp b/native/src/edf/gel_pl.cpp
index 53ae575..f0f9303 100644
--- a/native/src/edf/gel_pl.cpp
+++ b/native/src/edf/gel_pl.cpp
@@ -6,8 +6,10 @@
 #include <limits>
 #include <algorithm>
 #include <cmath>
+#include <iostream>
 
-static bool reversed_order(double first, double second) {
+static bool reversed_order(const fractional_t& first,
+                           const fractional_t& second) {
     return second < first;
 }
 
@@ -16,8 +18,8 @@ GELPl::GELPl(Scheduler sched, unsigned int num_processors, const TaskSet& ts,
 :no_cpus(num_processors), tasks(ts), rounds(num_rounds)
 {
     std::vector<unsigned long> pps;
-    double S;
-    std::vector<double> Y_ints;
+    fractional_t S = 0;
+    std::vector<fractional_t> Y_ints;
 
     int task_count = tasks.get_task_count();
     // Reserve capacity in all vectors to minimize allocation costs.
@@ -29,8 +31,8 @@ GELPl::GELPl(Scheduler sched, unsigned int num_processors, const TaskSet& ts,
     // For faster lookups
     utilizations.reserve(task_count);
     for (int i = 0; i < task_count; i++) {
-        utilizations.push_back(double(tasks[i].get_wcet())
-                               / double(tasks[i].get_period()));
+        utilizations.push_back(tasks[i].get_wcet());
+        utilizations[i] /= tasks[i].get_period();
     }
     
     unsigned long min_pp = std::numeric_limits<unsigned long>::max();
@@ -50,36 +52,60 @@ GELPl::GELPl(Scheduler sched, unsigned int num_processors, const TaskSet& ts,
 
     // Reduce to compute minimum.  Also compute Y intercepts, S_i values, and
     // S.
-    S = 0.0;
     for (int i = 0; i < task_count; i++) {
         pps[i] -= min_pp;
         const Task& task = tasks[i];
-        double wcet = double(task.get_wcet());
-        double period = double(task.get_period());
-        S_i[i] = std::max(0.0, wcet * (1.0 -  double(pps[i])/ period));
-        S += S_i[i];
-        Y_ints.push_back((0.0 - wcet/no_cpus) * (wcet / period)
-                         + task.get_wcet() - S_i[i]);
+        unsigned long wcet = task.get_wcet();
+        unsigned long period = task.get_period();
+        S_i.push_back(pps[i]);
+        fractional_t& S_i_i = S_i[i];
+        S_i_i *= -1;
+        S_i_i /= period;
+        S_i_i += 1;
+        S_i_i *= wcet;
+        if (S_i_i < 0) {
+            S_i_i = 0;
+        }
+        S += S_i_i;
+        Y_ints.push_back(wcet);
+        fractional_t& Y_ints_i = Y_ints[i];
+        Y_ints_i *= -1;
+        Y_ints_i /= no_cpus;
+        Y_ints_i *= utilizations[i];
+        Y_ints_i += wcet;
+        Y_ints_i -= S_i_i;
     }
 
-    double s;
+    fractional_t s;
     if (rounds == 0) {
-        s = compute_exact_s(S, Y_ints);
+        compute_exact_s(S, Y_ints, s);
     }
     else {
-        s = compute_binsearch_s(S, Y_ints);
+        compute_binsearch_s(S, Y_ints, s);
     }
 
     for (int i = 0; i < task_count; i++) {
+        fractional_t x_i = s;
+        fractional_t x_comp = tasks[i].get_wcet();
+        x_comp /= no_cpus;
+        x_i -= x_comp;
+        // Compute ceiling
+        integral_t xi_ceil = x_i.get_num();
+        mpz_cdiv_q(xi_ceil.get_mpz_t(),
+                   x_i.get_num().get_mpz_t(),
+                   x_i.get_den().get_mpz_t());
         bounds.push_back(pps[i]
                          + tasks[i].get_wcet()
-                         + (unsigned long)std::ceil(
-                         s - (double(tasks[i].get_wcet() / double(no_cpus)))));
-        G_i.push_back(Y_ints[i] + s * utilizations[i]);
+                         + xi_ceil.get_ui());
+        G_i.push_back(s);
+        G_i[i] *= utilizations[i];
+        G_i[i] += Y_ints[i];
     }
 }
 
-double GELPl::compute_exact_s(double S, const std::vector<double>& Y_ints) {
+void GELPl::compute_exact_s(const fractional_t& S,
+                            const std::vector<fractional_t>& Y_ints,
+                            fractional_t& s) {
     int task_count = tasks.get_task_count();
 
     std::vector<ReplacementType> replacements;
@@ -88,19 +114,22 @@ double GELPl::compute_exact_s(double S, const std::vector<double>& Y_ints) {
             // We can ignore parallel and identical lines - either don't
             // intersect or we don't care which is picked.
             if (utilizations[i] != utilizations[j]) {
-                double intersect = (Y_ints[j] - Y_ints[i])
-                                   / (utilizations[i] - utilizations[j]);
+                fractional_t intersect_den = utilizations[i];
+                intersect_den -= utilizations[j];