Area of the triangle generated by the hands of a clock

Question

Problem explanation:

Let’s suppose we have a clock that has an hour hand 3 units long, a minute hand 4 units long, and a second hand 5 units long. The hour hand moves once every hour, the minute hand moves once every minute, and the second hand moves once every second. Therefore, exactly every second, the triangle defined by the ends of the hands changes its area.

Which is the maximum area between two given times?

Input: A sequence of hours in sets of 2 in the format hh:mm:ss.

Output: The maximum area defined by the ends of the hands in square units.

For a more comprehensive explanation: https://jutge.org/problems/P17681_en/statement

My code:

I have written a function that takes as arguments the position of the hands in radians and returns its area (trigonometrically intensive).

double area(double secRad, double minRad, double houRad){
    double alfaSM = atan( (5-4*cos(abs(secRad-minRad)))/(4*sin(abs(secRad-minRad))) );
    double alfaSH = atan( (5-3*cos(abs(secRad-houRad)))/(3*sin(abs(secRad-houRad))) );
    double alfaMH = atan( (4-3*cos(abs(minRad-houRad)))/(3*sin(abs(minRad-houRad))) );

    double distSM = abs(5*sin(alfaSM)+4*sin(abs(secRad-minRad)-alfaSM));
    double distSH = abs(5*sin(alfaSH)+3*sin(abs(secRad-houRad)-alfaSH));
    double distMH = abs(4*sin(alfaMH)+3*sin(abs(minRad-houRad)-alfaMH));

    double s = (distSM+distSH+distMH)/2;

    return s*(s-distSM)*(s-distSH)*(s-distMH);
}

My main function reads the input from the user and calculates the area for each second in the range of time the user has entered, stores the maximum area and finally prints it.

int main(){
    char z;
    double hI,mI,sI, hF,mF,sF;
    while(cin >> hI >> z >> mI >> z >> sI >> hF >> z >> mF >> z >> sF){
        double maxArea = 0;
        for(int i = 3600*hI+60*mI+sI; i <= 3600*hF+60*mF+sF; i++){
            double cacheArea = area( 2*pi*((double)(i%60)/60) , 2*pi*(double)((i%3600)/60)/60, 2*pi*(double)((i%43200)/3600)/12);
            if(cacheArea > maxArea) maxArea = cacheArea; 
        }
        cout << fixed << setprecision(3) << sqrt(maxArea) << endl;
    }
}

Write this at the beginning so you can test the program:

#include <iostream>
#include <cmath>
#include <iomanip>
#define pi 3.14159265358979323846
using namespace std;

Where am I now and what is left?

The algorithm of the code above works as intended, but it isn't efficient enough. The online judge throws a time limit exceeded error.

I have managed to reduce the number of calculations of area needed by supposing that if the range of time is more than an hour, we just need to check all the combinations for the first hour, since the maximum area will be there and will be appearing again every hour(the relative positions repeat).

I have also implemented a variable that checks if the separation between hands of the clock in the time its checking is smaller than the separation of the hands in the maximum area found until that moment. If it is, the program skips the area calculation.

These changes seem to improve performance considerably (from 43,200 area calculations to ~1,500 in the complete spin of 12 hours), but the judge still doesn't admit it.

I seem to need a cleaner solution, a mathematically prettier one, but I'm struggling to find one.

Answer 1

The first thing to do is to get rid of trig functions which are slow (15-30 cycles) and replace them with precomputed tables. Since the hands have only 60 possible positions each this is a 60-slot table (plus a bit extra because who bothers with bounds checking, honestly?.... cough)

If you feel like cheating, you can also use a precomputed table in the source, but that would be a bit wimpy IMO.

Rotational Symmetry

Let's create a function A(h,m,s) which computes the area depending on hour, minute and second.

If we look at the clock face and rotate the entire thing (three clock hands included) by any angle, then the area will not change. Thus we can put the hours hand on "0" (the other hands also follow).

Therefore,

\$ A(h,m,s) = A(0, m-5h \mod 60, s-5h \mod 60)\$

Thus we can make a neat little precomputed table area[minute][second], filled with 3600 values, and use that! Populating this takes negligible time relative to the programs' libc initialization time, so I did not bother to optimize it further. It is just a bunch of MACs which modern cpus crunch at insane speed anyway.

Additionally, a per-minute max area is computed and stored in array marea. This corresponds to the max over h=0, m, s=[0..59] and it can be extended to any h:m by rotating as explained above.

Counting

The rest is very simple, We start, say at 00:01:02 and end at 00:35:15.

First we increment the seconds from 02 to 59 until we get to 00:02:00

Then we increment the minutes, until we reach 00:35:00

Then we increment the seconds again until we reach the final time.

At each step we pick the area from the precomputed table (or the per-minute precomputed max) and compare it to our running maximum.

It's pretty fast since it does basically nothing except a few int ops, a float cmp, and there should be a few CMOVs if the compiler does its job. No trig at all...

This ugly piece of code should thus do the deed...

#include <iostream>
#include <stdio.h>
#include <cmath>
#include <iomanip>
#define pi 3.14159265358979323846
using namespace std;

float sin_lut[76];
inline float sint( int x ) { return sin_lut[x]; }
inline float cost( int x ) { return sin_lut[15+x]; }

float sarea[60][60]; // area for 00:m:s
float marea[60];    // max area for 00:m and all seconds
float globalmax = 0.0;

// this is to avoid typing tons of if()
struct Maxi 
{
    float _max;
    Maxi( float x ) { _max=x; }

    // adds the area of h:m:s to the running max
    void hms( int h, int m, int s ) 
    { 
        m -= h*5;   // align hour hands on zero by rotating the dial
        s -= h*5;
        m = (m<0) ? (m+60) : m; // bounds check
        s = (s<0) ? (s+60) : s;
        float x = sarea[m][s];
        if(x>_max) _max=x; 
    };

    // adds the area of h:m to the running max (with all possible seconds)
    void hm( int h, int m )
    {
        m -= h*5;   // align hour hands on zero by rotating the dial
        m = (m<0) ? (m+60) : m; // bounds check
        float x = marea[m];
        if(x>_max) _max=x; 
    };
};


float get_max_area( int hI, int mI, int sI, int hF, int mF, int sF )
{
    Maxi maxa( sarea[mI][sI] ); // this is our running max area

    //~ printf( "\n\nSearch %02d:%02d:%02d - %02d:%02d:%02d\n\n", hI,mI,sI, hF,mF,sF );

    // if it makes one complete turn it will explore all positions...
    // thus return the global maximum.
    int tI = hI*3600 + mI*60 + sI;
    int tF = hF*3600 + mF*60 + sF;
    if( tF-tI >= 3600 )
        return globalmax;

    // increase time in clever steps...
    if( hI < hF || mI < mF )
    {
        // go to the end of the current minute.
        while( sI < 60 )
            maxa.hms( hI, mI, sI++ );
        sI = 0;
        mI++;
        if( mI>=60 )
        {
            mI=0;
            hI++;
        }
    }

    // Increase minutes...
    while( hI < hF || mI < mF )
    {
        maxa.hm( hI, mI++ );
        if( mI>=60 )
        {
            mI=0;
            hI++;
        }
    }

    // finish the last minute...
    while( sI <= sF )
    {
        maxa.hms( hI, mI, sI++ );
    }

    return maxa._max;
}


int main()
{
    /* Precalculate sine table because sin is slow ;) */
    for( int i=0; i<=15; i++ )
    {
        float s = sin( i* (pi/30.) );
        sin_lut[i]    = sin_lut[30-i] = sin_lut[i+60] =s;
        sin_lut[i+30] = sin_lut[60-i] = -s;
    }

    // for( int i=0; i<60; i++ ) printf( "%d %1.08f %1.08f\n", i, sint(i), cost(i) );

    /* Precompute area for all minutes and second positions
       for zero hour. We can use this to determine the area for
       any h:m:s by rotating the whole clock face to put the hour
       hand on zero, so this covers all cases. 

       Note: symmetry could be used here, also the constants could be optimized
       but this whole thing takes less than 1ms, much less time than libc startup/init 
       so... who cares? 
    */    
    for( int m=0; m<60; m++ )
    {
        float x2 = 4.0*sint(m);
        float y2 = 4.0*cost(m);
        float maxarea = 0.0;

        for( int s=0; s<60; s++ )
        {
            float x3 = 5.0*sint(s);
            float y3 = 5.0*cost(s);

            // cross product area would be: 0.5* (-x2*y1 + x3*y1 + x1*y2 - x3*y2 - x1*y3 + x2*y3
            // however, conveniently, hour is zero, so x1=0 and y1=3...
            float a = abs( 0.5* (-x2*3.0 + x3*3.0 - x3*y2 + x2*y3));
            sarea[m][s] = a;
            if( a>maxarea )
                maxarea = a;

            //~ printf( "00:%02d:%02d %2.08f\n", m,s,area[m][s] );
        }
        marea[m] = maxarea;
        //~ printf( "00:%02d:xx %2.08f\n", m,maxarea );
    }

    // and the max during one hour, which is also the global max...
    for( int i=0; i<60; i++ )
        if( marea[i]>globalmax )
            globalmax = marea[i];

    // process input
    char z;
    int  hI,mI,sI, hF,mF,sF;
    while(cin >> hI >> z >> mI >> z >> sI >> hF >> z >> mF >> z >> sF)
    {
        float m = get_max_area( hI, mI, sI, hF, mF, sF );

        //~ printf( "%02d:%02d:%02d - %02d:%02d:%02d - %.03f\n", hI,mI,sI, hF,mF,sF, m );
        printf( "%.03f\n", m );

    }

}

Oops I forgot the review.

The point here was to kinda reverse-engineer the judging process. For just a few input lines, even your implementation which uses lots of slow trig should be fast enough on a modern cpu to take negligible time compared to program init.

Therefore, if you get a time limit error, this means the input file must contain tons of lines.

Therefore, precalculated look-up tables are a good fit.

And using a non-constant time step (first second, then minutes) reduces the number of iterations substantially, but it does require a precalculated maximum area per each minute table. You could include it as text into the source code, but that would be lame... Thus it is best to precalculate it all and then rip through the input file.

I would not be surprised if cin>> was the slowest part of the program...

Answer 2

This is an addendum to Oscar's answer. Just doing some math.

Let t = 0 be 12 o'clock, when all hands are coinciding and the area is 0. We need to calculate the area as a function of t. Let s, m, h be the angles of second, minute and hour with respect to the trigonometric reference. Then $$s=\frac{\pi}2-\frac{\pi}{30}t\\ m=\frac{\pi}2-\frac{\pi}{1800}t\\ h=\frac{\pi}2-\frac{\pi}{21600}t$$ and let S, M, H be the lengths of clock's hands. Then the area is $$A=\frac12\left|SM\sin\alpha+MH\sin\beta+HS\sin\gamma\right|$$ where $$\alpha=s-m,\;\beta=m-h,\;\gamma=h-s$$ Therefore we have $$\alpha=-\frac{59}{60}\cdot\frac{\pi}{30}t\\ \beta=-\frac{11}{720}\cdot\frac{\pi}{30}t\\ \gamma=\frac{719}{720}\cdot\frac{\pi}{30}t$$ where t is a positive integer. As you see, doing some boring math can sometimes reduce the amount of codes needed for a job.

P.S. and here is a plot of A vs. t. It is periodic indeed.

Answer 3

So the first place to look is your area function. This currently uses 15 trig operations to calculate the area of a triangle (3 cos, 3 atan, and 9 sin). It is possible to instead do the same using only 3 sin functions and some other math that will take far less time. The basic idea of how to do this is that if you have two side lengths of a triangle, and the angle between them, the area of the tringle is the product of the side lengths divided by two times the sin of the angle between them, or sin(θ)*a*b/2. Thus, all you have to do is find the area of each piece of the three pieces of the triangle, and add them together.

Answer 4

Avoid using namespace std;

Importing all names of a namespace is a bad habit to get into, and can cause surprise when names like begin and size are in the global namespace. Get used to using the namespace prefix (std is intentionally very short), or importing just the names you need into the smallest reasonable scope.

The exceptions to this rule are namespaces explicitly intended to be imported wholesale, such as the std::literals namespaces.

Prefer named constants over preprocessor macros

The value of π can be declared as a double - and it can be calculated, rather than written:

static const double pi = 4 * std::atan(1);

Area calculation

The area() function is misnamed - it returns the square of the area.

The similar repeated expressions are very prone to simple typos that are easily missed.

Whenever I see std::atan() with a division, it's a sign that std::atan2() may be what's really meant.

Much of the trigonometry can be eliminated by using the cross product of any two of the sides as vectors (remember to divide by two to convert the parallelogram area to triangle area!).

Main program

It's inefficient to compute the triangle area for every single second in the interval. Observe the symmetry of the problem: if the time span is longer than 48 minutes, we can simply return a precomputed constant. For time spans of a minute or more, we can find the largest area by placing the second hand bisecting the large angle between minute and hour hands, and moving it toward the hour until the area starts decreasing. Only for the partial minute at each end of the span do we need to fall back to calculating all the areas.

With the changes from the above paragraph, there's no longer a need to compose hours, minutes and seconds into a single number and decompose them again. That's good, because it's complex and hard to understand; it also has a bug, because 12 hours is 43200 seconds, which is greater than the required range of (signed) int.

The program would be improved by separating the I/O from the calculation of maximum area for a set of start and end values. That would enable some unit testing, rather than relying on whole-program tests with input and output streams.

Answer 5

The area of the triangle given points A, B, C is half the magnitude of the cross product of vectors AB and AC. The points are defined in 3D, with a zero z component. The cross product will only be non-zero along the z-axis (it's perpendicular to the plane formed by the two vectors). Thus, you only need to compute the z component of the cross product. This z-component is the area of the parallelepiped spanned by the two vectors. This is a few multiplications, no trigonometry involved.

You do need trigonometry to compute the position of the end points of the hands. You do this only when hands change. Or you could even pre-compute sin and cos for all 60 possible angles. And note that you only have unique values over 8 of those angles (the first 45 degrees), the rest are trivially derived from those.

Answer 6

Since the problem is symmetric, over a period of one hour all possible areas are 'reached'. One optimization I can think of is to simply return the [pre-computed] maximal possible area whenever the input time range is longer than one hour.

In fact, the period required to ensure all areas are reached is not one hour, but merely 2844 seconds! This can be demonstrated using the following code:

double maxmax = 0;    
for(int i = 0; i <= 3600; i++)
{
    double cacheArea = area( 2*pi*((double)(i%60)/60) , 2*pi*(double)((i%3600)/60)/60, 2*pi*(double)((i%43200)/3600)/12);
    if (cacheArea > maxmax) 
        maxmax = cacheArea; 
}

for (int s = 0; s <= 23*60*60; s++)
{
    double maxArea = 0;
    for(int i = s; i <= s+2844; i++)
    {
        double cacheArea = area( 2*pi*((double)(i%60)/60) , 2*pi*(double)((i%3600)/60)/60, 2*pi*(double)((i%43200)/3600)/12);
        if (cacheArea > maxArea) 
            maxArea = cacheArea; 
    }

    if (fabs(maxArea - maxmax) > 1e-10)
        cout << fixed << setprecision(3) << sqrt(maxArea) << "   " << sqrt(maxmax) << endl;
}

Answer 7

Lior Kogan showed manually that any timespan strictly longer than 2844 seconds reaches a global maximum.

if timespan >= 2844 seconds:
    return maximum possible area

If the time span is not this long, we can break what is left into some number of whole minutes plus up to two partial minutes. The maximum over a whole minute is determined by only the hour and the minute hand; the second hand makes a complete turn so we do not depend on its value. Since the problem is rotationally invariant and symmetric, this only truly depends on the absolute difference between the hour and the minute, of which there are only 31 possibilities. These can be cached in a length-31 array.

for hour:minute in full_minutes(timespan):
    maximum = max(maximum, full_minute_cache[abs(minute - 5 * hour)])

Then you have at most two partial minutes left over, which means there are at most a further 118 specific timestamps. You can do a full check on these fairly quickly, but it is easy enough to skip it if the minute cache says the full minute never beats the current maximum.

first_hour:first_minute = timespan.start
if full_minute_cache[abs(first_minute - 5 * first_hour)] > maximum:
    for timestamp in first_minute(timespan):
        maximum = max(maximum, calculate_area(timestamp))

last_hour:last_minute = timespan.end
if full_minute_cache[abs(last_minute - 5 * last_hour)] > maximum:
    for timestamp in last_minute(timespan):
        maximum = max(maximum, calculate_area(timestamp))

This is a worst-case total of 118 calls to calculate_area and 49 array lookups, which is extremely fast even with an inefficient method to calculate the area.

---

We can do better if we cache more; for each of the 60 hour-minute differences, there is a given set of results for the areas of the triangle. We don't want to just store all of the many results, but we can cheat. Because the function is mostly smooth, this will have a small number of maxima; it turns out no more than three maxima total, no more than two internal maxima, and no less than one internal maxima, and when there are three at least one is at 0 or 59 seconds.

If we look at some subsection of this, \$x..y\$, we know either the global maxima of it is going to be one of the internal maxima strictly within that range, or the value of \$x\$ or \$y\$. There can only be two strictly within, so this leaves us with four calls to calculate_area and a cache of length 120.

Our total cost is now down to no more than eight calls to calculate_area, four for each of the two partial minutes, 47 array lookups for the minutes, a length-31 array of doubles and a length-120 array of bytes. The 47 lookups are probably no slower than eight calls to calculate_area, so this is about where optimisations stop being algorithmic, and where I shall leave this answer to rest.

Answer 8

There are few aspects to optimising this problem.

The first is recognising that there is a lot of redundancy in the problem space. With rotation, every hour of the day can actually be mapped to any other hour with an appropriate rotation. Equally, the first half of an hour is a reflection of the second. This means we can precompute a small segment of non-repetitive calculations and then map any other times on to that.

The second is that those pre-computations do not need to factor into your runtime; it makes a lot more sense to load them in and just do a lookup. It's not cheating if it's both fast and correct. As such, I'd wave off trying to optimize how you generate your values as unnecessary and just commit the results to memory.

If you look at the problem space carefully, we can go even further than that. polfosol's graph helps a lot here: we can see there are clear, regular local maxima roughly once per minute plus a clear maxima once per hour. If we already know those maxima and if any given time period contains one or more maxima the maximum is guaranteed to the the largest of those values. If we manage to select a time period that does not contain any maxima, the maximum value is guaranteed to be either the first or the last value (which, for 2 calculations, is cheap to calculate on the fly to be worth optimising away).

This makes the solution incredibly simple: for segments larger than one minute or so we're guaranteed to find a maxima with a simple search within boundaries; searching <30 maxima for the biggest is not computationally intensive. For segments less than a minute we still have to check the local maxima and, if none are contained, do two calculations for the bookend values: one of them is guaranteed to be the maximum in a no-maxima containing segment as the segment must be either a 'u' shape or a slope.

Sometimes optimisation is finding ways to do the minimum number of calculations rather than doing calculations themselves faster.