diff --git a/src/sage/graphs/hyperbolicity.pyx b/src/sage/graphs/hyperbolicity.pyx
index a34675f8f40..152084f8978 100644
--- a/src/sage/graphs/hyperbolicity.pyx
+++ b/src/sage/graphs/hyperbolicity.pyx
@@ -785,6 +785,181 @@ cdef tuple __hyperbolicity__(int N,
         return (h, certificate, h_UB if GOTO_RETURN else h)
 
 
+######################################################################
+# Compute the hyperbolicity using the algorithm of [BCC12]_
+######################################################################
+
+cdef tuple __hyper_new__(int N,
+                         unsigned short **  distances,
+                         unsigned short **  far_apart_pairs,
+                         int D,
+                         int h_LB,
+                         float approximation_factor,
+                         float additive_gap,
+                         verbose = False):
+    """
+    Return the hyperbolicity of a graph.
+
+    This method implements the exact and the approximate algorithms proposed in
+    [CCL12]_. See the module's documentation for more details.
+
+    This method assumes that the graph under consideration is connected.
+
+    INPUTS:
+
+    - ``N`` -- number of vertices of the graph
+
+    - ``distances`` -- path distance matrix
+
+    - ``far_apart_pairs`` -- 0/1 matrix of far-apart pairs. Pair ``(i,j)`` is
+      far-apart if ``far_apart_pairs[i][j]\neq 0``.
+
+    - ``D`` -- diameter of the graph
+
+    - ``h_LB`` -- lower bound on the hyperbolicity
+
+    - ``approximation_factor`` -- When the approximation factor is set to some
+      value larger than 1.0, the function stop computations as soon as the ratio
+      between the upper bound and the best found solution is less than the
+      approximation factor. When the approximation factor is 1.0, the problem is
+      solved optimaly.
+
+     - ``additive_gap`` -- When sets to a positive number, the function stop
+       computations as soon as the difference between the upper bound and the
+       best found solution is less than additive gap. When the gap is 0.0, the
+       problem is solved optimaly.
+
+    - ``verbose`` -- (default: ``False``) is boolean set to ``True`` to display
+      some information during execution
+
+    OUTPUTS:
+
+    This function returns a tuple ( h, certificate, h_UB ), where:
+
+    - ``h`` -- is an integer. When 4-tuples with hyperbolicity larger or equal
+     to `h_LB are found, h is the maximum computed value and so twice the
+     hyperbolicity of the graph. If no such 4-tuple is found, it returns -1.
+
+    - ``certificate`` -- is a list of vertices. When 4-tuples with hyperbolicity
+      larger that h_LB are found, certificate is the list of the 4 vertices for
+      which the maximum value (and so the hyperbolicity of the graph) has been
+      computed. If no such 4-tuple is found, it returns the empty list [].
+
+    - ``h_UB`` -- is an integer equal to the proven upper bound for `h`. When
+      ``h == h_UB``, the returned solution is optimal.
+    """
+    cdef int hh # can get negative value
+    cdef int a, b, c, d, h, h_UB
+    cdef int x, y, l1, l2, S1, S2, S3
+    cdef list certificate = []
+    cdef uint32_t *nb_p = <uint32_t *>sage_malloc(sizeof(uint32_t)) # The total number of pairs.
+
+    # Test if the distance matrix corresponds to a connected graph, i.e., if
+    # distances from node 0 are all less or equal to N-1.
+    for a from 0 <= a < N:
+        if distances[0][a]>=N:
+            raise ValueError("The input graph must be connected.")
+
+    # nb_pairs_of_length[d] is the number of pairs of vertices at distance d
+    cdef uint32_t * nb_pairs_of_length = <uint32_t *>sage_malloc((D+1) * sizeof(uint32_t))
+
+    if (nb_pairs_of_length == NULL):
+        raise MemoryError
+
+    cdef pair ** pairs_of_length = sort_pairs(N, D, distances, far_apart_pairs, nb_p, nb_pairs_of_length)
+    
+    if verbose:
+        print "Current 2 connected component has %d vertices and diameter %d" %(N,D)
+        if far_apart_pairs == NULL:
+            print "Number of pairs: %d" %(nb_p[0])
+            print "Repartition of pairs:", [ (i, nb_pairs_of_length[i]) for i in range(1, D+1) if nb_pairs_of_length[i]>0]
+        else:
+            print "Number of far-apart pairs: %d\t(%d pairs in total)" %(nb_p[0], binomial(N, 2))
+            print "Repartition of far-apart pairs:", [ (i, nb_pairs_of_length[i]) for i in range(1, D+1) if nb_pairs_of_length[i]>0]
+
+
+    approximation_factor = min(approximation_factor, D)
+    additive_gap = min(additive_gap, D)
+
+    # We use some short-cut variables for efficiency
+    cdef pair * pair_ab
+    cdef pair * pair_cd
+    cdef uint32_t ** mate
+    cdef uint32_t * cont_mate
+    h_UB = D
+    
+    for pair_ab from pairs_of_length[0] + 1 <= pair_ab < pairs_of_length[D]:
+        a = pair_ab.s
+        b = pair_ab.t
+        
+        dist_a = distances[a]
+        dist_b = distances[b]
+        h_UB = dist_a[b]
+        
+                # Termination if required approximation is found
+        if certificate and ( (h_UB <= h*approximation_factor) or (h_UB-h <= additive_gap) ):
+            GOTO_RETURN = 1
+            break
+
+        # If we cannot improve further, we stop
+        #
+        # See the module's documentation for a proof that this cut is
+        # valid. Remember that the triples are sorted in a specific order.
+        if h_UB <= h:
+            h_UB = h
+            break
+        
+        for pair_cd from pairs_of_length[0] <= pair_cd < pair_ab:
+            c = pair_cd.s
+            d = pair_cd.t
+            S2 = dist_a[c] + dist_b[d]
+            S3 = dist_a[d] + dist_b[c]
+            
+            if S2 > S3:
+                hh = S1 - S2
+            else:
+                hh = S1 - S3
+
+            if h < hh or not certificate:
+                # We update current bound on the hyperbolicity and the
+                # search space.
+                #
+                # Note that if hh==0, we first make sure that a,b,c,d are
+                # all distinct and are a valid certificate.
+                if hh > 0 or not (a == c or a == d or b == c or b == d):
+                    h = hh
+                    certificate = [a, b, c, d]
+
+                    if verbose:
+                        print "New lower bound:", ZZ(hh) / 2
+
+                    # If we cannot improve further, we stop
+                    if l2 <= h:
+                        GOTO_RETURN = 1
+                        h_UB = h
+                        break
+
+                    # Termination if required approximation is found
+                    if (h_UB <= h * approximation_factor) or (h_UB - h <= additive_gap):
+                        GOTO_RETURN = 1
+                        break
+
+
+
+    # We now free the memory
+    sage_free(nb_pairs_of_length)
+    sage_free(pairs_of_length[0])
+    sage_free(pairs_of_length)
+
+    # Last, we return the computed value and the certificate
+    if len(certificate) == 0:
+        return ( -1, [], h_UB )
+    else:
+        # When using far-apart pairs, the loops may end before improving the
+        # upper-bound
+        return (h, certificate, h_UB if GOTO_RETURN else h)
+
+
 def hyperbolicity(G, algorithm='CCL+FA', approximation_factor=None, additive_gap=None, verbose = False):
     r"""
     Returns the hyperbolicity of the graph or an approximation of this value.