|
Duplicate :
|
l.com/FolkThomasF/MATHTXT.HTM
/**
*main function calls testing routine
*/
public static final void main(String[] args)
{ cmplx.testrtn(); }
private static void testrtn()
{ cmplx c1=null,c2=null,c3=null,c4=null,c5=null;
double a=0,b=0,c=0,d=0;
cmplx[] res=null;
int ii=0,jj=0,nr=0;
System.out.println("testing complex manipulation");
c1=new cmplx(-117,44);
dispcmplx(c1);
c2=new cmplx(3,4);
dispcmplx(c2);
System.out.println("sum "+cmplx.add (c1,c2));
System.out.println("diff "+cmplx.subtract(c1,c2));
System.out.println("prod "+cmplx.multiply(c1,c2));
System.out.println("quot "+cmplx.divide (c1,c2));
dispcmplx(new cmplx(2,0));
dispcmplx(new cmplx(0,2));
dispcmplx(new cmplx(-2,0));
dispcmplx(new cmplx(0,-2));
dispcmplx(new cmplx(0,0));
dispcmplx(new cmplx(1,1));
dispcmplx(c2);
System.out.println("pow 2 "+cmplx.pow(c2,2));
System.out.println("pow 3 "+cmplx.pow(c2,3));
System.out.println("pow -1 "+cmplx.pow(c2,-1));
System.out.println("pow -2 "+cmplx.pow(c2,-2));
System.out.println("pow 0 "+cmplx.pow(c2,0));
System.out.println("zero pow 2 "+cmplx.pow
(cmplx.zero,2));
System.out.println("zero pow 0 "+cmplx.pow
(cmplx.zero,0));
dispcmplx(c1);
System.out.println("pow 1/3 "+cmplx.pow(c1,1.0/3.0));
System.out.println("zero");
dispcmplx(cmplx.zero);
System.out.println("one");
dispcmplx(cmplx.one);
System.out.println("prim cb root of 1");
dispcmplx(cmplx.primcbroot);
System.out.println("other prim cb root of 1");
dispcmplx(cmplx.primcbrootsq);
dispcmplx(c1);
System.out.println("conj "+cmplx.conjugate(c1));
res=new cmplx[4];
for (ii=0;ii<5;ii++)
{ switch(ii)
{ case 0: a=3; b=-3; c=-60; break;
case 1: a=0; b=5; c=4; break;
case 2: a=2; b=0; c=18; break;
case 3: a=0; b=6; c=9; break;
case 4: a=0; b=0; c=3; break; }
System.out.println("coefficients "+a+","+b+","+c);
nr=cmplx.solveQuad(a,b,c,res);
System.out.println("solvequad num roots "+nr);
for (jj=0;jj<nr;jj++) System.out.println(res
[jj]); }
for (ii=0;ii<22;ii++)
{ switch(ii)
{ case 0: a=0; b=0; c=0; d=0; break;
case 1: a=0; b=0; c=0; d=4; break;
case 2: a=0; b=0; c=2; d=3; break;
case 3: a=0; b=0; c=1; d=-8; break;
case 4: a=0; b=1; c=2; d=1; break;
case 5: a=0; b=1; c=4; d=4; break;
case 6: a=0; b=1; c=-4;d=4; break;
case 7: a=0; b=1; c=0; d=9; break;
case 8: a=0; b=1; c=0; d=-16; break;
case 9: a=0; b=1; c=-5;d=6; break;
case 10:a=0; b=4; c=0; d=-36; break;
case 11:a=0; b=10;c=3; d=-1; break;
case 12:a=1; b=0; c=0; d=27; break;
case 13:a=1; b=0; c=0; d=-64;break;
case 14:a=4; b=4; c=0; d=0; break;
case 15:a=2; b=0; c=18;d=0; break;
case 16:a=3; b=0; c=0; d=24; break;
case 17:a=5; b=15;c=15;d=5; break;
case 18:a=7; b=0; c=0; d=0; break;
case 19:a=1; b=2; c=-1;d=-2; break;
case 20:a=4; b=-36; c=0; d=0;break;
case 21:a=1; b=6; c=12;d=-19;break;
}
System.out.println(a+" *y*y*y + "+b+" *y*y + "+c+"
*y + "+d+" =0");
nr=cmplx.solveCubic(a,b,c,d,res);
System.out.println("solvecubic, num roots "+nr);
for (jj=0;jj<nr;jj++) System.out.println(res[jj]);
}
}
private static void dispcmplx(cmplx c)
{ System.out.println(c);
System.out.println(" magsquared "+cmplx.magsquared(c)
+" mag " +cmplx.mag(c)
+" arg " +cmplx.arg(c)
+" arg/pi " +cmplx.arg
(c)/Math.PI); } }
(Review ID: 146693)
======================================================================
Name: jk109818 Date: 08/01/2002
FULL PRODUCT VERSION :
\jdk1.3\bin\java cccont
\jdk1.3\bin\java cmplx
cccont.java:
...
cmplx.java:
...
A DESCRIPTION OF THE PROBLEM :
The CubicCurve2D method solveCubic misses some roots.
For example the equation
x*x*x - 4*x*x == 0
should return 0 and 4 as roots, but the solveCubic
method only returns 4 as a root.
A substitute method, that returns real and
imaginary parts of all three roots is enclosed.
This substitute is needed to process the "contains"
and other methods.
For example, suppose a test point
is located far from the cubic curve. Suppose also
that the point is such that a horizontal
line drawn from that test point to the curve intersects
the curve tangentially. This need not happen near the end
points. It could also happen near the midpoint.
The "contains"
method returns "true" in this case, even though the
answer should be "false" (the point is too far away).
To properly handle these cases, the solveCubic function must
return all roots.
Furthermore, the "inflect" array
is being used incorrectly. The program is currently
calculating the first derivative, or slope. Inflection
requires the second derivative, since inflection
refers to concavity. Concavity determines how tangential
crossings should be handled. A tangential crossing
at a point of inflection, or change of concavity,
behaves as a sigle crossing, whereas
other tangential crossings behave as 0 or 2 crossings.
The enclosed substitute for solveCubic contains
routines for multiplying complex numbers
and raising complex numbers to
given exponents. The "solveCubic" substitute
takes equation coefficients and an array
of complex numbers (to hold the roots) as parameters.
The "main" function includes examples
of how to call the function, and the function has, within
comment boxes, a detailed description of the technique
(it relies on a matrix determinant).
STEPS TO FOLLOW TO REPRODUCE THE PROBLEM :
1. Call the CubicCurve2D solveCubic method
such that the input equation has
three roots, with two of those roots coalescing.
For example: x*x*x -4*x*x == 0
2. Print out the number of roots returned.
3. Separately, create a CubicCurve2D.Double object so that
the line segment joining the end points is
not horizontal but the curve has a horizontal
tangent somewhere along its length and away from
the end points.
4. Call the "contains" method using a test
point far from the curve but such that
the point is on the horizontal line that
intersects the curve tangentially.
EXPECTED VERSUS ACTUAL BEHAVIOR :
The "solveCubic" method fails to return all of the roots,
and the "contains" method returns true for a point that
is far away. A substitute "solveCubic" method, that
finds complex roots, is enclosed.
This bug can be reproduced always.
---------- BEGIN SOURCE ----------
import java.lang.*;
import java.awt.*;
import java.awt.geom.*;
import java.awt.event.*;
import java.applet.*;
public class cccont extends Applet
{ private static boolean standalone=false;
private Frame win=null;
private TextArea ta=null;
private canvcls canvobj=null;
public static final void main(String[] args)
{ standalone=true;
(new cccont()).init(); }
public void init()
{ win=new Frame("cubics");
win.addWindowListener(new WindowAdapter()
{ public void windowClosing(WindowEvent we)
{ we.getWindow().dispose(); }
public void windowClosed(WindowEvent we)
{ if (standalone) System.exit(0); } });
win.setLayout(new FlowLayout(FlowLayout.CENTER,10,10));
ta=new TextArea("",10,60);
win.add(ta);
canvobj=new canvcls();
win.add(canvobj);
win.pack();
win.setSize(600,480);
win.show(); }
//see also http://members.aol.com/FolkThomasF/MATHTXT.HTM
private void dispmsg(String msg)
{ ta.append(msg+"\r\n");
if (standalone) System.out.println(msg); }
private class canvcls extends Canvas
{ private Dimension dim=null;
public canvcls() { dim=new Dimension(500,200); }
public Dimension getMinimumSize() { return dim; }
public Dimension getPreferredSize() { return dim; }
public void paint(Graphics g)
{ g.setColor(Color.green);
g.fillRect(0,0,dim.width,dim.height);
ta.setText("");
try { calculate(g); }
catch(Throwable te) { dispmsg(String.valueOf(te)); } }
private int mapx(double x)
{ return (int)Math.round(dim.width/2+20*x); }
private int mapy(double y)
{ return (int)Math.round(dim.height/3+20*y); }
private void calculate(Graphics g)
{ CubicCurve2D cc=null;
int ii=0,nr=0;
double tx=0,ty=0,xx=0,yy=0,tt=0;
double[] eqn=null,res=null;
boolean contained=false;
cc=new CubicCurve2D.Double(1,1,1,5,5,1,5,5);
tx=-10; ty=3;
dispmsg(" x1 "+cc.getX1()+" y1 "+cc.getY1()
+" x2 "+cc.getX2()+" y2 "+cc.getY2()
+" cpx1 "+cc.getCtrlX1()+" cpy1 "+cc.getCtrlY1()
+" cpx2 "+cc.getCtrlX2()+" cpy2 "+cc.getCtrlY2());
g.setColor(Color.blue);
g.drawLine(mapx(cc.getX1()),mapy(cc.getY1()),
mapx(cc.getX2()),mapy(cc.getY2()));
for (tt=0;tt<=1;tt+=0.02)
{ xx= cc.getX1() *(1-tt)*(1-tt)*(1-tt)
+ 3*cc.getCtrlX1()*(1-tt)*(1-tt)*tt
+ 3*cc.getCtrlX2()*(1-tt)*tt *tt
+ cc.getX2() *tt *tt *tt;
yy= cc.getY1() *(1-tt)*(1-tt)*(1-tt)
+ 3*cc.getCtrlY1()*(1-tt)*(1-tt)*tt
+ 3*cc.getCtrlY2()*(1-tt)*tt *tt
+ cc.getY2() *tt *tt *tt;
if (tt>=0 && tt<=1) g.setColor(Color.blue);
else g.setColor(Color.red);
g.fillOval(mapx(xx)-1,mapy(yy)-1,2,2); }
contained=cc.contains(tx,ty);
dispmsg("test point: "+tx+";"+ty+"; contains "+contained);
g.setColor(Color.blue);
g.fillOval(mapx(tx)-2,mapy(ty)-2,4,4);
g.drawString("contained: "+contained,mapx(0),mapy(-2));
eqn=new double[4];
res=new double[4];
eqn[3]=1; eqn[2]=-4; eqn[1]=0; eqn[0]=0;
nr=CubicCurve2D.solveCubic(eqn,res);
dispmsg("\r\n"+eqn[3]+" *x*x*x + "+eqn[2]+" *x*x + "
+eqn[1]+" *x + " + eqn[0]);
dispmsg("num roots "+nr);
for (ii=0;ii<nr;ii++) dispmsg(String.valueOf(res[ii]));
dispmsg("(misses root(s) at 0)"); } } }
---------- END SOURCE ----------
CUSTOMER WORKAROUND :
import java.lang.*;
/**
*cmplx represents a complex number.
*It has static methods for
*adding, subtracting, multiplying, dividing
*two complex numbers or for raising
*a complex number to a power.
*It also has a static method for solving
*linear, quadratic, or cubic equations.
*/
public class cmplx
{ /**
*real part
*/
public double re=0;
/**
*imaginary part
*/
public double im=0;
private static final double tol=1.0e-100;
private static final double bigtol=1.0e-50;
/**
*primitive cube root of unity
*/
public static final cmplx primcbroot
=new cmplx(- 1.0/2.0, Math.sqrt(3.0)/2 );
/**
*another primitive cube root of unity
*/
public static final cmplx primcbrootsq
=new cmplx(- 1.0/2.0, - Math.sqrt(3.0)/2 );
/**
*zero
*/
public static final cmplx zero =new cmplx(0,0);
/**
*one
*/
public static final cmplx one =new cmplx(1,0);
/**
*construct a complex number
*@param re double precision real part
*@param im double precision imaginary part
*/
public cmplx(double re,double im) { this.re=re;
this.im=im; }
/**
*add two complex numbers
*@param x cmplx first complex number
*@param y cmplx second complex number
*@return cmplx sum of complex numbers
*/
public static cmplx add(cmplx x,cmplx y)
{ return new cmplx(x.re+y.re,x.im+y.im); }
/**
*subtract two complex numbers
*@param x cmplx first complex number
*@param y cmplx second complex number
*@return cmplx difference of two complex numbers (x-y)
*/
public static cmplx subtract(cmplx x,cmplx y)
{ return new cmplx(x.re-y.re,x.im-y.im); }
/**
*multiply two complex numbers
*@param x cmplx first complex number
*@param y cmplx second complex number
*@return cmplx product of two complex numbers
*/
public static cmplx multiply(cmplx x,cmplx y)
{ return new cmplx(x.re * y.re - x.im * y.im,
x.re * y.im + x.im * y.re); }
/**
*square magnitude
*@param x cmplx complex number
*@return double precision magnitude squared
*/
public static double magsquared(cmplx x)
{ return x.re * x.re + x.im * x.im; }
/**
*magnitude
*@param x cmplx complex number
*@return double precision magnitude
*/
public static double mag(cmplx x)
{ return Math.sqrt(cmplx.magsquared(x)); }
/**
*argument or polar angle
*@param x cmplx complex number
*@return double precision argument or polar angle
*/
public static double arg(cmplx x)
{ return Math.atan2(x.im,x.re); }
/**
*complex conjugate
*@param x cmplx complex number
*@return cmplx complex conjugate
*/
public static cmplx conjugate(cmplx x)
{ return new cmplx(x.re, - x.im); }
/**
*divide complex numbers
*@param x cmplx first complex number (dividend)
*@param y cmplx second complex number (divisor)
*@return cmplx ratio of complex numbers (x/y)
*/
public static cmplx divide(cmplx x,cmplx y)
throws IllegalArgumentException
{ double d=0;
d=magsquared(y);
if (d<cmplx.tol)
throw new IllegalArgumentException("divisor too
small");
return new cmplx((x.re * y.re + x.im * y.im)/d,
(- x.re * y.im + x.im * y.re)/d); }
/**
*raise to a power
*@param base cmplx complex number
*@param exp power to raise to
*@return cmplx complex number raised to a power
*/
public static cmplx pow(cmplx base,double exp)
{ double r=0,a=0;
r=cmplx.mag(base);
a=cmplx.arg(base);
r=Math.pow(r,exp);
a=a*exp;
return new cmplx(r*Math.cos(a),r*Math.sin(a)); }
/**
*string representation of complex number
*@return String showing real and imaginary values
*/
public String toString()
{ return this.re + " + i * " + this.im; }
/**
*<PRE>
*Solve an equation of degree at least one (linear)
*and at most three (cubic).
*If the coefficient of the highest term
*is zero then call the solvequad function;
*otherwise divide the equation so that the
*coefficient of the highest term is 1:
*
* 0 == y*y*y + r *y*y + s *y + t
*
*Make a change of variables so that the quadratic term
vanishes.
*
* 0 == x*x*x + p *x + q
*
*where
*
* y=x-r/3
*
*The new coefficients are related to the original by:
*
* p=(s-r*r/3);
* q=(t-r*s/3+2*r*r*r/27);
*
*Set the result equal to the determinant of the matrix
*
* x,a,b
* b,x,a
* a,b,x
*
*The determinant = x*x*x-3*a*b*x+a*a*a+b*b*b, and is
zero for x=-a-b.
*Thus p=-3*a*b and q=a*a*a+b*b*b, so that
*x=-a-b=-a+p/(3*a)
*for a=CubeRoot((1/2)*(q PlusOrMinus SquareRoot(q*q+
(4/27)*p*p*p))).
*Given one cube root, other cube roots may be obtained by
*multiplying by cubic roots of unity.
*</PRE>
*
*@param cubecoef double precision coefficient of
* cubic term
*@param quadcoef double precision coefficient of
* quadratic term
*@param linecoef double precision coefficient of
* linear term
*@param constcoef double precision coefficient of
* constant term
*@param result array of cmplx (complex number)
* objects to return up to three roots of the
equation;
* this parameter may not be null and must
* be an array large enough to hold all the roots
returned
*@param int number of roots returned;
* this number is equal to the degree of the input
equation
*/
public static int solveCubic(double cubecoef,double
quadcoef,
double linecoef,double
constcoef,
cmplx[] result)
{ double r=0,s=0,t=0; /* 0==y*y*y+r*y*y+s*y+t; */
double p=0,q=0; /* 0==x*x*x +p*x+q; */
double dsc=0; /* (q*q+4*p*p*p/27)/4 */
boolean dbgsw=false;
cmplx c1=null,c2=null,
acubed=null,a=null,bcubed=null,b=null,
wa=null,wwa=null,wb=null,wwb=null,minuspdiv3=null,
res1=null,res2=null,res3=null,minusrdiv3=null;
int ii=0;
if (result==null)
throw new IllegalArgumentException("no array
allocated");
if (Math.abs(cubecoef)<tol)
{ return solveQuad
(quadcoef,linecoef,constcoef,result); }
r=quadcoef/cubecoef;
s=linecoef/cubecoef;
t=constcoef/cubecoef;
p=(s-r*r/3);
q=(t-r*s/3+2*r*r*r/27);
dsc=(q*q+4*p*p*p/27)/4;
c1=new cmplx(q/2,0);
c2=cmplx.pow(new cmplx(dsc,0),1.0/2.0);
acubed=cmplx.add(c1,c2);
bcubed=cmplx.subtract(c1,c2);
a=cmplx.pow(acubed,1.0/3.0);
b=cmplx.pow(bcubed,1.0/3.0);
minuspdiv3=new cmplx( - p/3.0, 0);
if (cmplx.magsquared(acubed)>bigtol)
{ b=cmplx.divide(minuspdiv3,a);
if (dbgsw) System.out.println("case one"); }
else
{ if (cmplx.magsquared(b)>bigtol)
{ a=cmplx.divide(minuspdiv3,b);
if (dbgsw) System.out.println("case
two"); }
else
{ a=cmplx.zero; b=cmplx.zero;
if (dbgsw) System.out.println("case
three"); } }
wa=cmplx.multiply(a,cmplx.primcbroot);
wwa=cmplx.multiply(a,cmplx.primcbrootsq);
wb=cmplx.multiply(b,cmplx.primcbroot);
wwb=cmplx.multiply(b,cmplx.primcbrootsq);
minusrdiv3=new cmplx( - r/3, 0);
res1=cmplx.subtract(minusrdiv3,a);
res1=cmplx.subtract(res1,b);
res2=cmplx.subtract(minusrdiv3,wa);
res2=cmplx.subtract(res2,wwb);
res3=cmplx.subtract(minusrdiv3,wb);
res3=cmplx.subtract(res3,wwa);
result[0]=res1;
result[1]=res2;
result[2]=res3;
return 3; }
/**
*solve a quadratic equation
*@param quadcoef double precision coefficient of
quadratic term
*@param linecoef double precision coefficient of linear
term
*@param constcoef double precision coefficient of
constant term
*@return int number of roots;
* this number is equal to the degree of the equation
*/
public static int solveQuad(double quadcoef,double
linecoef,
double constcoef,cmplx[]
result)
{ double a=0,b=0,c=0,dsc=0;
cmplx c1=null,c2=null;
if (Math.abs(quadcoef)<tol)
{ if (Math.abs(linecoef)<tol) return 0;
result[0]=new cmplx( - constcoef/linecoef,0);
return 1; }
a=quadcoef; b=linecoef; c=constcoef;
c1=new cmplx(- b/ (2*a), 0);
dsc= (b*b-4*a*c)/(4*a*a);
c2=cmplx.pow(new cmplx(dsc,0),1.0/2.0);
result[0]=cmplx.add(c1,c2);
result[1]=cmplx.subtract(c1,c2);
return 2; }
//see also http://members.ao