def getNSphereRadius( n, k ):
    if real_int( k ) < 3:
        raise ValueError( 'the number of dimensions must be at least 3' )

    if not isinstance( n, RPNMeasurement ):
        return RPNMeasurement( n, 'meter' )

    dimensions = n.getDimensions( )

    if dimensions == { 'length' : 1 }:
        return n
    elif dimensions == { 'length' : int( k - 1 ) }:
        m2 = n.convertValue( RPNMeasurement( 1, [ { 'meter' : int( k - 1 ) } ] ) )

        result = root( fdiv( fmul( m2, gamma( fdiv( k, 2 ) ) ),
                             fmul( 2, power( pi, fdiv( k, 2 ) ) ) ), fsub( k, 1 ) )

        return RPNMeasurement( result, [ { 'meter' : 1 } ] )
    elif dimensions == { 'length' : int( k ) }:
        m3 = n.convertValue( RPNMeasurement( 1, [ { 'meter' : int( k ) } ] ) )

        result = root( fmul( fdiv( gamma( fadd( fdiv( k, 2 ), 1 ) ),
                                   power( pi, fdiv( k, 2 ) ) ),
                             m3 ), k )

        return RPNMeasurement( result, [ { 'meter' : 1 } ] )
    else:
        raise ValueError( 'incompatible measurement type for computing the radius: ' +
                          str( dimensions ) )

def getEclipseTotality( body1, body2, location, date ):
    '''Returns the angular size of an astronomical object in radians.'''
    if isinstance( location, str ):
        location = getLocation( location )

    if not isinstance( body1, RPNAstronomicalObject ) or not isinstance( body2, RPNAstronomicalObject ) and \
       not isinstance( location, RPNLocation ) or not isinstance( date, RPNDateTime ):
        raise ValueError( 'expected two astronomical objects, a location and a date-time' )

    separation = body1.getAngularSeparation( body2, location, date ).value

    radius1 = body1.getAngularSize( ).value
    radius2 = body2.getAngularSize( ).value

    if separation > fadd( radius1, radius2 ):
        return 0

    distance1 = body1.getDistanceFromEarth( date )
    distance2 = body2.getDistanceFromEarth( date )

    area1 = fmul( pi, power( radius1, 2 ) )
    area2 = fmul( pi, power( radius2, 2 ) )

    area_of_intersection = fadd( getCircleIntersectionTerm( radius1, radius2, separation ),
                                 getCircleIntersectionTerm( radius2, radius1, separation ) )

    if distance1 > distance2:
        result = fdiv( area_of_intersection, area1 )
    else:
        result = fdiv( area_of_intersection, area2 )

    if result > 1:
        return 1
    else:
        return result

def getNthOctagonalTriangularNumber( n ):
    sign = power( -1, real( n ) )

    return nint( floor( fdiv( fmul( fsub( 7, fprod( [ 2, sqrt( 6 ), sign ] ) ),
                                    power( fadd( sqrt( 3 ), sqrt( 2 ) ),
                                           fsub( fmul( 4, real_int( n ) ), 2 ) ) ),
                              96 ) ) )

def mp_binom(k,n,p):
	"""
	Binomial function, returning the probability of k successes in n trials given the trial success probability p, and supporting multiprecission output.
	Parameters
	----------
	k : int, ndarray
		Successes.
	n : int, ndarray
		Trials.
	p : float,ndarray
		Trial (experiment) success probability.
	Returns
	-------
	val : float,ndarray
		Probability of k successes in n trials.
	Examples
	--------
	>>> k = 10
	>>> n = 10000
	>>> p = 0.9
	>>> mp_binom(k, n, p)
	9.56548769092821e-9958
	"""
	import mpmath as mp
	val = mp_comb(n,k) * mp.power(p,k) * mp.power(1-p,n-k)
	return val

def calculateWindChill( measurement1, measurement2 ):
    validUnitTypes = [
        [ 'velocity', 'temperature' ],
    ]

    arguments = matchUnitTypes( [ measurement1, measurement2 ], validUnitTypes )

    if not arguments:
        raise ValueError( '\'wind_chill\' requires velocity and temperature measurements' )

    wind_speed = arguments[ 'velocity' ].convert( 'miles/hour' ).value
    temperature = arguments[ 'temperature' ].convert( 'degrees_F' ).value

    if wind_speed < 3:
        raise ValueError( '\'wind_chill\' is not defined for wind speeds less than 3 mph' )

    if temperature > 50:
        raise ValueError( '\'wind_chill\' is not defined for temperatures over 50 degrees fahrenheit' )

    result = fsum( [ 35.74, fmul( temperature, 0.6215 ), fneg( fmul( 35.75, power( wind_speed, 0.16 ) ) ),
                   fprod( [ 0.4275, temperature, power( wind_speed, 0.16 ) ] ) ] )

    # in case someone puts in a silly velocity
    if result < -459.67:
        result = -459.67

    return RPNMeasurement( result, 'degrees_F' ).convert( arguments[ 'temperature' ].units )

def getNthAperyNumber( n ):
    result = 0

    for k in arange( 0, real( n ) + 1 ):
        result = fadd( result, fmul( power( binomial( n, k ), 2 ),
                                     power( binomial( fadd( n, k ), k ), 2 ) ) )

    return result

def getEulerPhi( n ):
    if real( n ) < 2:
        return n

    if g.ecm:
        return reduce( fmul, ( fmul( fsub( i[ 0 ], 1 ), power( i[ 0 ], fsub( i[ 1 ], 1 ) ) ) for i in getECMFactors( n ) ) )
    else:
        return reduce( fmul, ( fmul( fsub( i[ 0 ], 1 ), power( i[ 0 ], fsub( i[ 1 ], 1 ) ) ) for i in getFactors( n ) ) )

 def __call__(self, y):
     """
     This function should be used by the root finder.
     If the brentq root finder is used then the derivative does not help.
     """
     r = self.r2 / self.r1
     lhs = self.a1 - self.b1 * mp.root(y, 2)
     rhs = self.a2 * mp.power(y, r - 1) + self.b2 * mp.power(y, 1.5*r - 1)
     return rhs - lhs

def getNthNonagonalPentagonalNumber( n ):
    sqrt21 = sqrt( 21 )
    sign = power( -1, real_int( n ) )

    return nint( floor( fdiv( fprod( [ fadd( 25, fmul( 4, sqrt21 ) ),
                                       fsub( 5, fmul( sqrt21, sign ) ),
                                       power( fadd( fmul( 2, sqrt( 7 ) ), fmul( 3, sqrt( 3 ) ) ),
                                              fsub( fmul( 4, n ), 4 ) ) ] ),
                              336 ) ) )

def getNthDecagonalHeptagonalNumber( n ):
    sqrt10 = sqrt( 10 )

    return nint( floor( fdiv( fprod( [ fsub( 11,
                                             fmul( fmul( 2, sqrt10 ),
                                                   power( -1, real_int( n ) ) ) ),
                                       fadd( 1, sqrt10 ),
                                       power( fadd( 3, sqrt10 ),
                                              fsub( fmul( 4, n ), 3 ) ) ] ), 320 ) ) )

def get_prob_correct_sync(k, pc):
    """prob = \sum_{l=0}^{ceil(k/2)-1} (k_choose_l) * (pc)^l * (1 - pc)^(k-l)
    """
    retval = 0
    lmax = mpmath.ceil(float(k)/2) - 1
    for l in range(0, lmax + 1):
        retval += mpmath.binomial(k,l) * mpmath.power(pc, l) * \
            mpmath.power(1 - pc, k - l)
    return retval

def get_prob_correct_async(k, pc, rprime):
    """prob = \sum_{l=0}^{k-rprime} (k_choose_l) * (pc)^l * (1 - pc)^(k-l)
    """
    retval = 0
    lmax = k - rprime
    for l in range(0, lmax + 1):
        retval += mpmath.binomial(k,l) * mpmath.power(pc, l) * \
            mpmath.power(1 - pc, k - l)
    return retval

def getNthNonagonalTriangularNumber( n ):
    a = fmul( 3, sqrt( 7 ) )
    b = fadd( 8, a )
    c = fsub( 8, a )

    return nint( fsum( [ fdiv( 5, 14 ),
                         fmul( fdiv( 9, 28 ), fadd( power( b, real_int( n ) ), power( c, n ) ) ),
                         fprod( [ fdiv( 3, 28 ),
                                  sqrt( 7 ),
                                  fsub( power( b, n ), power( c, n ) ) ] ) ] ) )

def fibSum(n):    
    if n == 0 or n == 1:
        return 0
    elif n == 2:
        return 1
    else:
        n = n+1
        a = math.sqrt(5) ; b =1+a ; c = 1-a
        e = (mp.power(b,n)-mp.power(c,n))/(mp.power(2,n)*a)
        return e-1

def getNthSquareTriangularNumber( n ):
    neededPrecision = int( real_int( n ) * 3.5 )  # determined by experimentation

    if mp.dps < neededPrecision:
        setAccuracy( neededPrecision )

    sqrt2 = sqrt( 2 )

    return nint( power( fdiv( fsub( power( fadd( 1, sqrt2 ), fmul( 2, n ) ),
                                           power( fsub( 1, sqrt2 ), fmul( 2, n ) ) ),
                                    fmul( 4, sqrt2 ) ), 2 ) )

    def pdf_eval(self, x):
        if x < 0:
            return 0
        elif x < self.location:
            return 0
        else:
            x = mpf(x)
            a = self.shape/self.scale
            b = (x-self.location)/self.scale
            b = mpmath.power(b, self.shape-1)
            c = mpmath.exp(-mpmath.power(((x-self.location)/self.scale), self.shape))

            return a * b *c

def getNthDecagonalCenteredSquareNumber( n ):
    sqrt10 = sqrt( 10 )

    dps = 7 * int( real_int( n ) )

    if mp.dps < dps:
        mp.dps = dps

    return nint( floor( fsum( [ fdiv( 1, 8 ),
                              fmul( fdiv( 7, 16 ), power( fsub( 721, fmul( 228, sqrt10 ) ), fsub( n, 1 ) ) ),
                              fmul( fmul( fdiv( 1, 8 ), power( fsub( 721, fmul( 228, sqrt10 ) ), fsub( n, 1 ) ) ), sqrt10 ),
                              fmul( fmul( fdiv( 1, 8 ), power( fadd( 721, fmul( 228, sqrt10 ) ), fsub( n, 1 ) ) ), sqrt10 ),
                              fmul( fdiv( 7, 16 ), power( fadd( 721, fmul( 228, sqrt10 ) ), fsub( n, 1 ) ) ) ] ) ) )

def getNthNonagonalSquareNumber( n ):
    if real( n ) < 0:
        ValueError( '' )

    p = fsum( [ fmul( 8, sqrt( 7 ) ), fmul( 9, sqrt( 14 ) ), fmul( -7, sqrt( 2 ) ), -28 ] )
    q = fsum( [ fmul( 7, sqrt( 2 ) ), fmul( 9, sqrt( 14 ) ), fmul( -8, sqrt( 7 ) ), -28 ] )
    sign = power( -1, real_int( n ) )

    index = fdiv( fsub( fmul( fadd( p, fmul( q, sign ) ),
                              power( fadd( fmul( 2, sqrt( 2 ) ), sqrt( 7 ) ), n ) ),
                        fmul( fsub( p, fmul( q, sign ) ),
                              power( fsub( fmul( 2, sqrt( 2 ) ), sqrt( 7 ) ), fsub( n, 1 ) ) ) ), 112 )

    return nint( power( nint( index ), 2 ) )

def interval_bisection(equation, lower, upper) -> str:
    """ Calculate the root of the equation using
    `Interval Bisection` method.

    Examples:
    >>> interval_bisection([1, 0, -3, -4], 2, 3)[0]
    '2.19582335'
    >>> interval_bisection([1, 0, 1, -6], 1, 2)[0]
    '1.63436529'
    >>> interval_bisection([1, 0, 3, -12], 1, 2)[0]
    '1.85888907'
    """

    # Counts the number of iteration
    counter = 0

    equation = mpfiy(equation)
    lower, upper = mpf(lower), mpf(upper)
    index_length = len(equation)
    x = mean(lower, upper)

    previous_alpha = alpha = None
    delta = fsub(upper, lower)

    while delta > power(10, -10) or previous_alpha is None:
        ans = mpf(0)

        # Summing the answer
        for i in range(index_length):
            index_power = index_length - i - 1
            ans = fadd(ans, fmul(equation[i], power(x, index_power)))

        if ans > mpf(0):
            upper = x
        else:
            lower = x

        x = mean(lower, upper)
        previous_alpha, alpha = alpha, x

        if previous_alpha is None:
            continue
        elif previous_alpha == alpha:
            break

        delta = abs(fsub(alpha, previous_alpha))
        counter += 1

    return str(near(alpha, lower, upper)), counter

def get_prob_schedulable(msgs, mi, boot_time, sync = True):
    """Returns the probability that message m is successfully transmitted even
    in the presence of omission, commission, transmission faults. Here, although
    m represents a single message, we compute the probability of successful
    transmission collectively for all replicas of m, identified by a common
    'tid' field. """
    r = msgs.get_replication_factor(mi)
    rprime = int(mpmath.floor(r / 2.0) + 1)

    prob_msg_corrupted = 1 - br.get_prob_poisson(0, mi.deadline, msgs.po)
    prob_msg_omitted = 1 - br.get_prob_poisson(0, boot_time, msgs.po)

    # since pr(correct) is just a function of k <= r and pc,
    # we compute it beforehand for all values (pc is commission fault prob.)
    prob_correct = []
    for k in range(0, r + 1):
        if sync:
            prob_correct.append(get_prob_correct_sync(k, prob_msg_corrupted))
        else:
            prob_correct.append(get_prob_correct_async(k, prob_msg_corrupted, rprime))

    # need to iterate over all subsets of replica set of m
    replica_ids = []
    for mk in msgs:
        if mk.tid == mi.tid:
            replica_ids.append(mk.id)
    assert r == len(replica_ids)

    prob_success = 0
    for omitted_ids in powerset(replica_ids):
        for mk in msgs:
            if mk.id in omitted_ids:
                mk.omitted = True
            else:
                mk.omitted = False

        s = r - len(omitted_ids)
        prob_omitted = mpmath.power(prob_msg_omitted, r - s) * \
            mpmath.power(1 - prob_msg_omitted, s)

        prob_time_correct = 0
        for k in range(1, s + 1):
            prob_time_correct += get_prob_time_periodic(msgs, mi, k) * \
                                 prob_correct[k]

        prob_success += prob_omitted * prob_time_correct

    return min(prob_success, 1)