RELEASE 1.0

compiler/gdsMill/pyx/mathutils.py

@@ -0,0 +1,169 @@
+# -*- coding: ISO-8859-1 -*-
+#
+#
+# Copyright (C) 2005-2006 Michael Schindler <[email protected]>
+# Copyright (C) 2006 André Wobst <[email protected]>
+#
+# This file is part of PyX (http://pyx.sourceforge.net/).
+#
+# PyX is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# PyX is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with PyX; if not, write to the Free Software
+# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
+
+import math, types
+
+# try:
+#     import Numeric, LinearAlgebra
+#     _has_numeric = 1
+# except:
+#     _has_numeric = 0
+
+
+def sign(x):
+    """sign of x, i.e. +1 or -1; returns 1 for x == 0"""
+    if x >= 0:
+        return 1
+    return -1
+
+
+def asinh(x):
+    """Return the arc hyperbolic sine of x."""
+    return math.log(x+math.sqrt(x*x+1))
+
+
+def acosh(x):
+    """Return the arc hyperbolic cosine of x."""
+    return math.log(x+math.sqrt(x*x-1))
+
+
+def _realroots_quadratic(a1, a0):
+    """gives the real roots of x**2 + a1 * x + a0 = 0"""
+    D = a1*a1 - 4*a0
+    if D < 0:
+        return []
+    SD = math.sqrt(D)
+    return [0.5 * (-a1 + SD), 0.5 * (-a1 - SD)]
+
+
+def _realroots_cubic(a2, a1, a0):
+    """gives the real roots of x**3 + a2 * x**2 + a1 * x + a0 = 0"""
+    # see http://mathworld.wolfram.com/CubicFormula.html for details
+
+    Q = (3*a1 - a2*a2) / 9.0
+    R = (9*a2*a1 - 27*a0 - 2*a2*a2*a2) / 54.0
+    D = Q*Q*Q + R*R
+
+    if D > 0:         # one real and two complex roots
+        SD = math.sqrt(D)
+        if R + SD >= 0:
+            S = (R + SD)**(1/3.0)
+        else:
+            S = -(-R - SD)**(1/3.0)
+        if R - SD >= 0:
+            T = (R - SD)**(1/3.0)
+        else:
+            T = -(SD - R)**(1/3.0)
+        return [S + T - a2/3.0]
+    elif D == 0:
+        if Q == 0:    # one real root (R==0)
+            return [-a2/3.0]
+        else:         # two real roots (R>0, Q<0)
+            S = -math.sqrt(-Q)
+            return [2*S - a2/3.0, -S - a2/3.0]
+    else:             # three real roots (Q<0)
+        SQ = math.sqrt(-Q)
+        arg = R / (SQ**3)
+        if arg >= 1:
+            theta = 0
+        elif arg <= -1:
+            theta = math.pi
+        else:
+            theta = math.acos(R/(SQ**3))
+        return [2 * SQ * math.cos((theta + 2*2*i*math.pi)/3.0) - a2/3.0 for i in range(3)]
+
+
+def _realroots_quartic(a3, a2, a1, a0):
+    """gives the real roots of x**4 + a3 * x**3 + a2 * x**2 + a1 * x + a0 = 0"""
+    # see http://mathworld.wolfram.com/QuarticEquation.html for details
+    ys = _realroots_cubic(-a2, a1*a3 - 4*a0, 4*a0*a2 - a1*a1 - a0*a3*a3)
+    ys = [y for y in ys if a3*a3-4*a2+4*y >= 0 and  y*y-4*a0 >= 0]
+    if not ys:
+        return []
+    y1 = min(ys)
+    if a3*y1-2*a1 < 0:
+        return (_realroots_quadratic(0.5*(a3+math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1-math.sqrt(y1*y1-4*a0))) +
+                _realroots_quadratic(0.5*(a3-math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1+math.sqrt(y1*y1-4*a0))))
+    else:
+        return (_realroots_quadratic(0.5*(a3+math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1+math.sqrt(y1*y1-4*a0))) +
+                _realroots_quadratic(0.5*(a3-math.sqrt(a3*a3-4*a2+4*y1)), 0.5*(y1-math.sqrt(y1*y1-4*a0))))
+
+
+def realpolyroots(*cs):
+    """returns the roots of a polynom with given coefficients
+
+    polynomial with coefficients given in cs:
+      0 = \sum_i  cs[i] * x^(len(cs)-i-1)
+    """
+    if not cs:
+        return [0]
+    try:
+        f = 1.0/cs[0]
+        cs = [f*c for c in cs[1:]]
+    except ArithmeticError:
+        return realpolyroots(*cs[1:])
+    else:
+        n = len(cs)
+        if n == 0:
+            return []
+        elif n == 1:
+            return [-cs[0]]
+        elif n == 2:
+            return _realroots_quadratic(*cs)
+        elif n == 3:
+            return _realroots_cubic(*cs)
+        elif n == 4:
+            return _realroots_quartic(*cs)
+        else:
+            raise RuntimeError("realpolyroots solver currently limited to polynoms up to the power of 4")
+
+
+# def realpolyroots_eigenvalue(*cs):
+#     # as realpolyroots but using an equivalent eigenvalue problem
+#     # (this code is currently used for functional tests only)
+#     if not _has_numeric:
+#         raise RuntimeError("realpolyroots_eigenvalue depends on Numeric")
+#     if not cs:
+#         return [0]
+#     try:
+#         f = 1.0/cs[0]
+#         cs = [f*c for c in cs[1:]]
+#     except ArithmeticError:
+#         return realpolyroots_eigenvalue(*cs[1:])
+#     else:
+#         if not cs:
+#             return []
+#         n = len(cs)
+#         a = Numeric.zeros((n, n), Numeric.Float)
+#         for i in range(n-1):
+#             a[i+1][i] = 1
+#         for i in range(n):
+#             a[0][i] = -cs[i]
+#         rs = []
+#         for r in LinearAlgebra.eigenvalues(a):
+#             if type(r) == types.ComplexType:
+#                 if not r.imag:
+#                     rs.append(r.real)
+#             else:
+#                 rs.append(r)
+#         return rs
+#