It finds the fist natural solution of the diophantine equation
        a*x + b*y = c. Some lines of this code are taken from the project
        sympy.

        :param c: constant
        :param a: quotient of x
        :param b: quotient of y
        :return: the first natural solution of the diophatine equation

def diop_natural_solution_linear(c, a, b):
        """
        It finds the fist natural solution of the diophantine equation
        a*x + b*y = c. Some lines of this code are taken from the project
        sympy.

        :param c: constant
        :param a: quotient of x
        :param b: quotient of y
        :return: the first natural solution of the diophatine equation
        """
        def get_intersection(a, b, a_dir, b_dir):
            # Do the intersection between two
            # ranges.
            if (a_dir, b_dir) == (">=", ">="):
                lb = a if a > b else b
                ub = float('inf')
            elif (a_dir, b_dir) == ("<=", ">="):
                if a > b:
                    lb = b
                    ub = a
                else:
                    lb = None
                    ub = None
            elif (a_dir, b_dir) == (">=", "<="):
                if b > a:
                    lb = a
                    ub = b
                else:
                    lb = None
                    ub = None
            elif (a_dir, b_dir) == ("<=", "<="):
                ub = a if a < b else b
                lb = float('-inf')

            return lb, ub

        d = StridedInterval.igcd(a, StridedInterval.igcd(b, c))
        a = a // d
        b = b // d
        c = c // d

        if c == 0:
            return (0, 0)
        else:
            x0, y0, d = StridedInterval.extended_euclid(int(abs(a)), int(abs(b)))
            x0 = x0 * StridedInterval.sign(a)
            y0 = y0 * StridedInterval.sign(b)

            if c % d == 0:
                """
                Integer solutions are: (c*x0 + b*t, c*y0 - a*t)
                we have to get the first positive solution, which means
                that we have to solve the following disequations for t:
                c*x0 + b*t >= 0 and c*y0 - a*t >= 0.
                """
                assert b != 0
                assert a != 0

                t0 = (-c * x0) / float(b)
                t1 = (c * y0) / float(a)
                # direction of the disequation depends on b and a sign
                if b < 0:
                    t0_dir = '<='
                else:
                    t0_dir = '>='
                if a < 0:
                    t1_dir = '>='
                else:
                    t1_dir = '<='

                # calculate the intersection between the found
                # solution intervals to get the common solutions
                # for t.
                lb, ub = get_intersection(t0, t1, t0_dir, t1_dir)

                # Given that we are looking for the first value
                # which solve the diophantine equation, we have to
                # select the value of t closer to 0.
                if lb <= 0 and ub >= 0:
                    t = ub if abs(ub) < abs(lb) else lb
                elif lb == float('inf') or lb == float("-inf"):
                    t = ub
                elif ub == float('inf') or ub == float("-inf"):
                    t = lb
                else:
                    t = ub if abs(ub) < abs(lb) else lb
                # round the value of t
                if t == ub:
                    t = int(math.floor(t))
                else:
                    t = int(math.ceil(t))

                return (c*x0 + b*t, c*y0 - a*t)
            else:
                return (None, None)