# powm1 is based on the Boost implementation
# https://github.com/boostorg/math/blob/242e9d669198415d165bb0213f730c00985bee9a/include/boost/math/special_functions/powm1.hpp
#
# (C) Copyright John Maddock 2006.
#
# Boost Software License - Version 1.0 - August 17th, 2003
#
# Permission is hereby granted, free of charge, to any person or organization
# obtaining a copy of the software and accompanying documentation covered by
# this license (the "Software") to use, reproduce, display, distribute,
# execute, and transmit the Software, and to prepare derivative works of the
# Software, and to permit third-parties to whom the Software is furnished to
# do so, all subject to the following:
#
# The copyright notices in the Software and this entire statement, including
# the above license grant, this restriction and the following disclaimer,
# must be included in all copies of the Software, in whole or in part, and
# all derivative works of the Software, unless such copies or derivative
# works are solely in the form of machine-executable object code generated by
# a source language processor.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
# SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
# FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
# DEALINGS IN THE SOFTWARE.
"""
    powm1(x,y)

Compute `x^y-1` accurately when `x` is near 1. Based on Boost's implementation

"""
function powm1(x::Real,y::Real)
    if x < 0 && isinteger(y/2) # y is an even integer
        x = -x
    end
    if x > 0 && (abs(y * (x - 1)) < 0.5 || abs(y) < 0.2)
        l = y * log(x)
        if l < 0.5
            return expm1(l)
        end
    end
    x^y-1
end


"""
    poscubicroot(p3,p2,p1,0)

Find the smallest positive real root of the polynomial

    p3*x^3 + p2*x^2 + p1*x + p0

"""
function poscubicroot(p3,p2,p1,p0)
    a = p2/p3
    b = p1/p3
    C = p0/p3

    Q = (a^2 - 3*b) / 9
    r = (2*a^3 - 9*a*b + 27*C) / 54

    if r^2 - Q^3 >= 0
        # 1 real root
        A = -copysign(cbrt(abs(r) + sqrt(r^2 - Q^3)), r)
        B = A==0 ? A : Q/A
        return A + B - a/3

    else
        # 3 real roots
        θ = acos(r/(Q*sqrt(Q)))

        r1 = -2 * sqrt(Q) * cos(θ / 3) - a / 3
        r2 = -2 * sqrt(Q) * cos((θ + 2π)/3) - a / 3
        r3 = -2 * sqrt(Q) * cos((θ - 2π)/3) - a / 3

        r = Inf
        if 0 < r1
            r = r1
        end
        if 0 < r2 < r
            r = r2
        end
        if 0 < r3 < r
            r = r3
        end
        if r < Inf
            return r
        else
            throw(DomainError())
        end
    end
end
