#!/usr/bin/ruby
# Solutions to x for:
#    1/x = (k/x)^2 * (k + x^2) - k*x
define i = Num.i
func roots (k) {
    # Formulas from Wolfram|Alpha
    # https://www.wolframalpha.com/input/?i=1%2Fx+%3D+(k%2Fx)%5E2+*+(k%2Bx%5E2)++-+k*x
    var x1 = (
        ((2*k**6 + 27*k**5 - 9*k**3 + 3*sqrt(3)*sqrt(4*k**11 + 27*k**10 - 18*k**8 - k**6 + 4*k**3))**(1/3) / (3*2**(1/3) * k)) -
        ((2**(1/3) * (3*k - k**4)) / (3*(2*k**6 + 27*k**5 - 9*k**3 + 3*sqrt(3)*sqrt(4*k**11 + 27*k**10 - 18*k**8 - k**6 + 4*k**3))**(1/3) * k)) + k/3
    )
    var x2 = (
        -(((1 - i*sqrt(3)) * (2*k**6 + 27*k**5 - 9*k**3 + 3*sqrt(3)*sqrt(4*k**11 + 27*k**10 - 18*k**8 - k**6 + 4*k**3))**(1/3)) / (6*2**(1/3) * k)) +
         (((1 + i*sqrt(3)) * (3*k - k**4)) / (3*2**(2/3) * (2*k**6 + 27*k**5 - 9*k**3 + 3*sqrt(3)*sqrt(4*k**11 + 27*k**10 - 18*k**8 - k**6 + 4*k**3))**(1/3) * k)) + k/3
    )
    (x1, x2, x2.conj);
}
func S (n) {
    1..n -> sum_by {|k|
        var (x1, x2, x3) = roots(k)
        1..n -> sum_by {|p|
            (x1 + x2)**p * (x2 + x3)**p * (x3 + x1)**p
        }
    }
}
func S_int (n) {
    sum_by(2 .. n, {|k|
        var p = (k**2 - 1)
        p * ((-1)**n * p**n - 1) / (p + 1)
    })
}
assert_eq(S(4).round(-40), 51160)
assert_eq(S(5).round(-40), -8385346)
assert_eq(S(4).round(-40), S_int(4))
assert_eq(S(5).round(-40), S_int(5))
say "** Test passed!"