#!/usr/bin/env python3

import sys
import random
from fractions import Fraction
import re
from math import sqrt, ceil, floor

no_svg = False
try:
	import svgwrite
	from svgwrite import cm, mm
except ImportError:
    print("Warning: svgwrite is not installed. No picture will be created.")
    no_svg = True



filter_ET = False
filter_qx = False
check_pattern = False
use_frac = False
print_tree = False
print_grid = False
better = False
count_iter = False
best_iter = None

help = """
Usage: ./critical.py <n> <p> <x> [options]

n = number of jobs
p = either specific p value or a range in the form [10:20]. 
    a single number p will be interpreted as the range [0:p].
x = similar

Will produce a picture in a file named n{n}_p{p}_x{p}.svg. Every grid cell
contains to a specific (p,x) value in the range, and the color corresponds to
the equilibrium schedule of this range.

The equilibrium schedule is computed by exploring exhaustively the min-max
game tree. But this can be improved using the following option:

  -better     : uses the polynomial time algorithm to compute actual ratio
  
The computation is done using floating point numbers, unless the following
option is used. In that case p or x values have to be written as fractions
such as 10/4 unless they are integral.

  -fraction   : compute exactly using fractions

The output can be controlled using the following options.

  -filterET   : will show only the algorithm's decision in the schedules
  -filterqx   : will show only the adversary's decision in the schedules
  -pattern    : show schedules which do not follow the pattern (Tx)*(Tp)*(Ex)*(Ep)*.
  -count_iter : counts the number of iterations of the polynomial time algorithm

The full game tree can be shown for a specific p, x value pair with the
following option, producing a file name tree_n{n}_p{p}_x{x}.dot which can be
rendered into PDF using dot <file.dot> -T pdf -o <file.pdf>.

  -tree       : print the game tree for the given p, x values

Each step of the polynomial time algorithm can be shown for a specific p, x
value pair with the following option, producing a file name
grid_n{n}_p{p}_x{x}.html.

  -grid       : print the grid for specific p, x values
"""

def read_param(s):
    if "/" in s:
        assert use_frac
        return Fraction(s)
    elif "." in s:
        return float(s)
    else:
        return int(s)


def read_interval(s):
    assert s[0] == "[" and s[-1] == "]"
    t = s[1:-1].split(":")
    assert len(t) == 2
    return (read_param(t[0]), read_param(t[1]))


if __name__ == "__main__":
    pattern = re.compile('<(Tx)*(Tp)*(Ex)*(Ep)*x*>')

    if len(sys.argv) < 4:
        print(help)
        sys.exit(1)


    for a in sys.argv[4:]:
        filter_ET |= a == "-filterET"
        filter_qx |= a == "-filterqx"
        use_frac  |= a == "-fraction"
        check_pattern |= a == "-pattern"
        print_tree |= a == "-tree"
        print_grid |= a == "-grid"
        better |= a == "-better"
        count_iter |= a == "-count_iter"
        if not (filter_ET or filter_qx or use_frac or check_pattern or
                print_tree or print_grid or better):
            print("is_bad option")
            sys.exit(1)

    n = int(sys.argv[1])

    param_p = sys.argv[2]
    if param_p[0] == '[':
        min_p, max_p = read_interval(param_p)
    else:
        min_p = 0
        max_p = read_param(param_p)

    param_x = sys.argv[3]
    if param_x[0] == '[':
        min_x, max_x = read_interval(param_x)
    else:
        min_x = 0
        max_x = read_param(param_x)


def fmt_ratio(r):
    if use_frac:
        return str(r)
    else:
        return "%.6f" % r

# worst_ratio_ratio = 0

def optimal_ratio(n, p, x):
    # non adaptive competitive ratio
    if x < 2 + 1/p:
        return sqrt(1 + x / p)
    else:
        D = 8 * p * (x - 1) * x * x + (1 + p * x - x * x) ** 2
        return 1 + (x * x - p * x + sqrt(D)) / (2 * p * x * x)

def max_ratio(p, x, a, b):
    """a, b are pairs of ratio and corresponding schedule in string form
    assuming a is for p and b for x
    """
    # if a[0] == b[0] and a[0] > 1:
    #     print('arbitrary tie breaking in max, p=%f, x=%f, a=%f, b=%f' % (p, x, a, b))
    is_bad =  a[2] == BAD or b[2] == BAD
    if a[0] >= b[0]:
        return (a[0], a[1], BAD if is_bad else a[2])
    else:
        return (b[0], b[1], BAD if is_bad else b[2])


def min_ratio(p, x, a, b):
    # assuming a is for T and b for E 
    # if a[0] == b[0] and a[0] > 1:
    #     print('arbitrary tie breaking in min, p=%f, x=%f, a=%f, b=%f' % (p, x, a, b))
    is_bad =  a[2] == BAD or b[2] == BAD
    if a[0] <= b[0]:
        return (a[0], a[1], BAD if is_bad else GOOD_T)
    else:
        return (b[0], b[1], BAD if is_bad else b[2])

BAD = 0
GOOD_E = 1
GOOD_T = 2

def my_div(a, b):
    if b:
        return a / b
    elif a:
        return float('inf')
    else:
        return 1    # = 0 / 0


def cost(p, x, schedule):
    n = len(schedule) // 2
    rank = n 
    obj = p * 0      # will become a fraction if p is a fraction
    for i in range(n):
        a = schedule[2 * i: 2 * i + 2]
        if a == 'Tx':
            obj += rank
        elif a == 'Tp':
            obj += (1 + p) * rank
            rank -= 1
        elif a == 'Ex':
            obj += (p + x) * rank
            rank -= 1
        elif a == 'Ep':
            obj += p * rank
            rank -= 1
    while rank:
        obj += (p + x) * rank
        rank -= 1
    return obj


def compute_ratio(p, x, schedule):
    np = schedule.count("p")
    nx = schedule.count("x")
    opt = "Ep" * np + "Ex" * nx
    return my_div(cost(p, x, schedule), cost(p, x, opt))


def algorithm(p, x, schedule, f):
    """computes the equilibrium ratio from a node of the game tree where

    - the lengths of the short and long jobs are p and p+x resp.
    - the algorithm makes the first move
    - the schedule so far is schedule
    
    :returns: (ratio, schedule, label_of_node_in_game_tree)
    """
    if len(schedule) == 2 * n:
        postponed = schedule.count("Tx")
        nx = schedule.count("x")
        if filter_ET:
            S = "T{} E{}".format(schedule.count("T"), schedule.count("E"))
        elif filter_qx:
            np = n - nx
            S = "p{} x{}".format(np, nx)
        else:
            S = schedule + "x" * postponed
        r = compute_ratio(p, x, schedule)
        ratio_schedule = (r, S, GOOD_E)
        if f is not None:
            print('_%s [label="%s"];' % (schedule, fmt_ratio(ratio_schedule[0])), file=f)
        return ratio_schedule
    # internal node:
    assert len(schedule) < 2 * n
    alt_T = adversary(p, x, schedule + "T", f)
    alt_E = adversary(p, x, schedule + "E", f)
    best = min_ratio(p, x, alt_T, alt_E)
    if f is not None:
        if best[1] == alt_T[1]:
            print('_{} -> _{}T [label="T"];'.format(schedule, schedule), file=f)
            print('_{} -> _{}E [style=dotted,label="E"];'.format(schedule, schedule), file=f)
        else:
            print('_{} -> _{}T [style=dotted,label="T"];'.format(schedule, schedule), file=f)
            print('_{} -> _{}E [label="E"];'.format(schedule, schedule), file=f)
        print('_{} [label="{}"];'.format(schedule, fmt_ratio(best[0])), file=f)
    return best


def adversary(p, x, schedule, f):
    """ Now it is the turn of the adversary,
        the algorithm just decided to execute a job
    """
    alt_p = algorithm(p, x, schedule + "p", f)
    alt_x = algorithm(p, x, schedule + "x", f)
    best = max_ratio(p, x, alt_p, alt_x)
    if f is not None:
        if best[1] == alt_p[1]:
            print('_{} -> _{}p [label="p"];'.format(schedule, schedule), file=f)
            print('_{} -> _{}x [style=dotted,label="x"];'.format(schedule, schedule), file=f)
        else:
            print('_{} -> _{}p [style=dotted,label="p"];'.format(schedule, schedule), file=f)
            print('_{} -> _{}x [label="x"];'.format(schedule, schedule), file=f)
        print('_{} [label="{}"];'.format(schedule, fmt_ratio(best[0])), file=f)
    return best


# ---- helper functions for the dynamic program

fgrid = None

def print_the_grid(n, stop_ratio, direction, path, critical_cd):
    global fgrid
    if fgrid is None:
        fgrid = open(("grid_n%i_p%s_x%s.html" % (n, max_p, max_x)).replace("/","div"), "w")
        print("<style>" +
              "table {" +
              "  border: 1px solid black;" + 
              "  border-collapse: collapse;" +
              "}" +
              "td {"
              " text-align: center;" + 
              "}" +            
              "</style>" + 
              "Rows = #Tq, Columns = #Tp<br>" + 
              "Entries = stop ratio<br>" + 
              "Yellow = adversial strategy, Orange = algorithm's strategy" , file=fgrid)
    print("<p><table>", file=fgrid)
    for d in range(n + 1):
        print("<tr>", file=fgrid)
        for c in range(n - d + 1):
            if direction[c][d] == 'x':
                print("<td></td><td>↓</td>", file=fgrid)
            else:
                print("<td></td><td></td>", file=fgrid)
        print("</tr>", file=fgrid)
        print("<tr>", file=fgrid)
        for c in range(n - d + 1):
            if direction[c][d] == 'p':
                print("<td>→</td>", file=fgrid)
            else:
                print("<td></td>", file=fgrid)
            if stop_ratio[c][d]:
                if (c,d) in path:
                    if (c,d) == critical_cd:
                        color = "orange"
                    else:
                        color = "yellow"
                else:
                    color = "white"
                if use_frac:
                    print('<td width=80 bgcolor="{}">{}</td>'.format(color, 
                                    stop_ratio[c][d]), file=fgrid)
                else:
                    print('<td width=80 bgcolor="{}">{:.6f}</td>'.format(color, 
                                    stop_ratio[c][d]), file=fgrid)
            else:
                print("<td width=80 >0</td>", file=fgrid)
        print("</tr>", file=fgrid)
    print("</table>", file=fgrid)


def ratio_E_star(n, p, x, c, d, e):
    # try different alternatives for b and take the maximizer
    # see Mathematica file for details
    # we analyze the remaining tests
    if use_frac:
        return ratio_E_star_naive(n, p, x, c, d, e)
    alt = {0, n - c - d}
    root = ((-4 * e - 4 *n *p - 4 *n**2 *p - 4 *d *x - 4 *d**2 *x)**2 -
            4 *(-2 *x + 2 *c *x - 2 *d *x - 2 *n *x) *(-2 *e - 4 *d *e - 2 *c *n *p -
            2 *d *n *p + 2 *n**2 *p - 2 *c *n**2 *p - 2 *d *n**2 *p + 2 *n**3 *p - 2 *c *d *x -
            2 *d**2 *x - 2 *c *d**2 *x - 2 *d**3 *x + 2 *d *n *x + 2 *d**2 *n *x))
    A = e + n *p + n**2 *p + d *x + d**2 *x
    B = x * (-1 + c - d - n)
    if root > 0:
        alt.add(int(ceil( (A + sqrt(root)/4) / B)))
        alt.add(int(floor((A + sqrt(root)/4) / B)))
        alt.add(int(ceil( (A - sqrt(root)/4) / B)))
        alt.add(int(floor((A - sqrt(root)/4) / B)))
    best_R = 0
    best_b = None
    for b in alt:
        if 0 <= b <= n - c - d:
            ALG = p * n * (n + 1) + 2 * e + x *(n - c) * (n - c + 1) -  x *(n - b - c) *(n - b - c + 1) + x *d *(d + 1)
            OPT = p *n *(n + 1) + x *(b + d) *(b + d + 1)
            if use_frac:
                R = Fraction(ALG, OPT)
            else:
                R = ALG / OPT
            if R > best_R:
                best_R = R
                best_b = b
    assert best_b is not None
    return best_R, best_b

def ratio_E_star_naive(n, p, x, c, d, e):
    """This is just a sanity check for the previous procedure.
    Should give the same result.
    """
    best_R = 0
    best_b = None
    for b in range(n - c - d + 1):
        ALG = p * n * (n + 1) + 2 * e + x *(n - c) * (n - c + 1) -  x *(n - b - c) *(n - b - c + 1) + x *d *(d + 1)
        OPT = p *n *(n + 1) + x *(b + d) *(b + d + 1)
        if use_frac:
            R = Fraction(ALG, OPT)
        else:
            R = ALG / OPT
        if R > best_R:
            best_R = R
            best_b = b
    assert best_b is not None
    return best_R, best_b


def equilibrium_unfiltered(p, x, n, f):
    # computes equlibrium ratio and schedule
    # complexity O(n^2) is in fact O(n^3) because of string concatenation
    # could be made O(n^2) using list appends instead
    if not better:
        answer = algorithm(p, x, "", f)
        if answer[2] == BAD:
            print("Counter example for T*E* conjecture : n=%i p=%f x=%f" % (n, p, x))
            sys.exit(1)
        return answer
    else:
        return compute_ratio_using_grid(p, x, n)

def equilibrium(p, x, n, file=None):
    # global worst_ratio_ratio
    best = equilibrium_unfiltered(p, x, n, file)
    assert best[1] != ""
    if check_pattern:
        # worst_ratio_ratio = max(worst_ratio_ratio, best[0] / optimal_ratio(n, p, x) )
        if pattern.match("<" + best[1] + ">") is None:
            return best
        else:
            return (best[0], "")
    else:
        return best


def compute_grid_R(p, x, n, critical_ratio):
    """ c = #Tp, d = #Tq, e = sum of inverse ranks of tests
        for every cell (c,d) compute maximum stop ratio that 
        can be acheived after c Tp and d Tq such that the 
        intermediate stop ratios are strictly above the
        critical_ratio.
        stop_ratio[c,d] is this ratio
        E[c,d] is the corresponding e value of the path
        direction[c,d] indicates the direction from which the path is 
        reaching (c,d).
        path will be the path the adversary will choose
        and critical ratio the ratio obtained by an optimal algorithm
        against this path.
    """
    if use_frac:
        zero = Fraction(0, 1)
    else:
        zero = 0
    N = range(n + 1)
    stop_ratio = [[ zero  for d in N] for c in N]
    E = [[ 0  for d in N] for c in N]
    critical  = [[ 0  for d in N] for c in N]
    direction = [[' ' for d in N] for c in N]
    for c in N:
        for d in range(n - c + 1):
            if c == d == 0:
                E[c][d] = zero
                stop_ratio[c][d] = ratio_E_star(n, p, x, c, d, E[c][d])[0]
                critical[c][d] = stop_ratio[c][d]
                direction[c][d] = ' '
            else:
                if c > 0 and stop_ratio[c - 1][d] > critical_ratio:
                    eTp = E[c - 1][d] + n - c + 1  # the last Tp was done with value c - 1
                    rTp = ratio_E_star(n, p, x, c, d, eTp)[0]
                    cTp = min(critical[c - 1][d], rTp)
                else:
                    cTp = rTp = 0
                if d > 0 and stop_ratio[c][d - 1] > critical_ratio:
                    eTx = E[c][d - 1] + n - c
                    rTx = ratio_E_star(n, p, x, c, d, eTx)[0]
                    cTx = min(critical[c][d - 1], rTx)
                else:
                    cTx = rTx = 0
                if rTp > rTx and rTp > critical_ratio:
                    stop_ratio[c][d] = rTp
                    E[c][d] = eTp
                    critical[c][d] = cTp
                    direction[c][d] = 'p'
                elif rTx >= rTp and rTx > critical_ratio:
                    stop_ratio[c][d] = rTx
                    E[c][d] = eTx
                    critical[c][d] = cTx
                    direction[c][d] = 'x'
                else:
                    stop_ratio[c][d] = 0
                    E[c][d] = 0
                    critical[c][d] = 0
                    direction[c][d] = ' '
    # check if the adversary wins
    win = None
    for c in N:
        if stop_ratio[c][n-c] and (win is None or critical[win][n - win] < critical[c][n - c]):
            win = c
    if win is None:
        path = set()
        critical_ratio = 0
        critical_cd = (None, None)
    else:
        c = win
        d = n - c 
        path = {(c, d)}
        while (c, d) != (0, 0):
            if direction[c][d] == 'p':
                c -= 1
            else:
                d -= 1
            path.add((c, d))
        critical_ratio, critical_cd = min((stop_ratio[c][d], (c,d)) for (c,d) in path)
    # return stop_ratio, direction, path, critical_ratio
    return stop_ratio, E, direction, path, critical_ratio, critical_cd


def compute_ratio_using_grid(p, x, n):
    global best_iter
    R = 1
    nb_iter = 0
    while R:
        _stop_ratio, _E, _direction, _path, R, _critical_cd = compute_grid_R(p, x, n, R)
        nb_iter += 1
        if R:
            stop_ratio = _stop_ratio
            E = _E
            path = _path 
            ratio = R
            critical_cd = _critical_cd
    schedule = ""
    for c, d in sorted(path):
        if (c,d) != (0,0):
            if c > last_c:
                schedule += "Tp"
            else:
                schedule += "Tx"
        last_c = c
        if (c,d) == critical_cd:
            b = ratio_E_star(n, p, x, c, d, E[c][d])[1]
            schedule += "Ex" * b
            schedule += "Ep" * (n - b - c - d)
            schedule += "x" * d
            break
    assert not use_frac or type(ratio) == Fraction
    # best_iter = max(best_iter, (nb_iter, p, x))
    alt = (nb_iter, p, x)
    if best_iter is None or best_iter < alt:
        best_iter = alt
    if count_iter:
        return ratio, str(nb_iter)
    else:
        return ratio, schedule



def filter_compact(schedule):
    compact = ""
    last = ""
    count = 0
    count_last_x = 0
    if schedule == "" or schedule[0].isdigit():
        return schedule                            # don't filter special cases
    for i, z in enumerate(schedule + "$"):
        if count_last_x:
            count_last_x += 1
        elif i % 2 == 0:
            if z.isupper():
                code = z 
            else:
                compact += " (" + last + ")^" + str(count)
                count_last_x = 1
        else:
            code += z
            if code != last:
                if count:
                    compact += " (" + last + ")^" + str(count)
                last = code
                count = 1
            else:
                count += 1
    count_last_x -= 1 # ignore the terminating dollar sign
    if count_last_x:
        compact += " "+"x^" + str(count_last_x)
    return compact 


def axis_step(min_val, max_val):
    s = 1
    while 10 * s < max_val - min_val:
        s *= 10
    lo = int(min_val / s) * s
    hi = int(max_val / s) * s
    return range(lo, hi + 1, s)


colors = [
"grey",
"red",
"blue",
"cyan",
"magenta",
"orange",
"green",
"yellow",
"lightgreen",
"brown",
"lemonchiffon",
"salmon",
"violet",
"sienna",

"orchid",
"deeppink",
"slategray",
"lightsalmon",
"darkorange",
"beige",
"limegreen",
"yellowgreen",
"cadetblue",
"whitesmoke",
"fuchsia",
"darkkhaki",
"lightgoldenrodyellow",
"gold",
"ivory",
"blanchedalmond",
"darkgoldenrod",
"dimgrey",
"lime",
"darkmagenta",
"lavender",
"aquamarine",
"lightgrey",
"lightskyblue",
"indigo",
"maroon",
"dodgerblue",
"black",
"mediumturquoise",
"thistle",
"seagreen",
"snow",
"olive",
"greenyellow",
"silver",
"teal",
"mediumvioletred",
"olivedrab",
"tan",
"darkseagreen",
"darkblue",
"goldenrod",
"lightcoral",
"darkslategrey",
"mediumslateblue",
"lawngreen",
"chocolate",
"mintcream",
"deepskyblue",
"sandybrown",
"papayawhip",
"lightpink",
"skyblue",
"navajowhite",
"dimgray",
"crimson",
"darkorchid",
"forestgreen",
"orangered",
"darkgreen",
"peru",
"palegoldenrod",
"white",
"lightseagreen",
"antiquewhite",
"slateblue",
"blueviolet",
"seashell",
"moccasin",
"gray",
"bisque",
"darkslategray",
"powderblue",
"mediumseagreen",
"springgreen",
"mistyrose",
"midnightblue",
"mediumspringgreen",
"darksalmon",
"lightslategrey",
"mediumaquamarine",
"chartreuse",
"saddlebrown",
"lightgray",
"palevioletred",
"darkgray",
"indianred",
"azure",
"aqua",
"navy",
"darkturquoise",
"lightslategray",
"grey",
"coral",
"honeydew",
"peachpuff",
"mediumpurple",
"khaki",
"darkred",
"steelblue",
"linen",
"pink",
"rosybrown",
"darkviolet",
"darkgrey",
"lavenderblush",
"gainsboro",
"darkolivegreen",
"lightcyan",
"tomato",
"lightblue",
"darkcyan",
"palegreen",
"firebrick",
"royalblue",
"purple",
"cornsilk",
"mediumorchid",
"lightsteelblue",
"turquoise",
"mediumblue",
"lightyellow",
"slategrey",
"floralwhite",
"oldlace",
"plum",
"burlywood",
"aliceblue",
"cornflowerblue",
"paleturquoise",
"wheat",
"darkslateblue",
"hotpink"]
# random.seed(3)
# random.shuffle(colors)

if __name__ == "__main__":
    if print_tree:
        f = open(("tree_n%i_p%s_x%s.dot" % (n, param_p, param_x)).replace("/","div"), "w")
        print("digraph G {", file=f)
        print(equilibrium(max_p, max_x, n, f))
        print("}", file=f)
        f.close()
    elif print_grid:
        # print(equilibrium(max_p,  max_x, n))
        R = 1
        while R:
            stop_ratio, E, direction, path, R, critical_cd = compute_grid_R(max_p, max_x, n, R)
            print_the_grid(n, stop_ratio, direction, path, critical_cd)
    else:
        # f = open("n%i_p%g_x%g.txt" % (n, max_p, max_x), "w")
        # print("# x p ratio", file=f)
        grid = 128
        code = {}
        M_schedule = []
        M_ratio = []
        range_x = []
        range_p = []
        for ip in range(grid, 0, -1):
            M_schedule.append([])
            M_ratio.append([])
            for ix in range(1, grid + 1):
                if use_frac:
                    p = Fraction(min_p + ip * (max_p - min_p), grid)
                    x = Fraction(min_x + ix * (max_x - min_x), grid)
                else:
                    p = min_p + ip * (max_p - min_p) / grid
                    x = min_x + ix * (max_x - min_x) / grid
                if ix == 1:
                	range_p.append(p)
                if ip == grid:
                	range_x.append(x)
                ratio, schedule = equilibrium(p, x, n)[:2]
                if schedule not in code:
                    code[schedule] = len(code)
                M_schedule[-1].append(code[schedule])
                M_ratio[-1].append(ratio)
                # print(x, p, ratio, file=f)
            # print(file=f)
        # f.close()

        if no_svg:
	        name = ("n%i_p%s_x%s.py" % (n, param_p, param_x)).replace(":", "..").replace("/","div")
        	f = open(name, 'w')
        	print("n=", n, file=f)
        	print("code =", code, file=f)
        	print("M_schedule =", M_schedule, file=f)
        	print("M_ratio =", M_ratio, file=f)
        	print("range_x=", range_x, file=f)
        	print("range_p=", range_p, file=f)
        	f.close()
        else:
	        Z = len(colors)

	        name = ("n%i_p%s_x%s.svg" % (n, param_p, param_x)).replace(":", "..").replace("/","div")
	        dwg = svgwrite.Drawing(filename=name, debug=True)
	        for i, row in enumerate(M_schedule):
	            for j, z in enumerate(row):
	                # dwg.add(dwg.rect(insert=(j*mm, i*mm), size=(1*mm, 1*mm), fill=colors[z % Z], stroke_width=0))
	                dwg.add(dwg.rect(insert=(j*mm, i*mm), size=(1.1*mm, 1.1*mm), fill=colors[z % Z], stroke_width=0))

	        for i, schedule in enumerate(code):
	            z = code[schedule]
	            dwg.add(dwg.rect(insert=((grid +10 )* mm, (i * 5)*mm), size=(4*mm, 4*mm), fill=colors[z % Z], stroke='black', stroke_width=1))
	            dwg.add(dwg.text(filter_compact(schedule), ((grid + 15) * mm, (i * 5 + 4)*mm)))

	        C = 'white'
	        for x in axis_step(min_x, max_x):
	            i = int((x - min_x) / (max_x - min_x) * grid)
	            dwg.add(dwg.line(start=(i*mm, 0*mm), end=(i*mm, grid*mm), stroke=C))
	            dwg.add(dwg.text("%i"%x, (i * mm, (grid+5)*mm)))
	        dwg.add(dwg.text("x", ((grid/2) * mm, (grid+10)*mm)))

	        for p in axis_step(min_p, max_p):
	            j = int(grid - (p - min_p) / (max_p - min_p) * grid)
	            dwg.add(dwg.line(start=(0*mm, j*mm), end=(grid*mm, j*mm), stroke=C))
	            dwg.add(dwg.text("%i"%p, ((grid+1) * mm, j*mm)))
	        dwg.add(dwg.text("p", ((grid+5) * mm, (grid/2)*mm)))

	        dwg.save()

        print(code)
        # if check_pattern:
            # print("worst ratio ratio", worst_ratio_ratio)
    if best_iter is not None:
        print("largest number of iterations found, corresponding p, x = ", best_iter)