def pseudo_join(s, b, smart_join=True):
        """
        It two intervals in a way that the resulting SI is the one that has the least
        SI cardinality (i.e., which represents the least number of elements) possible if the smart_join flag is enabled,
        otherwise it just joins the SI according the order they are passed to the function.

        The pseudo-join operation is not associative in wrapping intervals (please refer to section 3.1 paper
        'Signedness-Agnostic Program Analysis: Precise Integer Bounds for Low-Level Code'), Therefore the join of three
        WI may  give us different results according on the order we join them. All of the results will be sound, though.

        Please use the function least_upper_bound as a stub.

        :param s:           The first SI
        :param b:           The other SI.
        :param smart_join:  Enable the smart join behavior. If this flag is set, this function joins the two SI in a way
                            that the resulting Si has least number of elements (more precise). If it is unset, this
                            function will join the two SI according on the order they are passed to the function.
        :return:            A new StridedInterval
        """

        assert s.bits == b.bits
        w = s.bits

        if s._reversed != b._reversed:
            logger.warning('Incoherent reversed flag between operands %s and %s', s, b)

        uninit_flag = s.uninitialized | b.uninitialized

        #
        # Trivial cases
        #

        if s.is_empty:
            return b
        if b.is_empty:
            return s

        if s.is_integer and b.is_integer:
            u = max(s.upper_bound, b.upper_bound) if smart_join else b.upper_bound
            l = min(s.lower_bound, b.lower_bound) if smart_join else s.lower_bound
            stride = abs(u - l)
            return StridedInterval(bits=w, stride=stride, lower_bound=l, upper_bound=u, uninitialized=uninit_flag)

        #
        # Other cases
        #

        if s._is_surrounded(b):
            # Containment: s <= b
            new_stride = StridedInterval.gcd(s.stride, b.stride) if s.is_interval else b.stride
            new_stride = StridedInterval.gcd(new_stride, s._modular_sub(s.lower_bound, b.lower_bound, w))
            return StridedInterval(bits=w, stride=new_stride, lower_bound=b.lower_bound,
                                   upper_bound=b.upper_bound, uninitialized=uninit_flag)

        elif b._is_surrounded(s):
            # Containment: b <= s
            # TODO: This case is missing in the original implementation. Is that a bug?
            new_stride = StridedInterval.gcd(s.stride, b.stride) if b.is_interval else s.stride
            new_stride = StridedInterval.gcd(new_stride, s._modular_sub(b.lower_bound, s.lower_bound, w))
            return StridedInterval(bits=w, stride=new_stride, lower_bound=s.lower_bound,
                                   upper_bound=s.upper_bound, uninitialized=uninit_flag)

        elif (s._surrounds_member(b.lower_bound) and s._surrounds_member(b.upper_bound) and
            b._surrounds_member(s.lower_bound) and b._surrounds_member(s.upper_bound)):
            # The union of them covers the entire sphere
            return StridedInterval.top(w, uninitialized=uninit_flag)

        elif s._surrounds_member(b.lower_bound):
            # Overlapping. Nor s or b are integer here.
            # We return the join with less values
            new_stride = StridedInterval.gcd(s.stride, b.stride)
            new_stride = StridedInterval.gcd(new_stride, s._modular_sub(b.lower_bound, s.lower_bound, w))
            return StridedInterval(bits=w, stride=new_stride, lower_bound=s.lower_bound,
                                    upper_bound=b.upper_bound, uninitialized=uninit_flag)
        elif b._surrounds_member(s.lower_bound):
            # Overlapping. Nor s or b are integer here.
            # We return the join with less values
            new_stride = StridedInterval.gcd(s.stride, b.stride)
            new_stride = StridedInterval.gcd(new_stride, s._modular_sub(s.lower_bound, b.lower_bound, w))
            return StridedInterval(bits=w, stride=new_stride, lower_bound=b.lower_bound,
                                   upper_bound=s.upper_bound, uninitialized=uninit_flag)
        else:
            # no overlapping.
            # we join the two intervals according on the order they are given
            if not smart_join:
                if s.is_integer:
                    new_stride = StridedInterval.gcd(b.stride, s._modular_sub(b.lower_bound, s.lower_bound, w))
                elif b.is_integer:
                    new_stride = StridedInterval.gcd(s.stride, s._modular_sub(b.lower_bound, s.lower_bound, w))
                else:
                    new_stride = StridedInterval.gcd(s.stride, b.stride)
                    new_stride = StridedInterval.gcd(new_stride, StridedInterval._wrapped_cardinality(s.lower_bound, b.lower_bound, w) - 1)
                return StridedInterval(bits=w, stride=new_stride, lower_bound=s.lower_bound, upper_bound=b.upper_bound, uninitialized=uninit_flag)

            # Else: smart join.
            # we return the join which produce an interval with the least number of values
            if s.is_integer:
                new_stride = b.stride
            elif b.is_integer:
                new_stride = s.stride
            else:
                new_stride = StridedInterval.gcd(s.stride, b.stride)

            # from b to s
            new_stride1 = StridedInterval.gcd(new_stride, StridedInterval._wrapped_cardinality(b.lower_bound, s.lower_bound, w) - 1)
            # from s to b
            new_stride2 = StridedInterval.gcd(new_stride, StridedInterval._wrapped_cardinality(s.lower_bound, b.lower_bound, w) - 1)
            si1 = StridedInterval(bits=w, stride=new_stride1, lower_bound=b.lower_bound,
                                   upper_bound=s.upper_bound, uninitialized=uninit_flag)
            si2 = StridedInterval(bits=w, stride=new_stride2, lower_bound=s.lower_bound,
                                   upper_bound=b.upper_bound, uninitialized=uninit_flag)

            if si1.n_values <= si2.n_values:
                return si1
            else:
                return si2