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 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 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 getNthKFibonacciNumber( n, k ):
    if real( n ) < 0:
        raise ValueError( 'non-negative argument expected' )

    if real( k ) < 2:
        raise ValueError( 'argument <= 2 expected' )

    if n < k - 1:
        return 0

    nth = int( n ) + 4

    precision = int( fdiv( fmul( n, k ), 8 ) )

    if ( mp.dps < precision ):
        mp.dps = precision

    poly = [ 1 ]
    poly.extend( [ -1 ] * int( k ) )

    roots = polyroots( poly )
    nthPoly = getNthFibonacciPolynomial( k )

    result = 0
    exponent = fsum( [ nth, fneg( k ), -2 ] )

    for i in range( 0, int( k ) ):
        result += fdiv( power( roots[ i ], exponent ), polyval( nthPoly, roots[ i ] ) )

    return floor( fadd( re( result ), fdiv( 1, 2 ) ) )

def getInvertedBits( n ):
    value = real_int( n )

    # determine how many groups of bits we will be looking at
    if value == 0:
        groupings = 1
    else:
        groupings = int( fadd( floor( fdiv( ( log( value, 2 ) ), g.bitwiseGroupSize ) ), 1 ) )

    placeValue = mpmathify( 1 << g.bitwiseGroupSize )
    multiplier = mpmathify( 1 )
    remaining = value

    result = mpmathify( 0 )

    for i in range( 0, groupings ):
        # Let's let Python do the actual inverting
        group = fmod( ~int( fmod( remaining, placeValue ) ), placeValue )

        result += fmul( group, multiplier )

        remaining = floor( fdiv( remaining, placeValue ) )
        multiplier = fmul( multiplier, placeValue )

    return result

def getSkyLocation( n, k ):
    if not isinstance( n, ephem.Body ) or not isinstance( k, RPNDateTime ):
        raise ValueError( '\'sky_location\' expects an astronomical object and a date-time' )

    n.compute( k.to( 'utc' ).format( ) )

    return [ fdiv( fmul( mpmathify( n.ra ), 180 ), pi ), fdiv( fmul( mpmathify( n.dec ), 180 ), pi ) ]

def calculateDistance( measurement1, measurement2 ):
    validUnitTypes = [
        [ 'length', 'time' ],
        [ 'velocity', 'time' ],
        [ 'acceleration', 'time' ],
        [ 'jerk', 'time' ],
        [ 'jounce', 'time' ]
    ]

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

    if not arguments:
        raise ValueError( '\'distance\' requires specific measurement types (see help)' )

    time = arguments[ 'time' ]

    if 'length' in arguments:
        distance = arguments[ 'length' ]
    elif 'acceleration' in arguments:
        # acceleration and time
        distance = getProduct( [ fdiv( 1, 2 ), arguments[ 'acceleration' ], time, time ] )
    elif 'jerk' in arguments:
        # jerk and time
        distance = calculateDistance( getProduct( [ fdiv( 1, 2 ), arguments[ 'jerk' ], time ] ), time )
    elif 'jounce' in arguments:
        # jounce and time
        distance = calculateDistance( getProduct( [ fdiv( 1, 2 ), arguments[ 'jounce' ], time ] ), time )
    else:
        # velocity and time
        distance = multiply( arguments[ 'velocity' ], time )

    return distance.convert( 'meter' )

def findCenteredPolygonalNumber( n, k ):
    if real_int( k ) < 3:
        raise ValueError( 'the number of sides of the polygon cannot be less than 3,' )

    s = fdiv( k, 2 )

    return nint( fdiv( fadd( sqrt( s ),
                       sqrt( fsum( [ fmul( 4, real( n ) ), s, -4 ] ) ) ), fmul( 2, sqrt( s ) ) ) )

def getNthMotzkinNumber( n ):
    result = 0

    for j in arange( 0, floor( fdiv( real( n ), 3 ) ) + 1 ):
        result = fadd( result, fprod( [ power( -1, j ), binomial( fadd( n, 1 ), j ),
                                      binomial( fsub( fmul( 2, n ), fmul( 3, j ) ), n ) ] ) )

    return fdiv( result, fadd( n, 1 ) )

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 getAngularSize( self, location=None, date=None ):
        if location and date:
            if isinstance( location, str ):
                location = getLocation( location )

            location.observer.date = date.to( 'utc' ).format( )
            self.object.compute( location.observer )

        # I have no idea why size seems to return the value in arcseconds... that
        # goes against the pyephem documentation that it always uses radians for angles.
        return RPNMeasurement( mpmathify( fdiv( fmul( fdiv( self.object.size, 3600 ), pi ), 180 ) ), 'radian' )

def expandDataUnits( ):
    # expand data measurements for all prefixes
    newConversions = { }

    for dataUnit in dataUnits:
        unitInfo = unitOperators[ dataUnit ]

        for prefix in dataPrefixes:
            newName = prefix[ 0 ] + dataUnit

            # constuct unit operator info
            helpText = '\n\'Using the standard SI prefixes, ' + newName + '\' is the equivalent\nof ' + \
                       '{:,}'.format( 10 ** prefix[ 2 ] ) + ' times the value of \'' + dataUnit + \
                       '\'.\n\nPlease see the help entry for \'' + dataUnit + '\' for more information.'

            if unitInfo.abbrev:
                newAbbrev = prefix[ 0 ] + unitInfo.abbrev
            else:
                newAbbrev = ''

            unitOperators[ newName ] = \
                RPNUnitInfo( unitInfo.unitType, prefix[ 0 ] + unitInfo.plural,
                             newAbbrev, [ ], unitInfo.categories, helpText, True )

            # create new conversions
            newConversion = power( 10, mpmathify( prefix[ 2 ] ) )
            newConversions[ ( newName, dataUnit ) ] = newConversion
            newConversion = fdiv( 1, newConversion )
            newConversions[ ( dataUnit, newName ) ] = newConversion

        for prefix in binaryPrefixes:
            newName = prefix[ 0 ] + dataUnit

            # constuct unit operator info
            helpText = '\n\'Using the binary data size prefixes, ' + newName + '\' is the equivalent\nof 2^' + \
                       str( prefix[ 2 ] ) + ' times the value of \'' + dataUnit + \
                       '\'.\n\nPlease see the help entry for \'' + dataUnit + '\' for more information.'

            if unitInfo.abbrev:
                newAbbrev = prefix[ 0 ] + unitInfo.abbrev
            else:
                newAbbrev = ''

            unitOperators[ newName ] = \
                RPNUnitInfo( unitInfo.unitType, prefix[ 0 ] + unitInfo.plural,
                             newAbbrev, [ ], unitInfo.categories, helpText, True )

            # create new conversions
            newConversion = power( 2, mpmathify( prefix[ 2 ] ) )
            newConversions[ ( newName, dataUnit ) ] = newConversion
            newConversion = fdiv( 1, newConversion )
            newConversions[ ( dataUnit, newName ) ] = newConversion

    return newConversions

def _crt( a, b, m, n ):
    d = getGCD( m, n )

    if fmod( fsub( a, b ), d ) != 0:
        return None

    x = floor( fdiv( m, d ) )
    y = floor( fdiv( n, d ) )
    z = floor( fdiv( fmul( m, n ), d ) )
    p, q, r = getExtendedGCD( x, y )

    return fmod( fadd( fprod( [ b, p, x ] ), fprod( [ a, q, y ] ) ), z )

def getAntiprismSurfaceArea( n, k ):
    if real( n ) < 3:
        raise ValueError( 'the number of sides of the prism cannot be less than 3,' )

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

    if k.getDimensions( ) != { 'length' : 1 }:
        raise ValueError( '\'antiprism_area\' argument 2 must be a length' )

    result = getProduct( [ fdiv( n, 2 ), fadd( cot( fdiv( pi, n ) ), sqrt( 3 ) ), getPower( k, 2 ) ] )
    return result.convert( 'meter^2' )

def getNthStern( n ):
    """Return the nth number of Stern's diatomic series recursively"""
    if real_int( n ) < 0:
        raise ValueError( 'non-negative, real integer expected' )

    if n in [ 0, 1 ]:
        return n
    elif n % 2 == 0: # even
        return getNthStern( floor( fdiv( n, 2 ) ) )
    else:
        return fadd( getNthStern( floor( fdiv( fsub( n, 1 ), 2 ) ) ),
                     getNthStern( floor( fdiv( fadd( n, 1 ), 2 ) ) ) )

def isFriendly( n ):
    first = True

    abundance = 0

    for i in n:
        if first:
            abundance = fdiv( getSigma( i ), i )
            first = False
        elif fdiv( getSigma( i ), i ) != abundance:
            return 0

    return 1

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 getRegularPolygonArea( n, k ):
    if real( n ) < 3:
        raise ValueError( 'the number of sides of the polygon cannot be less than 3,' )

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

    dimensions = k.getDimensions( )

    if dimensions != { 'length' : 1 }:
        raise ValueError( '\'polygon_area\' argument 2 must be a length' )

    return multiply( fdiv( n, fmul( 4, tan( fdiv( pi, n ) ) ) ), getPower( k, 2 ) ).convert( 'meter^2' )

def solveQuadraticPolynomial( a, b, c ):
    if a == 0:
        if b == 0:
            raise ValueError( 'invalid expression, no variable coefficients' )
        else:
            # linear equation, one root
            return [ fdiv( fneg( c ), b ) ]
    else:
        d = sqrt( fsub( power( b, 2 ), fmul( 4, fmul( a, c ) ) ) )

        x1 = fdiv( fadd( fneg( b ), d ), fmul( 2, a ) )
        x2 = fdiv( fsub( fneg( b ), d ), fmul( 2, a ) )

        return [ x1, x2 ]

def getPartitionNumber( n ):
    '''
    This version is, um, less recursive than the original, which I've kept.
    The strategy is to create a list of the smaller partition numbers we need
    to calculate and then start calling them recursively, starting with the
    smallest.  This will minimize the number of recursions necessary, and in
    combination with caching values, will calculate practically any integer
    partition without the risk of a stack overflow.

    I can't help but think this is still grossly inefficient compared to what's
    possible.  It seems that using this algorithm, calculating any integer
    partition ends up necessitating calculating the integer partitions of
    every integer smaller than the original argument.
    '''
    debugPrint( 'partition', int( n ) )

    if real_int( n ) < 0:
        raise ValueError( 'non-negative argument expected' )
    elif n in ( 0, 1 ):
        return 1

    sign = 1
    i = 1
    k = 1

    estimate = log10( fdiv( power( e, fmul( pi, sqrt( fdiv( fmul( 2, n ), 3 ) ) ) ),
                            fprod( [ 4, n, sqrt( 3 ) ] ) ) )
    if mp.dps < estimate + 5:
        mp.dps = estimate + 5

    partitionList = [ ]
    signList = [ ]

    while n - k >= 0:
        partitionList.append( ( fsub( n, k ), sign ) )
        i += 1

        if i % 2:
            sign *= -1

        k = getNthGeneralizedPolygonalNumber( i, 5 )

    partitionList = partitionList[ : : -1 ]

    total = 0

    for partition, sign in partitionList:
        total = fadd( total, fmul( sign, getPartitionNumber( partition ) ) )

    return total