Rectangle inscribed in another rectangle: solutions for all cases. ================================================================== Let X and Y be the sides of a rectangle XY exactly inscribed in another rectangle AB of sides A and B, i.e. each vertex of XY belongs to a side of AB. Let us call phi (0°<=phi<=90°) one of the angles formed by a side Y of XY and a side A of AB. It is graphically obvious that: (1) A = Y cos(phi) + X sin(phi) ( = a1 + a2) (2) B = Y sin(phi) + X cos(phi) ( = b1 + b2) Where we have put a1 = Y cos(phi) a2 = X sin(phi) = A - a1 b1 = Y sin(phi) b2 = X cos(phi) = B - b1 with the obvious meaning of distances between a vertex of AB and one of XY. We will treat three non trivial cases: 1). Given A,B,phi, find X and Y, or state that the problem has no solution. 2). Given A,X,Y, find B (and phi), or state that the problem has no solution. 2.a). Trigonometrical solution 2.b). Analytical solution 3). Given A,B,Y, find X (and phi), or state that the problem has no solution. 3.a) Graphical solution 3.b) Analytical solution 3.c) Numerical solution As a byproduct of solutions, useful for a graphical construction of inner rectangle, also a1,a2,b1 and b2 will be found. Remaining possible cases can be solved trivially: i.e., 0.1). Given X,Y,phi, find A and B: solution is given directly by (1) and (2). 0.2). Given A,Y,phi, find X and B: a1 = Y cos(phi) a2 = A-a1 (if negative, there is no solution) X = a2/sin(phi) B is given by (2). Case 1). Given A,B,phi, find X and Y ==================================== This problem is encountered in practice when a picture has been taken, holding the camera non perfectly right, and one wants to rotate the photo by an angle phi, so that the horizon becomes horizontal; and then, cut the photo's oblique sides, saving as much as possible of the image. This is the simplest case: it is enough to solve the system of equations (1) and (2) in unknowns X and Y. X = [B cos(phi)-A sin(phi)]/[cos(phi)^2-sin(phi)^2], Y = [A cos(phi)-B sin(phi)]/[cos(phi)^2-sin(phi)^2]. If both X and Y turn out non-negative, this is the required solution. And, of course, a1 = Y cos(phi) a2 = X sin(phi) b1 = Y sin(phi) b2 = X cos(phi). In the particular case in which the denominator [cos(phi)^2-sin(phi)^2] is zero, i.e. phi=45°, solutions exist only if A=B, and in this case X can assume any value between 0 and A*sqrt(2); Y is determined accordingly: Y = A*sqrt(2) - X. Case 2). Given A,X,Y, find B (and phi) ====================================== If BOTH X and Y are greater than A, the problem has no solution. If the diagonal of XY ( sqrt(X^2 + Y^2) ) is less than A, the problem has no solution. OTHERWISE: 2.a). Trigonometrical solution ------------------------------ If Y is less than or equal to A, the problem has at least one solution. Proof: Let us draw two parallel straight lines, r and s, at distance A. Let us insert rectangle XY between r and s, with sides X parallel to them. This can be done because Y<=A. Let us rotate XY, until both endpoints of one of its diagonals touch r and s. This can be done because its diagonal is greater than or equal to A. Let us draw two other straight lines, orthogonal to both r and s, passing throught the endpoints of the other diagonal of XY. Their distance is a value for B that satisfies (1) and (2). Now: Let us call theta the acute (or rect) angle between the first diagonal and r (or s): theta = arcsin( A / sqrt(X^2 + Y^2) ) Let us call alpha the angle between the same diagonal and one of the X sides: alpha = arctan(Y/X) It turns out phi = theta - alpha = arcsin(A/sqrt(X^2 + Y^2)) - arctan(Y/X) Once found phi, B can be computed from (2): B = Y sin(phi) + X cos(phi) and obviously a1 = Y cos(phi) a2 = X sin(phi) b1 = Y sin(phi) b2 = X cos(phi). BESIDES: Also if X is less than or equal to A, the problem has at least one solution. Proof proceeds as before: Let us draw two parallel straihgt lines, r and s, at distance A. Let us insert rectangle XY between r and s, this time with sides Y parallel to them. This can be done because X<=A. Let us rotate XY, until both endpoints of one of its diagonals touch r and s. This can be done because its diagonal is greater than or equal to A. Let us draw two other straight lines, orthogonal to both r and s, passing throught the endpoints of the other diagonal of XY. Their distance is a value for B that satisfies (1) and (2). In this case, alpha will be the angle between one diagonal of XY and one of its Y (not X) sides, i.e. the complementary of that seen above; and hence alpha = arctan(X/Y) phi = theta - alpha = arcsin(A/sqrt(X^2 + Y^2)) - arctan(X/Y). Once found another value for phi, another value for B can be computed from (2): B = Y sin(phi) + X cos(phi) and obviously, four other values for a1,a2,b1 and b2: a1 = Y cos(phi) a2 = X sin(phi) b1 = Y sin(phi) b2 = X cos(phi). Thus, if BOTH X and Y are <= A, there are TWO solutions for B (that coincide if X=Y, i.e. XY is a square). Of course, on can choose between them the smallest one. You can find, in the appendix below, a subroutine than implements the described algorithm and a test program that calls it. Both are written in Fortran 95 and can be compiled and run on a linux PC using f95. 2.b). Analytical solution ------------------------- putting cos(phi) = c in (1) (1) A = Y cos(phi) + X sin(phi) we obtain (1') A = Y c + X sqrt(1-c^2) A - Y c = X sqrt(1 - c^2) (A - Yc)^2 = X^2 (1 - c^2) (Y^2+X^2) c^2 - 2AY c +(A^2-C^2) = 0 This is a II degree equation in c that has two solutions: c_1 = {A Y + sqrt[A^2 Y^2 - (A^2-X^2) (X^2+Y^2)]} / (X^2+Y^2) c_2 = {A Y - sqrt[A^2 Y^2 - (A^2-X^2) (X^2+Y^2)]} / (X^2+Y^2) Acceptable solutions must be real (the term under square root must be non-negative) and within 0 and 1 included (phi=arcos(c) must be acute). Under these conditions, c1(X) e c2(X) coincide with cos(phi) found in two subcases of case 2.a) seen above. From both values of c, a value of phi can be found as arcos(c), and then B, and then a1,a2,b1,b2, as above. Caso 3). Given A,B,Y, find X (and phi) ====================================== This is the hardest case. It can be solved graphically, analitycally and numerically. 3.a) Graphical solution 3.b) Analytical solution 3.c) Numerical solution 3.a). Graphical solution ----------------------- If Y is greater than the diagonal of AB, Y > sqrt(A^2+B^2), the problem has no solution. Otherwise, solutions can be one or more. To find them, let us consider the solution found in case 2): once found phi, we obtained B = Y sin(phi) + X cos(phi) a1 = Y cos(phi) b1 = Y sin(phi) we can also write, straightfully, X = a2/sin(phi) = (A-a1)/sin(phi) = (A-Y*cos(phi))/sin(phi) and so B = Y sin(phi) + [(A-Y*cos(phi))/sin(phi)] cos(phi). Let us now consider phi as the independent variable, and plot curves B(phi), a1(phi), b1(phi), X(phi), and the horizontal straight line corresponding to the input value of B. Using gnuplot, we could write (or copy and paste): A = 1 # replace with your input value Y = .8 # replace with your input value Binput = 1 # replace with your input value phi(x)=x # phi = x, i.e. x-axis *is* angle phi X(x) = (A-Y*cos(phi(x)))/sin(phi(x)) B(x) = Y*sin(phi(x)) + X(x)*cos(phi(x)) a1(x)= Y*cos(phi(x)) b1(x) = Y*sin(phi(x)) set angles degrees set xrange [0:90] # phi must be acute set yrange [0:2*Binput] set samples 10000 # enough for a good resolution, # low enough to make zooms in less than 1 second plot B(x), Binput, X(x), a1(x), b1(x) Let's now locate intersections of B(x) with the horizontal straight line y=Binput: between 0° and 90°, there can be 0, 1, 2 or 3 of them (this is shown below, in 3.b)). For each intersection, zooming on it, let's locate corresponding angle (i.e. its x coordinate). With input data shown in the example given, A = 1, Binput = 1, Y = .8, we find phi = 17.1144°, 45°, 72.8856°. Zooming on vertical straight lines passing throught each intersection, we can locate corresponding values for X(phi): positive ones can be accepted; negative ones must be dropped, together with corresponding value of phi. In the example given, they turn out to be X = 0.8, 0.61421, 0.8 that, from a geometrical point of view, is quite obvious: all squares with sides Y < 1 can be inscribed in a square 1 by 1, but also a rectangle with two sides of given lenght Y, forming two isosceles triangles at 45° with two opposite verteces of the square 1 by 1 (second solution). The third solution is symmetric of the first one, with respect to both axes of the 1 by 1 square. Having obtained three values for phi, it is easy to find corresponding values of a1 and b1, for an easy drawing of three inner rectangles. 3.b). Analytical solution ------------------------ Any manipulation of (1) and (2), aimed at eliminating all unkowns but one, will lead to a quartic (fourth degree) equation in the remaining unkown. This appears somehow discouraging, because, while most people knows how to solve a quadratic equation with the well-known formula x=[-b +- sqrt(b^2-4ac)]/(2a), not everybody suspect that similar, thought more complex, formulas exist to compute exactly complex radices of a cubic and of a quartic polynomials. They have been found, respectively, by Girolamo Cardano and by his pupil Lodovico Ferrari, and published together in "Ars Magna" in 1545. Here they are:-) http://planetmath.org/encyclopedia/QuarticFormula.html Ferrari's algorithm looks actractive, because it is NOT iterative: however, it is NOT suggested for solutions of quartic equations: it has been said "beautiful mathematics but a bad algorithm", expensive, unstable and inaccurate. Use of Jenkins-Traub_algorithm is suggested instead ( http://en.wikipedia.org/wiki/Jenkins-Traub_algorithm ); or, create a companion matrix of the polynomial ( http://en.wikipedia.org/wiki/Companion_matrix ) and compute its eigenvalues with QR method ( http://en.wikipedia.org/wiki/QR_algorithm ). However, if implemented on a common PC, in this case Ferrari's algorithm works, following some precautions: - use complex double precision variables; - accept as valid solutions "real" numbers with imaginary part nonzero, negative lenghts, cosines greater than 1, if these differences vanish in conversion to single precision variables. Geometrically, rectangular triangles, cut by inner rectangle out of outer one, are all similar: hence (3) b1 / a1 = (A-a1) / (B-b1) and, for Pythagoras' theorem, (4) a1^2 + b1^2 = Y^2 (5) (A-a1)^2 + (B-b1)^2 = X^2 i.e. (3') b1(B-b1) = a1(A-a1) Bb1-b1^2 = Aa1-a1^2 Bb1 = Aa1 + b1^2 - a1^2 elevando al quadrato (3") B^2b1^2 = A^2a1^2 + b1^4 + a1^4 + 2Aa1b1^2 - 2Aa1^3 -2 b1^2a1^2 (4') b1^2 = Y^2 - a1^2 replacing b1^2 from (4') into (3") B^2(Y^2-a1^2) = A^2a1^2 + (Y^2-a1^2)^2 + a1^4 + 2Aa1(Y^2-a1^2) - 2Aa1^3 - 2(Y^2-a1^2)a1^2 B^2 Y^2 - B^2 a1^2 = A^2 a1^2 - 4 A a1^3 + 2 A a1 Y^2 + 4 a1^4 - 4 a1^2 Y^2 + Y^4 and, after rearrangement, 4a1^4 - 4Aa1^3 + (A^2+B^2-4Y^2)a1^2 + 2AY^2a1 + (Y^4 - B^2Y^2) = 0 that is a quartic in unkown a1 in canonical form: ax^4 + bx^3 + cx^2 + dx + e = 0 where a = 4 b = - 4A c = A^2+B^2-4Y^2 d = 2AY^2 e = Y^4 - B^2Y^2 An algorithm to compute step by step its solutions using Ferrari's method is well exposed here: http://en.wikipedia.org/wiki/Quartic_equation#Solving_a_quartic_equation "subroutine Lodovico_Ferrari", that can be found in appendix, follows strictly this method. It solves correctly quartic equations built starting from their solutions, e.g. (x-1)(x-7)(x+2/3)(x+50)=0 or (x-3)(x-15)(X^2+1)=0. It is written in fortan 95, together with a subroutine that computes a,b,c,d,e from input data, and a test program. Them all can be compiled with f95 and run on a Linux machine. 3.c). Numerical solution ------------------------ Follows same guidelines of given graphical solution: X(phi) = (A-Y*cos(phi))/sin(phi) B(phi) = Y*sin(phi) + X(phi)*cos(phi) and the problem is reduced to finding zeroes of B(phi)-Binput in phi interval ]0°,90°] (0 is excluded because, for phi=0°, B(phi) turns out +infinite, -infinite, or indeterminate.) Zeroes of B(phi)-Binput can be found iteratively with several algorithms, the simplest one being bisection. However, them all require knowledge of contiguous intervals [phi_1,phi_2], where B(phi)-Binput is defined and continuous and changes sign only once; i.e. phi_1 and phi_2 "bracket" a zero, and only one zero. Such intervals must be identified before. Let us study behaviour of function B(phi) in interval [0°,90°], considering A and Y parameters. It's (geometrically) obvious that B(90°) = Y. Limit of B(phi) for phi->0 has a geometrical meaning only if Y<=A; otherwise, phi is inferiorly limited by the need of inserting into rectangle AxB at least a rectangle Yx0 (i.e. a segment of lenght Y): it is phi_min = arcos(A/Y), e B(phi_min) = Y sin(phi_min). case 3.c.1): Y >= A In this case, if Y>A, B(phi) has always a maximum between phi_min and 90°. If Y=A, B(0°) turns out indetermined: in facts, for any Binput, choosing X=Binput suffices to obtain a solution (outer rectangle inscribes itself). However, for phi->0, B(phi) tends to 0, and, in semi-open interval ]0°,90°], B(phi) and also X(phi) are always positive. If Y0(B(phi))=+infinite, and, for what concerns behaviour of B(phi), one must distinguish two subcases: case 3.c.2): Y =< A/sqrt(2) In this case, between 0° and 90° B(phi) is monotone and decreasing. In the particular case Y = A/sqrt(2), B(phi) has an horizontal inflection point for phi=45°, where B(phi)=A; anyway, between 0° e 90° B(phi) is decreasing. Hence, if Y =< A/sqrt(2), for any Binput>=Y one and only one solution exists; for Binput Y > A/sqrt(2) In this case, B(phi) has a minimum for phi between 0° and 45°, and a maximum for phi between 45° and 90°. For any Binput>=Y at least one solution exists (but there can be 3 of them); for Binput0; so minima and maxima can fall only where (a) y*cos(x)*(1+2*sin(x)**2) = 1. The value of cos(x)*(1+2*sin(x)**2) is 1 for x=0, 0 for x=pi/2 (phi=90°), and its maximum value is sqrt(2) where x=pi/4 (phi=45°). So, for y = sqrt(1/2), equation (a) can be verified straightforwardly for x=pi/4 (phi=45°); but B(x) has there an horizontal inflection point. For any other value of x in the interval, the left member of (a) is <1. A fortiori, if y < sqrt(1/2), the left member of (a) is always <1, and so B(x) is monotone. For y > sqrt(1/2), but < 1, the value of left member of (a) is y (i.e. < 1) for x=0, 0 for x=pi/2 (phi=90°), ad for x=pi/4 (phi=45°) it has a maximum y*sqrt(2) (that is > 1). So, it equals 1 for two values of x, one in the interval 01, the value of left member of (a) is > 1 for x=0, and (a) holds only in one point, in the interval pi/445°, it's enough to show geometrically that, for phi=phi_min, B(phi) is increasing: in facts, having inserted exactly a rectangle Yx0 in AxB as its diagonal, any rotation of one of its Y sides requires B to increase. ------------------------------------------------------------ Appendix ======== Subroutines and test programs (Fortran 95): ------------------------------------------- Follows on next post. This margin is too narrow to contain them :-)