def test():
    # Basic test of the Weibull functions
    w = Weibull(shape=mpf(2.0), scale=mpf(12), location=6)

    print "Weibull(%s,%s,%s): " % (w.shape, w.scale, w.location)

    sum = 0

    for i in range(10000):
        sum += w.draw_truncated(6)

    print "MEAN: ", mpf(sum) / 10000.
    
    print w.draw_inverse_transform(0)
    print w.draw_inverse_transform(0)
    print w.draw_inverse_transform(0)
    print w.draw_inverse_transform(0)

    print "Max hazard rate is %e\n" % w.get_max_hazard_rate(100)
    
    for i in range(0,200,5):
        print "CDF at time %d is %f\n" % (i, w.cdf_eval(i))
    
    w = Weibull(shape=mpf(1.0), scale=mpf(120000))

    print "Bunch of draws:"
    for i in range(10):
        print w.draw_inverse_transform(1000000)

    print "Weibull(%s,%s,%s): " % (w.shape, w.scale, w.location)

    print "Max hazard rate is %e\n" % w.get_max_hazard_rate(1000)
    
    for i in range(0,1000,100):
        print "Hazard rate at time %d is %e\n" % (i, w.hazard_rate(i))

def test_mpmath_lambda():
    mpmath.mp.dps = 50
    sin02 = mpmath.mpf("0.19866933079506121545941262711838975037020672954020")
    f = lambdify(x, sin(x), "mpmath")
    prec = 1e-49  # mpmath precision is around 50 decimal places
    assert -prec < f(mpmath.mpf("0.2")) - sin02 < prec
    raises(TypeError, lambda: f(x))

def make_gamma_vals():
    from mpmath import gamma

    x = [-mpf('0.5'), -mpf('0.01'), mpf('0.01')] + linspace(mpf('0.1'), 10, 100) + [mpf(20), mpf(30)]
    ga = [gamma(val) for val in x]

    return make_special_vals('gamma_vals', ('x', x), ('ga', ga))

	def Pprime(self,z):
		# A+S 18.9.
		from mpmath import ellipfun, sqrt, cos, sin, mpc, mpf
		Delta = self.Delta
		e1, e2, e3 = self.__roots
		if self.__ng3:
			z = mpc(0,1) * z
		if Delta > 0:
			zs = sqrt(e1 - e3) * z
			m = (e2 - e3) / (e1 - e3)
			retval = -2 * sqrt((e1 - e3)**3) * ellipfun('cn',u=zs,m=m) * ellipfun('dn',u=zs,m=m) / (ellipfun('sn',u=zs,m=m)**3)
		elif Delta < 0:
			H2 = (sqrt((e2 - e3) * (e2 - e1))).real
			assert(H2 > 0)
			m = mpf(1) / mpf(2) - 3 * e2 / (4 * H2)
			zp = 2 * z * sqrt(H2)
			retval = -4 * sqrt(H2**3) * ellipfun('sn',u=zp,m=m) * ellipfun('dn',u=zp,m=m) / ((1 - ellipfun('cn',u=zp,m=m))**2)
		else:
			g2, g3 = self.__invariants
			if g2 == 0 and g3 == 0:
				retval = -2 / (z**3)
			else:
				c = e1 / 2
				A = sqrt(3 * c)
				retval = -6 * c * A * cos(A * z) / (sin(A * z))**3
		if self.__ng3:
			return mpc(0,-1) * retval
		else:
			return retval

def findPoissonChangePoint( data, factorial ):
	# data is a list of counts in each time period, uniformly spaced

	# the denominator (including both P(D|H1) and constant parts of P(D|H2) )
	C = data.sum()
	N = mpf(len(data))
	denominator = factorial[C-1] * pi / ( 2 * N**C )

	# the numerator (trickier)
	# this needs to be averaged over the possible change points 
	weights = zeros(N,dtype=object)
	CA = 0
	CB = C
	for i in range(1,N) :
		# points up through i are in data set A; the rest are in B
		datapoint = data[i-1]	
		NA = mpf(i)   ; CA += datapoint
		NB = mpf(N-i) ; CB -= datapoint
		
		fraction_num = factorial[CA] * factorial[CB] 
		fraction_den = NA**(CA+1) * NB**(CB+1) * ( (CA/NA)**2 + (CB/NB)**2 )
		#weights.append( fraction_num/fraction_den )
		weights[i-1] = mpf(fraction_num)/fraction_den

	numerator = weights.mean()
	lognum= inv_log10 * log( numerator )
	logden= inv_log10 * log( denominator )
	logodds = lognum - logden
	print "num:",numerator, "log num:", lognum, "| denom:", denominator, "log denom:", logden, "|| log odds:", logodds 

	# If there is a change point, then logodds will be greater than 0
	if logodds < 0 : return None
	return ( weights.argmax(), logodds ) 

def _f_cdf(dfn, dfd, x):
    if x < 0:
        return mpmath.mpf(0)
    dfn, dfd, x = mpmath.mpf(dfn), mpmath.mpf(dfd), mpmath.mpf(x)
    ub = dfn*x/(dfn*x + dfd)
    res = mpmath.betainc(dfn/2, dfd/2, x2=ub, regularized=True)
    return res

def lagrange_inversion(a):
    """Given a series

    f(x) = a[1]*x + a[2]*x**2 + ... + a[n-1]*x**(n - 1),

    use the Lagrange inversion formula to compute a series

    g(x) = b[1]*x + b[2]*x**2 + ... + b[n-1]*x**(n - 1)

    so that f(g(x)) = g(f(x)) = x mod x**n. We must have a[0] = 0, so
    necessarily b[0] = 0 too.

    The algorithm is naive and could be improved, but speed isn't an
    issue here and it's easy to read.

    """
    n = len(a)
    f = sum(a[i]*x**i for i in range(len(a)))
    h = (x/f).series(x, 0, n).removeO()
    hpower = [h**0]
    for k in range(n):
        hpower.append((hpower[-1]*h).expand())
    b = [mp.mpf(0)]
    for k in range(1, n):
        b.append(hpower[k].coeff(x, k - 1)/k)
    b = map(lambda x: mp.mpf(x), b)
    return b

def testing_kbes(Rt,Xt):
    [R0,R1,NR]=Rt
    [X0,X1,NX]=Xt
    NRr=mpmath.mpf(NR)
    NXr=mpmath.mpf(NX)
    for j in range(1,NR):
        rj=mpmath.mpf(j)
        R=R0+R1*rj/NRr
        iR=mpmath.mpc(0,R)
        for k in range(1,NX):
            rk=mpmath.mpf(k)
            x=X0+X1*rk/NXr
            print "r,x=",R,x
            if(x>R):
                print "kbes_asymp="
                timeit( "kbes_asymp(R,x)",repeat=1)
            else:
                print "kbes_rec="
                timeit( "kbes_rec(R,x)",repeat=1)
            print "mpmath.besselk="
            timeit("mpmath.besselk(iR,x)",repeat=1)
            

            #print "t1(",R,x,")=",t1
            #print "t2(",R,x,")=",t2
            if(R<15.0):
                if(x<0.3 *R):
                    print "Case 1"
                elif(x<=max(10.0 +1.2*R,2 *R)):
                    print "Case 2"
            elif(R>20  and x>4 *R):
                print "Case 3"
            else:
                print "Case 4"

def dedekind(tau, floatpre):
    """
    Algorithm 22 (Dedekind eta)
    Input : tau in the upper half-plane, k in N
    Output : eta(tau)
    """
    a = 2 * mpmath.pi / mpmath.mpf(24)
    b = mpmath.exp(mpmath.mpc(0, a))
    p = 1
    m = 0
    while m <= 0.999:
        n = nearest_integer(tau.real)
        if n != 0:
            tau -= n
            p *= b ** n
        m = tau.real * tau.real + tau.imag * tau.imag
        if m <= 0.999:
            ro = mpmath.sqrt(mpmath.power(tau, -1) * 1j)
            if ro.real < 0:
                ro = -ro
            p = p * ro
            tau = (-p.real + p.imag * 1j) / m
    q1 = mpmath.exp(a * tau * 1j)
    q = q1 ** 24
    s = 1
    qs = mpmath.mpc(1, 0)
    qn = 1
    des = mpmath.mpf(10) ** (-floatpre)
    while abs(qs) > des:
        t = -q * qn * qn * qs
        qn = qn * q
        qs = qn * t
        s += t + qs
    return p * q1 * s

def findGaussianChangePoint( data ):
	
	# the denominator. This is the easy part.
	N = len( data )

	if N<6 : return None # can't find a cp in data this small

	# set up gamma function table
	#for i in range(N):
		

	s2 = mpf(data.var())
	gpart = gamma( mpf(N)/2.0 - 1 )
	denom = (pi**1.5) * mpf((N*s2))**( -N/2.0 + 0.5 ) * gpart

	# the numerator. A little trickier.
	# calc_twostate_weights() already deals with ts<3 and ts>N-2.
	weights=calc_twostate_weights( data )
	if weights is None: return None

	num = 2.0**2.5 * abs(data.mean()) * weights.mean()

	logodds = log( num ) - log( denom ) 	

	print "num:", num, "log num:", log(num), "| denom:", denom, "log denom:", log(denom), "|| log odds:", logodds 
	
	# If there is a change point, then logodds will be greater than 0
	if logodds < 0 : 
		return None
	
	return ( weights.argmax(), logodds ) 

def nodes(n):
    left  = mp.mpf(0)
    right = mp.mpf(n+(n-1)*sqrt(n)*1.01)

    i = 2
    factor = 2

    while True:
        l = [ (x,mp.laguerre(n,0,x)) for x in mp.linspace(left,right,n*factor**i)]

        intervals = []
        for j in range(len(l)-1):
            prod = l[j][1]*l[j+1][1]
            if prod < 0:
                intervals.append([l[j][0], l[j+1][0]])

        if len(intervals) == n:
            break
        i += 1

    roots = []
    f = lambda x: mp.laguerre(n,0,x)
    for ab in intervals:
        a,b = ab
        try:
            z = mp.findroot(f, (a, b), tol=1e-50, solver='bisect')
        except:
            z = bisect(f, a, b, tol=1e-50)
        roots.append( z )

    return roots

    def locate(self,times,exclude):
        """locate(t,exclude)
        input: list-of-times in microseconds
        and what to exclude in calculations
        output (x, y) in meters
        """
        sos = 340.2;
        
        times = times[:];
        pts = self.sensors[:];
        del pts[exclude];
        del times[exclude];
        x,y,t = sympy.symbols('x y t');
        a,b = pts[0];
        c,d = pts[1];
        e,f = pts[2];
        #print(a,b,c,d,e,f);
        t_1 = mpf(times[0])/1e6;
        t_2 = mpf(times[1])/1e6;
        t_3 = mpf(times[2])/1e6;
        #t_1,t_2,t_3 = 0,t_2-t_1,t_3-t_1;
        
        f_1 = (x-a)**2 + (y-b)**2 - ((t_1-t)*sos)**2;
        f_2 =  (x-c)**2 + (y-d)**2 - ((t_2-t)*sos)**2;
        f_3 = (x-e)**2 + (y-f)**2 - ((t_3-t)*sos)**2;

        try:
            res = sympy.nsolve((f_1,f_2,f_3),(x,y,t),(0,0,min(t_1,t_2,t_3)));
            return ( float(res[0]), float(res[1]) );
        except(ValueError,ZeroDivisionError):
            return ( float('NaN'), float('NaN') );

 def hazard_rate(self, x):
     if x < self.location:
         return 0
     elif self.shape == 1:
         return mpf(1) / self.scale
     else:
         return abs(self.pdf_eval(x) / (mpf(1) - self.cdf_eval(x)))

    def cdf_eval(self, x):
        x = mpf(x)

        if x < self.location:
            return 0

        return mpf(1) - mpmath.exp(-mpmath.power(((x-self.location)/self.scale), self.shape))

 def Test_fix():
     f = mpFormat()
     x = mpf("1.2345678901234567890")
     f.digits( 0); assert(f.fix(x)  == " 1.")
     f.digits( 1); assert(f.fix(x)  == " 1.2")
     f.digits( 2); assert(f.fix(x)  == " 1.23")
     f.digits( 3); assert(f.fix(x)  == " 1.235")
     f.digits(20); assert(f.fix(x)  == " 1.23456789012345678900")
     f.digits(20); assert(f.fix(-x) == "-1.23456789012345678900")
     x *= mpf("1e9")
     f.digits( 0); assert(f.fix(x) == " 1234567890.")
     f.digits( 1); assert(f.fix(x) == " 1234567890.1")
     f.digits( 2); assert(f.fix(x) == " 1234567890.12")
     f.digits( 3); assert(f.fix(x) == " 1234567890.123")
     f.digits( 4); assert(f.fix(x) == " 1234567890.1235")
     f.digits(10); assert(f.fix(x) == " 1234567890.1234567890")
     x /= mpf("1e12")
     f.digits( 0); assert(f.fix(x)  == " 0.")
     f.digits( 0); assert(f.fix(-x) == "-0.")
     f.digits( 1); assert(f.fix(x)  == " 0.0")
     f.digits( 1); assert(f.fix(-x) == "-0.0")
     f.digits( 2); assert(f.fix(x)  == " 0.00")
     f.digits( 2); assert(f.fix(-x) == "-0.00")
     f.digits( 3); assert(f.fix(x)  == " 0.001")
     f.digits( 3); assert(f.fix(-x) == "-0.001")
     f.digits( 4); assert(f.fix(x)  == " 0.0012")
     f.digits( 5); assert(f.fix(x)  == " 0.00123")
     f.digits( 7); assert(f.fix(x)  == " 0.0012346")
     f.digits(25); assert(f.fix(x)  == " 0.0012345678901234567890000")
     f.digits(25); assert(f.fix(-x) == "-0.0012345678901234567890000")

 def __call__(self, s):
     assert len(s) > 0
     suffix = 1
     if s != "now" and s != "today":
         if len(s) > 1 and s[:2] != "0x":
             if s[-1] in suffixes_ln:
                 exponent = suffixes_ln[s[-1]]
                 if exponent >= 0:
                     suffix = Zn(10**exponent)
                 else:
                     suffix = mpf("1e" + str(exponent))
                 s = s[:-1]
     for func in (self.j, self.i, self.q, self.v, self.r, self.c):
         x = func(s)
         if x != None: 
             if suffix == 1:
                 return x
             if isint(x):
                 if isint(suffix): return Zn(suffix*x)
                 else:             return suffix*mpf(int(x))
             elif isinstance(x, Rational):
                 if isint(suffix): return Rational(suffix*x.n, x.d)
                 else:             return suffix*x.mpf()
             elif isinstance(x, mpf) or \
                  isinstance(x, mpc) or \
                  isinstance(x, mpi):
                 if isint(suffix): return mpf(int(suffix))*x
                 else:             return suffix*x
             else:
                 return None
     return None

 def Test_sci():
     f = mpFormat()
     mpFormat.show_zero_exponent = True
     x = mpf("1.2345678901234567890")
     f.digits( 0); assert(f.sci( x) == " 1.e+0")
     f.digits( 0); assert(f.sci(-x) == "-1.e+0")
     f.digits( 1); assert(f.sci(x)  == " 1.e+0")
     f.digits( 1); assert(f.sci(-x) == "-1.e+0")
     f.digits( 2); assert(f.sci(x)  == " 1.2e+0")
     f.digits( 3); assert(f.sci(x)  == " 1.23e+0")
     f.digits( 3); assert(f.sci(-x) == "-1.23e+0")
     f.digits(20); assert(f.sci(x)  == " 1.2345678901234567890e+0")
     mpFormat.show_zero_exponent = False
     f.digits( 0); assert(f.sci( x) == " 1.")
     f.digits( 0); assert(f.sci(-x) == "-1.")
     f.digits( 1); assert(f.sci(x)  == " 1.")
     f.digits( 1); assert(f.sci(-x) == "-1.")
     f.digits( 2); assert(f.sci(x)  == " 1.2")
     f.digits( 3); assert(f.sci(x)  == " 1.23")
     f.digits( 3); assert(f.sci(-x) == "-1.23")
     f.digits(20); assert(f.sci(x)  == " 1.2345678901234567890")
     x *= mpf("1e9")
     f.digits( 0); assert(f.sci(x)  == " 1.e+9")
     f.digits( 0); assert(f.sci(-x) == "-1.e+9")
     f.digits( 1); assert(f.sci(x)  == " 1.e+9")
     f.digits( 2); assert(f.sci(x)  == " 1.2e+9")
     f.digits( 2); assert(f.sci(-x) == "-1.2e+9")
     x /= mpf("1e18")
     f.digits( 0); assert(f.sci(x)  == " 1.e-9")
     f.digits( 1); assert(f.sci(x)  == " 1.e-9")
     f.digits( 2); assert(f.sci(x)  == " 1.2e-9")
     f.digits(20); assert(f.sci( x) == " 1.2345678901234567890e-9")
     f.digits(20); assert(f.sci(-x) == "-1.2345678901234567890e-9")

def solve():
    # mpmath.mp = 100
    fn1 = [0]*9
    fn2 = [0]*9
    fn1[8] = 1
    fn2[8] = 1
    k = 3
    vn1 = mpmath.mpf(1.000000000e0)
    vn2 = mpmath.mpf(1.000000000e0)

    while True:
        # if k % 1000 == 0:
        #     print k

        fn = add_last_digits(fn1,fn2)
        vn = vn1+vn2
        if last_digits_are_pandigital(fn):
            if first_digits_are_pandigital(vn):
                mpmath.nprint(vn,n=10)
                return k
      
        vn2 = vn1
        vn1 = vn

        fn2 = fn1
        fn1 = fn
        k += 1

def _noncentral_chi_cdf(x, df, nc, dps=None):
    if dps is None:
        dps = mpmath.mp.dps
    x, df, nc = mpmath.mpf(x), mpmath.mpf(df), mpmath.mpf(nc)
    with mpmath.workdps(dps):
        res = mpmath.quad(lambda t: _noncentral_chi_pdf(t, df, nc), [0, x])
        return res

def multiprec_pdf(df, mu, x):
    df = mp.mpf(str(df))
    mu = mp.mpf(str(mu))
    x  = mp.mpf(str(x))
    pdf = (mp.exp(-mu**2 / 2) * 2**(-df) * df**(df/2) * mp.gamma(df+1) * \
           ((mp.sqrt(2) * x * mu * (x**2+df)**(-df/2-1) * mp.hyp1f1(df/2+1,1.5,(mu**2 * x**2)/(2 * (x**2+df))))/(mp.gamma((df+1)/2))+((x**2+df)**(-df/2-0.5) * mp.hyp1f1((df+1)/2,0.5,(mu**2 * x**2)/(2 * (x**2+df))))/(mp.gamma(df/2+1))))/(mp.gamma(df/2))
    return pdf