common:Math
|
|||||||||
| PREV CLASS NEXT CLASS | FRAMES NO FRAMES | ||||||||
| SUMMARY: NESTED | FIELD | CONSTR | METHOD | DETAIL: FIELD | CONSTR | METHOD | ||||||||
Object
class
Mathematical tools library.
This class contains a number of useful functions for performing root-solving and other high-level math functions. They are implemented as static methods so that you do not have to maintain an instance of this class; you need only reference the class itself.
For example, to solve 4x^2 - 4x + 1 = 0:
myroots = new ComplexArray(2) Math.SolveRealQuadraticPolynomial(4,-4,1,0,myroots)
myroots.m_Elements[0] will be 0.5 and myroots.m_Elements[1] will be -0.5. Note that you don't have to create a Math object first, you just use the function directly.
If you are constructing polynomials inside a loop, you should create your roots array object just once, outside the loop, and re-use it each time. This is much faster that constructing a new roots object for each polynomial.
Root solvers are provided for quadratic, cubic, and quartic functions, with either real or complex coefficients. Using the real version (if you know your coefficients are only real values) may be faster in some cases; note that real coefficients may still produce complex roots.
In general, the pflags parameter controls some aspects of root-finding. It is the sum of these values:
1: Return real roots only (only available in the real-coefficient versions).
2: If a root occurs more than once, only return it once.
4: Use iterative improvement for roots (only available on cubic and quartic).
8: Emphasize the accuracy of smaller-magnitude roots, at the expense of accuracy in larger-magnitude roots (only available on cubic and quartic).
16: Force real roots by setting discriminant to zero if it's negative (quadratic only; used internall by the quartic solver).
For example, pflags = 5 will return real roots only and use iterative improvement.
Because of round-off errors, real roots (especially duplicate roots) may be found to have small imaginary parts. If real roots are specified, imaginary parts less than a threshold (in magnitude) are set to zero. This threshold is currently hard-coded inside these formulas and is generally acceptable, but may not be sufficient for your specific use.
Similarly, if duplicate roots are to be removed and two roots differ by less than the threshold, they will be considered equal.
Consider
a*x^n + b*x^(n-1) + c*x^(n-2) + ... = 0 <h5>Order n = 2, quadratic.</h5>
This solution is well known.
x = (-b +- sqrt(b*b-4*a*c))/(2*a)
But there can be cancellation and loss of accuracy. Use this for the root that doesn't lose accuracy. Since x1*x2=c/a, use x2 = c/(a*x1) for the other root.
For real coefficients, the two roots are either both real or a complex-conjugate pair. <h5>Order n = 3, cubic.</h5>
Divide by a, then use the substitution x = t-b/(3*a) to reduce to the "depressed" form,
t^3 + aa*t + bb = 0
For real coefficients, the three roots are either all real or one real and a complex-conjugate pair. <h5>Order n = 4, quartic.</h5>
Divide by a, then use the substitution x = t-b/(4*a) to reduce to the "depressed" form,
t^4 + aa*t^2 + bb*t + cc = 0 There are several ways (at least 5) to proceed. All of them require the solution of a cubic polynomial.
For real coefficients, the four roots are in pairs. In each pair, there are either two reals or a complex-conjugate pair.
If there are multiple roots, round-off may cause the roots to have a small imaginary part. Use the "real" flag setting if you want only real roots.
Quadratic: no other problems.
Cubic and Quartic: If there is a large variation in the magnitude of the roots, the smaller roots may be found less accurately. Some of this is inherent: the polynomial is ill-conditioned (small changes in the coefficients lead to larger changes in the roots). Sometimes it helps to solve for the roots of w = 1/x.
Iterative refinement also sometimes helps. In solving p(x) = 0, suppose we have a "close enough" guess at a root, x0. Then Newton's method can be used to improve it, i.e.,
x1 = x0 - p(x0)/p'(x0)
which converges quadratically if the root is a simple one. Quadratically means that the number of incorrect digits is halved at each iteration.
However, Newton's method converges only linearly to a multiple root. (For a double root, the error is cut in half, approximately, at each iteration.)
Instead we use a formula that includes the second derivative as well:
x1 = x0 - p(x0)*p'(x0)/[ p'(x0)^2 - k*p(x0)*p"(x0)] with k = 1
This converges quadratically even with multiple roots. Note that setting k = 1/2 gives Halley's method, which converges cubically for simple roots, but only linearly with multiple roots.
NOTE: If higher accuracy is needed, use the "Additional Precision" setting on the Formula tab in Layer Properties.
Currently all thresholds are set to 1e-6, a "reasonable" value.
These approximations to the Bessel functions J_0(x) and J_1(x)
(for real x) are accurate to 15 or 16 digits. They come from
Computer Approximations
J. F. Hart, et al.
Wiley, 1968
Spot values have been checked from the text and as calculated by
Mathematica.
These "Bessel" functions J_0(z) and J_1(z) (for complex z) are not quite the actual Bessel functions. Let R = cabs(z). For R <= R1, these are the actual Bessel functions, computed by power series, to full precision except for round-off. For R > R1, these are similar to the actual Bessel functions, with appropriate asymptotic behavior for large R, and matched at R = R1. The actual Bessel functions cannot be calculated with full precision for large R without going to higher-precision calculations, and this is (a) not possible in UF and (b) not really necessary for making fractals.
R1 = 8 is used below.
Jim Blue, June 2008
class Math {
; Mathematical tools library.
;
; <p>
; This class contains a number of useful functions for
; performing root-solving and other high-level math functions.
; They are implemented as static
; methods so that you do not have to maintain an instance
; of this class; you need only reference the class itself.
; <p>
; For example, to solve 4x^2 - 4x + 1 = 0:
; <p>
; myroots = new ComplexArray(2)
; Math.SolveRealQuadraticPolynomial(4,-4,1,0,myroots)
; <p>
; myroots.m_Elements[0] will be 0.5 and myroots.m_Elements[1]
; will be -0.5. Note that you don't have to create a Math
; object first, you just use the function directly.
; <h4>Root Solver General Notes</h4>
; <p>
; If you are constructing polynomials inside a loop, you
; should create your roots array object just once, outside
; the loop, and re-use it each time. This is much faster
; that constructing a new roots object for each polynomial.
; <p>
; Root solvers are provided for quadratic, cubic, and quartic
; functions, with either real or complex coefficients. Using the
; real version (if you know your coefficients are only real
; values) may be faster in some cases; note that real coefficients
; may still produce complex roots.
;
; <p>
; In general, the pflags parameter controls some aspects of root-finding.
; It is the sum of these values:
; <p>
; 1: Return real roots only (only available in the real-coefficient versions).<br>
; 2: If a root occurs more than once, only return it once.<br>
; 4: Use iterative improvement for roots (only available on cubic and quartic).<br>
; 8: Emphasize the accuracy of smaller-magnitude roots, at the expense of accuracy in larger-magnitude roots (only available on cubic and quartic).<br>
; 16: Force real roots by setting discriminant to zero if it's negative (quadratic only; used internall by the quartic solver).
; <p>
; For example, pflags = 5 will return real roots only and use iterative improvement.
; <p>
; Because of round-off errors, real roots (especially duplicate roots) may be
; found to have small imaginary parts. If real roots are specified, imaginary
; parts less than a threshold (in magnitude) are set to zero. This threshold is
; currently hard-coded inside these formulas and is generally acceptable, but
; may not be sufficient for your specific use.
; <p>
; Similarly, if duplicate roots are to be removed and two roots differ by less
; than the threshold, they will be considered equal.
;
; <h4>Specific Solver Notes</h4>
; <p>
; Consider
; <p>
; a*x^n + b*x^(n-1) + c*x^(n-2) + ... = 0
; <h5>Order n = 2, quadratic.</h5>
; <p>
; This solution is well known.
; <p>
; x = (-b +- sqrt(b*b-4*a*c))/(2*a)
; <p>
; But there can be cancellation and loss of accuracy. Use this for the root that
; doesn't lose accuracy. Since x1*x2=c/a, use x2 = c/(a*x1) for the other root.
; <p>
; For real coefficients, the two roots are either both real or a complex-conjugate
; pair.
;
; <h5>Order n = 3, cubic.</h5>
; <p>
; Divide by a, then use the substitution x = t-b/(3*a) to reduce to the "depressed"
; form,
; <p>
; t^3 + aa*t + bb = 0
; <p>
; For real coefficients, the three roots are either all real or one real and a
; complex-conjugate pair.
;
; <h5>Order n = 4, quartic.</h5>
; <p>
; Divide by a, then use the substitution x = t-b/(4*a) to reduce to the "depressed"
; form,
; <p>
; t^4 + aa*t^2 + bb*t + cc = 0
; There are several ways (at least 5) to proceed. All of them require the solution
; of a cubic polynomial.
; <p>
; For real coefficients, the four roots are in pairs. In each pair, there are
; either two reals or a complex-conjugate pair.
;
; <h4>Accuracy</h4>
; <p>
; If there are multiple roots, round-off may cause the roots to have a small
; imaginary part. Use the "real" flag setting if you want only real roots.
; <p>
; Quadratic: no other problems.
; <p>
; Cubic and Quartic: If there is a large variation in the magnitude of the roots,
; the smaller roots may be found less accurately. Some of this is inherent: the
; polynomial is ill-conditioned (small changes in the coefficients lead to larger
; changes in the roots). Sometimes it helps to solve for the roots of w = 1/x.
; <p>
; Iterative refinement also sometimes helps. In solving p(x) = 0, suppose we have
; a "close enough" guess at a root, x0. Then Newton's method can be used to
; improve it, i.e.,
; <p>
; x1 = x0 - p(x0)/p'(x0)
; <p>
; which converges quadratically if the root is a simple one. Quadratically means
; that the number of incorrect digits is halved at each iteration.
; <p>
; However, Newton's method converges only linearly to a multiple root. (For a
; double root, the error is cut in half, approximately, at each iteration.)
; <p>
; Instead we use a formula that includes the second derivative as well:
; <p>
; x1 = x0 - p(x0)*p'(x0)/[ p'(x0)^2 - k*p(x0)*p"(x0)] with k = 1
; <p>
; This converges quadratically even with multiple roots. Note that setting k = 1/2
; gives Halley's method, which converges cubically for simple roots, but only
; linearly with multiple roots.
; <p>
; NOTE: If higher accuracy is needed, use the "Additional Precision" setting on
; the Formula tab in Layer Properties.
; <p>
; Currently all thresholds are set to 1e-6, a "reasonable" value.
; <p>
;
; <h4>Bessel Function Notes</h4>
; <p>
; These approximations to the Bessel functions J_0(x) and J_1(x)
; (for <strong>real</strong> x) are accurate to 15 or 16 digits. They come from <br>
; Computer Approximations<br>
; J. F. Hart, et al.<br>
; Wiley, 1968<br>
; Spot values have been checked from the text and as calculated by
; Mathematica.
; <p>
; These "Bessel" functions J_0(z) and J_1(z)
; (for <strong>complex</strong> z) are not quite the actual Bessel functions.
; Let R = cabs(z). For R <= R1, these are the actual Bessel functions,
; computed by power series, to full precision except for round-off.
; For R > R1, these are
; similar to the actual Bessel functions, with appropriate asymptotic
; behavior for large R, and matched at R = R1. The actual Bessel functions
; cannot be calculated with full precision for large R without going to
; higher-precision calculations, and this is (a) not possible in UF and
; (b) not really necessary for making fractals.
; <p>
; R1 = 8 is used below.
; <p>
; Jim Blue, June 2008
public:
; Constructor
; <p>
; You should never need to actually construct a Math object, as all of its
; methods are static and can be called without an object (just use the class
; name).
func Math()
endfunc
; Returns the sign of a number.
;
; @param px the number to test
; @return the sign; negative numbers give -1, positive numbers (including zero) return 1
static float func Sgn(const float px)
if (px >= 0.0)
return 1.0
else
return -1.0
endif
endfunc
; Determine the maximum of two numbers.
;
; @param px first number
; @param py second number
; @return the largest of the two numbers
static float func MaxF(const float px, const float py)
if (px >= py)
return px
else
return py
endif
endfunc
; Determine the minimum of two numbers.
;
; @param px first number
; @param py second number
; @return the smallest of the two numbers
static float func MinF(const float px, const float py)
if (px <= py)
return px
else
return py
endif
endfunc
; Solve (real) quadratic equation.
; <p>
; See the class overview for details on the solvers.
; This solves a generic quadratic equation of the form
; <p>
; (a1)(x^2) + (a2)(x) + (a3) = 0
;
; @param pa1 a1 term
; @param pa2 a2 term
; @param pa3 a3 term
; @param pflags options
; @param proots reference to ComplexArray object to hold roots; this should already be sized to at least two elements (regardless of whether return values are limited)
; @return the number of roots found (will be 2 unless return values are limited)
static int func SolveRealQuadraticPolynomial( \
const float pa1, const float pa2, const float pa3, \
const int pflags, ComplexArray &proots)
bool forcerealroots = (floor(pflags/16)%2 == 1)
;bool smallrootsbest = (floor(pflags/8)%2 == 1)
;bool iterimprove = (floor(pflags/4)%2 == 1)
bool removeduplicates = (floor(pflags/2)%2 == 1)
bool realrootsonly = ( (pflags )%2 == 1)
; Degenerate case: not a quadratic
if (pa1 == 0)
if (pa2 == 0)
return 0
else
proots.m_Elements[0] = -pa3/pa2
return 1
endif
endif
; Degenerate case: zero root
if (pa3 == 0)
proots.m_Elements[0] = 0
if (removeduplicates && pa2 == 0)
return 1
endif
proots.m_Elements[1] = -pa2/pa1
return 2
endif
float RealThreshold = 1.e-6
float d = pa2*pa2 - 4*pa1*pa3 ; discriminant
if (d < 0 && forcerealroots)
d = 0
endif
if (d < 0)
int nroots = 2
float i2pa1 = 0.5 / pa1
float im = sqrt(-d) * i2pa1
if (realrootsonly || removeduplicates)
if (abs(im) <= RealThreshold)
im = 0
if (removeduplicates)
nroots = 1
endif
else
return 0
endif
endif
float re = -pa2 * i2pa1
proots.m_Elements[0] = re + flip(im)
proots.m_Elements[1] = re - flip(im)
return nroots
else ; d >= 0 so do carefully to avoid possible cancellation
float t = -0.5*(pa2 + sgn(pa2)*sqrt(d))
proots.m_Elements[0] = t/pa1 ; this root is larger in magnitude
proots.m_Elements[1] = pa3/t
if (removeduplicates && \
abs(real(proots.m_Elements[0]-proots.m_Elements[1])) < RealThreshold)
return 1
else
return 2
endif
endif
endfunc ; end of SolveRealQuadraticPolynomial
; Solve (complex) quadratic equation.
; <p>
; See the class overview for details on the solvers.
; This solves a generic quadratic equation of the form
; <p>
; (a1)(x^2) + (a2)(x) + (a3) = 0
;
; @param pa1 a1 term
; @param pa2 a2 term
; @param pa3 a3 term
; @param pflags options
; @param proots reference to ComplexArray object to hold roots; this should already be sized to at least two elements (regardless of whether return values are limited)
; @return the number of roots found (will be 2 unless return values are limited)
static int func SolveComplexQuadraticPolynomial( \
const complex pa1, const complex pa2, const complex pa3, \
const int pflags, ComplexArray &proots)
bool removeduplicates = (floor(pflags/2)%2 == 1)
; Degenerate case: not a quadratic
if (pa1 == 0)
if (pa2 == 0)
return 0
else
proots.m_Elements[0] = -pa3/pa2
return 1
endif
endif
; Degenerate case: zero root
if (pa3 == 0)
proots.m_Elements[0] = 0
if (removeduplicates && pa2 == 0)
return 1
endif
proots.m_Elements[1] = -pa2/pa1
return 2
endif
float Thresholdsq = sqr(1.e-6)
complex sqrtd = sqrt(pa2*pa2 - 4*pa1*pa3)
; do carefully to avoid possible cancellation
; (this gives the same result as the real case when sqrtd is real)
complex t1 = (pa2 + sqrtd)
complex t2 = (pa2 - sqrtd)
complex t
if (|t1| >= |t2|)
t = -0.5*t1
else
t = -0.5*t2
endif
proots.m_elements[0] = t/pa1
proots.m_elements[1] = pa3/t
if (removeduplicates && \
|proots.m_elements[0] - proots.m_elements[1]| < Thresholdsq)
return 1
else
return 2
endif
endfunc ; end of SolveComplexQuadraticPolynomial
; Solve (real) cubic equation.
; <p>
; See the class overview for details on the solvers.
; This solves a generic cubic equation of the form
; <p>
; (a1)(x^3) + (a2)(x^2) + (a3)(x) + (a4) = 0
;
; @param pa1 a1 term
; @param pa2 a2 term
; @param pa3 a3 term
; @param pa4 a4 term
; @param pflags options
; @param proots reference to ComplexArray object to hold roots; this should already be sized to at least three elements (regardless of whether return values are limited)
; @return the number of roots found (will be 3 unless return values are limited)
static int func SolveRealCubicPolynomial( \
const float pa1, const float pa2, const float pa3, const float pa4, \
const int pflags, ComplexArray &proots)
if (pa1 == 0) ; degenerate-not a cubic
return SolveRealQuadraticPolynomial(pa2, pa3, pa4, pflags, proots)
endif
if (pa4 == 0) ; one root of zero
int n = SolveRealQuadraticPolynomial(pa1, pa2, pa3, pflags, proots)
proots.m_Elements[n] = 0
return n+1
endif
bool smallrootsbest = (floor(pflags/8)%2 == 1)
bool iterimprove = (floor(pflags/4)%2 == 1)
bool removeduplicates = (floor(pflags/2)%2 == 1)
bool realrootsonly = ( (pflags )%2 == 1)
; scale out leading coefficient
float a
float b
float c
if (smallrootsbest) ; w = 1/x, solve for roots of w
a = pa3/pa4
b = pa2/pa4
c = pa1/pa4
else
a = pa2/pa1
b = pa3/pa1
c = pa4/pa1
endif
; Now have x^3 + aa*x + bb = 0
; change of variable t = x - a/3 removes the quadratic term, giving
; t^3 - 3*Q*t + 2*R = 0
float a3 = a/3
float aa = b - 3*a3*a3
float bb = c - a3*(b - 2*a3*a3)
float Q = -aa/3 ; (a*a-3*b)/9
float R = bb/2 ; (a*(2*a*a-9*b)+27*c)/54
float D = Q*Q*Q-R*R
; Could here do special cases Q=0, R=0, D=0
; If two equal roots, D may be slightly negative, making them into a
; c.c. pair with small imaginary part.
int nroots = 3
if (D >= 0) ; three real roots
; Note - theta usually done as acos, but atan2 is better
float theta = atan2(R +flip(sqrt(D)))/3
; float theta = acos(R/sqrt(Q*Q*Q))/3 ; is the usual version
float cc = -2*sqrt(Q) ; D > 0, so Q > 0
proots.m_Elements[0] = cc*cos(theta ) -a3
proots.m_Elements[1] = cc*cos(theta + 2*#pi/3) -a3
proots.m_Elements[2] = cc*cos(theta - 2*#pi/3) -a3
else ; one real root, two complex-conjugate roots
; no duplicates possible
float Ar = -Sgn(R) * ((abs(R) + sqrt(-D)))^(1/3)
float Br = 0
if (Ar != 0) ; always Ar!=0 with real coefficients
Br = Q/Ar
endif
proots.m_Elements[0] = Ar+Br - a3
float im = 0.5*sqrt(3)*(Ar-Br)
proots.m_Elements[1] = -0.5*(Ar+Br) -a3 + flip(im)
proots.m_Elements[2] = conj(proots.m_Elements[1])
endif
if (smallrootsbest)
int n = 0
while (n < nroots)
proots.m_Elements[n] = 1/proots.m_Elements[n]
n = n + 1
endwhile
endif
; iterative refinement
if (iterimprove)
int n = 0
while (n < nroots)
complex x = proots.m_Elements[n]
complex v0 = pa4 + x*(pa3 + x*(pa2 + pa1*x))
complex v1 = pa3 + x*(2*pa2 + 3*pa1*x)
if (v0*v1 != 0)
complex v2 = 2*pa2 + 6*pa1*x
; complex xnew = x - v0*v1/(v1*v1-v0*v2)
; x = xnew
x = x - v0*v1/(v1*v1-v0*v2)
endif
proots.m_Elements[n] = x
n = n + 1
endwhile
endif
if (realrootsonly)
nroots = RealOnly(nroots, proots)
endif
if (removeduplicates)
nroots = SingleOnly(nroots, proots)
endif
return nroots
endfunc ; end of SolveRealCubicPolynomial
; Solve (complex) cubic equation.
; <p>
; See the class overview for details on the solvers.
; This solves a generic cubic equation of the form
; <p>
; (a1)(x^3) + (a2)(x^2) + (a3)(x) + (a4) = 0
;
; @param pa1 a1 term
; @param pa2 a2 term
; @param pa3 a3 term
; @param pa4 a4 term
; @param pflags options
; @param proots reference to ComplexArray object to hold roots; this should already be sized to at least three elements (regardless of whether return values are limited)
; @return the number of roots found (will be 3 unless return values are limited)
static int func SolveComplexCubicPolynomial( \
const complex pa1, const complex pa2, const complex pa3, const complex pa4, \
const int pflags, ComplexArray &proots)
if (pa1 == 0) ; degenerate-not a cubic
return SolveComplexQuadraticPolynomial(pa2, pa3, pa4, pflags, proots)
endif
if (pa4 == 0) ; one root of zero
int n = SolveComplexQuadraticPolynomial(pa1, pa2, pa3, pflags, proots)
proots.m_Elements[n] = 0
return n+1
endif
bool smallrootsbest = (floor(pflags/8)%2 == 1)
bool iterimprove = (floor(pflags/4)%2 == 1)
bool removeduplicates = (floor(pflags/2)%2 == 1)
;bool realrootsonly = ( (pflags )%2 == 1)
; scale out leading coefficient
complex a
complex b
complex c
if (smallrootsbest) ; w = 1/x, solve for roots of w
a = pa3/pa4
b = pa2/pa4
c = pa1/pa4
else
a = pa2/pa1
b = pa3/pa1
c = pa4/pa1
endif
; Now have x^3 + a*x^2 + b*x + c = 0
; change of variable t = x - a/3 removes the quadratic term, giving
; t^3 - 3*Q*t + 2*R = 0
complex a3 = a/3
complex aa = b - 3*a3*a3
complex bb = c - a3*(b - 2*a3*a3)
complex Q = -aa/3 ; (a*a-3*b)/9
complex R = bb/2 ; (a*(2*a*a-9*b)+27*c)/54
complex D = Q*Q*Q-R*R
complex sm = sqrt(-D)
complex t1 = r + sm
complex t2 = r - sm
complex Ar
if (|t1| >= |t2|)
Ar = -(t1^(1/3))
else
Ar = -(t2^(1/3))
endif
complex Br = 0
if (Ar != 0)
Br = Q/Ar
endif
proots.m_Elements[0] = Ar + Br - a3
proots.m_Elements[1] = -0.5*(Ar+Br) -a3 + (0,0.5)*sqrt(3)*(Ar-Br)
proots.m_Elements[2] = -0.5*(Ar+Br) -a3 - (0,0.5)*sqrt(3)*(Ar-Br)
int nroots = 3
if (smallrootsbest)
int n = 0
while (n < nroots)
proots.m_Elements[n] = 1/proots.m_Elements[n]
n = n + 1
endwhile
endif
; iterative refinement
if (iterimprove)
int n = 0
while (n < nroots)
complex x = proots.m_Elements[n]
complex v0 = pa4 + x*(pa3 + x*(pa2 + pa1*x))
complex v1 = pa3 + x*(2*pa2 + 3*pa1*x)
if (v0*v1 != 0)
complex v2 = 2*pa2 + 6*pa1*x
; complex xnew = x - v0*v1/(v1*v1-v0*v2)
; x = xnew
x = x - v0*v1/(v1*v1-v0*v2)
endif
proots.m_Elements[n] = x
n = n+1
endwhile
endif
if (removeduplicates)
nroots = SingleOnly(nroots, proots)
endif
return nroots
endfunc ; end of SolveComplexCubicPolynomial
; Solve (real) quartic equation.
; <p>
; See the class overview for details on the solvers.
; This solves a generic quartic equation of the form
; <p>
; (a1)(x^4) + (a2)(x^3) + (a3)(x^2) + (a4)(x) + (a5) = 0
;
; @param pa1 a1 term
; @param pa2 a2 term
; @param pa3 a3 term
; @param pa4 a4 term
; @param pa5 a5 term
; @param pflags options
; @param proots reference to ComplexArray object to hold roots; this should already be sized to at least four (regardless of whether return values are limited)
; @return the number of roots found (will be 4 unless return values are limited)
static int func SolveRealQuarticPolynomial( \
const float pa1, const float pa2, const float pa3, const float pa4, \
const float pa5, const int pflags, ComplexArray &proots)
if (pa1 == 0) ; degenerate: not a quartic
return SolveRealCubicPolynomial(pa2, pa3, pa4, pa5, pflags, proots)
endif
if (pa5 == 0) ; one root of zero
int n = SolveRealCubicPolynomial(pa1, pa2, pa3, pa4, pflags, proots)
proots.m_Elements[n] = 0
return n+1
endif
bool smallrootsbest = (floor(pflags/8)%2 == 1)
bool iterimprove = (floor(pflags/4)%2 == 1)
bool removeduplicates = (floor(pflags/2)%2 == 1)
bool realrootsonly = ( (pflags )%2 == 1)
float a
float b
float c
float d
; scale out leading coefficient
if (smallrootsbest) ; w = 1/x, solve for roots of w
a = pa4/pa5
b = pa3/pa5
c = pa2/pa5
d = pa1/pa5
else
a = pa2/pa1
b = pa3/pa1
c = pa4/pa1
d = pa5/pa1
endif
; Now have x^4 + a*x^3 + b*x^2 + c*x + d = 0
; change of variable t = x - a/4 removes the cubic term, giving
; the "depressed quartic"
; t^4 + aa*t^2 + bb*t + cc = 0
float a4 = a/4
float aa = b - 6*a4*a4
float bb = c - 2*a4*(b - 4*a4*a4)
float cc = d - a4*(c-a4*(b-3*a4*a4))
; Could here do any special cases
; The next substitution is from
; http://en.wikipedia.org/wiki/Quartic_equation
; It gives the same result as Faucette, who is quoted on Wolfram's site,
; http://mathworld.wolfram.com/QuarticEquation.html
;
; Factor the depressed quartic into two quadratic terms
; (t^2 +p*t + q)*(t^2 + r*t + s) = 0
; Equating terms gives
; 0 = p + r so r = -p
; aa = q - p^2 + s
; bb = p * (s - q)
; cc = q * s
; Some algebra gives a cubic in p^2
; p^6 + 2*aa*p^4 + (aa^2-4*cc)*p^2 - b^2 = 0
; Solve this and choose a real root for p^2.
; if p^2 > 0, one sign of the square root gives
; s and the other sign gives q,
; then 2s = aa +-b/p +p^2
; then 2q = aa -+b/p +p^2
; but we choose the better equation for s and then q = cc/s.
; if p^2 = 0, use the aa and cc equations to get a quadratic for s and q.
float c1 = 2*aa
float c2 = aa*aa-4*cc
float c3 = -bb*bb
int n1 = SolveRealCubicPolynomial(1, c1, c2, c3, 1, proots)
; Since -p is the sum of 2 roots of the quartic, at least one of the 3 roots
; of the cubic is real and non-negative (within round-off).
; If all roots of the quartic are real, all 3 roots of the cubic are real
; and non-negative. Also, n1 is 1 or 3.
float p2
if (n1 == 1)
p2 = real(proots.m_Elements[0])
else
p2 = Maxf(real(proots.m_Elements[0]), real(proots.m_Elements[1]))
p2 = Maxf(p2, real(proots.m_Elements[2]))
; minimum positive root maybe gives most accuracy
; float p2_0, float p2_1, float p2_2
; p2_0 = real(proots.m_Elements[0])
; p2_1 = real(proots.m_Elements[1])
; p2_2 = real(proots.m_Elements[2])
; Sort3f(p2_0, p2_1, p2_2)
; p2 = p2_0
; if (p2 + RealThreshold < 0)
; p2 = p2_1
; if (p2 + RealThreshold < 0)
; p2 = p2_2
; endif
; endif
endif
if (p2 < 0) ; only because of roundoff
p2 = 0
endif
float s
float q
int n3 = SolveRealQuadraticPolynomial(1, -(aa+p2), cc, 16, proots)
if (n3 == 0) ; can't happen -- means p2 = 0, but round-off...
return 0
endif
s = real(proots.m_Elements[0])
q = real(proots.m_Elements[1])
float p = sgn(bb)*sgn(s-q)*sqrt(p2)
; Now solve the two quadratics; since the leading term is 1, each has two
; roots, either both real or c.c.
int n2 = SolveRealQuadraticPolynomial(1, p, q, 0, proots)
proots.m_Elements[2] = proots.m_Elements[0] - a4
proots.m_Elements[3] = proots.m_Elements[1] - a4
int n0 = SolveRealQuadraticPolynomial(1, -p, s, 0, proots)
proots.m_Elements[0] = proots.m_Elements[0] - a4
proots.m_Elements[1] = proots.m_Elements[1] - a4
int nroots = n0 + n2
if (nroots == 0) ; shouldn't happen
return 0
endif
if (smallrootsbest)
int n = 0
while (n < nroots)
proots.m_Elements[n] = 1/proots.m_Elements[n]
n = n + 1
endwhile
endif
; iterative refinement
if (iterimprove)
int n = 0
while (n < nroots)
complex x = proots.m_Elements[n]
complex v0 = pa5 + x*(pa4 + x*(pa3 + x*(pa2 + pa1*x)))
complex v1 = pa4 + x*(2*pa3 + x*(3*pa2 + 4*pa1*x))
if (v0*v1 != 0)
; x = x - v0/v1
complex v2 = 2*pa3 + x*(6*pa2 + 12*pa1*x)
x = x - v0*v1/(v1*v1-v0*v2)
endif
proots.m_Elements[n] = x
n = n + 1
endwhile
endif
if (realrootsonly)
nroots = RealOnly(nroots, proots)
endif
if (removeduplicates)
nroots = SingleOnly(nroots, proots)
endif
return nroots
endfunc ; end of SolveRealQuartic Polynomial
; Solve (complex) quartic equation.
; <p>
; See the class overview for details on the solvers.
; This solves a generic quartic equation of the form
; <p>
; (a1)(x^4) + (a2)(x^3) + (a3)(x^2) + (a4)(x) + (a5) = 0
;
; @param pa1 a1 term
; @param pa2 a2 term
; @param pa3 a3 term
; @param pa4 a4 term
; @param pa5 a5 term
; @param pflags options
; @param proots reference to ComplexArray object to hold roots; this should already be sized to at least four (regardless of whether return values are limited)
; @return the number of roots found (will be 4 unless return values are limited)
static int func SolveComplexQuarticPolynomial( \
const complex pa1, const complex pa2, const complex pa3, const complex pa4, \
const complex pa5, const int pflags, ComplexArray &proots)
if (pa1 == 0) ; degenerate: not a quartic
return SolveComplexCubicPolynomial(pa2, pa3, pa4, pa5, pflags, proots)
endif
if (pa5 == 0) ; one root of zero
int n = SolveComplexCubicPolynomial(pa1, pa2, pa3, pa4, pflags, proots)
proots.m_Elements[n] = 0
return n+1
endif
bool smallrootsbest = (floor(pflags/8)%2 == 1)
bool iterimprove = (floor(pflags/4)%2 == 1)
bool removeduplicates = (floor(pflags/2)%2 == 1)
;bool realrootsonly = ( (pflags )%2 == 1)
complex a
complex b
complex c
complex d
; scale out leading coefficient
if (smallrootsbest) ; w = 1/x, solve for roots of w
a = pa4/pa5
b = pa3/pa5
c = pa2/pa5
d = pa1/pa5
else
a = pa2/pa1
b = pa3/pa1
c = pa4/pa1
d = pa5/pa1
endif
; Now have x^4 + a*x^3 + b*x^2 + c*x + d = 0
; change of variable t = x - a/4 removes the cubic term, giving
; the "depressed quartic"
; t^4 + aa*t^2 + bb*t + cc = 0
complex a4 = a/4
complex aa = b - 6*a4*a4
complex bb = c - 2*a4*(b - 4*a4*a4)
complex cc = d - a4*(c-a4*(b-3*a4*a4))
; Could here do any special cases
; The next substitution is from
; http://en.wikipedia.org/wiki/Quartic_equation
; It gives the same result as Faucette, who is quoted on Wolfram's site,
; http://mathworld.wolfram.com/QuarticEquation.html
;
; Factor the depressed quartic into two quadratic terms
; (t^2 +p*t + q)*(t^2 + r*t + s) = 0
; Equating terms gives
; 0 = p + r so r = -p
; aa = q - p^2 + s
; bb = p * (s - q)
; cc = q * s
; Some algebra gives a cubic in p^2
; p^6 + 2*aa*p^4 + (aa^2-4*cc)*p^2 - b^2 = 0
; Solve this and choose a real root for p^2.
; if p^2 > 0, one sign of the square root gives
; s and the other sign gives q,
; then 2s = aa +-b/p +p^2
; then 2q = aa -+b/p +p^2
; but we choose the better equation for s and then q = cc/s.
; if p^2 = 0, use the aa and cc equations to get a quadratic for s and q.
complex c1 = 2*aa
complex c2 = aa*aa-4*cc
complex c3 = -bb*bb
SolveComplexCubicPolynomial(1, c1, c2, c3, 0, proots)
; Since -p is the sum of 2 roots of the quartic, at least one of the 3 roots
; of the cubic is real and non-negative (within round-off).
; If all roots of the quartic are real, all 3 roots of the cubic are real
; and non-negative. Also, n1 is 1 or 3.
complex p2 = proots.m_Elements[0]
if (|p2| < |proots.m_Elements[1]|)
p2 = proots.m_Elements[1]
endif
if (|p2| < |proots.m_Elements[2]|)
p2 = proots.m_Elements[2]
endif
complex s
complex q
int n3 = SolveComplexQuadraticPolynomial(1, -(aa+p2), cc, 0, proots)
if (n3 == 0) ; can't happen -- means p2 = 0, but round-off...
return 0
endif
s = proots.m_Elements[0]
q = proots.m_Elements[1]
complex p = sqrt(p2)
complex t1 = bb - p * (s-q)
complex t2 = bb = p * (s-q)
if (|t1| > |t2|)
p = -p
endif
; Now solve the two quadratics
int n2 = SolveComplexQuadraticPolynomial(1, p, q, 0, proots)
proots.m_Elements[2] = proots.m_Elements[0] - a4
proots.m_Elements[3] = proots.m_Elements[1] - a4
int n0 = SolveComplexQuadraticPolynomial(1, -p, s, 0, proots)
proots.m_Elements[0] = proots.m_Elements[0] - a4
proots.m_Elements[1] = proots.m_Elements[1] - a4
int nroots = n0 + n2
if (smallrootsbest)
int n = 0
while (n < nroots)
proots.m_Elements[n] = 1/proots.m_Elements[n]
n = n + 1
endwhile
endif
; iterative refinement
if (iterimprove)
int n = 0
while (n < nroots)
complex x = proots.m_Elements[n]
complex v0 = pa5 + x*(pa4 + x*(pa3 + x*(pa2 + pa1*x)))
complex v1 = pa4 + x*(2*pa3 + x*(3*pa2 + 4*pa1*x))
if (v0*v1 != 0)
; x = x - v0/v1
complex v2 = 2*pa3 + x*(6*pa2 + 12*pa1*x)
x = x - v0*v1/(v1*v1-v0*v2)
endif
proots.m_Elements[n] = x
n = n + 1
endwhile
endif
if (removeduplicates)
nroots = SingleOnly(nroots, proots)
endif
return nroots
endfunc ; end of SolveRealQuartic Polynomial
; Complex Gamma() David Makin June 2008
; @param pz complex argument
; @return Gamma function of argument
static complex func ComplexGamma(complex pz)
complex w = pz
; float gplus = 5.2421875
if pz==0.0
w = recip(0)
else
if real(pz)<0.0
w = -pz
if imag(pz)==0.0
if (real(pz)%1)==0.0
w = recip(0)
endif
endif
endif
if real(w)<1e300
w = (sqrt(2.0*#pi)* \
(0.99999999999999709182 \
+ 57.156235665862923517/(w+1.0) \
- 59.597960355475491248/(w+2.0) \
+ 14.136097974741747174/(w+3.0) \
- 0.49191381609762019978/(w+4.0) \
+ .33994649984811888699E-4/(w+5.0) \
+ .46523628927048575665E-4/(w+6.0) \
- .98374475304879564677e-4/(w+7.0) \
+ .15808870322491248884E-3/(w+8.0) \
- .21026444172410488319e-3/(w+9.0) \
+ .21743961811521264320e-3/(w+10.0) \
- .16431810653676389022e-3/(w+11.0) \
+ .84418223983852743293e-4/(w+12.0) \
- .26190838401581408670e-4/(w+13.0) \
+ .36899182659531622704e-5/(w+14.0)) \
/w) * exp(-w-5.2421875)*(w+5.2421875)^(w+0.5)
if real(pz)<0.0
w = pz*w*sin(#pi*pz)
if (w==0.0)
w = recip(0)
else
w = -#pi/w
endif
endif
endif
endif
return w
endfunc
; Real Gamma() David Makin June 2008
; @param x float argument
; @return Gamma function of argument
static float func RealGamma(float x)
float w = x
; float gplus = 5.2421875
if x==0.0
w = recip(0)
else
if x<0.0
w = -x
if (x%1)==0.0
w = recip(0)
endif
endif
if w<1e300
w = (sqrt(2.0*#pi)* \
(0.99999999999999709182 \
+ 57.156235665862923517/(w+1.0) \
- 59.597960355475491248/(w+2.0) \
+ 14.136097974741747174/(w+3.0) \
- 0.49191381609762019978/(w+4.0) \
+ .33994649984811888699E-4/(w+5.0) \
+ .46523628927048575665E-4/(w+6.0) \
- .98374475304879564677e-4/(w+7.0) \
+ .15808870322491248884E-3/(w+8.0) \
- .21026444172410488319e-3/(w+9.0) \
+ .21743961811521264320e-3/(w+10.0) \
- .16431810653676389022e-3/(w+11.0) \
+ .84418223983852743293e-4/(w+12.0) \
- .26190838401581408670e-4/(w+13.0) \
+ .36899182659531622704e-5/(w+14.0)) \
/w) * exp(-w-5.2421875)*(w+5.2421875)^(w+0.5)
if x<0.0
w = x*w*sin(#pi*x)
if (w==0.0)
w = recip(0)
else
w = -#pi/w
endif
endif
endif
endif
return w
endfunc
; Complex Beta() David Makin June 2008
; @param a complex argument
; @param b complex argument
; @return Beta function of arguments
static complex func ComplexBeta(complex a, complex b)
return ComplexGamma(a)*ComplexGamma(b)/ComplexGamma(a+b)
endfunc
; Real Beta() David Makin June 2008
; @param a float argument
; @param b float argument
; @return Beta function of arguments
static float func RealBeta(float a, float b)
return RealGamma(a)*RealGamma(b)/RealGamma(a+b)
endfunc
; Real Bessel J0.
; @param x float argument
; @return Bessel function of order 0
static float func RealBesselJ0(const float x)
; JZERO 5846
float ax = abs(x)
if (ax <= 8.0)
float xx = x*x
return ( \
+.181083899214164308956057747e9 + xx*( \
-.4446475956150241214875200307e8 + xx*( \
+.2629297609723391264477929426e7 + xx*( \
-.6635120296335057670760789901e5 + xx*( \
+.9000010866385970012308555062e3 + xx*( \
-.7411779995119348858162812772e1 + xx*( \
+.39771954892980603699828143e-1 + xx*( \
-.1446351863454593575005971793e-3 + xx*( \
+.3612785949562383677336235392e-6 + xx*( \
-.6089358750808789415483968656e-9 + xx*( \
+.6424736784979100876734782112e-12 + xx*( \
-.3313360272836170910319585771e-15)))))))))))) / ( \
+.1810838992141643042977083347e9 + xx*( \
+.806215242038697916042541338e6 + xx*( \
+.1415495011707417263126447751e4 + xx)))
else
float z = sqr(8.0/ax)
; PZERO 6546
float P0 = ( \
+.8554822541506661710252074e4 + z*( \
+.8894437532960619440762804e4 + z*( \
+.2204501043965180428995069e4 + z*( \
+.128677585748714193298851e3 + z*( \
+.9004793474802880316384e0))))) / ( \
+.8554822541506662842462151e4 + z*( \
+.8903836141709595355210343e4 + z*( \
+.2214048851914710419825683e4 + z*( \
+.130884900499923882835109e3 + z ))))
; QZERO 6946
float Q0 = (8.0/x) * ( \
-.37510534954957111594836e2 + z*( \
-.46093826814625174976377e2 + z*( \
-.13990976865960680088016e2 + z*( \
-.1049732798234554833126e1 + z*( \
-.93525953294031893049e-2)))))/ ( \
+.2400674237117267479318819e4 + z*( \
+.2971983745208491990065486e4 + z*( \
+.921566975526530895082307e3 + z*( \
+.74428389741411178824152e2 + z ))))
z = ax - #pi/4
return sqrt(2/#pi/ax)*(P0*cos(z)-Q0*sin(z))
endif
endfunc
; Real Bessel J1.
; @param x float argument
; @return Bessel function of order 1
static float func RealBesselJ1(const float x)
; JONE 6045
float ax = abs(x)
if (ax <= 8.0)
float xx = x*x
return x * ( \
+.695364226329838502166085207e8 + xx*( \
-.8356785487348914291918495672e7 + xx*( \
+.3209027468853947029888682298e6 + xx*( \
-.58787877666568200462726094e4 + xx*( \
+.6121876997356943874446879769e2 + xx*( \
-.3983107983952332023421699105e0 + xx*( \
+.1705769264349617107854016566e-2 + xx*( \
-.4910599276555129440130592573e-5 + xx*( \
+.9382193365140744507653268479e-8 + xx*( \
-.1107352224453730633782671362e-10 + xx*( \
+.63194310317443161294700346e-14 ))))))))))) / ( \
+.139072845265967685120764336e9 + xx*( \
+.6705346835482299302199750802e6 + xx*( \
+.1284593453966301898121332163e4 + xx)))
else
float z = sqr(8.0/ax)
; PONE 6747
float P1 = ( \
+.1290918471896188077350689e5 + z*( \
+.1309042051103506486292571e5 + z*( \
+.313275295635506951011069e4 + z*( \
+.17431379748379024599685e3 + z*( \
+.122850537643590432633e1 )))))/ ( \
+.1290918471896187879332737e5 + z*( \
+.1306678308784401036110575e5 + z*( \
+.310928141677002883350924e4 + z*( \
+.16904721775008609992033e3 + z ))))
; QONE 7147
float Q1 = (8.0/ax) * ( \
+.14465282874995208675225e3 + z*( \
+.1744291689092425885102e3 + z*( \
+.5173653281836591636536e2 + z*( \
+.379944537969806734901e1 + z*( \
+.36363466476034710809e-1 )))))/ ( \
+.308592701333231723110639e4 + z*( \
+.373434010601630179517765e4 + z*( \
+.11191098527047487025919e4 + z*( \
+.8522392064341340397334e2 + z ))))
z = ax - #pi*3/4
z = sqrt(2/#pi/ax)*(P1*cos(z)-Q1*sin(z))
if (x < 0)
return -z
else
return z
endif
endif
endfunc
; Complex Bessel J0.
; @param z complex argument
; @return "Bessel" function of order 0
static complex func ComplexBesselJ0(const complex z)
float R = cabs(z)
if (R <= 8)
return ComplexBesselJ0_Inner(z)
else
complex zz = z * 8/R
return ComplexBesselJ0_Outer(z) * (ComplexBesselJ0_Inner(zz)/ComplexBesselJ0_Outer(zz))
endif
endfunc
; Complex Bessel J1.
; @param z complex argument
; @return "Bessel" function of order 0
static complex func ComplexBesselJ1(const complex z)
float R = cabs(z)
if (R <= 8)
return ComplexBesselJ1_Inner(z)
else
complex zz = z * 8/R
return ComplexBesselJ1_Outer(z) * (ComplexBesselJ1_Inner(zz)/ComplexBesselJ1_Outer(zz))
endif
endfunc
private:
; Helper function to eliminate duplicate roots.
static int func SingleOnly(int nroots, ComplexArray &proots)
float DuplThreshold = 1.e-12
int n = nroots
while (n > 0)
n = n - 1
complex xx = proots.m_Elements[n]
bool dupl = false
int k = 0
while (k < n && !dupl) ; check earlier roots against nth
if (|xx -proots.m_Elements[k]| < DuplThreshold) ; nth is duplicate
nroots = nroots - 1
dupl = true
if (n < nroots) ; need to move down
int m = n
while (m < nroots)
proots.m_Elements[m] = proots.m_Elements[m+1]
m = m + 1
endwhile
endif
endif
k = k+1
endwhile
endwhile
return nroots
endfunc
; Helper function to eliminate complex roots.
static int func RealOnly(int nroots, ComplexArray &proots)
float RealThreshold = 1.e-6
int n = nroots
while (n > 0)
n = n - 1
if (abs(imag(proots.m_Elements[n])) <= RealThreshold) ; real
proots.m_Elements[n] = real(proots.m_Elements[n])
else ; not real
nroots = nroots - 1
if (n < nroots) ; need to move down
int m = n
while (m < nroots)
proots.m_Elements[m] = proots.m_Elements[m+1]
m = m + 1
endwhile
endif
endif
endwhile
return nroots
endfunc
; Power series for Bessel Function J0.
; <p>
; See Abramowitz & Stegun, Handbook of Mathematical Functions, p. 360.
; <p>
; Number of terms to use determined experimentally.
;
; @param z complex argument
; @return Bessel function of order 0 via power series for small argument
static complex func ComplexBesselJ0_Inner(const complex z)
int Nmax = 3 + round(6*sqrt(cabs(z)))
complex zz = sqr(z/2)
complex term = 1
complex sum = term
int n = 1
while (n <= Nmax)
term = -term * zz / sqr(n)
sum = sum + term
n = n + 1
endwhile
return sum
endfunc
; Function with appropriate asymptotic behavior for Bessel Function J0.
; <p>
; See Abramowitz & Stegun, Handbook of Mathematical Functions, p. 364.
; <p>
; Two terms of P0 and Q0 are used, rather arbitrarily.
;
; @param z complex argument
; @return "Bessel" function of order 0 via asymptotic series for large argument
static complex func ComplexBesselJ0_Outer(const complex z)
complex z1 = 1/(8*z)
complex z2 = sqr(z1)
complex arg = z - #pi/4
return sqrt(2/#pi/z)*((1-z2*9/2)*cos(arg)+z1*(1-z2*9*25/6)*sin(arg))
endfunc
; Power series for Bessel Function J1.
; <p>
; See Abramowitz & Stegun, Handbook of Mathematical Functions, p. 360.
; <p>
; Number of terms to use determined experimentally.
;
; @param z complex argument
; @return Bessel function of order 1 via power series for small argument
static complex func ComplexBesselJ1_Inner(const complex z)
int Nmax = 3 + round(6*sqrt(cabs(z)))
complex zz = sqr(z/2)
complex term = z/2
complex sum = term
int n = 1
while (n <= Nmax)
term = -term * zz / (n*(n+1))
sum = sum + term
n = n + 1
endwhile
return sum
endfunc
; Function with appropriate asymptotic behavior for Bessel Function J1.
; <p>
; See Abramowitz & Stegun, Handbook of Mathematical Functions, p. 364.
; <p>
; Two terms of P1 and Q1 are used, rather arbitrarily.
;
; @param z complex argument
; @return "Bessel" function of order 1 via asymptotic series for large argument
static complex func ComplexBesselJ1_Outer(const complex z)
complex z1 = 1/(8*z)
complex z2 = sqr(z1)
complex arg = z - #pi*3/4
return sqrt(2/#pi/z)*((1-z2*3*5/2)*cos(arg)+3*z1*(1-z2*5*21/6)*sin(arg))
endfunc
}
| Constructor Summary | |
|---|---|
Math()
Constructor |
|
| Method Summary | |
|---|---|
static complex |
ComplexBesselJ0(complex z)
Complex Bessel J0. |
static complex |
ComplexBesselJ1(complex z)
Complex Bessel J1. |
static complex |
ComplexBeta(complex a,
complex b)
Complex Beta() David Makin June 2008 |
static complex |
ComplexGamma(complex pz)
Complex Gamma() David Makin June 2008 |
static float |
MaxF(float px,
float py)
Determine the maximum of two numbers. |
static float |
MinF(float px,
float py)
Determine the minimum of two numbers. |
static float |
RealBesselJ0(float x)
Real Bessel J0. |
static float |
RealBesselJ1(float x)
Real Bessel J1. |
static float |
RealBeta(float a,
float b)
Real Beta() David Makin June 2008 |
static float |
RealGamma(float x)
Real Gamma() David Makin June 2008 |
static float |
Sgn(float px)
Returns the sign of a number. |
static int |
SolveComplexCubicPolynomial(complex pa1,
complex pa2,
complex pa3,
complex pa4,
int pflags,
ComplexArray proots)
Solve (complex) cubic equation. |
static int |
SolveComplexQuadraticPolynomial(complex pa1,
complex pa2,
complex pa3,
int pflags,
ComplexArray proots)
Solve (complex) quadratic equation. |
static int |
SolveComplexQuarticPolynomial(complex pa1,
complex pa2,
complex pa3,
complex pa4,
complex pa5,
int pflags,
ComplexArray proots)
Solve (complex) quartic equation. |
static int |
SolveRealCubicPolynomial(float pa1,
float pa2,
float pa3,
float pa4,
int pflags,
ComplexArray proots)
Solve (real) cubic equation. |
static int |
SolveRealQuadraticPolynomial(float pa1,
float pa2,
float pa3,
int pflags,
ComplexArray proots)
Solve (real) quadratic equation. |
static int |
SolveRealQuarticPolynomial(float pa1,
float pa2,
float pa3,
float pa4,
float pa5,
int pflags,
ComplexArray proots)
Solve (real) quartic equation. |
| Methods inherited from class Object |
|---|
|
| Constructor Detail |
|---|
public Math()
You should never need to actually construct a Math object, as all of its methods are static and can be called without an object (just use the class name).
| Method Detail |
|---|
public static float Sgn(float px)
px - the number to test
public static float MaxF(float px,
float py)
px - first numberpy - second number
public static float MinF(float px,
float py)
px - first numberpy - second number
public static int SolveRealQuadraticPolynomial(float pa1,
float pa2,
float pa3,
int pflags,
ComplexArray proots)
See the class overview for details on the solvers. This solves a generic quadratic equation of the form
(a1)(x^2) + (a2)(x) + (a3) = 0
pa1 - a1 termpa2 - a2 termpa3 - a3 termpflags - optionsproots - reference to ComplexArray object to hold roots; this should already be sized to at least two elements (regardless of whether return values are limited)
public static int SolveComplexQuadraticPolynomial(complex pa1,
complex pa2,
complex pa3,
int pflags,
ComplexArray proots)
See the class overview for details on the solvers. This solves a generic quadratic equation of the form
(a1)(x^2) + (a2)(x) + (a3) = 0
pa1 - a1 termpa2 - a2 termpa3 - a3 termpflags - optionsproots - reference to ComplexArray object to hold roots; this should already be sized to at least two elements (regardless of whether return values are limited)
public static int SolveRealCubicPolynomial(float pa1,
float pa2,
float pa3,
float pa4,
int pflags,
ComplexArray proots)
See the class overview for details on the solvers. This solves a generic cubic equation of the form
(a1)(x^3) + (a2)(x^2) + (a3)(x) + (a4) = 0
pa1 - a1 termpa2 - a2 termpa3 - a3 termpa4 - a4 termpflags - optionsproots - reference to ComplexArray object to hold roots; this should already be sized to at least three elements (regardless of whether return values are limited)
public static int SolveComplexCubicPolynomial(complex pa1,
complex pa2,
complex pa3,
complex pa4,
int pflags,
ComplexArray proots)
See the class overview for details on the solvers. This solves a generic cubic equation of the form
(a1)(x^3) + (a2)(x^2) + (a3)(x) + (a4) = 0
pa1 - a1 termpa2 - a2 termpa3 - a3 termpa4 - a4 termpflags - optionsproots - reference to ComplexArray object to hold roots; this should already be sized to at least three elements (regardless of whether return values are limited)
public static int SolveRealQuarticPolynomial(float pa1,
float pa2,
float pa3,
float pa4,
float pa5,
int pflags,
ComplexArray proots)
See the class overview for details on the solvers. This solves a generic quartic equation of the form
(a1)(x^4) + (a2)(x^3) + (a3)(x^2) + (a4)(x) + (a5) = 0
pa1 - a1 termpa2 - a2 termpa3 - a3 termpa4 - a4 termpa5 - a5 termpflags - optionsproots - reference to ComplexArray object to hold roots; this should already be sized to at least four (regardless of whether return values are limited)
public static int SolveComplexQuarticPolynomial(complex pa1,
complex pa2,
complex pa3,
complex pa4,
complex pa5,
int pflags,
ComplexArray proots)
See the class overview for details on the solvers. This solves a generic quartic equation of the form
(a1)(x^4) + (a2)(x^3) + (a3)(x^2) + (a4)(x) + (a5) = 0
pa1 - a1 termpa2 - a2 termpa3 - a3 termpa4 - a4 termpa5 - a5 termpflags - optionsproots - reference to ComplexArray object to hold roots; this should already be sized to at least four (regardless of whether return values are limited)
public static complex ComplexGamma(complex pz)
pz - complex argument
public static float RealGamma(float x)
x - float argument
public static complex ComplexBeta(complex a,
complex b)
a - complex argumentb - complex argument
public static float RealBeta(float a,
float b)
a - float argumentb - float argument
public static float RealBesselJ0(float x)
x - float argument
public static float RealBesselJ1(float x)
x - float argument
public static complex ComplexBesselJ0(complex z)
z - complex argument
public static complex ComplexBesselJ1(complex z)
z - complex argument
|
|||||||||
| PREV CLASS NEXT CLASS | FRAMES NO FRAMES | ||||||||
| SUMMARY: NESTED | FIELD | CONSTR | METHOD | DETAIL: FIELD | CONSTR | METHOD | ||||||||