from math import *
from numpy import *

f = lambda x1, x2, x3: - x1 * x2 * x3

# >= 0
ge = [
    lambda x1, x2, x3: x1,
    lambda x1, x2, x3: x2,
    lambda x1, x2, x3: x3,
]

# <= 0
le = [
    lambda x1, x2, x3: x1 - 42,
    lambda x1, x2, x3: x2 - 42,
    lambda x1, x2, x3: x3 - 42,
    lambda x1, x2, x3: x1 + 2 * x2 + 2 * x3 - 72
]

# == 0
eq = [

]

eps_barrier = 1e-2

x0 = array([1.0, 1.0, 1.0])
eps = 1e-4
eps_min = 1e-7
dx = 1e-7

###
r0 = 10.0
z = 4.0
phi = lambda x: min(0.0, x) ** 2  # -min(0.0, x)
phi_eq = lambda x: x ** 2  # abs(x)


def min_range(f, x_0, h):
    x_k, f_k, a, b = [0] * 4
    k = 2

    a, x_k, f_k = x_0, x_0 + h, f(x_0 + h)

    while True:
        x_k_prev, f_k_prev = x_k, f_k
        x_k = x_0 + 2 ** (k - 1) * h
        f_k = f(x_k)
        if isnan(f_k):
            return a, x_k_prev
        if f_k_prev <= f_k:
            return a, x_k
        else:
            a = x_k_prev
            k += 1


def dih(f, a, b, eps, delta):
    while (b - a) / 2.0 >= eps:
        x_1 = (a + b) / 2.0 - delta
        x_2 = (a + b) / 2.0 + delta
        if f(x_1) <= f(x_2):
            b = x_2
        else:
            a = x_1

    return (a + b) / 2.0

def grad(f, x):
    res = zeros(x.size)
    for i in range(x.size):
        t = x[:]
        t[i] += dx
        ft = f(t)
        t[i] -= 2 * dx
        res[i] = (ft - f(t)) / (2 * dx)

    return res


def m3(f, xk, H):
    k = 0
    g = None
    beta = 0.0
    while True:
        print("\nшаг 1: k =", k)
        gm1 = g
        g = grad(f, xk)
        norm = sqrt(g.dot(g))
        print("\tшаг 2: градиент =", g, ", |g| = ", norm)
        if norm <= eps:
            print("\tшаг 3: ||grad|| < eps, {", norm, "<", eps, "} остановка")
            print("\tшаг 12:\n\tx* =", xk)
            print("\tf(x*) =", f(xk))
            return xk
        if k != 0:
            sk = xk - xkm1
            print("\tшаг 5: sk = ", sk)
            yk = g - gm1
            print("\tшаг 6: yk = ", yk)
            H = H - H.dot(outer(yk, yk)).dot(H) / yk.dot(H).dot(yk) + outer(sk, sk) / yk.dot(sk)
            print("\tшаг 7: H = \n", H)
        p = -H.dot(g)
        print("\tшаг 8: p = ", p)
        tf = lambda alpha: f((xk + alpha * p))
        a, b = min_range(tf, 0.0, eps_min)
        print("\tшаг 9: одномерная минимизация:\n\t\tотрезок локализации минимума: [", a, ';', b, ']')
        alpha = dih(tf, a, b, eps_min, eps_min / 10.0)
        print("\t\tak =", alpha)
        xkm1 = xk
        xk = xk + alpha * p
        print("\tшаг 10: x(k+1) = xk + ak * g = ", xk)
        print("\tf(xk) = ", f(xk))

        # костыль
        if alpha <= eps_min:
            break

        k += 1
        print("\tшаг 11: k = k + 1 = ", k)

    return xk


# m3(x0, H)

P = lambda x, rk: rk * (
    sum(
        list(map(lambda t: phi(t(*x)), ge))
    ) +
    sum(
        list(map(lambda t: phi(-t(*x)), le))
    ) +
    sum(
        list(map(lambda t: phi_eq(t(*x)), eq))
    ))


def alg(xk, r, z):
    k = 1
    while True:
        F = lambda x: f(*x) + P(x, r ** k)
        H = identity(xk.size)
        xrk = m3(F, xk.copy(), H)

        p = P(xrk, r ** k)

        print("P: ", p)

        if p < eps_barrier:
            return xrk
        else:
            r *= z
            xk = xrk
            k += 1


x = alg(x0, r0, z)
print(*list(map(lambda x: "%.8f" % x, x)))
print(f(*x))