def rgb_distance(r1, g1, b1, r2, g2, b2):
    """Calculates numerical distance between two colors in RGB color space.

    The distance is calculated by CIE94 formula.

    :params: Two colors with ``r, g, b`` values in ``0..1`` range
    :returns: A number in ``0..100`` range.  The lesser - the
              closer colors are.
    """

    # Formulae from wikipedia article re CIE94

    L1, A1, B1 = xyz_to_lab(*rgb_to_xyz(r1, b1, g1))
    L2, A2, B2 = xyz_to_lab(*rgb_to_xyz(r2, b2, g2))

    dL = L1 - L2
    C1 = _sqrt(A1 * A1 + B1 * B1)
    C2 = _sqrt(A2 * A2 + B2 * B2)
    dCab = C1 - C2
    dA = A1 - A2
    dB = B1 - B2

    dEab = _sqrt(dL ** 2 + dA ** 2 + dB ** 2)

    dHab = _sqrt(max(dEab ** 2 - dL ** 2 - dCab ** 2, 0.0))

    dE = _sqrt((dL ** 2) + ((dCab / (1 + 0.045 * C1)) ** 2) + (
        dHab / (1 + 0.015 * C1)) ** 2)
    return dE

def test_isqrt():
    from math import sqrt as _sqrt
    limit = 17984395633462800708566937239551
    assert int(_sqrt(limit)) == integer_nthroot(limit, 2)[0]
    assert int(_sqrt(limit + 1)) != integer_nthroot(limit + 1, 2)[0]
    assert isqrt(limit + 1) == integer_nthroot(limit + 1, 2)[0]
    assert isqrt(limit + 1 + S.Half) == integer_nthroot(limit + 1, 2)[0]

 def vonmisesvariate(self, mu, kappa):
     random = self.random
     if kappa <= 9.9999999999999995e-007:
         return TWOPI * random()
     
     a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa)
     b = (a - _sqrt(2.0 * a)) / 2.0 * kappa
     r = (1.0 + b * b) / 2.0 * b
     while 1:
         u1 = random()
         z = _cos(_pi * u1)
         f = (1.0 + r * z) / (r + z)
         c = kappa * (r - f)
         u2 = random()
         if u2 >= c * (2.0 - c):
             pass
         if not (u2 > c * _exp(1.0 - c)):
             break
         
     u3 = random()
     if u3 > 0.5:
         theta = mu % TWOPI + _acos(f)
     else:
         theta = mu % TWOPI - _acos(f)
     return theta

    def vonmisesvariate(self, mu, kappa):
        """Circular data distribution.
        
        mu is the mean angle, expressed in radians between 0 and 2*pi, and
        kappa is the concentration parameter, which must be greater than or
        equal to zero.  If kappa is equal to zero, this distribution reduces
        to a uniform random angle over the range 0 to 2*pi.
        
        """
        random = self.random
        if kappa <= 1e-06:
            return TWOPI * random()
        a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa)
        b = (a - _sqrt(2.0 * a)) / (2.0 * kappa)
        r = (1.0 + b * b) / (2.0 * b)
        while 1:
            u1 = random()
            z = _cos(_pi * u1)
            f = (1.0 + r * z) / (r + z)
            c = kappa * (r - f)
            u2 = random()
            if u2 < c * (2.0 - c) or u2 <= c * _exp(1.0 - c):
                break

        u3 = random()
        if u3 > 0.5:
            theta = mu % TWOPI + _acos(f)
        else:
            theta = mu % TWOPI - _acos(f)
        return theta

def arcball(x, y):
    h2 = x*x+y*y
    if h2 > 1.:
        h = _sqrt(h2)
        v = x/h, y/h, 0.
    else:
        v = x, y, _sqrt(1.-h2)
    return 0., v

def __parse_word(word):
    """
    Read an string and determines if it is a word or a number
    or a repetitive list os values such as 3*4.5=[4.5 , 4.5, 4.5]

    Args:
        word: An string that should be parsed

    Returns:
       result:
          Value extracted (word, number or list)

       kind:
          'word' such as 'ntime', 'eV', '*1', 'angstrom', etc
          'int' such as 1, 2, 3
          'float' such as 4.5, 6.7, etc
          'list' such as [4.5 , 4.5, 4.5]

    """
    result = None
    kind = None
    if word[0].isalpha() and word[:4] != 'sqrt' and word[:5] != '-sqrt':
        result = word
        kind = 'word'

    elif word[:4] == 'sqrt':
        result = _sqrt(float(word[5:-1]))
        kind = 'float'

    elif word[:5] == '-sqrt':
        result = -_sqrt(float(word[6:-1]))
        kind = 'float'

    elif word[0] == '*':
        result = word
        kind = 'word'

    elif word.isdigit():
        result = int(word)
        kind = "int"

    elif '*' in word:
        splt = word.split('*')

        if splt[0].isdigit():
            mult = int(splt[0])
            number, kind = string2number(splt[1])
            if number is not None:
                result = mult * [number]
                kind = 'list'
            else:
                result = None
                kind = None
    else:
        result, kind = string2number(word)

    return result, kind

    def _facmod(self, n, q):
        res, N = 1, int(_sqrt(n))

        # Exponent of prime p in n! is e_p(n) = [n/p] + [n/p**2] + ...
        # for p > sqrt(n), e_p(n) < sqrt(n), the primes with [n/p] = m,
        # occur consecutively and are grouped together in pw[m] for
        # simultaneous exponentiation at a later stage
        pw = [1]*N

        m = 2 # to initialize the if condition below
        for prime in sieve.primerange(2, n + 1):
            if m > 1:
                m, y = 0, n // prime
                while y:
                    m += y
                    y //= prime
            if m < N:
                pw[m] = pw[m]*prime % q
            else:
                res = res*pow(prime, m, q) % q

        for ex, bs in enumerate(pw):
            if ex == 0 or bs == 1:
                continue
            if bs == 0:
                return 0
            res = res*pow(bs, ex, q) % q

        return res

    def gauss(self, mu, sigma):
        # """Gaussian distribution.
        # mu is the mean, and sigma is the standard deviation.  This is
        # slightly faster than the normalvariate() function.
        # Not thread-safe without a lock around calls.
        # """

        # When x and y are two variables from [0, 1), uniformly
        # distributed, then
        #
        #    cos(2*pi*x)*sqrt(-2*log(1-y))
        #    sin(2*pi*x)*sqrt(-2*log(1-y))
        #
        # are two *independent* variables with normal distribution
        # (mu = 0, sigma = 1).
        # (Lambert Meertens)
        # (corrected version; bug discovered by Mike Miller, fixed by LM)

        # Multithreading note: When two threads call this function
        # simultaneously, it is possible that they will receive the
        # same return value.  The window is very small though.  To
        # avoid this, you have to use a lock around all calls.  (I
        # didn't want to slow this down in the serial case by using a
        # lock here.)

        __random = self.random
        z = self.gauss_next
        self.gauss_next = None
        if z is None:
            x2pi = __random() * TWOPI
            g2rad = _sqrt(-2.0 * _log(1.0 - __random()))
            z = _cos(x2pi) * g2rad
            self.gauss_next = _sin(x2pi) * g2rad

        return mu + z*sigma

    def vonmisesvariate(self, mu, kappa):
        """Circular data distribution.
        
        mu is the mean angle, expressed in radians between 0 and 2*pi, and
        kappa is the concentration parameter, which must be greater than or
        equal to zero.  If kappa is equal to zero, this distribution reduces
        to a uniform random angle over the range 0 to 2*pi.
        
        """
        random = self.random
        if kappa <= 1e-06:
            return TWOPI * random()
        s = 0.5 / kappa
        r = s + _sqrt(1.0 + s * s)
        while 1:
            u1 = random()
            z = _cos(_pi * u1)
            d = z / (r + z)
            u2 = random()
            if u2 < 1.0 - d * d or u2 <= (1.0 - d) * _exp(d):
                break

        q = 1.0 / r
        f = (q + z) / (1.0 + q * z)
        u3 = random()
        if u3 > 0.5:
            theta = (mu + _acos(f)) % TWOPI
        else:
            theta = (mu - _acos(f)) % TWOPI
        return theta

    def _swing(cls, n):
        if n < 33:
            return cls._small_swing[n]
        else:
            N, primes = int(_sqrt(n)), []

            for prime in sieve.primerange(3, N + 1):
                p, q = 1, n

                while True:
                    q //= prime

                    if q > 0:
                        if q & 1 == 1:
                            p *= prime
                    else:
                        break

                if p > 1:
                    primes.append(p)

            for prime in sieve.primerange(N + 1, n//3 + 1):
                if (n // prime) & 1 == 1:
                    primes.append(prime)

            L_product = R_product = 1

            for prime in sieve.primerange(n//2 + 1, n + 1):
                L_product *= prime

            for prime in primes:
                R_product *= prime

            return L_product*R_product

def get_sphere_mapping(x, y, width, height):
    x = min([max([x, 0]), width])
    y = min([max([y, 0]), height])

    sr = _sqrt((width/2)**2 + (height/2)**2)
    sx = ((x - width / 2) / sr)
    sy = ((y - height / 2) / sr)

    sz = 1.0 - sx**2 - sy**2

    if sz > 0.0:
        sz = _sqrt(sz)
        return (sx, sy, sz)
    else:
        sz = 0
        return norm((sx, sy, sz))

    def vonmisesvariate(self, mu, kappa):
        """Circular data distribution.

        mu is the mean angle, expressed in radians between 0 and 2*pi, and
        kappa is the concentration parameter, which must be greater than or
        equal to zero.  If kappa is equal to zero, this distribution reduces
        to a uniform random angle over the range 0 to 2*pi.

        """
        # mu:    mean angle (in radians between 0 and 2*pi)
        # kappa: concentration parameter kappa (>= 0)
        # if kappa = 0 generate uniform random angle

        # Based upon an algorithm published in: Fisher, N.I.,
        # "Statistical Analysis of Circular Data", Cambridge
        # University Press, 1993.

        # Thanks to Magnus Kessler for a correction to the
        # implementation of step 4.

        random = self.random
        if kappa <= 1e-6:
            return TWOPI * random()

        a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa)
        b = (a - _sqrt(2.0 * a)) / (2.0 * kappa)
        r = (1.0 + b * b) / (2.0 * b)

        while 1:
            u1 = random()

            z = _cos(_pi * u1)
            f = (1.0 + r * z) / (r + z)
            c = kappa * (r - f)

            u2 = random()

            if not (u2 >= c * (2.0 - c) and u2 > c * _exp(1.0 - c)):
                break

        u3 = random()
        if u3 > 0.5:
            theta = (mu % TWOPI) + _acos(f)
        else:
            theta = (mu % TWOPI) - _acos(f)

        return theta

    def eval(cls, n, k):
        n, k = map(sympify, (n, k))

        if k.is_Number:
            if k.is_Integer:
                if k < 0:
                    return S.Zero
                elif k == 0 or n == k:
                    return S.One
                elif n.is_Integer and n >= 0:
                    n, k = int(n), int(k)

                    if k > n:
                        return S.Zero
                    elif k > n // 2:
                        k = n - k

                    M, result = int(_sqrt(n)), 1

                    for prime in sieve.primerange(2, n+1):
                        if prime > n - k:
                            result *= prime
                        elif prime > n // 2:
                            continue
                        elif prime > M:
                            if n % prime < k % prime:
                                result *= prime
                        else:
                            N, K = n, k
                            exp = a = 0

                            while N > 0:
                                a = int((N % prime) < (K % prime + a))
                                N, K = N // prime, K // prime
                                exp = a + exp

                            if exp > 0:
                                result *= prime**exp

                    return C.Integer(result)
                else:
                    result = n - k + 1

                    for i in xrange(2, k+1):
                        result *= n-k+i
                        result /= i

                    return result
        elif k.is_negative:
            return S.Zero
        elif (n - k).simplify().is_negative:
            return S.Zero
        else:
            d = n - k

            if d.is_Integer:
                return cls.eval(n, d)

def is_prime3(num):
    '''Tests if a given number is prime. Written with reduce.'''
    if num == 2:
        return True
    elif num % 2 == 0 or num <= 1:
        return False
    root = _ceil(_sqrt(num))
    return _reduce(lambda acc, d: False if not acc or num % d == 0 else True,
                   range(3, root+1, 2), True)

def is_prime2(num):
    '''Tests if a given number is prime. Written with a map.'''
    if num == 2:
        return True
    elif num % 2 == 0 or num <= 1:
        return False
    root = _ceil(_sqrt(num))
    return all(map(lambda div: False if num % div == 0 else True, 
                   range(3, root+1, 2)))

def is_prime(integer):
    """Returns True if ``integer`` is a prime, otherwise False."""
    assert integer < primes[-1] ** 2
    integer = -integer if integer < 0 else integer
    limit = _floor(_sqrt(integer)) + 1
    for i in _takewhile(lambda elem: elem < limit, primes):
        if integer % i == 0:
            return False
    return integer > 1

 def gauss(self, mu, sigma):
     random = self.random
     z = self.gauss_next
     self.gauss_next = None
     if z is None:
         x2pi = random() * TWOPI
         g2rad = _sqrt(-2.0 * _log(1.0 - random()))
         z = _cos(x2pi) * g2rad
         self.gauss_next = _sin(x2pi) * g2rad
     return mu + z * sigma

def sqrt(S):
    """Convenience function for taking square roots of PowerSeries.
    
    This can also replace the ``math.sqrt`` function, extending it to
    take a PowerSeries as an argument.
    """
    from math import sqrt as _sqrt
    if isinstance(S, PowerSeries):
        return S.squareroot()
    return _sqrt(S)

    def gammavariate(self, alpha, beta):
        """Gamma distribution.  Not the gamma function!
        
        Conditions on the parameters are alpha > 0 and beta > 0.
        
        The probability distribution function is:
        
                    x ** (alpha - 1) * math.exp(-x / beta)
          pdf(x) =  --------------------------------------
                      math.gamma(alpha) * beta ** alpha
        
        """
        if alpha <= 0.0 or beta <= 0.0:
            raise ValueError, "gammavariate: alpha and beta must be > 0.0"
        random = self.random
        if alpha > 1.0:
            ainv = _sqrt(2.0 * alpha - 1.0)
            bbb = alpha - LOG4
            ccc = alpha + ainv
            while 1:
                u1 = random()
                if not 1e-07 < u1 < 0.9999999:
                    continue
                u2 = 1.0 - random()
                v = _log(u1 / (1.0 - u1)) / ainv
                x = alpha * _exp(v)
                z = u1 * u1 * u2
                r = bbb + ccc * v - x
                if r + SG_MAGICCONST - 4.5 * z >= 0.0 or r >= _log(z):
                    return x * beta

        else:
            if alpha == 1.0:
                u = random()
                while u <= 1e-07:
                    u = random()

                return -_log(u) * beta
            while 1:
                u = random()
                b = (_e + alpha) / _e
                p = b * u
                if p <= 1.0:
                    x = p ** (1.0 / alpha)
                else:
                    x = -_log((b - p) / alpha)
                u1 = random()
                if p > 1.0:
                    if u1 <= x ** (alpha - 1.0):
                        break
                elif u1 <= _exp(-x):
                    break

            return x * beta

    def _eval(self, n, k):
        # n.is_Number and k.is_Integer and k != 1 and n != k
        from sympy.functions.elementary.exponential import log
        from sympy.core import N

        if k.is_Integer:
            if n.is_Integer and n >= 0:
                n, k = int(n), int(k)

                if k > n:
                    return S.Zero
                elif k > n // 2:
                    k = n - k

                if HAS_GMPY:
                    from sympy.core.compatibility import gmpy
                    return Integer(gmpy.bincoef(n, k))

                prime_count_estimate = N(n / log(n))

                # if the number of primes less than n is less than k, use prime sieve method
                # otherwise it is more memory efficient to compute factorials explicitly
                if prime_count_estimate < k:
                    M, result = int(_sqrt(n)), 1
                    for prime in sieve.primerange(2, n + 1):
                        if prime > n - k:
                            result *= prime
                        elif prime > n // 2:
                            continue
                        elif prime > M:
                            if n % prime < k % prime:
                                result *= prime
                        else:
                            N, K = n, k
                            exp = a = 0

                            while N > 0:
                                a = int((N % prime) < (K % prime + a))
                                N, K = N // prime, K // prime
                                exp = a + exp

                            if exp > 0:
                                result *= prime**exp
                else:
                    result = ff(n, k) / factorial(k)
                return Integer(result)
            else:
                d = result = n - k + 1
                for i in range(2, k + 1):
                    d += 1
                    result *= d
                    result /= i
                return result