def tables(f, r, v, xi, symmetric, verbose=False):
	"""Calculate k[i], w[i], and f[i]."""
	ki = [None] * len(xi)
	wi = ki[:]
	fi = ki[:]
	if symmetric:
		im = 2**31
	else:
		im = 2**32
	for i, x in enumerate(xi):
		if verbose and i & 7 == 0:
			print('\r{0}/{1}'.format(i, len(xi)), end='', file=sys.stderr)
		if i == 0:
			ki[0] = mpmath.floor(im * r*f(r)/v)
			wi[0] = v / f(r) / im
		else:
			ki[i] = mpmath.floor(im * xi[i-1]/x)
			wi[i] = x / im
		fi[i] = f(x)
	if verbose:
		print('\r{0}/{0}'.format(len(xi)), file=sys.stderr)
	assert all(v is not None for v in ki)
	assert all(v is not None for v in wi)
	assert all(v is not None for v in fi)
	return ki, wi, fi

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 convertFromEphemDate( ephem_date ):
        dateValues = list( ephem_date.tuple( ) )

        dateValues.append( int( fmul( fsub( dateValues[ 5 ], floor( dateValues[ 5 ] ) ), 1000000 ) ) )
        dateValues[ 5 ] = int( floor( dateValues[ 5 ] ) )

        return RPNDateTime( *dateValues )

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 _fmod(x, y):
    """
    returns x - n * y, where n is the quotient of x / y, rounded
    towards zero to an integer.
    """
    fquot = mpmath.mpf(x) / y
    if fquot < 0:
        n = -mpmath.floor(-fquot)
    else:
        n = mpmath.floor(fquot)
    return x - n * y

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 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 reduce_to_fpp(self,z):
		# TODO: need to check what happens here when periods are infinite.
		from mpmath import floor, matrix, lu_solve
		T1, T2 = self.periods
		R1, R2 = T1.real, T2.real
		I1, I2 = T1.imag, T2.imag
		A = matrix([[R1,R2],[I1,I2]])
		b = matrix([z.real,z.imag])
		x = lu_solve(A,b)
		N = int(floor(x[0]))
		M = int(floor(x[1]))
		alpha = x[0] - N
		beta = x[1] - M
		assert(alpha >= 0 and beta >= 0)
		return alpha,beta,N,M

def getNthLucasNumber( n ):
    if real( n ) == 0:
        return 2
    elif n == 1:
        return 1
    else:
        return floor( fadd( power( phi, n ), 0.5 ) )

    def getRoot( self, operand ):
        if ( floor( operand ) != operand ):
            raise ValueError( 'cannot take a fractional root of a measurement' )

        newUnits = RPNUnits( self.units )

        for unit, exponent in newUnits.items( ):
            if fmod( exponent, operand ) != 0:
                if operand == 2:
                    name = 'square'
                elif operand == 3:
                    name = 'cube'
                else:
                    name = getOrdinalName( operand )

                baseUnits = self.convertToBaseUnits( )

                if ( baseUnits != self ):
                    return baseUnits.getRoot( operand )
                else:
                    raise ValueError( 'cannot take the ' + name + ' root of this measurement: ', self.units )

            newUnits[ unit ] /= operand

        value = root( self.value, operand )

        return RPNMeasurement( value, newUnits ).normalizeUnits( )

 def __call__(self,t):
 
   with mp.extradps(self.prec):
     t = mpf(t)
     if t == 0:
       print "ERROR:   Inverse transform can not be calculated for t=0"
       return ("Error");
           
     N = 2*self.N
     # Initiate the stepsize (mit aktueller Präsision)
     h = 2*pi/N
  
   # The for loop is evaluating the Laplace inversion at each point theta i
   #   which is based on the trapezoidal rule
     ans =  0.0
     for k in range(self.N):
       theta = -pi + (k+0.5)*h
       z = self.shift + N/t*(Talbot.c1*theta/tan(Talbot.c2*theta) - Talbot.c3 + Talbot.c4*theta)
       dz = N/t * (-Talbot.c1*Talbot.c2*theta/sin(Talbot.c2*theta)**2 + Talbot.c1/tan(Talbot.c2*theta)+Talbot.c4)
       v1 = exp(z*t)*dz
       prec = floor(max(log10(abs(v1)),0))
       with mp.extradps(prec):
         value = self.F(z)
       ans += v1*value
           
     return ((h/pi)*ans).imag

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 convertToBaseN( value, base, outputBaseDigits, numerals ):
    if outputBaseDigits:
        if ( base < 2 ):
            raise ValueError( 'base must be greater than 1' )
    else:
        if not ( 2 <= base <= len( numerals ) ):
            raise ValueError( 'base must be from 2 to {0}'.format( len( numerals ) ) )

    if value == 0:
        return 0

    if value < 0:
        return '-' + convertToBaseN( fneg( value ), base, outputBaseDigits, numerals )

    if base == 10:
        return str( value )

    if outputBaseDigits:
        result = [ ]
    else:
        result = ''

    leftDigits = mpmathify( value )

    while leftDigits > 0:
        modulo = fmod( leftDigits, base )

        if outputBaseDigits:
            result.insert( 0, int( modulo ) )
        else:
            result = numerals[ int( modulo ) ] + result

        leftDigits = floor( fdiv( leftDigits, base ) )

    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 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 convertUnitList( self, other ):
        if not isinstance( other, list ):
            raise ValueError( 'convertUnitList expects a list argument' )

        result = [ ]

        nonIntegral = False

        for i in range( 1, len( other ) ):
            conversion = g.unitConversionMatrix[ ( other[ i - 1 ].getUnitName( ), other[ i ].getUnitName( ) ) ]

            if conversion != floor( conversion ):
                nonIntegral = True

        if nonIntegral:
            source = self

            for count, measurement in enumerate( other ):
                with extradps( 2 ):
                    conversion = source.convertValue( measurement )

                    if count < len( other ) - 1:
                        result.append( RPNMeasurement( floor( conversion ), measurement.units ) )
                        source = RPNMeasurement( chop( frac( conversion ) ), measurement.units )
                    else:
                        result.append( RPNMeasurement( conversion, measurement.units ) )

            return result
        else:
            source = self.convert( other[ -2 ] )

            with extradps( 2 ):
                result.append( source.getModulo( other[ -2 ] ).convert( other[ -1 ] ) )

            source = source.subtract( result[ -1 ] )

            for i in range( len( other ) - 2, 0, -1 ):
                source = source.convert( other[ i - 1 ] )

                with extradps( 2 ):
                    result.append( source.getModulo( other[ i - 1 ] ).convert( other[ i ] ) )

                source = source.subtract( result[ -1 ] )

            result.append( source )

            return result[ : : -1 ]

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 real_int( n ):
    if im( n ) != 0:
        raise ValueError( 'real argument expected ({})'.format( n ) )

    if n != floor( n ):
        raise ValueError( 'integer argument expected ({})'.format( n ) )

    return int( n )

def duplicateOperation( valueList ):
    if g.duplicateOperations > 0:
        raise ValueError( "'dupop' must be followed by another operation" )

    if isinstance( valueList[ -1 ], list ):
        raise ValueError( "'dupop' cannot accept a list argument" )

    g.duplicateOperations = nint( floor( valueList.pop( ) ) )

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