"Fossies" - the Fresh Open Source Software Archive

Member "NZMATH-1.2.0/nzmath/equation.py" (19 Nov 2012, 9601 Bytes) of package /linux/misc/old/NZMATH-1.2.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "equation.py" see the Fossies "Dox" file reference documentation.

    1 """
    2 equation -- methods to solve algebraic equations.
    3 
    4 If you would like to solve an equation in the complex field and need
    5 more precision, then you can import and use plugins.SETPRECISION.
    6 
    7 e.g:
    8 from nzmath.plugins import SETPRECISION
    9 SETPRECISION(200)
   10 
   11 Then, the following computations will be done in such precision, if
   12 plugins.PRECISION_CHANGEABLE is true.
   13 """
   14 
   15 from __future__ import division
   16 import logging
   17 
   18 import nzmath.arith1 as arith1
   19 import nzmath.bigrandom as bigrandom
   20 import nzmath.bigrange as bigrange
   21 import nzmath.finitefield as finitefield
   22 import nzmath.imaginary as imaginary
   23 import nzmath.poly.uniutil as uniutil
   24 
   25 from nzmath.plugins import MATHMODULE as math, CMATHMODULE as cmath, \
   26      FLOATTYPE as Float, COMPLEXTYPE as Complex
   27 
   28 _log = logging.getLogger('nzmath.equation')
   29 _log.setLevel(logging.DEBUG)
   30 
   31 # x is (list,tuple)
   32 # t is variable
   33 def e1(x):
   34     """
   35     0 = x[0] + x[1]*t
   36     """
   37     if x[1] == 0:
   38         raise ZeroDivisionError("No Solution")
   39     else:
   40         return -x[0]/x[1]
   41 
   42 def e1_ZnZ(x, n):
   43     """
   44     Return the solution of x[0] + x[1]*t = 0 (mod n).
   45     x[0], x[1] and n must be positive integers.
   46     """
   47     try:
   48         return (-x[0] * arith1.inverse(x[1], n)) % n
   49     except ZeroDivisionError:
   50         raise ValueError("No Solution")
   51 
   52 def e2(x):
   53     """
   54     0 = x[0] + x[1]*t + x[2]*t**2
   55     """
   56     c, b, a = x
   57     d = b**2 - 4*a*c
   58     if d >= 0:
   59         sqrtd = math.sqrt(d)
   60     else:
   61         sqrtd = cmath.sqrt(d)
   62     return ((-b + sqrtd)/(2*a), (-b - sqrtd)/(2*a))
   63 
   64 def e2_Fp(x, p):
   65     """
   66     p is prime
   67     f = x[0] + x[1]*t + x[2]*t**2
   68     """
   69     c, b, a = [_x % p for _x in x]
   70     if a == 0:
   71         return [e1_ZnZ([c, b], p)]
   72     if p == 2:
   73         solutions = []
   74         if x[0] & 1 == 0:
   75             solutions.append(0)
   76         if (x[0] + x[1] + x[2]) & 1 == 0:
   77             solutions.append(1)
   78         if len(solutions) == 1:
   79             return solutions * 2
   80         return solutions
   81     d = b**2 - 4*a*c
   82     if arith1.legendre(d, p) == -1:
   83         return []
   84     sqrtd = arith1.modsqrt(d, p)
   85     a = arith1.inverse(2*a, p)
   86     return [((-b+sqrtd)*a)%p, ((-b-sqrtd)*a)%p]
   87 
   88 def e3(x):
   89     """
   90     0 = x[0] + x[1]*t + x[2]*t**2 + x[3]*t**3
   91     """
   92     x3 = Float(x[3])
   93     a = x[2] / x3
   94     b = x[1] / x3
   95     c = x[0] / x3
   96     p = b - (a**2)/3
   97     q = 2*(a**3)/27 - a*b/3 + c
   98     w = (-1 + cmath.sqrt(-3)) / 2
   99     W = (1, w, w.conjugate())
  100     k = -q/2 + cmath.sqrt((q**2)/4 + (p**3)/27)
  101     l = -q/2 - cmath.sqrt((q**2)/4 + (p**3)/27)
  102     m = k ** (1/3)
  103     n = l ** (1/3)
  104 
  105     # choose n*W[i] by which m*n*W[i] be the nearest to -p/3
  106     checker = [abs(3*m*n*z + p) for z in W]
  107     n = n * W[checker.index(min(checker))]
  108 
  109     sol = []
  110     for i in range(3):
  111         sol.append(W[i]*m + W[-i]*n - a/3)
  112     return sol
  113 
  114 def e3_Fp(x, p):
  115     """
  116     p is prime
  117     0 = x[0] + x[1]*t + x[2]*t**2 + x[3]*t**3
  118     """
  119     x.reverse()
  120     lc_inv = finitefield.FinitePrimeFieldElement(x[0], p).inverse()
  121     coeff = []
  122     for c in x[1:]:
  123         coeff.append((c * lc_inv).n)
  124     sol = []
  125     for i in bigrange.range(p):
  126         if (i**3 + coeff[0]*i**2 + coeff[1]*i + coeff[2]) % p == 0:
  127             sol.append(i)
  128             break
  129     if len(sol) == 0:
  130         return sol
  131     X = e2_Fp([coeff[1] + (coeff[0] + sol[0])*sol[0], coeff[0] + sol[0], 1], p)
  132     if len(X) != 0:
  133         sol.extend(X)
  134     return sol
  135 
  136 def Newton(f, initial=1, repeat=250):
  137     """
  138     Compute x s.t. 0 = f[0] + f[1] * x + ... + f[n] * x ** n
  139     """
  140     length = len(f)
  141     df = []
  142     for i in range(1, length):
  143         df.append(i * f[i])
  144     l = initial
  145     for k in range(repeat):
  146         coeff = Float(0)
  147         dfcoeff = Float(0)
  148         for i in range(1, length):
  149             coeff = coeff * l + f[-i]
  150             dfcoeff += dfcoeff * l + df[-i]
  151         coeff = coeff * l + f[0]
  152         if coeff == 0:
  153             return l
  154         elif dfcoeff == 0: # Note coeff != 0
  155             raise ValueError("There is not solution or Choose different initial")
  156         else:
  157             l = l - coeff / dfcoeff
  158     return l
  159 
  160 
  161 def SimMethod(f, NewtonInitial=1, repeat=250):
  162     """
  163     Return zeros of a polynomial given as a list.
  164     """
  165     if NewtonInitial != 1:
  166         ni = NewtonInitial
  167     else:
  168         ni = None
  169     return _SimMethod(f, newtoninitial=ni, repeat=repeat)
  170 
  171 
  172 def _SimMethod(g, initials=None, newtoninitial=None, repeat=250):
  173     """
  174     Return zeros of a polynomial given as a list.
  175 
  176     - g is the list of the polynomial coefficient in ascending order.
  177     - initial (optional) is a list of initial approximations of zeros.
  178     - newtoninitial (optional) is an initial value for Newton method to
  179       obtain an initial approximations of zeros if 'initial' is not given.
  180     - repeat (optional) is the number of iteration. The default is 250.
  181     """
  182     if initials is None:
  183         if newtoninitial:
  184             z = _initialize(g, newtoninitial)
  185         else:
  186             z = _initialize(g)
  187     else:
  188         z = initials
  189 
  190     f = uniutil.polynomial(enumerate(g), imaginary.theComplexField)
  191     deg = f.degree()
  192     df = f.differentiate()
  193 
  194     value_list = [f(z[i]) for i in range(deg)]
  195     for loop in range(repeat):
  196         sigma_list = [0] * deg
  197         for i in range(deg):
  198             if not value_list[i]:
  199                 continue
  200             sigma = 0
  201             for j in range(i):
  202                 sigma += 1 / (z[i] - z[j])
  203             for j in range(i+1, deg):
  204                 sigma += 1 / (z[i] - z[j])
  205             sigma_list[i] = sigma
  206 
  207         for i in range(deg):
  208             if not value_list[i]:
  209                 continue
  210             ratio = value_list[i] / df(z[i])
  211             z[i] -= ratio / (1 - ratio*sigma_list[i])
  212             value_list[i] = f(z[i])
  213 
  214     return z
  215 
  216 def _initialize(g, newtoninitial=None):
  217     """
  218     create initial values of equation given as a list g.
  219     """
  220     g = [float(c) for c in g[:]]
  221     q = [-abs(c) for c in g[:-1]]
  222     q.append(abs(g[-1]))
  223     if newtoninitial is None:
  224         r = Newton(q, _upper_bound_of_roots(q))
  225     else:
  226         r = Newton(q, newtoninitial)
  227 
  228     deg = len(g) - 1
  229     center = -g[-2]/(deg*g[-1])
  230     about_two_pi = 6
  231     angular_step = cmath.exp(Complex(0, 1) * about_two_pi / deg)
  232     angular_move = r
  233     z = []
  234     for i in range(deg):
  235         z.append(center + angular_move)
  236         angular_move *= angular_step
  237 
  238     return z
  239 
  240 def _upper_bound_of_roots(g):
  241     """
  242     Return an upper bound of roots.
  243     """
  244     weight = len(filter(None, g))
  245     assert g[-1]
  246     return max([pow(weight*abs(c/g[-1]), 1/len(g)) for c in g])
  247 
  248 def root_Fp(g, p, flag=True):
  249     """
  250     Return a root over F_p of nonzero polynomial g.
  251     p must be prime.
  252     If flag = False, return a root randomly
  253     """
  254     if isinstance(g, list):
  255         if not isinstance(g[0], tuple):
  256             g = zip(range(len(g)), g)
  257     Fp = finitefield.FinitePrimeField(p)
  258     g = uniutil.FinitePrimeFieldPolynomial(g, Fp)
  259     h = uniutil.FinitePrimeFieldPolynomial({1:-1, p:1}, Fp)
  260     g = g.gcd(h)
  261     deg_g = g.degree()
  262     if g[0] == 0:
  263         deg_g = deg_g - 1
  264         g = g.shift_degree_to(deg_g)
  265     while True:
  266         if deg_g == 0:
  267             return None
  268         if deg_g == 1:
  269             return (-g[0]/g[1]).toInteger()
  270         elif deg_g == 2:
  271             d = g[1]*g[1] - 4*g[0]
  272             e = arith1.modsqrt(d.toInteger(), p)
  273             return ((-g[1]-e)/(2*g[2])).toInteger()
  274         deg_h = 0
  275         x = uniutil.FinitePrimeFieldPolynomial({0:-1, (p-1)>>1:1}, Fp)
  276         if flag:
  277             a = 0
  278             while deg_h == 0 or deg_h == deg_g:
  279                 b = uniutil.FinitePrimeFieldPolynomial({0:-a, 1:1}, Fp)
  280                 v = g(b)
  281                 h = x.gcd(v)
  282                 a = a + 1
  283                 deg_h = h.degree()
  284                 b = uniutil.FinitePrimeFieldPolynomial({0:a-1, 1:1}, Fp)
  285         else:
  286             while deg_h == 0 or deg_h == deg_g:
  287                 a = bigrandom.randrange(p)
  288                 b = uniutil.FinitePrimeFieldPolynomial({0:-a, 1:1}, Fp)
  289                 v = g(b)
  290                 h = x.gcd(v)
  291                 deg_h = h.degree()
  292                 b = uniutil.FinitePrimeFieldPolynomial({0:a, 1:1}, Fp)
  293         g = h(b)
  294         deg_g = deg_h
  295 
  296 def allroots_Fp(g, p):
  297     """
  298     Return roots over F_p of nonzero polynomial g.
  299     p must be prime.
  300     """
  301     if isinstance(g, list):
  302         if not isinstance(g[0], tuple):
  303             g = zip(range(len(g)), g)
  304     Fp = finitefield.FinitePrimeField(p)
  305     g = uniutil.FinitePrimeFieldPolynomial(g, Fp)
  306     h = uniutil.FinitePrimeFieldPolynomial({1:-1, p:1}, Fp)
  307     g = g.gcd(h)
  308     deg_g = g.degree()
  309     roots = []
  310     if g[0] == 0:
  311         roots.append(0)
  312         deg_g = deg_g - 1
  313         g = g.shift_degree_to(deg_g)
  314     return roots + roots_loop(g, deg_g, p, Fp)
  315 
  316 def roots_loop(g, deg_g, p, Fp):
  317     if deg_g == 0:
  318         return []
  319     if deg_g == 1:
  320         return [(-g[0]/g[1]).toInteger()]
  321     elif deg_g == 2:
  322         d = g[1]*g[1]-4*g[0]
  323         e = arith1.modsqrt(d.toInteger(), p)
  324         g1 = -g[1]
  325         g2 = 2*g[2]
  326         return [((g1 - e) / g2).toInteger(), ((g1 + e) / g2).toInteger()]
  327     deg_h = 0
  328     x = uniutil.FinitePrimeFieldPolynomial({0:-1, (p-1)>>1:1}, Fp)
  329     a = 0
  330     while deg_h == 0 or deg_h == deg_g:
  331         b = uniutil.FinitePrimeFieldPolynomial({0:-a, 1:1}, Fp)
  332         v = g(b)
  333         h = x.gcd(v)
  334         a = a + 1
  335         deg_h = h.degree()
  336     b = uniutil.FinitePrimeFieldPolynomial({0:a-1, 1:1}, Fp)
  337     s = h(b)
  338     deg_s = deg_h
  339     t = g.exact_division(s)
  340     deg_t = t.degree()
  341     return roots_loop(s, deg_s, p, Fp) + roots_loop(t, deg_t, p, Fp)