sidef-scripts/Math/power_integers.sf

#!/usr/bin/ruby

# Author: Trizen
# Date: 10 May 2022
# https://github.com/trizen

# Implementation of "PowerInteger" numbers: a + b*w^e

class PowerInteger(a, b, w, e=1) {     # represents = a + b*w^e

    method to_s {
        "PowerInteger(#{a}, #{b}, #{w}, #{e.as_rat})"
    }

    method stringify {
        return a.stringify if (b == 0)
        var im = "(#{b.stringify})*(#{w.stringify})^(#{e.as_rat})"
        return im if (a == 0)
        "#{a.stringify} + #{im}"
    }

    method to_n {
        var v = (a + b*w**e)
        if (v.class != :Number) {
            return v.to_n
        }
        return v
    }

    method neg {
        __CLASS__(-a, -b, w, e)
    }

    method ==(Number n) {
        self.to_n == n
    }

    method ==(__CLASS__ z) {
        (a == z.a) && (b == z.b) && (w == z.w) && (e == z.e)
    }

    method +(Number n) {
        __CLASS__(a + n, b, w, e)
    }

    method +(__CLASS__ z) {
        var (x,y) = (z.a, z.b)

        if ((w == z.w) && (e == z.e)) {
            return __CLASS__(a+x, b+y, w, e)
        }

        __CLASS__(z + a, b, w, e)
    }

    __CLASS__.alias_method('+', 'add')

    method -(Number n) {
        self + -n
    }

    method -(__CLASS__ z) {
        self + -z
    }

    __CLASS__.alias_method('-', 'sub')

    method *(Number n) {
        __CLASS__(a*n, b*n, w, e)
    }

    method *(__CLASS__ z) {

        var (x,y) = (z.a, z.b)

        # (a + b*w^e) * (x + y*w^e) = (a*y + b*x)*w^e + a*x + b*y*w^(2*e)
        if ((z.w == w) && (z.e == e)) {
            return (__CLASS__(0, a*y + b*x, w, e) + __CLASS__(a*x, b*y, w, 2*e))
        }

        # (a + b*w^e) * (x + y*w^r) = y w^r (a + b w^e) + x (a + b w^e)
        # (a + b*w^e) * (x + y*w^r) = a y w^r + a x + b y w^(r + e) + b w^e x
        if (z.w == w) {
            var r = z.e
            return (__CLASS__(0, a*y, w, r) + __CLASS__(a*x, b*y, w, r+e) + __CLASS__(0, b*x, w, e))
        }

        __CLASS__(z * a, z * b, w, e)
    }

    __CLASS__.alias_method('*', 'mul')

    method inv {
        Fraction(1, self)
    }

    method **(Number n) {

        if (!n.is_int) {
            return self.to_n**n
        }

        if (n < 0) {
            return self**n.abs
        }

        var c = __CLASS__(1, 0, w, e)
        var x = self

        for bit in (n.as_bin.chars.flip) {
             c *= x if bit
             x *= x
        }

        return c
    }

    __CLASS__.alias_method('**', 'pow')

    method /(Number n) {
        __CLASS__(a/n, b/n, w, e)
    }

    method /(__CLASS__ z) {
        self * z.inv
    }
}

class Number {
    method *(PowerInteger z) {
        z * self
    }

    __CLASS__.alias_method('*', 'mul')

    method +(PowerInteger z) {
        z + self
    }

    __CLASS__.alias_method('+', 'add')

    method -(PowerInteger z) {
        PowerInteger(self, 0, z.w, z.e) - z
    }

    __CLASS__.alias_method('-', 'sub')

    method /(PowerInteger z) {
        Fraction(self, z)
    }

    __CLASS__.alias_method('/', 'div')
}

func sqrtP(n) {
    PowerInteger(0, 1, n, 1/2)
}

func cbrtP(n) {
    PowerInteger(0, 1, n, 1/3)
}

func solve_cubic_equation(a,b,c,d) {

    var D0 = (b*b - 3*a*c)
    var D1 = (2*b**3 - 9*a*b*c + 27*a*a*d)

    var roots = []
    var z = (-1 + sqrtP(-3))/2

    var C = cbrtP((D1 - (sgn(D0)||-1)*sqrtP(D1*D1 - 4*D0**3))/2)

    0..2 -> each {|k|
        var t = (C * z**k)
        var x = -((b + t + D0/t))/(3*a)
        roots << x
    }

    return roots
}

say ":: Solutions to: x^3 + 5*x^2 + 2*x - 8 = 0"
solve_cubic_equation(1, 5, 2, -8).each_kv{|k,v| say ("x_#{k+1} = ", v.stringify) }

say "\n:: Solutions to: x^3 + 4*x^2 + 7*x + 6 = 0"
solve_cubic_equation(1, 4, 7, 6).each_kv{|k,v| say ("x_#{k+1} = ", v.stringify) }

say "\n:: Solutions to: -36*x^3 + 8*x^2 - 82*x + 2850986 = 0:"
solve_cubic_equation(-36, 8, -82, 2850986).each_kv{|k,v| say ("x_#{k+1} = ", v.stringify) }

say "\n:: Solutions to: 15*x^3 - 22*x^2 + 8*x - 7520940423059310542039581 = 0:"
solve_cubic_equation(15, -22, 8, -7520940423059310542039581).each_kv{|k,v| say ("x_#{k+1} = ", v.stringify) }

# Run some tests

assert_eq(
    PowerInteger(3, 4, 5, 1/3)*PowerInteger(12, 4, 5, 1/3),
    ((3 + 4*cbrt(5)) * (12 + 4*cbrt(5)))
)

assert_eq(
    PowerInteger(3, 4, 5, 6)*PowerInteger(12, 4, 5, 2),
    ((3 + 4*5**6) * (12 + 4*5**2))
)

assert_eq(
    PowerInteger(3, 4, 5, 1/3)*PowerInteger(12, 4, 13, 1/3),
    ((3 + 4*cbrt(5)) * (12 + 4*cbrt(13)))
)

assert_eq(
    PowerInteger(3, 4, 5, 1/3)*PowerInteger(12, 4, 13, 2),
    ((3 + 4*cbrt(5)) * (12 + 4*13**2))
)

__END__
:: Solutions to: x^3 + 5*x^2 + 2*x - 8 = 0
x_1 = (-19 + (-5)*(-28 + (-1/2)*(-24300)^(1/2))^(1/3))/((3)*(-28 + (-1/2)*(-24300)^(1/2))^(1/3)) + (-1/3)*(-28 + (-1/2)*(-24300)^(1/2))^(1/3)
x_2 = (-19 + (5/2 + (-5/2)*(-3)^(1/2))*(-28 + (-1/2)*(-24300)^(1/2))^(1/3))/((-3/2 + (3/2)*(-3)^(1/2))*(-28 + (-1/2)*(-24300)^(1/2))^(1/3)) + (1/6 + (-1/6)*(-3)^(1/2))*(-28 + (-1/2)*(-24300)^(1/2))^(1/3)
x_3 = (-19 + (-5/4 + (-5/4)*(-3)^(1) + (5/2)*(-3)^(1/2))*(-28 + (-1/2)*(-24300)^(1/2))^(1/3))/((3/4 + (3/4)*(-3)^(1) + (-3/2)*(-3)^(1/2))*(-28 + (-1/2)*(-24300)^(1/2))^(1/3)) + (-1/12 + (-1/12)*(-3)^(1) + (1/6)*(-3)^(1/2))*(-28 + (-1/2)*(-24300)^(1/2))^(1/3)

:: Solutions to: x^3 + 4*x^2 + 7*x + 6 = 0
x_1 = (5 + (-4)*(19 + (1/2)*(1944)^(1/2))^(1/3))/((3)*(19 + (1/2)*(1944)^(1/2))^(1/3)) + (-1/3)*(19 + (1/2)*(1944)^(1/2))^(1/3)
x_2 = (5 + (2 + (-2)*(-3)^(1/2))*(19 + (1/2)*(1944)^(1/2))^(1/3))/((-3/2 + (3/2)*(-3)^(1/2))*(19 + (1/2)*(1944)^(1/2))^(1/3)) + (1/6 + (-1/6)*(-3)^(1/2))*(19 + (1/2)*(1944)^(1/2))^(1/3)
x_3 = (5 + (-1 + (-1)*(-3)^(1) + (2)*(-3)^(1/2))*(19 + (1/2)*(1944)^(1/2))^(1/3))/((3/4 + (3/4)*(-3)^(1) + (-3/2)*(-3)^(1/2))*(19 + (1/2)*(1944)^(1/2))^(1/3)) + (-1/12 + (-1/12)*(-3)^(1) + (1/6)*(-3)^(1/2))*(19 + (1/2)*(1944)^(1/2))^(1/3)

:: Solutions to: -36*x^3 + 8*x^2 - 82*x + 2850986 = 0:
x_1 = (8792 + (-8)*(49880745296 + (1/2)*(9952355007856165026816)^(1/2))^(1/3))/((-108)*(49880745296 + (1/2)*(9952355007856165026816)^(1/2))^(1/3)) + (1/108)*(49880745296 + (1/2)*(9952355007856165026816)^(1/2))^(1/3)
x_2 = (8792 + (4 + (-4)*(-3)^(1/2))*(49880745296 + (1/2)*(9952355007856165026816)^(1/2))^(1/3))/((54 + (-54)*(-3)^(1/2))*(49880745296 + (1/2)*(9952355007856165026816)^(1/2))^(1/3)) + (-1/216 + (1/216)*(-3)^(1/2))*(49880745296 + (1/2)*(9952355007856165026816)^(1/2))^(1/3)
x_3 = (8792 + (-2 + (-2)*(-3)^(1) + (4)*(-3)^(1/2))*(49880745296 + (1/2)*(9952355007856165026816)^(1/2))^(1/3))/((-27 + (-27)*(-3)^(1) + (54)*(-3)^(1/2))*(49880745296 + (1/2)*(9952355007856165026816)^(1/2))^(1/3)) + (1/432 + (1/432)*(-3)^(1) + (-1/216)*(-3)^(1/2))*(49880745296 + (1/2)*(9952355007856165026816)^(1/2))^(1/3)

:: Solutions to: 15*x^3 - 22*x^2 + 8*x - 7520940423059310542039581 = 0:
x_1 = (-124 + (22)*(-45689713070085311542890452111/2 + (-1/2)*(2087549880426724544732454788847681930990126035285976729825)^(1/2))^(1/3))/((45)*(-45689713070085311542890452111/2 + (-1/2)*(2087549880426724544732454788847681930990126035285976729825)^(1/2))^(1/3)) + (-1/45)*(-45689713070085311542890452111/2 + (-1/2)*(2087549880426724544732454788847681930990126035285976729825)^(1/2))^(1/3)
x_2 = (-124 + (-11 + (11)*(-3)^(1/2))*(-45689713070085311542890452111/2 + (-1/2)*(2087549880426724544732454788847681930990126035285976729825)^(1/2))^(1/3))/((-45/2 + (45/2)*(-3)^(1/2))*(-45689713070085311542890452111/2 + (-1/2)*(2087549880426724544732454788847681930990126035285976729825)^(1/2))^(1/3)) + (1/90 + (-1/90)*(-3)^(1/2))*(-45689713070085311542890452111/2 + (-1/2)*(2087549880426724544732454788847681930990126035285976729825)^(1/2))^(1/3)
x_3 = (-124 + (11/2 + (11/2)*(-3)^(1) + (-11)*(-3)^(1/2))*(-45689713070085311542890452111/2 + (-1/2)*(2087549880426724544732454788847681930990126035285976729825)^(1/2))^(1/3))/((45/4 + (45/4)*(-3)^(1) + (-45/2)*(-3)^(1/2))*(-45689713070085311542890452111/2 + (-1/2)*(2087549880426724544732454788847681930990126035285976729825)^(1/2))^(1/3)) + (-1/180 + (-1/180)*(-3)^(1) + (1/90)*(-3)^(1/2))*(-45689713070085311542890452111/2 + (-1/2)*(2087549880426724544732454788847681930990126035285976729825)^(1/2))^(1/3)