Python mpmath.power函数代码示例

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 ) )