Python mpmath.fsub函数代码示例

def getNthNonagonalOctagonalNumber( n ):
    sqrt6 = sqrt( 6 )
    sqrt7 = sqrt( 7 )

    return nint( floor( fdiv( fmul( fsub( fmul( 11, sqrt7 ), fmul( 9, sqrt6 ) ),
                                    power( fadd( sqrt6, sqrt7 ), fsub( fmul( 8, real_int( n ) ), 5 ) ) ),
                              672 ) ) )

def OLDgetPartitionNumber( n ):
    if n < 0:
        return 0

    if n < 2:
        return 1

    result = mpmathify( 0 )

    for k in arange( 1, n + 1 ):
        #n1 = n - k * ( 3 * k - 1 ) / 2
        n1 = fsub( n, fdiv( fmul( k, fsub( fmul( 3, k ), 1 ) ), 2 ) )
        #n2 = n - k * ( 3 * k + 1 ) / 2
        n2 = fsub( n, fdiv( fmul( k, fadd( fmul( 3, k ), 1 ) ), 2 ) )

        result = fadd( result, fmul( power( -1, fadd( k, 1 ) ), fadd( getPartitionNumber( n1 ), getPartitionNumber( n2 ) ) ) )

        if n1 <= 0:
            break

    #old = NOT_QUITE_AS_OLDgetPartitionNumber( n )
    #
    #if ( old != result ):
    #    raise ValueError( "It's broke." )

    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 cont_frac_expansion_sqrt(n):
	"""
	n is NOT square
	e.g. 2 --> (1,2) (2 repeats)
	"""
	if is_square(n):
		return 0
	seq = []
	r = mp.sqrt(n,prec=1000) # DOESNT MATTER?
	a = floor(r)
	fls = [r]
	seq.append(int(a))
	r = mp.fdiv(1.,mp.fsub(r,a,prec=1000),prec=1000)
	a = floor(r)
	fls.append(r)
	seq.append(int(a))
	r = mp.fdiv(1.,mp.fsub(r,a,prec=1000),prec=1000) #THESE TWO MATTER!!!
	a = floor(r)
	fls.append(r)
	seq.append(int(a))
	while not close(r, fls[1]):
		r = mp.fdiv(1.,mp.fsub(r,a,prec=1000),prec=1000) #THESE TWO MATTER!!!
		a = floor(r)
		fls.append(r)
		seq.append(int(a))
	# print seq
	seq.pop()
	return seq

 def _findLeaving(self, enteringPos):
   increase = mp.fsub(mp.inf,'1')
   idx = mp.inf
   pos = -1
   
   # find the variable with the smaller upper bound to the increase in the entering variable
   # if there are multiple choices, set the one with the smaller index on the list of basic indexes
   for i in range(self.m):
     upperBound = self._getUpperBound(self.b[i], self.A[i, enteringPos])
     
     # if this variable impose more constraint on the increase of the entering variable
     # than the current leaving variable, then this is the new leaving variable
     if upperBound <= mp.fsub(increase, self.tolerance):
       idx = self.basicIdx[i]
       pos = i
       increase = upperBound
     # if this variable impose the same constraint on the increase of the entering variable
     # than the current leaving variable  (considering the tolerance) but has a lower index, 
     # then this is the new leaving variable
     elif mp.almosteq(upperBound, increase, self.tolerance) and self.basicIdx[i] < idx:
       idx = self.basicIdx[i]
       pos = i
   
   if pos >= 0:
     return idx, pos
   else:
     return Dictionary.UNBOUNDEDCODE, Dictionary.UNBOUNDED

def getSigma( target ):
    '''
    Returns the sum of the divisors of n, including 1 and n.

    http://math.stackexchange.com/questions/22721/is-there-a-formula-to-calculate-the-sum-of-all-proper-divisors-of-a-number
    '''
    n = floor( target )

    if real( n ) == 0:
        return 0
    elif n == 1:
        return 1

    factors = getECMFactors( n ) if g.ecm else getFactors( n )

    result = 1

    for factor in factors:
        numerator = fsub( power( factor[ 0 ], fadd( factor[ 1 ], 1 ) ), 1 )
        denominator = fsub( factor[ 0 ], 1 )
        #debugPrint( 'sigma', numerator, denominator )
        result = fmul( result, fdiv( numerator, denominator ) )

        if result != floor( result ):
            raise ValueError( 'insufficient precision for \'sigma\', increase precision (-p))' )

    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 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 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 subtract( self, other ):
        if isinstance( other, RPNMeasurement ):
            if self.units == other.units:
                return RPNMeasurement( fsub( self.value, other.value ), self.units )
            else:
                return RPNMeasurement( fsub( self.value, other.convertValue( self ) ), self.units )

        else:
            return RPNMeasurement( fsub( self.value, other ), self.units )

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 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 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 getStandardDeviation( args ):
    if isinstance( args, RPNGenerator ):
        return getStandardDeviation( list( args ) )
    elif isinstance( args[ 0 ], ( list, RPNGenerator ) ):
        return [ getStandardDeviation( arg ) for arg in args ]

    if len( args ) < 2:
        return 0

    mean = fsum( args ) / len( args )

    dev = [ power( fsub( i, mean ), 2 ) for i in args ]
    return sqrt( fdiv( fsum( dev ), fsub( len( dev ), 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 ) )