def generateRealArgument( range = [ 0, 10 ], allowNegative = True ):
    factor = 1

    if allowNegative and getRandomInteger( 2 ) == 1:
        factor = -1

    return nstr( fmul( exp( fmul( getRandomNumber( ), random.uniform( *range ) ) ), factor ) )

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

    return nint( fdiv( fsum( [ sqrt( fsum( [ power( k, 2 ), fprod( [ 8, k, real( n ) ] ),
                                             fneg( fmul( 8, k ) ), fneg( fmul( 16, n ) ), 16 ] ) ),
                               k, -4 ] ), fmul( 2, fsub( k, 2 ) ) ) )

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 get_prob_poisson(events, length, rate):
    """ P(k, lambda = t * rate) = """
    avg_events = mpmath.fmul(rate, length) # lambda
    prob = mpmath.fmul((-1), avg_events)
    for i in range(1, events + 1):
        prob = mpmath.fadd(prob, mpmath.log(mpmath.fdiv(avg_events, i)))
    prob = mpmath.exp(prob)
    return prob

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 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 getNthDecagonalNonagonalNumber( n ):
    dps = 8 * int( real_int( n ) )

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

    return nint( floor( fdiv( fmul( fadd( 15, fmul( 2, sqrt( 14 ) ) ),
                                    power( fadd( fmul( 2, sqrt( 2 ) ), sqrt( 7 ) ),
                                           fsub( fmul( 8, n ), 6 ) ) ), 448 ) ) )

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 joinNumber( digits, base ):
    place = 1
    result = 0

    for digit in digits:
        result = fadd( result, fmul( digit, place ) )
        place = fmul( place, base )

    return result

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