import math


def pretty(values, classes=5, start=None, end=None):
    """
    R's pretty algorithm implemented in Python
    Code based on R implementation from 'labeling' R package.
    
    Compute a sequence of about 'n+1' equally spaced 'round' values
    which cover the range of the values in 'values'.  The values are chosen
    so that they are 1, 2 or 5 times a power of 10.

    Returns a number of breaks not necessarily equal to 'classes' using 
    rpretty, but likely to be legible.

    Parameters:
        values : list of input values
        classes     : number of class intervals
    """

    # if too few values, just return breakpoints for each unique value, ignoring classes
    if values and len(values) == 1:
        return values * 2

    # auto detect value range from values if not custom specified
    if start == None:
        dmin = min(values)
    else:
        dmin = start
    if end == None:
        dmax = max(values)
    else:
        dmax = end
        
    # begin
    n = classes
    min_n = int(n / 3) # Nonnegative integer giving the minimal number of intervals
    shrink_sml = 0.75  # Positive numeric by a which a default scale is shrunk 
                                         # in the case when range(x) is very small (usually 0).
    high_u_bias = 1.5  # Non-negative numeric, typically > 1. The interval unit 
                                         # is determined as {1,2,5,10} times b, a power of 10.
                                         # Larger high.u.bias values favor larger units
    u5_bias = 0.5 + 1.5 * high_u_bias
                                         # Non-negative numeric multiplier favoring 
                                         # factor 5 over 2.
    h = high_u_bias
    h5 = u5_bias
    ndiv = n

    dx = dmax - dmin
    if dx == 0 and dmax == 0:
        cell = 1.0
        i_small = True
        U = 1
    else:
        cell = max(abs(dmin), abs(dmax))
        if h5 >= 1.5 * h + 0.5:
            U = 1 + (1.0/(1 + h))
        else:
            U = 1 + (1.5 / (1 + h5))
        i_small = dx < (cell * U * max(1.0, ndiv) * 1e-07 * 3.0)

    if i_small:
        if cell > 10:
            cell = 9 + cell / 10
            cell = cell * shrink_sml
        if min_n > 1:
            cell = cell / min_n
    else:
        cell = dx
        if ndiv > 1:
            cell = cell / ndiv
    if cell < 20 * 1e-07:
        cell = 20 * 1e-07
    
    base = 10.0**math.floor(math.log10(cell))
    unit = base
    if (2 * base) - cell < h * (cell - unit):
        unit = 2.0 * base
        if (5 * base) - cell < h5 * (cell - unit):
            unit = 5.0 * base
            if (10 * base) - cell < h * (cell - unit):
                unit = 10.0 * base
    # Maybe used to correct for the epsilon here??
    ns = math.floor(dmin / unit + 1e-07)
    nu = math.ceil(dmax / unit - 1e-07)

    # Extend the range out beyond the data. Does this ever happen??
    while ns * unit > dmin + (1e-07 * unit):
        ns = ns-1
    while nu * unit < dmax - (1e-07 * unit):
        nu = nu+1
    # If we don't have quite enough labels, extend the range out to make more (these labels are beyond the data :( )
    k = math.floor(0.5 + nu-ns)
    if k < min_n:
        k = min_n - k
        if ns >= 0:
            nu = nu + k / 2
            ns = ns - k / 2 + k % 2
        else:
            ns = ns - k / 2
            nu = nu + k / 2 + k % 2
        ndiv = min_n
    else:
        ndiv = k
    graphmin = ns * unit
    graphmax = nu * unit
    count = int(math.ceil(graphmax - graphmin)/unit)
    res = [graphmin + k*unit for k in range(count+1)]
    if res[0] < dmin:
        res[0] = dmin
    if res[-1] > dmax:
        res[-1] = dmax
    return res