Вычисление треугольного корня с помощью сложения, вычитания и деления пополам

Выполните поиск пополам, но убедитесь, что y - x всегда является степенью двойки. (Это не увеличивает асимптотическое время выполнения.) Затем T((x + y) / 2) = T(x) + T(h) + x * h, где h — степень двойки, поэтому x * h можно вычислить со сдвигом.

Вот доказательство концепции Python (написано наспех, более или менее неоптимизировано, но позволяет избежать дорогостоящих операций).

def tri(n):
    return ((n * (n + 1)) >> 1)

def triroot(t):
    y = 1
    ty = 1

    # Find a starting point for bisection search by doubling y using
    # the identity T(2*y) = 4*T(y) - y. Stop when T(y) exceeds t.
    # At the end, x = 2*y, tx = T(x), and ty = T(y).
    while (ty <= t):
        assert (ty == tri(y))
        tx = ty
        ty += ty
        ty += ty
        ty -= y
        x = y
        y += y

    # Now do bisection search on the interval [x .. x + h),
    # using these identities:
    # T(x + h) = T(x) + T(h) + x*h
    # T(h/2) = (T(h) + h/2)/4
    th = tx
    h = x
    x_times_h = ((tx + tx) - x)
    while True:
        assert (tx == tri(x))
        assert (x_times_h == (x * h))

        # Divide h by 2
        h >>= 1
        x_times_h >>= 1
        if (not h):
            break
        th += h
        th >>= 1
        th >>= 1

        # Calculate the midpoint of the search interval
        tz = ((tx + th) + x_times_h)
        z = (x + h)
        assert (tz == tri(z))

        # If the midpoint is below the target, move the lower bound
        # of the search interval up to the midpoint
        if (t >= tz):
            tx = tz
            x = z
            x_times_h += ((th + th) - h)
    return x
for q in range(1, 100):
    p = triroot(q)
    assert (tri(p) <= q < tri((p + 1)))
    print(q, p)

Как видно на связанной странице на math.stackexchange.com, существует прямая формула для решения этой проблемы, и если x = n*(n+1)/2, то наоборот:

n = (sqrt(1+8*x) - 1)/2

Теперь есть квадратный корень и другие вещи, но я бы предложил использовать эту прямую формулу с реализацией, подобной следующей:

tmp  = x + x;   '2*x
tmp += tmp;     '4*x
tmp += tmp + 1; '8*x + 1
n = 0;
n2 = 0;
while(n2 <= tmp){
    n2 += n + n + 1; 'remember that (n+1)^2 - n^2 = 2*n + 1
    n++;
}
'here after the loops  n = floor(sqrt(8*x+1)) + 1
n -= 2;         'floor(sqrt(8*x+1)) - 1
n /= 2;         '(floor(sqrt(8*x+1)) - 1) / 2

Конечно, это можно улучшить для повышения производительности, если необходимо, например, учитывая, что целые значения floor(sqrt(8*x+1)) + 1 равны, поэтому n можно увеличивать с шагом 2 (соответственно переписав вычисление n2: n2 += n + n + n + n + 4, которое само по себе может быть записано лучше, чем это).