A Simplified Algorithm for Solution Classification of the

Document Sample

```					MM Research Preprints,135–145
No. 17, Dec. 1998. Beijing                                                                         135

A Simpliﬁed Algorithm for Solution Classiﬁcation
of the Perspective-three-Point Problem1)

Lu Yang
Chengdu Institute of Computer Applications

Abstract.
Given the distance between every pair of 3 control points, and given the angle to every
pair of the control points from an additional point called the Center of Perspectivity
(P ), ﬁnd the lengths of the segments joining the P to each of the control points. We call
this the “perspective-three-point (p3p) problem”. The corresponding algebraic equation
system is called the p3p equation system. So-called a solution classiﬁcation of the p3p
equation system is to give explicit conditions under which the system has none, one,
two, ... real physical solutions, respectively. This problem had been open for many years
[?, 1-8]. An eﬃcient algorithm is presented here for such a classiﬁcation under some
non-degenerate conditions that gives:
• The condition under which the p3p equation system has a unique real physical
solution.
• The conditions under which the p3p equation system has 2, 3, ... real physical
solutions, respectively, provided the system has any real physical solution.
The above result gives the complete classiﬁcaton (under some non-degenerate conditions)
in practice because at least one real solution exists in that case. A program written in
Maple was implemented on PC computers eﬃciently.

The problem is originated from camera calibration [?]. Given the distance between every
pair of 3 control points, and given the angle to every pair of the control points from an
additional point called the Center of Perspectivity (P ), ﬁnd the lengths of the segments
joining the P to each of the control points. We call this the “perspective-three-point (p3p)
problem”. The corresponding algebraic equation system
 2   2          2
y + z − 2pyz − a = 0
z 2 + x2 − 2 q z x − b2 = 0                                 (1)
x2 + y 2 − 2 r x y − c2 = 0


is called the p3p equation system, where a, b, c denote the distances between the control
points A, B, C, parameters p, q, r denote the cosines of the angles

α = BP C,        β = CP A,      γ = AP B,
1)
This work is supported in part by CAS and CNRS under a cooperative project between CICA and
LEIBNIZ, and in part by the National Climbing Project and “863” Project of China.
136                                              Lu Yang

and x, y, z denote the lengths of the segments P A, P B, P C, respectively.
So-called a solution classiﬁcation of the p3p equation system is to give explicit conditions
under which the system (1) has none, one, two, ... real physical solutions, respectively.
Espesially, one regards when the system (1) has a unique real physical solution. This problem
had been open for many years [?, ?].
An eﬃcient algorithm is presented here for such a classiﬁcation (with some non-degenerate
conditions) that gives:
• The condition under which the p3p equation system has a unique real physical solution.
• The conditions under which the p3p equation system has 2, 3, ... real physical solutions,
respectively, provided the system has any real physical solution.
The above result gives the complete classiﬁcaton (with some non-degenerate conditions)
for the camera calibration in practice because at least one real solution exists in that case.
It was pointed out that the classiﬁcation would be given by employing the Sturm-Habicht
sequence or the complete discrimination system for parametric polynomials [?], but one has
to meet with very large polynomials in that way. I guess this is why such a classiﬁcation was
not given before.
Now, by variable changes
√
¯
x = z x,          ¯
y = z y,    z = 1/ z ,
¯
we transform the system (1) into the equivalent one:
¯2         2¯

1 + y − 2py − a z = 0
¯
1 + x2 − 2 q x − b2 z = 0
¯         ¯      ¯                                (2)
x2 + y 2 − 2 r x y − c2 z = 0.

¯     ¯        ¯¯       ¯
Using Wu’s method [?], we ﬁnd the characteristic set of (2):
CS = { h1 = 0, h2 = 0, h3 = 0 }                            (3)
where
h1 = (−b4 − c4 − a4 + 2 b2 c2 + 4 a2 r2 b2 − 2 b2 a2 + 2 a2 c2 ) x4
¯
+ (−4 b2 c2 q + 4 b4 p r − 4 a2 p b2 r + 4 b2 a2 q − 4 c2 p b2 r
+ 4 c4 q + 4 a4 q − 8 a2 c2 q − 8 a2 q r2 b2 ) x3 + (−4 b4 r2
¯
− 4 c4 q 2 − 4 b4 p2 − 4 a4 q 2 + 4 b2 p2 c2 − 2 a4 − 2 c4
+ 8 a2 q p b2 r + 2 b4 + 8 c2 q p b2 r + 4 a2 c2 + 8 a2 q 2 c2
+ 4 a2 r2 b2 ) x2 + (−4 c2 p b2 r + 4 a4 q + 4 b4 p r − 8 a2 c2 q
¯
+ 4 b2 c2 q − 8 b2 p2 c2 q − 4 b2 a2 q − 4 a2 p b2 r + 4 c4 q) x
¯
− a4 − b4 − c4 + 2 b2 a2 + 2 a2 c2 − 2 b2 c2 + 4 b2 p2 c2 ,
h2 = −2 b2 (r x − p) y + (a2 − c2 + b2 ) x2 − 2 q (a2 − c2 ) x
¯      ¯                   ¯                   ¯
− b2 + a2 − c2 ,
h3 = 1 + x2 − 2 q x − b2 z .
¯        ¯      ¯
P3P Problem                                         137

¯ ¯ ¯
We will investigate the solutions to system (3) instead of system (1), where x, y , z are
regarded as unknowns while a, b, c, p, q, r as parameters. Since h3 = 0 always has a positive
¯           ¯
solution for z whenever x takes any real value, we need only investigate how many values of
¯
x solving equation h1 = 0 make the equation h2 = 0 has positive solutions for the unknown
¯
y.
The discriminant of h1 , denoted by D4 , is a polynomial in parameters a, b, c, p, q, r, which
is used to determine whether h1 has repeated root(s) or not. By computing we have

D4 = ((a2 + b2 − c2 ) p2 − 2 (a2 − c2 ) p q r − (b2 + c2 − a2 ) r2 )2 · DD

where DD is a polynomial of 617 terms in a, b, c, p, q, r.
We assume b = max{a, b, c} without loss of generality, and set the following non-degenerate
conditions (n.d.c.):

• | cos α| = | cos A|, | cos β| = | cos B|, | cos γ| = | cos C|.

• Amongest p, q, r, at least two are non-vanishing.

• D4 = 0.

It is easy to see that | cos α| = | cos A| means

−a4 − b4 − c4 + 2 b2 a2 + 2 a2 c2 − 2 b2 c2 + 4 b2 p2 c2 = 0,

and | cos γ| = | cos C| means

−b4 − c4 − a4 + 2 b2 c2 + 4 a2 r2 b2 − 2 b2 a2 + 2 a2 c2 = 0,

i.e., the leading coeﬃcient of h1 is non-vanishing.              The n.d.c will be supposed holding
throughout this talk.
¯
Solving h2 = 0 for y , we have

(a2 − c2 + b2 ) x2 − 2 q (a2 − c2 ) x − b2 + a2 − c2
¯                   ¯
¯
y=                         2 (r x − p)
,
2b      ¯

¯
so, to make clear when this equation has positive solution(s) for y , we need consider the
roots of both equations r x − p = 0 and
¯

(a2 − c2 + b2 ) x2 − 2 q (a2 − c2 ) x − b2 + a2 − c2 = 0.
¯                   ¯                                       (4)
p
The former has a root x1 = , and equation (4) has two real roots, one positive and one
r
negative, because b = max{a, b, c}. The unique positive root of (4) is

q (a2 − c2 ) +     b4 − (1 − q 2 ) (a2 − c2 )2
x2 =                                                  .
a2 − c2 + b2
Let n0 , n1 , n2 denote the numbers of the real roots of h1 in intervals (0, ∞), (x1 , ∞), (x2 , ∞),
respectively, and τ the number of the positive solutions of system (3), i.e., the number of
138                                                   Lu Yang

the real physical solutions of system (1) with n.d.c. Clearly τ can be denoted in terms of
n0 , n1 and n2 . For example, if p > 0, r > 0, then

τ = n0 − |n1 − n2 |.

It is obvious that n1 and n2 equal the numbers of the positive roots of polynomials
∆                        ∆
x        x               x        x
f1 (¯) = h1 (¯ + x1) and f2 (¯) = h1 (¯ + x2), respectively.
x       x       x
The coeﬃcient lists of h1 (¯), f1 (¯), f2 (¯) are denoted by

[u0 , u1 , u2 , u3 , u4 ], [v0 , v1 , v2 , v3 , v4 ], [w0 , w1 , w2 , w3 , w4 ],

respectively, where

¯                      ¯                      ¯
ui = coeﬀ (h1 , x, i), vi = coeﬀ (f1 , x, i), wi = coeﬀ (f2 , x, i),

for i = 0, 1, 2, 3, 4. And let uu, vv, ww denote the numbers of sign changes of these coeﬃcient
lists, respectively. Then, we have
Theorem 1. Assume system (3) has positive solution(s), the n.d.c hold and DD > 0.
Deﬁne a new function ν as follows:

• ν = uu − |vv − ww|, if p > 0 and r > 0;

• ν = ww, if p ≤ 0 and r ≥ 0;

• ν = uu − ww, if p ≥ 0 and r ≤ 0;

• ν = |vv − ww|, if p < 0 and r < 0.

Then ν = τ , the number of the positive solutions of (3).
The proof is not diﬃcult. DD > 0 implies D4 > 0, by a well-known result, that means
x                                    x          x
h1 (¯) has 4 real roots, and so have f1 (¯) and f2 (¯). Using the Descartes’s rule of signs, we
have
n0 = uu, n1 = vv, n2 = ww,                                    (5)
from which the conclusion can be inferred easily.
Theorem 2. Assume system (3) has positive solution(s), the n.d.c hold and DD < 0.
Deﬁne a new function µ as follows:

• µ = −sign(u4 ) sign(u0 ) sign(w0 ), if p > 0 and r > 0;

• µ = sign(w4 ) sign(w0 ), if p ≤ 0 and r ≥ 0;

• µ = sign(u0 ) sign(w0 ), if p ≥ 0 and r ≤ 0;

• µ = −sign(w0 ), if p < 0 and r < 0.
P3P Problem                                           139

Then, the system (3) has a unique positive solution if µ < 0, otherwise, it has two positive
solutions.
Since DD < 0, there should be one or two solution(s) in this case. It only depends on the
parities of uu, vv, ww, and these parities are determined by the signs of u0 , u4 , v0 , v4 , w0 , w4 .

Corollary 1. Assume the n.d.c. hold. The system (1) has a unique real physical solution
if and only if
((DD > 0) ∧ (ν = 1)) ∨ ((DD < 0) ∧ (µ = −1)).                       (6)
Corollary 2. Assume the n.d.c. hold. The system (1) has 3 real physical solution if and
only if
(DD > 0) ∧ (ν = 3).                                    (7)
One may wonder whether the conditions given in the above theorems and corollaries
(except Theorem 2) are “explicit” enough, because some of them are in terms of uu, vv, ww,
the numbers of sign changes of the coeﬃcient lists of polynomials h1 , f1 , f2 . In fact, these
conditions can be easily decomposed into atomic formulas. It may be rather lengthy in that
way, of course.
Example 1. Assume p ≤ 0, r ≥ 0, and the n.d.c. hold. The system (1) has a unique
real physical solution if and only if

[DD > 0, w0 > 0, w1 ≥ 0, w2 ≥ 0, w3 ≥ 0, w4 < 0 ]
∨ [DD > 0, w0 > 0, w1 ≥ 0, w2 ≥ 0, w3 ≤ 0, w4 < 0 ]
∨ [DD > 0, w0 > 0, w1 ≥ 0, w2 ≤ 0, w3 ≤ 0, w4 < 0 ]
∨ [DD > 0, w0 > 0, w1 ≤ 0, w2 ≤ 0, w3 ≤ 0, w4 < 0 ]
∨ [DD > 0, w0 < 0, w1 ≤ 0, w2 ≤ 0, w3 ≤ 0, w4 > 0 ]
∨ [DD > 0, w0 < 0, w1 ≤ 0, w2 ≤ 0, w3 ≥ 0, w4 > 0 ]
∨ [DD > 0, w0 < 0, w1 ≤ 0, w2 ≥ 0, w3 ≥ 0, w4 > 0 ]
∨ [DD > 0, w0 < 0, w1 ≥ 0, w2 ≥ 0, w3 ≥ 0, w4 > 0 ]
∨ [DD < 0, w0 w4 < 0 ],

x
where wi is the coeﬃcent of the i-th term of f2 (¯).
A progam called “P3P” was written according to the two theorems. When the parameters
a, b, c, p, q, r are speciﬁed as constants, the program gives the number of the real physical
solutions promptly.

References

[1] Abidi, M. A., and Chandra, T., A New Eﬃcient and Direct Solution for Pose Estimation Using
Quadrangular Targets: Algorithm and Evaluation, IEEE Transaction on Pattern Analysis and
Machine Intelligence, Vol.17, No.5, pp.534-538, May 1995.
[2] DeMenthon, D., and Davis, L. S., Exact and Approximate Solutions of the Perspective-Three-
Point, IEEE Transaction on Pattern Analysis and Machine Intelligence, Vol.14, No.11, pp.1100-
1105, November 1992.
140                                           Lu Yang

[3] Fishler, M. A., and Bolles, R. C., Random Sample Consensus: A Paradigm for Model Fitting
with Applications to Image Analysis and Automated Cartography, Communications of the ACM,
Vol.24, No.6, pp.381-395, June 1981.
[4] X. S. Gao and H. F. Cheng, On the Solution Classiﬁcation of the “P3P” Problem, Proceedings
of the Third Asian Symposium on Computer Mathematics, eds Z. B. Li, pp. 185-200, LanZhou
University Press, 1998.
[5] Haralick., R. M., Chung-nan Lee, Karsten Ottenberg and Michael Nolle, Analysis and Solutions
of The Three Point Perspective Pose Estimation Problem, Proc. of the Int. Conf. on Computer
Vision and Pattern Recognition, 1991.
[6] Cheng Su, Yingqing Xu, Hua Li, Shenquan Liu, Wu’s Method’s Application in Computer Ani-
mation, The Fifth Int. Conf. on CAD/CG, Vol.1, pp.211-215, Shenzhen, China, 1997.
[7] Wolfe, W. J., and Jones, K., Camera calibration using the perspective view of a triangle, Proc.
SPIE Conf. Automated Inspection Measurement, Vol. 730, Cambridge, MA, Oct.28-30, 1986.
[8] Wolfe, W. J., Mathis, D., Weber, C., and Magee, M., The Perspective View of Three Points,
IEEE Transaction on Pattern Analysis and Machine Intelligence, Vol.13, No.1, pp.66-73, January
1991.
[9] Wu, Wentsun, Basic principles of mechanical theorem proving in elementary geometries, J. Sys.
Sci. and Math. Sci. 4(3) (1984), pp.207-235.
[10] Yang, L., Hou, X. R., and Zeng, Z. B., A complete discrimination system for polynomials, Science
in China, Series E, 39:6(1996), 628-646.
[11] Joseph S. -C. Yuan, A General Photogrammetric Method for Determining Object Position and
Orientation, IEEE Transactions on Robotics and Automation, Vol.5, No.2, pp.129-142, April
1989.
P3P Problem                            141

Appendix. an online program in Maple

p3p:=proc(plist)
local pp,a,b,c,p,q,r,dd,DD,w,cc,cd,i,j,ss,vv,nu:
pp:=plist:
if max(pp[1],pp[2],pp[3])=pp[1]
then pp:=subsop(1=pp[2],2=pp[1],4=pp[5],5=pp[4],pp)
elif max(pp[1],pp[2],pp[3])=pp[3]
then pp:=subsop(2=pp[3],3=pp[2],5=pp[6],6=pp[5],pp)
fi:
a:=pp[1]: b:=pp[2]: c:=pp[3]: p:=pp[4]: q:=pp[5]: r:=pp[6]:
print(a,b,c,p,q,r):
if 4*b^2*c^2*p^2-(b^2+c^2-a^2)^2=0 or
4*a^2*b^2*r^2-(a^2+b^2-c^2)^2=0 or
4*a^2*c^2*q^2-(a^2+c^2-b^2)^2=0 or
(a^2+b^2-c^2)*p^2+(2*c^2-2*a^2)*r*q*p+(a^2-c^2-b^2)*r^2=0 or
(p=0 and r=0)
then print(‘degenerate case, please use another program‘):
RETURN()
fi:
DD:=
4*q^2*b^2*c^10-2*r^4*a^12*q^2-4*p^6*b^2*c^10+4*q^2*a^10*b^2-32*q^2*a^6*c^6-4*p^2*
b^6*a^6+16*a^8*q^6*c^4-2*a^10*r^4*c^2+p^4*b^12-a^8*b^4-a^4*b^8+p^4*c^12+r^4*a^12-
a^8*c^4-a^4*c^8+q^4*c^12-b^4*c^8-b^8*c^4+q^4*a^12-24*q^6*a^6*c^6+b^12+c^12+a^12+q
^4*a^8*r^4*c^4+r^4*b^8*p^4*c^4+q^4*p^4*b^4*c^8+r^4*b^12+a^8*q^4*b^4+q^4*p^4*c^12+
r^4*b^8*c^4+q^4*a^12*r^4+p^4*a^4*c^8+r^4*a^8*c^4+r^4*b^12*p^4+p^4*b^8*a^4-24*r^6*
b^6*a^6+16*r^6*b^8*a^4+6*r^4*b^10*a^2+6*a^10*r^4*b^2-4*a^6*q^2*b^6-4*r^6*b^2*a^10
-2*r^2*b^12*p^4-2*r^2*q^4*a^12-2*a^10*q^4*b^2+2*a^4*q^2*b^8-32*b^6*p^2*c^6+6*b^2*
p^4*c^10-33*p^4*b^4*c^8+4*p^2*q^2*c^12-2*q^4*p^2*c^12-33*p^4*b^8*c^4+r^4*b^8*p^4*
a^4+q^4*b^4*c^8+q^4*a^4*p^4*c^8-8*a^8*b^4*r^2*p^2-4*r^4*b^10*a^2*p^2-4*r^6*b^8*a^
4*q^2-4*r^5*b^6*p^3*a^6*q+20*r^2*b^6*p^2*a^6-6*r^4*b^4*p^2*a^8*q^2-6*q^2*a^4*b^8*
r^4*p^2+8*r^3*b^6*p^3*a^6*q+24*a^6*b^6*q*r*p-8*r^2*b^8*a^4*q^2-12*r^6*b^8*a^4*p^2
-16*a^8*b^4*q*r*p+10*r^4*b^8*a^4*q^2+20*a^6*q^2*b^6*r^2+24*a^4*r^4*b^8*p^2-4*a^10
*r^4*b^2*q^2-4*a^10*q^3*b^2*r^5*p+4*r^5*b^2*a^10*p*q+4*a^10*q^2*b^2*r^6+4*r^6*b^
10*a^2*p^2+52*p^4*b^6*c^6+16*p^6*b^8*c^4+4*r^2*q^2*a^12-2*q^2*p^4*c^12-2*r^2*a^12
-2*p^2*c^12+4*b^6*c^6-2*b^10*c^2-2*q^2*c^12-2*a^2*c^10-2*b^2*c^10-2*a^10*c^2+4*a^
6*c^6-2*b^2*a^10+4*b^6*a^6-2*b^10*a^2-2*q^2*a^12-2*p^2*b^12+q^4*a^8*b^4*r^4+12*q^
2*a^6*b^6*r^4*p^2-20*a^6*r^5*b^4*p*q*c^2-20*q^3*p^5*b^2*a^2*r*c^8+4*r^3*b^6*p^5*a
^2*q*c^4-28*q*p^3*b^2*a^2*r*c^8+44*q*p^3*b^2*a^4*r*c^6+16*q^3*p^5*b^6*r*a^2*c^4+
20*r^5*b^2*a^8*q*p*c^2+12*r^2*b^2*p^2*a^8*c^2+4*r^5*b^6*p^3*a^4*q*c^2-12*r^2*b^2*
p^2*a^6*c^4-28*r^4*b^8*a^2*p^2*c^2-56*r^2*b^4*a^6*q^2*c^2-84*q^2*a^4*r^2*b^4*p^2*
c^4+32*a^2*b^6*q*p^3*r*c^4-8*q^4*b^4*a^4*r^4*c^4+78*p^2*b^4*q^2*a^4*c^4+32*r^3*b^
4*a^6*q*p*c^2+6*a^2*b^8*c^2-4*a^6*b^2*c^4-4*a^2*b^4*c^6+6*a^2*b^2*c^8-4*b^6*a^2*c
^4+6*b^2*p^4*a^2*c^8+4*a^6*q^5*b^2*r^3*p*c^4+22*q^4*a^4*b^4*c^4+40*r^2*b^2*a^8*q^
2*c^2+32*a^4*r^3*b^6*q*p*c^2+44*a^4*q*p^3*b^6*r*c^2-8*a^2*q^2*b^8*c^2+8*p^2*b^2*q
^2*a^6*c^4-12*q^2*p^4*b^2*r^2*a^2*c^8+4*a^4*q^2*b^6*c^2+8*a^2*q*p*b^2*r*c^8+16*q^
2*a^6*b^4*c^2-28*q^2*p^4*b^2*a^2*c^8+12*q^2*a^8*p^2*b^2*c^2+48*p^3*b^8*a^2*r^3*q*
c^2+40*a^2*b^8*r^2*p^2*c^2+4*a^2*q^2*b^6*c^4-8*r^2*b^2*a^2*c^8-12*a^4*q*p*b^6*r*c
^2-12*a^4*b^6*p^4*c^2-104*a^4*b^4*r*p^3*q*c^4+40*q^2*p^4*b^4*a^2*c^6-12*a^2*b^8*p
^2*c^2-12*a^2*b^8*r^2*c^2+8*a^2*q*p*b^8*r*c^2-4*a^6*q^4*b^2*c^4-4*r^4*b^4*a^6*c^2
-12*r^2*b^2*a^8*c^2-4*r^4*b^6*a^4*c^2+22*p^4*b^4*a^4*c^4-12*p^2*b^2*a^2*c^8-8*a^8
*p^2*b^2*c^2+4*p^2*b^2*a^6*c^4+4*r^2*b^2*a^4*c^6+16*r^2*b^2*a^6*c^4+8*q^2*a^6*b^2
*c^4-12*r^4*b^6*a^4*q^2*c^2+8*p^2*b^6*a^2*c^4+16*a^2*q^2*b^4*c^6+8*q^2*a^4*b^2*c^
6+6*q^4*a^8*b^2*c^2+8*r^2*b^4*a^6*c^2+16*a^4*b^6*p^2*c^2+32*r^4*b^6*p^2*a^2*c^4-
48*p^3*b^6*r^3*a^2*q*c^4+8*r^2*b^6*a^4*q^2*c^2-12*a^2*q^2*b^6*r^2*c^4+6*r^4*b^8*a
^2*c^2-12*a^2*q^2*b^4*r^2*c^6-4*p^4*b^4*a^2*c^6-8*a^6*q*p^3*b^4*r*c^2+4*a^6*b^4*p
^2*c^2-12*q^4*a^2*r^2*b^2*p^2*c^8-12*r^4*b^6*a^2*c^4-28*r^3*b^8*a^2*q*p*c^2+8*a^6
*b^4*r^2*p^2*c^2-56*a^4*b^6*r^2*p^2*c^2-4*a^2*b^6*p^4*c^4+16*a^4*p^2*b^2*c^6+16*r
^2*b^6*a^2*c^4-32*a^4*q^2*b^4*c^4-16*a^4*q^3*p^3*b^6*r^5*c^2-22*b^4*p^4*a^4*q^2*c
^4-12*q^2*a^2*b^6*p^4*c^4-32*r^2*b^4*a^4*c^4+40*r^4*b^6*a^4*p^2*c^2-12*q^2*a^6*b^
4*p^2*c^2+12*q^2*a^2*b^8*p^2*c^2-12*q^2*a^8*b^2*c^2-8*r^2*b^2*p^4*a^4*c^6+12*q^2*
p^4*b^6*r^2*a^2*c^4-48*a^4*q*p^3*b^6*r^3*c^2+8*p^2*b^6*a^2*q^2*c^4+20*b^8*p^5*a^2
142                                            Lu Yang

*r*q*c^2+24*a^2*q^3*r^3*b^4*p*c^6-12*q^2*a^4*b^6*p^2*c^2+12*q^2*a^4*b^6*r^4*p^2*c
^2+32*q^4*b^4*a^2*r^2*p^4*c^6+32*p^4*b^6*a^4*r^4*q^2*c^2-12*a^6*q*p*b^4*r*c^2+22*
r^4*b^8*p^4*a^2*c^2-8*q^4*b^4*a^4*r^4*p^2*c^4+32*q^2*p^4*b^6*a^2*r^4*c^4+12*r^2*b
^2*p^2*a^2*c^8+6*p^4*b^8*a^2*c^2+40*q^2*p^2*b^2*a^2*c^8-56*q^2*p^2*b^2*a^4*c^6+8*
r^2*b^6*a^4*c^2-12*a^4*q^4*r^2*b^2*c^6+32*q^2*a^6*r^4*b^2*c^4-8*r^3*b^4*p^3*a^2*q
*c^6-4*a^4*b^2*c^6-4*a^4*b^6*c^2+10*a^4*b^4*c^4-33*r^4*b^4*a^8+2*q^2*b^8*c^4-4*q^
2*b^6*c^6-24*p^6*b^6*c^6+8*a^2*b^4*p^2*c^6+44*q*a^6*r^3*b^2*p*c^4+44*r^3*b^6*p*a^
2*q*c^4-56*q^2*p^2*b^4*a^2*c^6+40*q^4*a^6*b^2*r^2*c^4-28*r^2*b^8*p^4*a^2*c^2+22*r
^4*b^4*a^4*c^4-20*p^5*b^8*a^2*r^3*q*c^2+20*q*p^5*b^2*r*a^2*c^8+32*r*b^4*p^3*q*a^2
*c^6+40*p^4*b^6*r^2*a^2*c^4+4*r^2*b^4*a^2*c^6-8*a^2*q^4*b^4*r^2*c^6+24*r^3*b^2*p^
3*a^4*q*c^6-20*q^4*a^4*p^4*b^2*c^6-12*q^2*a^2*b^2*c^8+6*a^8*r^4*b^2*c^2+16*q*p^5*
b^4*r^3*a^2*c^6-4*a^2*q^2*b^8*r^4*c^2-22*q^4*b^4*a^4*r^2*c^4+88*q^3*a^4*r^3*b^4*p
*c^4+78*r^2*b^4*p^2*a^4*c^4+88*a^4*q*p^3*b^4*r^3*c^4-4*r^3*b^6*p*q*c^6-6*q^2*p^4*
b^8*r^2*c^4+54*q^4*b^4*a^4*r^2*p^2*c^4+32*q^4*a^6*b^4*r^2*c^2-28*q^4*a^8*b^2*r^2*
c^2-48*q^3*a^6*b^4*r^3*p*c^2-32*a^4*b^4*p^2*c^4+16*a^6*q^5*r*b^2*p^3*c^4-8*b^10*r
^3*p*a^2*q-12*r^4*b^8*p^2*a^2*q^2*c^2-56*r^2*b^6*p^2*a^2*c^4+16*a^4*q^5*b^2*r^3*p
*c^6-12*b^4*p^4*a^2*r^2*c^6+8*a^4*q^2*b^2*r^2*c^6+78*a^4*q^2*b^4*r^2*c^4+24*p^3*b
^2*r^3*a^6*q*c^4-24*q^2*a^6*r^4*b^2*p^2*c^4-8*r^4*b^2*p^2*a^6*c^4-20*a^4*q^5*r*b^
2*p*c^6-22*p^4*b^4*r^2*a^4*c^4+22*q^4*p^4*b^2*a^2*c^8-20*p^5*b^6*r*a^2*q*c^4+4*a^
2*q^3*p^5*b^4*r*c^6-8*a^4*q^3*p^3*b^6*r^3*c^2+4*a^4*q^5*r*b^2*p^3*c^6+20*r*a^8*q^
5*b^2*p*c^2+32*r^2*b^6*p^4*a^4*c^2-20*p^4*b^6*r^4*a^2*c^4-22*a^4*q^2*b^4*r^4*c^4-
22*r^4*b^4*p^2*a^4*c^4+20*a^2*q^5*r*b^2*p*c^8-16*a^6*q^3*p^3*b^4*r^5*c^2-20*r^5*b
^6*a^4*p*q*c^2+40*r^4*b^4*a^6*q^2*c^2-16*a^2*q^3*p^5*b^4*r^3*c^6-8*r^4*b^4*p^4*a^
4*q^2*c^4+48*a^2*q^3*r*b^2*p^3*c^8+32*q^2*p^4*b^2*a^4*c^6-8*a^4*q^3*r^3*b^2*p*c^6
-12*b^2*p^4*a^4*c^6+12*q^2*a^2*b^8*r^2*c^2-8*r^3*b^4*p*a^2*q*c^6+54*a^4*q^2*p^2*b
^4*r^4*c^4-48*q^3*a^6*r^3*b^2*p*c^4-4*a^4*q^4*b^2*c^6-24*q^2*b^6*a^2*r^4*p^2*c^4-
12*q^4*a^2*b^4*c^6-8*a^2*q^3*r^3*b^4*p^3*c^6-8*q^3*b^6*a^2*r^3*p^3*c^4-12*a^6*r^4
*b^2*c^4-12*q*a^4*b^2*r*p*c^6-12*a^2*q*b^4*r*p*c^6-24*q^4*p^2*b^4*a^2*r^2*c^6+6*q
^4*a^2*b^2*c^8-4*q^4*b^2*a^2*r^2*c^8-12*q*a^6*b^2*r*p*c^4+12*a^2*q^2*b^2*r^2*c^8-
12*a^2*q*b^6*r*p*c^4+56*a^4*q*b^4*r*p*c^4-20*a^8*q^5*b^2*p*r^3*c^2-20*a^6*q^5*r*b
^2*p*c^4+12*q^2*p^4*b^4*a^2*r^2*c^6-16*q^3*p^5*b^6*r^3*a^2*c^4+12*r^4*b^4*p^2*a^6
*q^2*c^2+12*q^4*a^4*r^2*b^2*p^2*c^6-56*a^6*q^2*b^2*r^2*c^4-8*q^4*p^4*b^4*a^4*r^2*
c^4+8*r^2*b^4*p^2*a^2*c^6+4*a^6*q^3*b^4*r^5*p*c^2+16*q^4*p^4*b^2*a^2*r^2*c^8-136*
q^3*b^4*a^4*r^3*p^3*c^4-8*a^6*b^2*p^3*q*r*c^4-8*q^3*b^2*a^4*r^3*p^3*c^6-20*q^3*b^
2*a^8*r^5*p*c^2-4*b^10*r^2*p^2*c^2+24*q^2*p^4*b^4*c^8+4*q*p^5*b^2*r*c^10+20*q^2*p
^2*b^6*c^6+10*q^2*p^4*b^8*c^4-4*q^2*p^4*b^2*c^10-8*q*p^3*b^10*r*c^2-8*q*p^3*b^2*r
*c^10-8*q^2*b^8*p^2*c^4-4*p^6*b^10*c^2-2*r^4*b^10*c^2+18*b^8*p^2*c^4+6*b^2*a^8*c^
2-2*b^2*q^4*c^10-4*b^6*r^2*c^6+16*b^4*p^6*c^8+4*b^10*r^2*c^2+6*p^4*b^10*c^2-4*a^6
*b^4*c^2+16*a^4*q^6*c^8+4*a^2*p^2*c^10+2*r^2*a^4*c^8+6*q^4*a^10*c^2-33*a^4*q^4*c^
8-4*a^2*q^6*c^10+32*p^3*b^8*r*q*c^4-16*b^8*r*p*q*c^4+4*r*b^2*p*q*c^10+4*r*b^10*p*
q*c^2-12*a^6*r^4*b^4*p^2*c^2-48*q^3*a^4*p^3*b^2*r*c^6-4*p^4*b^2*r^2*a^2*c^8+4*p^5
*b^10*r*q*c^2-48*r*b^6*p^3*q*c^6+10*r^2*b^4*p^4*c^8-4*q^3*p^5*b^2*r*c^10+32*r*b^4
*p^3*q*c^8-16*q*p^5*b^4*r*c^8-8*r^2*b^4*p^2*c^8+8*r*b^4*p^5*q^3*c^8-28*r^2*b^6*p^
4*c^6+4*p^6*b^2*q^2*c^10-4*q^2*p^2*b^2*c^10-12*q^2*p^2*b^4*c^8-12*r^2*b^8*p^2*c^4
-4*r^2*b^4*p^6*c^8+24*q*p^5*b^6*r*c^6-16*p^5*b^8*r*q*c^4-2*q^4*p^4*b^2*c^10-12*q^
2*p^6*b^4*c^8-4*q*p^5*b^10*r^3*c^2+20*r^2*b^6*p^2*c^6+24*p^4*b^8*r^2*c^4-4*r^2*b^
10*p^4*c^2-16*r*b^4*p*q*c^8+12*r^2*b^6*p^6*c^6+4*r^2*b^4*p^6*q^2*c^8-4*q*p^5*b^6*
r^3*c^6+8*r*b^2*p^3*q^3*c^10+4*q^4*p^2*b^2*c^10+24*r*b^6*p*q*c^6-12*r^2*b^8*p^6*c
^4+4*p^6*b^10*r^2*c^2+4*q^2*p^6*b^8*r^2*c^4+4*r^4*b^10*p^2*c^2+2*q^2*b^8*r^2*c^4+
12*p^6*b^6*q^2*c^6+8*q*p^5*b^8*r^3*c^4-4*q^2*b^6*r^2*c^6-6*r^2*b^4*p^4*q^2*c^8-2*
q^4*p^2*b^4*c^8+12*r^2*b^6*p^4*q^2*c^6+8*q^3*r*b^4*p*c^8-4*q^3*b^6*r*p*c^6-16*r^3
*b^8*p^3*q*c^4-4*r*b^2*p*q^3*c^10-16*q^3*r*b^4*p^3*c^8+8*q^3*p^3*b^6*r*c^6-4*r^3*
b^10*p*q*c^2-2*r^4*b^8*p^2*c^4+8*r^3*b^6*p^3*q*c^6-4*q^3*p^5*b^6*r*c^6-4*p^6*b^8*
q^2*c^4+8*r^3*b^10*p^3*q*c^2+8*r^3*b^8*p*q*c^4-2*r^4*b^10*p^4*c^2+8*a^8*q*p*b^2*r
*c^2-28*a^2*q*p^3*b^8*r*c^2-28*p^4*b^6*q^2*c^6+16*q^4*p^4*b^4*a^4*r^4*c^4-8*q^2*p
^6*b^6*r^2*c^6+2*r^2*b^4*q^2*c^8+8*a^8*q^5*r^3*p*c^4+8*q^3*a^10*r^3*p*c^2-4*q*p^3
*r*a^6*c^6+4*q^2*p^4*a^2*c^10+4*r^4*a^10*q^2*c^2-8*q^2*a^8*p^2*c^4+24*q^4*a^4*p^2
*c^8+20*q^2*p^2*a^6*c^6-4*q^4*p^2*a^2*c^10-4*r^2*q^2*a^10*c^2+8*r^3*p*a^8*q*c^4-6
*q^4*a^8*r^2*p^2*c^4-2*q^4*a^10*r^4*c^2+12*a^6*q^6*p^2*c^6+4*a^8*q^6*p^2*r^2*c^4-
4*r^3*p*a^6*q*c^6+8*r*p^3*q*a^4*c^8+4*a^2*q^6*p^2*c^10-2*q^4*a^2*p^4*c^10-8*a^6*q
^6*r^2*p^2*c^6-4*q^6*a^10*c^2+2*p^2*a^8*c^4+52*q^4*a^6*c^6-33*q^4*a^8*c^4+6*q^4*a
^2*c^10+18*p^2*b^4*c^8+2*r^2*b^4*c^8+18*q^2*a^8*c^4-4*p^2*a^6*c^6-2*p^4*a^2*c^10+
4*r^2*a^10*c^2+18*q^2*a^4*c^8-4*r^2*a^6*c^6+8*q^3*a^6*p^3*r*c^6-4*r^2*p^2*a^6*c^6
P3P Problem                            143

+12*q^4*a^6*r^2*p^2*c^6-4*a^10*q*r^3*p*c^2+2*r^2*p^2*a^8*c^4-6*q^4*a^4*r^2*p^2*c^
8-4*a^8*q^6*p^2*c^4-4*q^5*a^10*r^3*p*c^2-4*q^5*a^2*r*p^3*c^10-4*a^2*b^10*r^2*p^2+
24*r*p*a^6*q*c^6-16*r*p*a^8*q*c^4+4*r*a^2*q^5*p*c^10-8*r^2*q^2*a^4*c^8-48*r*a^6*q
^3*p*c^6-16*r*a^8*q^5*p*c^4+24*r*a^6*q^5*p*c^6+32*r*q^3*a^8*p*c^4-8*r*a^2*q^3*p*c
^10+4*r*q^5*a^10*p*c^2+32*r*q^3*a^4*p*c^8+4*r*a^10*q*p*c^2-8*r*q^3*a^10*p*c^2-16*
q^3*a^4*r*p^3*c^8+8*q^3*a^2*r*p^3*c^10-28*q^4*a^6*p^2*c^6-4*q^5*a^6*p^3*r*c^6+8*q
^5*a^4*r*p^3*c^8+10*q^4*a^8*p^2*c^4-4*q^2*a^2*b^8*p^4*c^2-24*q^2*a^4*b^6*p^4*r^2*
c^2-4*a^6*q^5*r^3*p*c^6-28*r^2*q^4*a^6*c^6+8*a^6*q^3*r^3*p*c^6+4*a^4*q^6*p^2*r^2*
c^8-16*r*a^4*q^5*p*c^8-16*q^3*a^8*r^3*p*c^4+4*r*q*a^2*p*c^10-16*r*p*a^4*q*c^8-4*r
^2*q^4*a^10*c^2+12*r^2*q^6*a^6*c^6+10*r^2*a^4*q^4*c^8-2*a^8*q^2*r^4*c^4-12*q^2*a^
4*p^2*c^8-12*r^2*a^8*q^6*c^4-2*q^2*p^4*a^4*c^8+20*r^2*q^2*a^6*c^6-12*r^2*q^2*a^8*
c^4-4*r*p^3*q*a^2*c^10-28*a^8*r^3*b^2*q*p*c^2+4*r^2*q^6*a^10*c^2-20*r^5*b^8*p^3*a
^2*q*c^2+24*q^3*b^6*a^2*r^3*p*c^4-24*p^4*b^2*a^4*q^2*r^2*c^6-8*b^4*p^4*a^4*r^4*c^
4-104*r^3*b^4*p*a^4*q*c^4-8*q^3*a^4*b^6*r^3*p*c^2-20*a^2*q^5*p^3*b^2*r*c^8+16*a^2
*q^2*p^4*b^8*r^4*c^2-12*a^4*q^6*p^2*c^8+2*r^2*p^2*a^4*c^8-4*a^2*q^2*p^2*c^10+24*r
^2*a^8*q^4*c^4-4*r^2*a^4*q^6*c^8-8*r^4*b^6*a^2*q^2*c^4-12*a^6*q^4*b^2*p^2*c^4+40*
q^4*a^4*p^2*b^2*c^6-4*a^8*b^2*r^4*p^2*c^2-8*q^2*a^4*b^6*p^4*c^2+48*q^3*a^8*r^3*b^
2*p*c^2-20*q*p^5*b^4*r*a^2*c^6+12*q^4*p^2*b^2*a^6*r^2*c^4+16*q^4*p^2*b^2*a^8*r^4*
c^2+20*r^5*b^8*p*a^2*q*c^2+32*q^4*b^2*a^4*r^2*p^4*c^6+24*q^3*p^3*b^6*a^4*r*c^2+24
*a^6*q^3*b^4*r*p^3*c^2-8*q^3*a^2*b^6*r*p^3*c^4-8*q*a^4*r^3*b^2*p*c^6+22*q^4*a^8*r
^4*b^2*c^2+88*q^3*a^4*b^4*p^3*r*c^4-48*q^3*a^2*r*b^4*p^3*c^6-12*q^2*a^2*b^8*r^2*p
^4*c^2-4*q^4*a^8*p^2*b^2*c^2-28*a^2*q^4*b^2*p^2*c^8+32*q^4*b^2*a^6*r^4*p^2*c^4+32
*q^4*a^6*r^4*b^4*p^2*c^2-28*a^8*q^3*p*b^2*r*c^2-12*r^2*b^2*p^2*a^4*c^6-8*q^3*b^2*
a^6*r^3*p^3*c^4-12*q^2*b^2*a^8*r^4*p^2*c^2+32*a^2*q^4*b^4*p^2*c^6-20*q^4*a^6*r^4*
b^2*c^4-8*q^3*p^3*b^2*a^6*r*c^4-24*q^4*a^6*b^4*r^2*p^2*c^2-28*r*b^2*p*q^3*a^2*c^8
-28*a^8*r^4*b^2*q^2*c^2-8*a^4*q^3*p*b^6*r*c^2-8*q^3*a^2*b^6*r*p*c^4+32*a^6*q^3*b^
2*p*r*c^4-22*q^4*a^4*b^4*p^2*c^4+44*a^2*q^3*r*b^4*p*c^6-16*r^3*b^2*p^3*q^5*a^4*c^
6-16*r^3*b^2*p^3*q^5*a^6*c^4-104*a^4*q^3*r*b^4*p*c^4+16*r^5*b^4*p^3*a^6*q*c^2-20*
q^4*a^6*b^4*r^4*c^2-8*q^4*a^6*b^4*p^2*c^2+32*q^3*a^4*r*b^2*p*c^6-8*q^3*a^6*b^4*r^
3*p^3*c^2-20*q^4*a^2*b^4*p^4*c^6-8*q^4*a^4*b^4*p^4*c^4+44*a^6*q^3*p*b^4*r*c^2-20*
r^4*b^6*p^4*a^4*c^2-12*q^4*a^8*r^2*b^2*p^2*c^2-8*p^3*b^4*a^6*r^3*q*c^2+54*r^2*b^4
*p^4*a^4*q^2*c^4-2*r^2*b^12+16*a^4*q^3*b^6*r^5*p*c^2-12*a^6*q^4*b^4*c^2-32*r^2*b^
6*a^6+16*r^6*b^4*a^8+52*r^4*b^6*a^6-33*r^4*b^8*a^4+18*r^2*b^8*a^4+18*a^8*b^4*r^2+
4*r^2*b^12*p^2+4*p^2*b^10*a^2-2*p^4*b^10*a^2-2*r^4*b^12*p^2+2*a^8*b^4*p^2-4*r^6*b
^10*a^2+4*a^10*q*p*b^2*r-2*r^4*b^10*p^4*a^2-16*a^4*b^8*q*r*p-28*a^6*b^6*r^4*p^2+8
*r^3*b^10*p^3*a^2*q-16*r^3*b^8*p^3*a^4*q+2*p^2*b^8*q^2*a^4+12*r^6*b^6*a^6*p^2+12*
r^6*b^6*a^6*q^2-4*a^8*r^6*b^4*p^2-2*q^4*a^10*r^4*b^2+24*a^6*r^5*b^6*p*q-4*p^3*b^6
*r*a^6*q-4*q^3*a^6*b^6*r^5*p+4*q^2*p^2*b^4*a^8*r^6-12*q^2*a^8*b^4*r^2+8*q^3*a^8*r
^5*b^4*p-16*r^5*b^8*a^4*p*q+8*q^3*a^6*b^6*r^3*p-2*q^4*a^8*b^4*r^2-8*r^6*b^6*p^2*a
^6*q^2+24*a^8*b^4*q^2*r^4-28*a^6*b^6*q^2*r^4+2*p^2*b^4*q^2*a^8+4*r^6*b^8*p^2*a^4*
q^2+8*r^5*b^8*p^3*a^4*q+10*r^4*b^4*p^2*a^8-12*r^2*b^8*p^2*a^4-4*r^5*b^10*p^3*a^2*
q-16*r^5*b^4*a^8*q*p-4*r^2*b^2*a^10*q^2-4*p^3*b^10*r*a^2*q+32*r^3*b^4*a^8*q*p-4*p
^2*b^6*a^6*q^2+8*q^3*a^8*b^4*r*p-4*q^3*a^10*b^2*r*p+8*q^3*a^10*r^3*b^2*p-16*q^3*a
^8*b^4*r^3*p-4*q^3*a^6*b^6*r*p-2*p^4*b^8*r^2*a^4+4*p^4*b^10*r^2*a^2-8*a^10*r^3*b^
2*q*p-12*r^6*b^4*a^8*q^2+4*r^5*b^10*a^2*p*q+8*p^3*b^8*r*a^4*q-48*a^6*r^3*b^6*q*p+
32*r^3*b^8*a^4*q*p+4*b^10*a^2*q*r*p+4*q^4*a^10*b^2*r^2:

if DD=0
then print(‘exists repeated root, please use another program‘):
RETURN()
fi:

cc.0:= [4*b^2*c^2*p^2-(b^2+c^2-a^2)^2, 4*q*c^4+4*q*a^4-8*q*p^2*b^2*
c^2-4*a^2*q*b^2-4*r*b^2*p*a^2+4*q*b^2*c^2-8*q*a^2*c^2-4*r*b^2*p*c^2+4*r*b^4*p, 8*
q^2*a^2*c^2-4*p^2*b^4-2*a^4+4*a^2*c^2-4*q^2*a^4+4*p^2*b^2*c^2-2*c^4-4*b^4*r^2+2*b
^4+8*a^2*q*p*b^2*r+4*a^2*r^2*b^2+8*r*q*b^2*p*c^2-4*q^2*c^4, -4*q*b^2*c^2-8*q*a^2*
r^2*b^2+4*a^2*q*b^2+4*q*c^4-4*r*b^2*p*c^2+4*q*a^4-4*r*b^2*p*a^2-8*q*a^2*c^2+4*r*b
^4*p, 4*a^2*b^2*r^2-(a^2+b^2-c^2)^2]:

cc.1:=
[-1, -r*(2*p*r*q*c^2-2*p*r*q*a^2+a^2*p^2+p^2*b^2-c^2*p^2+r^2*a^2-r^2*b^2-r^2*c^2)
*(r*q*c^2-r*q*a^2+a^2*p+b^2*p-p*c^2-b^2*p*r^2), 2*a^2*r^2*c^2+b^4*r^2-3*p^2*b^4-6
*a^2*p^2*b^2+6*a^2*q*p*b^2*r-2*b^4*r^4-2*r^2*q^2*c^4-3*a^4*p^2-3*p^2*c^4-r^2*c^4+
144                                            Lu Yang

6*r*q*a^4*p+6*a^2*p^2*c^2-r^2*a^4+6*p^2*b^2*c^2+2*b^2*a^2*r^4+4*r^2*b^4*p^2-12*r*
q*a^2*p*c^2+6*r*q*p*c^4-8*r^3*a^2*q*b^2*p+4*r^3*q*b^2*p*c^2+6*b^2*a^2*r^2*p^2-4*b
^2*r^2*p^2*c^2-6*r*q*b^2*p*c^2+4*r^2*q^2*a^2*c^2-2*r^2*q^2*a^4, -r*(2*r^3*a^2*q*b
^2+r^2*c^2*b^2*p-3*r^2*a^2*b^2*p-r^2*b^4*p-b^2*r*q*a^2-c^4*r*q+2*a^2*r*q*c^2-a^4*
r*q+b^2*r*q*c^2+p*c^4+a^4*p-2*b^2*p*c^2-2*a^2*p*c^2+b^4*p+2*a^2*b^2*p),
4*a^2*b^2*r^2-(a^2+b^2-c^2)^2]:

cc.2 := [2*c^6*a^2*p^2-r^2*b^8-p^2*b^8-a^8*r^2+r^2*c^8-6*r*b^2*p*q*c^6-2*q*p*b^6*r
*c^2+6*a^2*p^2*b^4*c^2+6*a^2*r*b^2*p*q*c^4-2*b^2*a^4*r^2*c^2-6*a^6*r*b^2*p*q+6*a^
4*q*p*b^2*r*c^2+4*b^2*a^2*r^2*c^4-6*p^2*b^4*c^4+4*p^2*b^6*c^2+6*q*p*b^4*r*c^4-12*
a^2*r*b^4*p*q*c^2-8*a^2*q^2*r^2*b^2*c^4+6*a^2*r*b^6*p*q+10*a^2*q^2*r^2*b^4*c^2-8*
a^4*q^2*r^2*b^4+4*a^4*q^2*r^2*b^2*c^2+6*a^4*r*b^4*p*q-6*a^2*p^2*b^2*c^4+4*p^2*b^2
*c^6+2*a^4*r^2*b^4-2*a^2*p^2*b^6+2*r^2*b^6*c^2+2*c^2*a^6*r^2-2*c^6*a^2*r^2+a^8*p^
2-c^8*p^2-2*c^2*a^6*p^2+8*q^2*a^8*r^2+2*b^2*a^6*p^2-2*c^6*r^2*b^2-6*q*a^8*r*p+18*
c^4*q^2*a^4*r^2+2*c^8*r*p*q+16*c^2*a^6*q*p*r-22*c^2*q^2*a^6*r^2-2*c^6*q^2*a^2*r^2
-2*c^2*a^2*r^2*b^4+6*q^2*a^6*p^2*c^2-4*b^2*q^2*a^2*p^2*c^4+8*b^2*q^2*a^4*p^2*c^2+
8*b^2*q^3*a^2*r*p*c^4-16*b^2*q^3*a^4*r*p*c^2-24*q^4*a^4*r^2*c^4+24*q^4*a^6*r^2*c^
2+24*q^3*a^4*p*r*c^4-24*q^3*a^6*p*r*c^2-8*q^3*a^2*p*r*c^6+8*q^4*a^2*r^2*c^6+2*b^4
*q^2*a^2*p^2*c^2-2*q^2*r^2*c^8-2*q^2*a^8*p^2-8*q^4*a^8*r^2-12*a^4*q*p*r*c^4+4*b^2
*q^2*r^2*c^6-2*b^4*q^2*r^2*c^4+2*q^2*a^2*p^2*c^6-4*b^2*q^2*a^6*p^2+8*b^2*q^3*a^6*
r*p-2*b^4*q^2*a^4*p^2+8*q^3*a^8*p*r-6*a^4*q^2*p^2*c^4-2*(-r*b^6*p+a^4*r*b^2*p+2*a
^2*q*r^2*b^4-a^2*r*b^4*p+a^2*p^2*b^4*q-4*q^2*a^4*r*b^2*p+2*a^4*p^2*b^2*q+q*a^6*p^
2+r*p*a^6-q*r^2*c^6+r*p*c^6-a^2*r*p*c^4-q*b^4*r^2*c^2+a^2*p^2*q*c^4-a^4*r*p*c^2+4
*q^3*a^6*r^2-2*q*a^6*r^2-3*r*b^2*p*c^4+3*r*b^4*p*c^2-2*a^2*q*r^2*b^2*c^2+2*q*r^2*
b^2*c^4+2*a^2*r*b^2*p*c^2-4*q^2*a^6*p*r-2*a^2*b^2*p^2*q*c^2-2*a^4*p^2*q*c^2+4*q^3
*a^2*r^2*c^4-8*q^3*a^4*r^2*c^2+3*a^4*r^2*q*c^2+4*q^2*a^2*r*b^2*p*c^2+8*a^4*q^2*r*
p*c^2-4*q^2*a^2*r*p*c^4)*(b^4-(1-q^2)*(-c^2+a^2)^2)^(1/2), -r*p*c^8+a^8*r*p-2*a^8
*r^2*q+2*b^8*r*p+2*a^8*q^3*r^2+b^4*q^2*p*r*c^4+q^2*p*r*c^8+b^6*q*r^2*c^2+b^2*q*r^
2*c^6+b^2*r*p*c^6-a^6*r*b^2*p-a^6*b^2*p^2*q-a^2*b^2*q*p^2*c^4-a^2*b^6*p^2*q+a^2*r
*b^6*p-a^8*q^2*r*p+3*r*b^4*p*c^4-5*b^6*r*p*c^2+4*a^6*b^2*r^2*q+2*a^6*q^2*p*r*c^2-
2*a^6*r*p*c^2-3*a^2*r*b^2*p*c^4+2*a^2*q*r^2*c^6-2*b^2*q^2*r*p*c^6-2*a^2*q^3*r^2*c
^6-2*a^2*q^2*r*p*c^6-4*a^2*b^6*q*r^2+2*a^2*r*p*c^6+3*a^4*b^2*r*p*c^2-10*a^4*b^2*q
^2*p*r*c^2-2*a^4*b^4*p^2*q+2*a^4*b^2*q*p^2*c^2-3*a^4*b^4*r*p+5*a^4*b^4*q^2*r*p+12
*a^4*b^2*q^3*r^2*c^2+6*a^6*q*r^2*c^2-6*a^6*b^2*q^3*r^2-6*a^6*q^3*r^2*c^2+2*a^4*b^
4*r^2*q-6*a^4*q*r^2*c^4-7*a^4*b^2*q*r^2*c^2+6*a^4*q^3*r^2*c^4+4*a^6*b^2*q^2*r*p-6
*a^2*b^4*q^2*p*r*c^2+2*a^2*b^4*p^2*q*c^2+8*a^2*b^2*q^2*p*r*c^4-6*a^2*b^2*q^3*r^2*
c^4+2*a^2*b^2*q*r^2*c^4-2*b^4*q*r^2*c^4-(b^6*r^2+p^2*b^6-a^2*r^2*b^4-a^4*r^2*b^2+
3*p^2*b^2*c^4+a^4*p^2*b^2+2*a^2*p^2*b^4+6*q^2*a^4*r^2*b^2-4*r*b^2*p*a^4*q-5*a^2*r
*b^4*p*q-c^6*p^2-a^4*c^2*p^2+a^2*r^2*c^4+b^2*r^2*c^4+2*a^2*c^4*p^2-2*a^4*c^2*r^2-
2*a^6*q^2*r^2-2*c^2*r^2*b^4-3*p^2*b^4*c^2+c^6*q*p*r-a^2*c^4*r*p*q+a^6*r^2-a^4*c^2
*q*p*r+a^6*q*r*p+r*b^4*p*q*c^2+6*a^2*q*p*b^2*r*c^2-4*a^2*p^2*b^2*c^2-2*a^2*c^4*q^
2*r^2-6*a^2*q^2*r^2*b^2*c^2+4*a^4*c^2*q^2*r^2-2*r*b^2*p*q*c^4)*(b^4-(1-q^2)*(-c^2
+a^2)^2)^(1/2), c^8+a^8-b^8-r^2*b^8-a^8*q^2-p^2*b^8-q^2*c^8-4*a^2*c^6+6*a^4*c^4-4
*a^6*c^2-2*b^2*c^6+2*b^6*c^2+2*b^2*a^6-2*b^6*a^2-r*b^2*p*q*c^6-q*p*b^6*r*c^2+4*a^
2*p^2*b^4*c^2+a^2*r*b^2*p*q*c^4+10*b^2*a^4*r^2*c^2-a^6*r*b^2*p*q+a^4*p^2*b^2*c^2+
a^4*q*p*b^2*r*c^2-2*a^6*b^2*q^2+4*a^2*q^2*c^6-5*b^2*a^2*r^2*c^4+2*b^2*q^2*c^6-3*p
^2*b^4*c^4+3*p^2*b^6*c^2+4*a^6*q^2*c^2-6*b^2*a^4*c^2+6*b^2*a^2*c^4+2*q*p*b^4*r*c^
4-6*a^2*r*b^4*p*q*c^2+6*a^2*q^2*r^2*b^2*c^4+5*a^2*r*b^6*p*q+6*a^2*q^2*r^2*b^4*c^2
+2*a^2*b^4*q^2*c^2+6*a^6*q^2*r^2*b^2-6*a^4*q^2*r^2*b^4+6*a^4*b^2*q^2*c^2-12*a^4*q
^2*r^2*b^2*c^2+4*a^4*r*b^4*p*q-2*a^2*p^2*b^2*c^4-6*a^2*b^2*q^2*c^4-r^2*b^4*c^4+p^
2*b^2*c^6-b^4*q^2*c^4-a^4*b^4*q^2-a^4*p^2*b^4+a^4*r^2*b^4-2*a^2*p^2*b^6-6*a^4*q^2
*c^4-5*b^2*a^6*r^2+5*b^6*a^2*r^2+2*r^2*b^6*c^2+(-3*r*b^2*(a^2+b^2-c^2)*(a^2+c^2-b
^2)*p+6*b^2*a^2*(a^2-c^2-b^2)*r^2*q)*(b^4-(1-q^2)*(-c^2+a^2)^2)^(1/2), -(2*a*r*b-
c^2+b^2+a^2)*(-2*a*r*b-c^2+b^2+a^2)*(b^4-(1-q^2)*(-c^2+a^2)^2)^(1/2)-r*b^2*(a^2+b
^2-c^2)*(a^2+c^2-b^2)*p+2*b^2*a^2*(a^2-c^2-b^2)*r^2*q, 4*a^2*b^2*r^2-(a^2+b^2-c^2
)^2]:

if DD<0
then if p>0 and r>0
then cd := -sign(cc.0[5])*sign(cc.0[1])*sign(evalf(cc.2[1]))
elif p<=0 and r>=0
then cd := sign(cc.2[5])*sign(evalf(cc.2[1]))
elif p>=0 and r<=0
P3P Problem   145

then cd := sign(cc.0[1])*sign(evalf(cc.2[1]))
elif p<0 and r<0
then cd := -sign(evalf(cc.2[1]))
fi:
if cd<0
then print(‘there is 1 real physical solution‘)
else print(‘there are 2 real physical solutions‘)
fi
else for i from 0 to 2 do
ss.i := map(sign,evalf(subs(0=NULL,cc.i))):
vv.i := 0:
for j to nops(ss.i)-1 do
if ss.i[j]+ss.i[j+1]=0
then vv.i := vv.i+1
fi
od:
print(vv.i):
od:
if p>0 and r>0
then nu := vv0-abs(vv1-vv2)
elif p<=0 and r>=0
then nu := vv2
elif p>=0 and r<=0
then nu := vv0-vv2
elif p<0 and r<0
then nu := abs(vv1-vv2)
fi:
print(‘there be ‘.nu.‘ real physical solution(s)‘)
fi
end:
lprint(‘THIS PROGRAM IS VALID TO PRACTICAL DATA‘):