NTHCIR - Editorial

PROBLEM LINK:

Practice
Contest

Author: Adury Surya Kiran
Tester: Mugurel Ionut Andreica
Editorialist: Lalit Kundu

DIFFICULTY:

Medium-Hard

PREREQUISITES:

advanced maths, geometry, recurrences

PROBLEM:

We are given radius of circles C_1, C_2, C_3 and C_4 in the following image:

Now, we are going to draw circles C_5, C_6 and so on in such a way that C_n touches circles C_1, C_2 and C_{n-1}. Like as shown in the following image:

Given upto 10^7 values of n, find sum of radius of circle C_n for all n.

QUICK EXPLANATION:

========================
Using Descartesā€™ Theorem, if we have 3 existing circles of curvature(inverse of radius) a_1, a_2, a_3, we can draw a 4^{th} circle of curvature a_4, touching all three circles such that a_4 = a_1 + a_2 + a_3 \pm 2 \sqrt{a_1*a_2 + a_2*a_3 + a_1*a_3}, where a_4 is negative if it represents a circle that circumscribes the first three.

Now, n^{th} circle in our problem is formed by finding the 4^{th} circle touching circles C_1, C_2 and C_{n-1}. So, we can define our recurrence, solving which we get general term for the radius of n^{th} circle.

EXPLANATION:

================

One of the ways to solve this problem is via Descartesā€™ Theorem. First, lets define what curvature of a circle is.

The curvature of a circle with radius r is defined as \pm \frac{1}{r}. The plus sign here applies to a circle that is externally tangent to the other circles, like the red circle is externally tangent to the three black circles in the image. For an internally tangent circle like the big red circle, that circumscribes the other circles, the minus sign applies.

Descartesā€™ Theorem states that if four circles are tangent to each other at six distinct points, and the circles have curvatures a_i (for i = 1, ..., 4), then:

(a_1+a_2+a_3+a_4)^2=2(a_1^2+a_2^2+a_3^2+a_4^2).

When trying to find the radius of a fourth circle tangent to three given kissing circles(circles that touch each other mutually), the equation is best rewritten as:

a_4 = a_1 + a_2 + a_3 \pm 2 \sqrt{a_1*a_2 + a_2*a_3 + a_1*a_3}

The \pm sign reflects the fact that there are in general two solutions. One solution is positive and the other is either positive or negative; if negative, it represents a circle that circumscribes the first three (as shown in the diagram above).

Now, letā€™s see how we employ this theorem to find solution to the problem.

We know that C_n(for n \ge 4) is formed as a 4^{th} circle tangent to three kissing circles C_1, C_2 and C_{n-1}, where C_1 is being touched by all circles internally. And also, the circle C_n will touch all other circles externally, so we are interested in positive curvature of C_n. And, curvature of C_1 is -\frac{1}{r_1}.

So, we can write using Descartesā€™ Theorem that a_n = a_1 + a_2 + a_{n-1} + 2 \sqrt{a_1*a_2 + (a_2 + a_1)*a_{n-1}}. Now, we can employ this recurrence to find n^{th} circle, in O(n). This will be enough to pass subtask 1.

For subtask 2, we need a constant time solution for each n, that means we need to arrive at a closed form of the n^{th} term of recurrence. Letā€™s see how we can solve this recurrence. Basically, our recurrence is:

a_n = c_1 + a_{n-1} + 2 \sqrt{c_2 + c_1 a_{n-1}},
where c_1 = a_1 + a_2 and
c_2 = a_1 * a_2

Here, we substitute x_n(a new recurrence) as \sqrt{c_2 + c_1 a_{n}}, which implies that
a_n = \frac{x_n^2 - c_2}{c_1} --------Eqn (1)

So, our recurrence is now a_n = c_1 + a_{n-1} + 2 x_{n-1}. In this equation, substitute a_{n-1} and a_n using equation 1 to get,

x_n^2 = x_{n-1}^2 + 2c_1x_{n-1} + c_1^2, which is basically x_n^2 = (x_{n-1}+ c_1)^2. Now, this has greatly been simplified to x_n = x_{n-1} + c_1, whose general term can be written as x_n = nc_1 + C, where C is an arbitrary constant, which depends on x_1.

Now, to find a_n we substitute general term for x_n in equation 1, which yields a_n = \frac{(nc_1 + C)^2 - c_2}{c_1}.

Now, we know the value of a_3 as \frac{1}{r_3}, so we can find value of C by substituting n as 3 in the general form of a_n. As you might note, it will give us two values for C because its a quadratic in C, so, we take the value which gives us correct value for a_4 whose value we already know as \frac{1}{r_4}.

Now, our solution for a single n is constant time.

ALTERNATE SOLUTION:

========================

There is another solution using the concept of inversion of a plane. Letā€™s see this concept in detail.
We are going to transform all points in an x-y plane(say plane 1) to another x-y plane(say plane 2). For each point (x, y) in plane 1, we map it to \frac{x}{r^2}, \frac{y}{r^2}, where r is distance of point (x, y) from origin i.e. \sqrt{x^2 + y^2}. If we consider point (x, y) in polar coordinates as r e^{i \theta}, then we map it to \frac{1}{r} e^{i \theta} in the plane 2. Transformation from plane 2 to plane 1 is simple, we just have to apply the same transformation again for all points in plane 2.

Now, letā€™s see how it affects simple figures: circle and lines. Here I am going to summarise the results, you can easily derive these results by considering equations of lines and circles(doing in polar coordinates will be easier).

  • All the circles which pass through origin in plane 1 will be mapped to straight lines in plane 2. A circle with center at (r, 0) and radius r will be mapped to line x = \frac{1}{2 r} in plane 2.

  • All the straight lines which do not pass through origin in the first plane will be mapped to circles which pass through origin in the second.

  • All the circles which do not pass through origin in the first plane will be mapped to some other circles in the second plane.

  • All the straight lines which pass through origin in the first plane will be mapped to the same lines in the second. Also, note that if there are two points on such a line where one point is at a distance r_1 from origin and another at distance r_2, where r_2 \gt r_1, then distance between them in plane 1 is r_2 - r_1 and distance between these points in transformed plane will be \frac{1}{r_1} - \frac{1}{r_2}.

Now, letā€™s see how we employe this transformation to solve our problem. First, weā€™ll choose a suitable origin in plane 1(plane where all our circles are present). Now, we see that all circles are tangent to circles C_1 and C_2 and we remember that we can transform circles passing through origin to straight lines, so, we choose the point where C_1 and C_2 touch as the origin. Now, we also assume that their centers lie on y-axis at points R_1 and R_2 respectively.

Now if we transform the plane, then circles C_1 and C_2 will become straight lines which are parallel to each other. All the circles C_3, C_4, C_5 \cdots will become some other circles in the transformed plane. But the fact that they all are touching C_1 and C_2 before transformation means that they all touch the two straight lines which are parallel to each other after transformation, which makes all their radii equal in the transformed plane(this radius is equal to \frac{1}{4R_1} - \frac{1}{4R_2}). See the following figure.

So, the equations of lines C_1^{\prime} and C_2^{\prime} are y = \frac{1}{2 R_1} and y = \frac{1}{2 R_2} respectively.

Now, to find the radius of any n^{th} circle in plane 1, we first find coordinates(in plane 2) of two points on the circle and we transform these points back to plane 2 and get the distance between them as the diameter.

We already know the y-coordinate of the centre of circles C_3^{\prime}, C_4^{\prime} and so on. We need to find their x coordinate now. If we have x-coordinate of the center of circle C_3^{\prime}, then we can easily find centre of any circle and we already know the radius of circles in plane 2, so we can also find the coordinates of a point on the circumference, which weā€™ll transform back.

Letā€™s see how we find x-coordinate of center of C_3^{\prime} in plane 2. Keep referring to image below while reading this for more clarity. We know that a line passing through origin in plane 1 remains same in plane 2, so we draw a line Z_1 in plane 1 passing from origin and center of circle C_3. Now, this line in plane represented as Z_1^{\prime} will also pass through center of circle C_3^{\prime} and origin of plane 2. If we know distance of center of C_3^{\prime} from origin somehow, we can calculate its x-coordinate also, since we already have y-coordinate.

We now, try to relate distances in plane 2 with distances in plane 1. Say, the distance d = distance between origin and center of C_3^{\prime}. So, circle C_3^{\prime} cuts line Z_1^{\prime} at two points, distance between these points in plane 1 must be the diameter of C_3 i.e. 2R_3 because the same line intersects circle C_3 at same two points. Distance of these two points from origin in plane 2 can be written as d - r_3^{\prime} and d + r_3^{\prime}. Now, if we transform them to plane 1 and see distance between them, we can say that:

2 r_3 = \frac{1}{d - r_3^{\prime}} - \frac{1}{d + r_3^{\prime}}

In above equation, we can find the d as all other variables are known. And using value of d we can find x-coordinate of center of x_3^{\prime} = C_3^{\prime} as \pm \sqrt{d^2 - y^2}. So, center of n^{th} circle will be x_3^{\prime} + 2(n-3)R_3 or x_3^{\prime} - 2(n-3)R_3, we can now check which one to use by finding radius of C_4 and match with given value in input.

Now, our problem is solved. You can do the calculations yourself to find the closed formula.

COMPLEXITY:

Constant per query.

AUTHORā€™S, TESTERā€™S SOLUTIONS:

setter
tester

7 Likes

There is one more recurrence relation which can be formedā€¦

as we know ā€¦C(N)=C(1)+C(2)+C(N-1) +/- 2*SQRT(REMAINING PART :D)

we can also see the complement circle for n th circle will be (n-2)th circleā€¦

so if C(N)=C(1)+C(2)+C(N-1) + 2*SQRT(REMAINING PART :D)

THEN C(N-2)=C(1)+C(2)+C(N-1) - 2*SQRT(REMAINING PART :D)

ADDING BOTH EQUATIONS WE GET : C(N)+C(N-2)=2C(1)+2C(2)+2C(N-1)ā€¦

4 Likes

Why dont we square on both sides or in other words use (k_n + k_n-1 + k_1 +k_2)^2 = 2*((k_n)^2 + (k_n-1)^2 + (k_1)^2 + (k_2)^2)ā€¦and then solve recursively.
By the way this question was like the ā€œMother Of Simplificationā€. I couldnt believe my own eyes when it boiled down to such a small formula which is kn=k4 + (n-4)((k4-k3) + (n-3)(k1+k2)) where n>4.

Instead of solving the recurrence, what I noticed was that k[i] (curvature of i-th circle) had a sequence such that the difference was in arithmetic progression. In this way too, we can calculate the curvature of n-th circle in O(1).

Just curious - was it expected to be solved by ourselves, without using Google? :slight_smile: I canā€™t understand the logic behind giving tasks of ā€œenter statement to Google and copy-paste formula from thereā€ kind. It isnā€™t even useful for ACM ICPC skills - I have strong doubt that solving such task does not affect problem solving skills at all.

16 Likes

i got wa in the 2nd task of 2nd subtask ā€¦probably due to precision ā€¦i used the formula kn=(n-3)^2 * (k1+k2) + k3 + 2*(n-3)sqrt(k1k2+k2k3+k1k3) in order ā€¦what precision error could possibly pop out ?

Even in my JEE I never came across even such terms ā€œApollonian Gasketsā€,ā€œSoddy Circlesā€,ā€œDescartes theoremā€,ā€œInversionsā€!! :confused:

2 Likes

What a beautiful questionā€¦ got the same recurrence but was not able to prove itā€¦

Also, Codechef please give the 2nd test case (1st subtask) values and solution as many got WA in this (including me)ā€¦

Later found matrix while surfing net, which can also help to solve for nth radius using Matrix exponentiation.

What do you mean by complement circle? Also how to do you prove that complement of nth circle is (n-2)th circle

if u see by descartes theoremā€¦for example initially u are given 4 circlesā€¦
1 outer C1 and three inner circles C2 C3 C4ā€¦now to make a new circle C5ā€¦u will put curvatures of C1,C2,and C4 in descartes equation which will provideā€¦two solutionsā€¦as it is quadratic in natureā€¦therefore one solution is the new circleā€¦what will be the other solutionā€¦the one which was already touching C2 C1 and C4ā€¦and what was that? C3 rytā€¦therefore ā€¦C(n) and C(n-2) are complement circlesā€¦u can see this by drawing :smiley: ā€¦though the recurrence relation u guys made was easier to solve than this one.

Yes, this is the recurrence I used in the testerā€™s solution. Then I simply found the general term of the linear recurrence.

1 Like

I have a doubt. The recursion is completely fine and theoritically we will get the right answer using it. But as test cases can be as large as 10^7 and the radius can be a real number as large as 10^12, the final answer can contain more than 20 digits (including decimal digits). So how are you guys getting the accuracy???..

Consider the test case:
10000000
1 2 1 1
99999999999.999999 r2 r3 r4

(since ā€˜mā€™ is 1 and B is 1 ā€˜nā€™ will always be 1 and so values of r2, r3 and r4 doesnā€™t affect the answer and so we can give any values for now)

I ran the setters solution on this test case and it gives wrong answer. The correct answer should be 999999999999999990 because we have to add r1 10^7 times.

Please clarify my doubt.

1 Like

I was wondering why this code didnā€™t work for the first subtask? It uses Descartes formula, yet it still gave WA. Can anyone tell me where is the mistake?

LL n, p, m, b;
double r1, r2, r3, r4, t1, t2, res;
vector R;

int main() {
//freopen(ā€œtest.inā€,ā€œrā€,stdin);
//freopen(ā€œtest.outā€,ā€œwā€,stdout);

int T;

scanf("%d", &T);

scanf("%lld%lld%lld%lld", &n, &p, &m, &b);
scanf("%lf%lf%lf%lf", &r1, &r2, &r3, &r4);
R.PB(r1); R.PB(r2); R.PB(r3); R.PB(r4);

r1 = 1.0 / r1; r2 = 1.0 / r2;
rb(i, 4, 2000) {
    r3 = r4; r3 = 1.0 / r3;
    t1 = r2 + r3 - r1;
    t2 = (r2 * r3) - r1 * (r2 + r3);
    r4 = t1 + 2.0 * sqrt(t2);
    r4 = 1.0 / r4;
    R.PB(r4);
}

res = 0.0;
while(T--) {
    //scanf("%d%d", &n, &q);
    n = (n * p) % m + b;
    res += R[n - 1];

    //printf("%lld\n", res + K[n][k]);
}

printf("%.9lf\n", res);


return 0;

}

I also have some doubts about precision :slight_smile: I was lazy to check it by myself. I had different solutions passing different tests - and I am not really sure that one which passed all tests is most precise :smiley:

I think the point was to see if we could figure out Subtask 2 using our algebra skills AND our knowledge of Descartes formula from Google.

2 Likes

You have to be very careful with your floating point calculations. I donā€™t think there was any test case in the original problem as brutal with precision as that (although I didnā€™t pass Subtask 2, so I donā€™t really know), but the second test case of Subtask 1 was hard with precision. I had to try multiple times to get my floating point calculations as precise as possible before I passed that test case.

Second test case of first subtask was probably wrong due to floating point errors. I also got WA on this test case for a long time until I finally fixed my floating point calculations and made them as precise as possible.

It was funny when one of my solutions was passing ONLY second case of subtask 2 but failing other 3 :smiley:

What does rb(i, 4, 2000) do? If it loops i from 4 to 2000, that was probably the problem because for Subtask 1, I think N_i could be as big as 3000.

please tell me why this did not pass even first sub task?

#include <iostream>
#include <cmath>
#include <cstdio>
 
typedef long long int ll;
typedef long double ldf;
 
using namespace std;
 
ll mulmod(ll a, ll b, ll c)
{
	ll x = 0;
	ll y = a%c;
 
	while(b>0)
	{
		if(b%2 == 1)
			x = (x+y)%c;
		y = (y*2)%c;
		b/=2;
	}
	return x%c;
}
 
int main()
{
	ll t;
	scanf("%lld",&t);
 
	ll n,p,m,b;
	scanf("%lld%lld%lld%lld",&n,&p,&m,&b);
 
	ll r1,r2,r3,r4;
	scanf("%lld%lld%lld%lld",&r1,&r2,&r3,&r4);
 
	ldf k1,k2,k3,k4;
 
	k1 = 1.0/r1;
	k2 = 1.0/r2;
	k3 = 1.0/r3;
	k4 = 1.0/r4;
 
	ldf ans = 0;
 
	for(ll x=1;x<=t;x++)
	{
		n = mulmod(p,n,m) + b;
		
		if(n==1)
		{
			ans += 1.0/k1;
			continue;
		}
		if(n==2)
		{
			ans += 1.0/k2;
			continue;
		}
		if(n==3)
		{
			ans += 1.0/k3;
			continue;
		}
		if(n==4)
		{
			ans += 1.0/k4;
			continue;
		}
		
		ldf z = k4;
		for(ll i=5;i<=n;i++)
			z = k2 - k1 + z + 2.0*sqrt(k2*z - k1*(k2 + z));
		ans += z;
	}
 
	ll ans1 = ans*10000000;
	ans = ans1/10000000.0;
	printf("%.6Lf",ans);
	return 0;
}