JMK no matter what

ACM ICPC 2008 Aizu

Monday, the local practice session at Yonsei for the ACM competition used the problemset from Asia Regional Aizu 2008. I ended up solving six before leaving early; and I solved the rest of the set yesterday.

Summary: Tons of (possibly) interesting geometry problems.. but four geometry out of ten? And there are problems which require crazy long implementations; No wonder nobody solved more than eight. :-p I could have solved eight if I could complete competing, though. Beware: long post

A: Grey Area

A simple implementation problem.

lang:cpp
#include<iostream>
#include<algorithm>
#include<vector>
#include<cstdio>
using namespace std;

#define all(x) (x).begin(),(x).end()
#define CLEAR(x,with) memset(x,with,sizeof(x))

int main()
{
    int n, w;
    while(scanf("%d %d", &n, &w) == 2)
    {
        if(n == 0 && w == 0) break;
        vector<int> cnt(101);
        for(int i = 0; i < n; ++i)
        {
            int v;
            scanf("%d", &v);
            cnt[v/w]++;
        }
        while(cnt.size() > 1 && cnt.back() == 0) cnt.pop_back();
        double ret = 0.01;
        int mv = *max_element(all(cnt));
        for(int i = 0; i < cnt.size(); ++i)
        {
            double shade = (cnt.size() - 1 - i) * 1.0 / (cnt.size() - 1);
            double height = cnt[i] * 1.0 / mv;
            ret += shade * height;
        }
        printf("%.8lf\n", ret);
    }
}

B: Expected Allowance

Given m n-sided dices m\times n \le 100,000), calculate E(max(1, sum_of_faces - W))?

A standard dynamic programming problem - however, since the state space has a size of m\times n \times m, using a quadratic array is just not feasible. We maintain a single row of the whole matrix to reduce the space compexity to O(mn).

lang:cpp
#include<iostream>
#include<cstring>
#include<algorithm>
#include<vector>
#include<cstdio>
using namespace std;

#define all(x) (x).begin(),(x).end()
#define CLEAR(x,with) memset(x,with,sizeof(x))

int main()
{
    int n, m, k;
    while(scanf("%d %d %d", &n, &m, &k) == 3)
    {
        if(n == 0) break;
        vector<double> prob(n*m+10);
        prob[0] = 1;
        for(int i = 0; i < n; ++i)
        {
            double sum = 0;
            for(int j = (i+1)*m; j >= 0; --j)
            {
                sum = sum - prob[j]; if(j-m >= 0) sum += prob[j-m];
                prob[j] = sum / m;
            }
        }        
        double ret = 0;
        for(int i = 0; i <= n*m; ++i)
            ret += prob[i] * max(1, (i-k));
        printf("%.10lf\n", ret);
    }
}

C: Stopped Watches

We are given n analog watches with their faces erased, and hands in equal length. Therefore, a watch can be interpreted into more than one possible time; we need to find a smallest interval where all watches is possibly interpreted.

The solution of this problem is twofold:

  1. Find all the times which can be interpreted from a given watch.
  2. From the list of times, find the smallest interval which contains at least one of each clock's interpretation.

The first subproblem usually boils down to finding a normalizing scheme so that instances in each classes will reduce to a same canonical representation. We want to find a normalizing scheme where times that have equal representations ignoring 1. rotation, 2. interchange of hands. I chose to select the lexicographically smallest representation. After we generate the signatures (canonized form) of all the given clocks, we can iterate through all possible times and see if their canonized form is among the given clocks.

Part 2 is brute-forcible; however, checking all O((246060)^2) intervals are not feasible. This is solved by iterating only through the times that might be represented by a clock.

lang:cpp
#include<iostream>
#include<algorithm>
#include<vector>
#include<cstdio>
#include<set>
using namespace std;

#define all(x) (x).begin(),(x).end()
#define CLEAR(x,with) memset(x,with,sizeof(x))
typedef vector<int> VI;
typedef vector<VI> VVI;

VI normalize(const VI& original)
{
    VI ret;
    for(int i = 0; i < 3; ++i)
    {
        VI rep;
        for(int j = 0; j < 3; ++j)
            rep.push_back((original[(i+j)%3] + 60 - original[i]) % 60);
        sort(all(rep));
        if(ret.empty() || rep < ret) ret = rep;
    }
    return ret;
}

int main()
{
    int n;
    while(scanf("%d", &n) == 1)
    {
        if(n == 0) break;
        map<VI,int> matched;
        for(int i = 0; i < n; ++i)
        {
            VI hand(3);
            for(int j = 0; j < 3; ++j)
                scanf("%d", &hand[j]);
            matched[normalize(hand)] += (1<<i);
        }
        VI match(12*60*60);
        set<int> allTimes;
        for(int time = 0; time < 12*60*60; ++time)
        {
            int hour = time / 720;
            int min = (time / 60) % 60;
            int sec = time % 60;
            vector<int> hands(3);
            hands[0] = hour;
            hands[1] = min;
            hands[2] = sec;
            hands = normalize(hands);
            if(matched.count(hands))
            {
                match[time] = matched[hands];
                allTimes.insert(time);                
            }
        }
        VI times(all(allTimes));
        int sBegin = 0, sDuration = 43199;
        for(int begin = 0; begin < times.size(); ++begin)
        {
            int m = 0;
            for(int span = 0; span < times.size(); ++span)
            {
                m |= match[times[(begin + span) % times.size()]];
                if(m == (1<<n)-1)
                {
                    int dur = (times[(begin+span) % times.size()] + 43200 - times[begin]) % 43200;
                    if(dur < sDuration)
                    {
                        sDuration = dur;
                        sBegin = times[begin];
                        break;
                    }
                }
            }
        }
        int sEnd = (sBegin + sDuration) % 43200;
        printf("%.2d:%.2d:%.2ld %.2d:%.2d:%.2d\n", 
                sBegin/3600,(sBegin/60)%60,sBegin%60,
                sEnd/3600,(sEnd/60)%60,sEnd%60);
    }
}

D: Digits On The Floor

A simple 2D geometry problem, actually a bit similar to a past problem appeared in Algospot Contest. Given many line segments on a 2D plane, we are asked to classify them. The normal forms of line segments are given; the given segments always preserve the connection between segments, only their lengths can be changed.

Using the number of two types of joints (end-to-end, end-to-middle) and the number of segments as a deciding factor, we can pretty much uniquely classify most digits. But 2 and 5 pose a problem; they are mostly same but only direction. Thus, we use cross products. Simple as pie.

lang:cpp
#include<iostream>
#include<algorithm>
#include<vector>
#include<cmath>
#include<cstdio>
#include<queue>
#include<cassert>
using namespace std;

#define all(x) (x).begin(),(x).end()
#define CLEAR(x,with) memset(x,with,sizeof(x))
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef long long ll;

struct point
{
    int y, x;
    point(int y_ = 0, int x_ = 0) : y(y_), x(x_) {}
    bool operator == (const point& rhs) const { return y == rhs.y && x == rhs.x; }
    int cross(const point& rhs) const { return y * rhs.x - x * rhs.y; }
    point operator - (const point& rhs) const { return point(y - rhs.y, x - rhs.x); }
};

int n;
point a[1000], b[1000];
bool seen[1000];

bool onSegment(const point& a, const point& b, const point p)
{
    point ba = b - a;
    point pa = p - a;
    if(ba.cross(pa)) return false;
    return 
        min(a.x, b.x) <= p.x && max(a.x, b.x) >= p.x &&
        min(a.y, b.y) <= p.y && max(a.y, b.y) >= p.y;
}

int meet(int x, int y)
{
    // meet in the end
    if(a[x] == a[y] || a[x] == b[y] || b[x] == a[y] || b[x] == b[y]) return 1;
    // on a segment
    if(onSegment(a[x], b[x], a[y])) return 2;
    if(onSegment(a[x], b[x], b[y])) return 2;
    if(onSegment(a[y], b[y], a[x])) return 2;
    if(onSegment(a[y], b[y], b[x])) return 2;
    return 0;
}

bool turnRight(int x, int y)
{
    point A, B, C;
    if(a[x] == a[y])
    {
        A = b[x];
        B = a[x];
        C = b[y];
    }
    else if(a[x] == b[y])
    {
        A = b[x];
        B = a[x];
        C = a[y];
    }
    else if(b[x] == a[y])
    {
        A = a[x];
        B = b[x];
        C = b[y];
    }
    else if(b[x] == b[y])
    {
        A = a[x];
        B = b[x];
        C = a[y];
    }
    else
    {
        assert(0);
    }
    return (B-A).cross(C-A) > 0;
}

int solve(int here)
{
    seen[here] = true;
    VI sticks(1, here);
    queue<int> q;
    q.push(here);
    while(!q.empty())
    {
        int here = q.front(); q.pop();
        for(int there = 0; there < n; ++there) if(!seen[there] && meet(here, there))
        {
            seen[there] = true;
            q.push(there);
            sticks.push_back(there);
        }
    }

    VI degree(sticks.size(), 0);
    int ones = 0, twos = 0;
    for(int i = 0; i < sticks.size(); ++i)
        for(int j = 0; j < i; ++j)
        {
            int judge = meet(sticks[i], sticks[j]);
            if(judge)
            {
                degree[i] ++;
                degree[j] ++;
                if(judge == 1) ++ones;
                if(judge == 2) ++twos;
            }
        }
    // edge count
    // 1: 1
    // 2: 
    // 3: 4, 7
    // 4: 0, 9, 3
    // 5: 2, 5, 6, 8
    int m = sticks.size();
    if(m == 1) return 1;
    if(m == 3) return (twos ? 4 : 7);
    if(m == 4) 
    {
        if(twos == 0) return 0;
        if(count(all(degree), 3) > 0) return 3;
        return 9;
    }
    if(m == 5)
    {
        if(twos == 1) return 6;
        if(twos == 2) return 8;
        assert(twos == 0);
        assert(ones == 4);
        for(int i = 0; i < m; ++i)
            if(degree[i] == 1)
            {
                for(int j = 0; j < m; ++j) if(i != j && meet(sticks[i], sticks[j]))
                {
                    if(turnRight(sticks[i], sticks[j]))
                        return 2;
                    return 5;
                }
            }
        assert(0);
    }
    return -1;
}

int main()
{    
    while(scanf("%d", &n) == 1)
    {
        if(n == 0) break;
        for(int i = 0; i < n; ++i)
        {
            scanf("%d %d", &a[i].x, &a[i].y);
            scanf("%d %d", &b[i].x, &b[i].y);
        }
        CLEAR(seen,0);
        VI cnt(10);
        for(int i = 0; i < n; ++i)
            if(!seen[i])
                cnt[solve(i)]++;
        for(int i = 0; i < 10; ++i)
        {
            printf("%d ", cnt[i]);
        }
        printf("\n");
    }
}

E: Spherical Mirrors

My favorite problem in this set; ray tracing in 3D space with spheres. Simple linear algebra gives pretty much all the formula needed to solve this problem, including reflection on a sphere.

lang:cpp
#include<iostream>
#include<algorithm>
#include<vector>
#include<cmath>
#include<cstdio>
#include<cassert>
using namespace std;

const double eps = 1e-6;

inline int sgn(double x)
{
    if(fabs(x) < eps) return 0;
    if(x < 0) return -1;
    return 1;
}

struct vector3
{
    double x, y, z;
    vector3(double x_ = 0, double y_ = 0, double z_ = 0) : x(x_), y(y_), z(z_) {}
    vector3 operator + (const vector3& other) const
    {
        return vector3(x + other.x, y + other.y, z + other.z);
    }
    vector3 operator - (const vector3& other) const
    {
        return vector3(x - other.x, y - other.y, z - other.z);
    }
    vector3 operator * (double p) const
    {
        return vector3(x * p, y * p, z * p);
    }
    // dot product
    double operator * (const vector3& b) const
    {
        const vector3& a = *this;
        return a.x * b.x + a.y * b.y + a.z * b.z;
    }
    // cross product
    vector3 operator ^ (const vector3& b) const
    {
        const vector3& a = *this;
        return vector3(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
    }
    double length() const
    {
        return sqrt(x*x + y*y + z*z);
    }
    vector3 normalized() const
    {   
        double l = length();
        return vector3(x / l, y / l, z / l);
    }
    // returns the projection of vector v on this
    vector3 projected(const vector3& v) const
    {
        vector3 n = normalized();
        return n * (n * v);
    }
};
vector3 operator * (double a, const vector3& b) { return b * a; }

int n;
vector3 center[100];
double radius[100];

vector<double> solve2(double a, double b, double c)
{
    double d = b*b - 4*a*c;
    if(d < -eps) return vector<double>();
    if(d < eps) return vector<double>(1, -b / (2*a));
    vector<double> ret;
    ret.push_back((-b + sqrt(d)) / (2*a));
    ret.push_back((-b - sqrt(d)) / (2*a));
    return ret;
}

bool simulate(vector3& here, vector3& direction)
{
    double minT = 1e200;
    int minSphere = -1;
    for(int i = 0; i < n; ++i)
    {
        vector3 v = center[i] - here;
//      (here + t * direction - center)^2 = radius^2        
        vector<double> t = solve2(
                direction * direction,
                2 * (direction * here - direction * center[i]),
                here * here + center[i] * center[i] - 2 * (here * center[i]) - radius[i] * radius[i]);
        double positiveMin = 1e201;
        for(int j = 0; j < t.size(); ++j)
            if(t[j] > eps && t[j] < positiveMin)
                positiveMin = t[j];
        if(positiveMin < minT)
        {
            vector3 p = here + direction * positiveMin;
            assert(sgn((center[i] - p).length() - radius[i]) == 0);
            minT = positiveMin;
            minSphere = i;
        }
    }
    if(minSphere == -1) return false;
    vector3 point_of_contact = here + minT * direction;
    vector3 from_center = point_of_contact - center[minSphere];
    vector3 displacement = here - point_of_contact;
    vector3 there = here - 2 * (displacement - from_center.projected(displacement));
    direction = (there - point_of_contact).normalized();
    here = point_of_contact;
    return true;
}

int main()
{
    while(cin >> n)
    {
        if(n == 0) break;
        vector3 here(0, 0, 0);
        double u, v, w;
        cin >> u >> v >> w;
        vector3 dir(u, v, w);
        for(int i = 0; i < n; ++i)
        {
            double x, y, z;
            cin >> x >> y >> z >> radius[i];
            center[i] = vector3(x, y, z);
        }
        while(simulate(here, dir));
        printf("%.3lf %.3lf %.3lf\n", here.x, here.y, here.z);

    }
}

F: Traveling Cubes

Simulation of cube rolling on a grid; the rest is pretty much BFS. However, its implementation is a bit messy, and very easy to get wrong.

lang:cpp
#include<iostream>
#include<cstring>
#include<algorithm>
#include<string>
#include<vector>
#include<cstdio>
#include<queue>
using namespace std;

#define all(x) (x).begin(),(x).end()
#define CLEAR(x,with) memset(x,with,sizeof(x))
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef long long ll;

#define checkmin(a,b) { if((a)==-1 || (b)<(a)) a=(b); }

const string COLORNAME("rcgmbywk#");
enum COLOR { RED, CYAN, GREEN, MAGENTA, BLUE, YELLOW };
enum SIDE { TOP, BOTTOM, NORTH, SOUTH, EAST, WEST };

const int dx[4] = { 1, -1, 0, 0 };
const int dy[4] = { 0, 0, 1, -1 };
// roll[x][i] replaces roll[x][(i+1)%4]
const SIDE rollReplace[4][4] = { 
    { TOP, EAST, BOTTOM, WEST },
    { TOP, WEST, BOTTOM, EAST },
    { TOP, SOUTH, BOTTOM, NORTH },
    { TOP, NORTH, BOTTOM, SOUTH },
};

int width, height;
int best[31][31][6][6];
int c[31][31][6][6];
string m[31];

COLOR norm[6][6][6];

void roll(COLOR cube[6], const SIDE replace[4])
{
    COLOR last = cube[replace[3]];
    for(int i = 2; i >= 0; --i)
        cube[replace[i+1]] = cube[replace[i]];
    cube[replace[0]] = last;
}

struct state
{
    int y, x, top, north;
    state(int y, int x, int top, int north) : y(y), x(x), top(top), north(north) {}
};

void traverse(char begin, char end)
{
//    printf("================== %c to %c ============================\n", begin, end);
    for(int y = 0; y < height; ++y)
        for(int x = 0; x < width; ++x)
            if(m[y][x] == begin)
            {
                for(int top = 0; top < 6; ++top)
                    for(int north = 0; north < 6; ++north)
                    {
                        if(best[y][x][top][north] == -1) continue;
                        int add = best[y][x][top][north];
    //                    printf("begin (%d,%d) color %c add %d top %c north %c\n", y, x, begin, add, COLORNAME[top], COLORNAME[north]);
/*                        {
                                COLOR cube[6];
                                for(int i = 0; i < 6; ++i)
                                    cube[i] = norm[top][north][i];
                                printf("top %c east %c bottom %c west %c\n", 
                                        COLORNAME[cube[TOP]], 
                                        COLORNAME[cube[EAST]], 
                                        COLORNAME[cube[BOTTOM]], 
                                        COLORNAME[cube[WEST]]);
                        }*/
                        CLEAR(c,-1);
                        queue<state> q;
                        q.push(state(y, x, top, north));

                        c[y][x][top][north] = 0;
                        while(!q.empty())
                        {
                            int north = q.front().north;
                            int top = q.front().top;
                            int x = q.front().x; 
                            int y = q.front().y;
                            q.pop();
/*                            printf("\t(%d,%d) top %c north %c\n", y, x, COLORNAME[top], COLORNAME[north]);
                        {
                                COLOR cube[6];
                                for(int i = 0; i < 6; ++i)
                                    cube[i] = norm[top][north][i];
                                printf("\t\ttop %c east %c bottom %c west %c\n", 
                                        COLORNAME[cube[TOP]], 
                                        COLORNAME[cube[EAST]], 
                                        COLORNAME[cube[BOTTOM]], 
                                        COLORNAME[cube[WEST]]);
                        }*/

                            for(int k = 0; k < 4; ++k)
                            {
                                int ny = y + dy[k], nx = x + dx[k];
                                if(ny < 0 || nx < 0 || ny >= height || nx >= width) continue;
                                if(m[ny][nx] != '#' && m[ny][nx] != 'w' && m[ny][nx] != end) continue;
                                COLOR cube[6];
                                for(int i = 0; i < 6; ++i)
                                    cube[i] = norm[top][north][i];
                                roll(cube, rollReplace[k]);
                                int ntop = cube[TOP];
                                int nnorth = cube[NORTH];

                                if(m[ny][nx] == end && COLORNAME[ntop] != end) continue;

                                int& thereCost = c[ny][nx][ntop][nnorth];
                                if(thereCost == -1)
                                {
                                    thereCost = c[y][x][top][north] + 1;
                                    if(m[ny][nx] == end)
                                    {
                                //        printf("arrived at (%d,%d) with (%d,%d) => %d\n", ny, nx, ntop, nnorth, thereCost + add);
                                        checkmin(best[ny][nx][ntop][nnorth], thereCost + add);
                                    }
                                    q.push(state(ny, nx, ntop, nnorth));
                                }
                            }
                        }                        
                    }
            }
}

int main()
{
    COLOR cube[6] = { RED, CYAN, GREEN, MAGENTA, BLUE, YELLOW };
    SIDE rtype[3][4] = { { TOP, EAST, BOTTOM, WEST}, {TOP, NORTH, BOTTOM, SOUTH}, {NORTH, EAST, SOUTH, WEST}};
    for(int i = 0; i < 4; ++i)
    {
        for(int j = 0; j < 4; ++j)
        {
            for(int k = 0; k < 4; ++k)
            {
                for(int f = 0; f < 6; ++f)
                    norm[ cube[TOP] ][ cube[NORTH] ][f] = cube[f];
                roll(cube, rtype[2]);
            }
            roll(cube, rtype[1]);
        }
        roll(cube, rtype[0]);
    }
    while(cin >> width >> height)
    {
        if(width == 0) break;
        CLEAR(best,-1);
        for(int i = 0; i < height; ++i)
        {
            cin >> m[i];
            for(int j = 0; j < width; ++j)
                if(m[i][j] == '#')
                    best[i][j][RED][GREEN] = 0;
        }
        string order;
        cin >> order;
        order = string("#") + order;
        for(int i = 0; i+1 < order.size(); ++i)
            traverse(order[i], order[i+1]);
        const int INF = 987654321;
        int ret = INF;
        for(int i = 0; i < height; ++i)
            for(int j = 0; j < width; ++j)
                if(m[i][j] == order[order.size()-1])
                {
                    for(int top = 0; top < 6; ++top)
                        for(int north = 0; north < 6; ++north)
                        {
                            if(best[i][j][top][north] != -1)
                                ret = min(best[i][j][top][north], ret);
                        }
                }
        if(ret == INF)
            printf("unreachable\n");
        else
            printf("%d\n", ret);
    }
}

G: Search of Concatenated Strings

Given a text and a number of tokens, return the number of substrings of text which are concatenations of all tokens.

I used Rabin-Karp to find out if there are any matches of each token in the text, and did a pruned backtracking to see if a given substring can be generated by concatenating the given tokens.

lang:cpp
#include<iostream>
#include<algorithm>
#include<string>
#include<vector>
#include<cstdio>
using namespace std;

#define all(x) (x).begin(),(x).end()
#define CLEAR(x,with) memset(x,with,sizeof(x))
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef long long ll;

int n, m;
string element[12], text;
int match[5000];
unsigned char seen[5000][(1<<12)>>3];

const long long M = 100000007;

void scanThrough(const string& text, const string& word, int record[], int id)
{
    if(text.size() < word.size()) return;
    long long key = 0, here = 0;
    long long mod = 1;
    for(int i = 0; i < word.size()-1; ++i)
        mod = (mod * 26) % M;
    for(int i = 0; i < word.size(); ++i)
    {
        key = (key * 26 + word[i] - 'a') % M;
        here = (here * 26 + text[i] - 'a') % M;
    }
    for(int start = 0; start + word.size() <= text.size(); ++start)
    {
        if(key == here)        
            if(strncmp(text.c_str() + start, word.c_str(), word.size()) == 0)
                match[start] |= id;
        here = (((here - (mod * (text[start] - 'a')) % M + M) % M) * 26 + text[start+word.size()] - 'a') % M;
    }
}

bool go(int here, int mask)
{
    if(mask == (1<<n)-1) return true;
    if(seen[here][mask>>3] & (1<<(mask & 7))) return false;
    seen[here][mask>>3] |= (1<<(mask & 7));
    for(int next = 0; next < n; ++next)
        if((mask & (1<<next)) == 0 && (match[here] & (1<<next)))
            if(go(here + element[next].size(), (1<<next) | mask))
                return true;
    return false;
}

int main()
{
    while(cin >> n >> m)
    {
        if(n == 0) break;
        for(int i = 0; i < n; ++i)
            cin >> element[i];
        text = "";
        for(int i = 0; i < m; ++i)
        {
            string tok;
            cin >> tok;
            text += tok;
        }
        CLEAR(match,0);
        for(int i = 0; i < n; ++i)
            scanThrough(text, element[i], match, 1<<i);
        CLEAR(seen,0);
        int tl = 0;
        for(int i = 0; i < n; ++i)
            tl += element[i].size();
        int ret = 0;
        for(int i = 0; i + tl <= text.size(); ++i)
            if(go(i, 0))
            {
                ++ret;
            }
        printf("%d\n", ret);
    }
}

H: Top Spinning

Another geometry; given a closed shape consisting of arcs and line segments, determine its center of gravity.

I was, at first, tempted to solve it by cutting the polygons piece-wise and using Simpson's Rule to approximate the center. However, the problem statement gives some formula for calculating the center of gravity given a circular segment, so I just used it.

The implementation could have been messy, however the input is well normalized to keep it nice. Still, a long-liner with 185 lines.

lang:cpp
#include<iostream>
#include<algorithm>
#include<vector>
#include<sstream>
#include<cmath>
#include<cstdio>
using namespace std;

#define all(x) (x).begin(),(x).end()
#define CLEAR(x,with) memset(x,with,sizeof(x))
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef long long ll;

struct point
{
    double y, x;
    point(double y = 0, double x = 0) : y(y), x(x) {}
    double cross(const point& rhs) const { return y * rhs.x - x * rhs.y; }
    point operator + (const point& rhs) const { return point(y + rhs.y, x + rhs.x); }
    point operator - (const point& rhs) const { return point(y - rhs.y, x - rhs.x); }
    point operator * (double f) const { return point(y * f, x * f); }
    point normalized() const { double d = norm(); return point(y / d, x / d); }
    double norm() const { return hypot(y, x); }
};

int n;
point p[110]; // number of points
double sideRadius[110]; 

double otherSide(double c, double a)
{
    return sqrt(c*c - a*a);
}

point getCenter()
{
    point center(0, 0);
    double totalArea = 0;
    for(int i = 0; i < n; ++i)
    {
        const point& a = p[i];
        const point& b = p[(i+1)%n];
        // triangle
        double area = -a.cross(b) / 2;
        point c = (a+b) * (1.0 / 3.0);
        center = center + c * area;
        totalArea += area;

        if(fabs(sideRadius[i]) > 1e-7)
        {

            double d = (b-a).norm();
            double r = sideRadius[i];
            double theta = acos( (2*r*r - d*d) / (2*r*r) );
            double area = r * r * (theta - sin(theta)) / 2;

            point direction = (b-a).normalized();
            point right(-direction.x, direction.y);
            double y = fabs(2*pow(r,3)*pow(sin(theta/2),3)/(3*area));
            double over = y - otherSide(r, d/2);
            assert(over > 0);

            if(r > 0) over = -over;
            point c = ((a+b) * 0.5) + right * over;

            if(r > 0) area = -area;

            center = center + c * area;
            totalArea += area;
        }
    }
    return center * (1.0 / totalArea);
}

point getArcCenter(const point& a, const point& b, double r)
{
    double d = (b-a).norm();
    double h = otherSide(r, d/2);
    point direction = (b-a).normalized();
    point right(-direction.x, direction.y);
    if(r < 0) h = -h;
    return ((a + b) * 0.5) + right * h;
}

bool intersect(const point& a, const point& b, const point& c, const point& d, point& p)
{
    // y = (dy/dx) * x + (y1 - (dy/dx)*x1)
    double A1 = b.y - a.y, B1 = a.x - b.x, C1 = a.y * (b.x - a.x) - a.x * (b.y - a.y);
    double A2 = d.y - c.y, B2 = c.x - d.x, C2 = c.y * (d.x - c.x) - c.x * (d.y - c.y);
    if(fabs(A1*B2 - A2*B1) < 1e-8) return false;
    // A1 x + B1 y + C1 = 0
    // A2 x + B2 y + C2 = 0
    // (A1*B2 - A2*B1)x + C1*B2 - C2*B1= 0
    p.x = (C2*B1 - C1*B2) / (A1*B2 - A2*B1);
    p.y = (C2*A1 - C1*A2) / (A2*B1 - A1*B2);
    return true;
}

int sgn(double x)
{
    if(x > 1e-8) return 1;
    if(x < -1e-8) return -1;
    return 0;
}

bool onSegment(const point& a, const point& b, const point& p)
{
    return min(a.x, b.x) <= p.x && p.x <= max(a.x, b.x) &&
        min(a.y, b.y) <= p.y && p.y <= max(a.y, b.y);
}

bool inside(const point& P)
{
    // see if it is inside a negative circular segment (then, true)
    for(int i = 0; i < n; ++i) if(sideRadius[i] < -1e-8)
    {
        double r = sideRadius[i];
        const point& a = p[i];
        const point& b = p[(i+1)%n];
        if((b-a).cross(P-a) <= 0) continue;
        point c = getArcCenter(a, b, sideRadius[i]);
        if((P-c).norm() > fabs(r)) continue;
        return true;
    }
    // see if it is inside a positive circular segment (then, false)
    for(int i = 0; i < n; ++i) if(sideRadius[i] > 1e-8)
    {
        double r = sideRadius[i];
        const point& a = p[i];
        const point& b = p[(i+1)%n];
        if((b-a).cross(P-a) >= 0) continue;
        point c = getArcCenter(a, b, sideRadius[i]);
        if((P-c).norm() > fabs(r)) continue;
        return false;
    }

    // see if it is inside the polygon (then, true)

    int cross = 0;
    point Q(P.y + 1111, P.x + 1000), x;
    for(int i = 0; i < n; ++i)
    {
        const point& a = p[i];
        const point& b = p[(i+1)%n];
        if(intersect(a, b, P, Q, x))
        {
            if(onSegment(a, b, x) && onSegment(P, Q, x))
                ++cross;
        }
    }

    return cross % 2 > 0;
}

int main()
{
    string command;
    while(cin >> command)
    {
        if(command == "end") break;
        n = 0;
        CLEAR(sideRadius,0);
        double y, x;
        cin >> x >> y;
        p[n++] = point(y, x);
        while(cin >> command)
        {
            if(command == "close") break;
            double y, x, r = 0;
            cin >> x >> y;
            if(command == "arc")
                cin >> r;
            p[n] = point(y, x);

            sideRadius[n-1] = r;
            n++;
        }
        point c = getCenter();
        if(inside(c))
            printf("%.4lf %.4lf +\n", c.x, c.y);
        else
            printf("%.4lf %.4lf -\n", c.x, c.y);

    }
}

I: Common Polynomials

Given two polynomials, find the GCD of them. It's a bit tricky, we have to parse the string representation, do the GCD, and everything. And finding GCDs of polynomials are not that intuitive.

lang:cpp
#include<iostream>
#include<cstring>
#include<algorithm>
#include<sstream>
#include<string>
#include<cmath>
#include<cstdio>
#include<deque>
using namespace std;

#define all(x) (x).begin(),(x).end()
#define CLEAR(x,with) memset(x,with,sizeof(x))

struct poly
{
    deque<int> coef;

    bool operator == (const poly& rhs) const { return coef == rhs.coef; }
    bool operator < (const poly& rhs) const 
    {
        if(coef.size() != rhs.coef.size()) return coef.size() < rhs.coef.size();
        for(int i = coef.size()-1; i >= 0; --i)
            if(coef[i] != rhs.coef[i])
                return coef[i] < rhs.coef[i];
        return false;
    }

    poly normalized() const 
    {
        int gcd = abs(coef[0]);
        for(int i = 1; i < coef.size(); ++i)
            gcd = __gcd(gcd, abs(coef[i]));
        deque<int> newCoef = coef;
        for(int i = 0; i < newCoef.size(); ++i)
            newCoef[i] /= gcd;
        return newCoef;
    }

    poly(const deque<int>& coef_ = deque<int>(1, 0)) : coef(coef_)
    {
        while(coef.size() > 1 && coef.back() == 0) coef.pop_back();
    }

    poly operator * (int x) const 
    {
        deque<int> newCoef = coef;
        for(int i = 0; i < newCoef.size(); ++i)
            newCoef[i] *= x;
        return newCoef;
    }

    poly operator * (const poly& rhs) const
    {
        int n = coef.size(), m = rhs.coef.size();
        deque<int> newCoef(n + m + 2);
        for(int i = 0; i < n; ++i)
            for(int j = 0; j < m; ++j)
                newCoef[i+j] += coef[i] * rhs.coef[j];
        return newCoef;
    }

    poly operator + (const poly& rhs) const
    {
        deque<int> newCoef(max(coef.size(), rhs.coef.size()));
        for(int i = 0; i < newCoef.size(); ++i)
        {
            if(i < coef.size()) newCoef[i] += coef[i];
            if(i < rhs.coef.size()) newCoef[i] += rhs.coef[i];
        }
        return newCoef;
    }
    poly operator - (const poly& rhs) const
    {
        deque<int> newCoef(max(coef.size(), rhs.coef.size()));
        for(int i = 0; i < newCoef.size(); ++i)
        {
            if(i < coef.size()) newCoef[i] += coef[i];
            if(i < rhs.coef.size()) newCoef[i] -= rhs.coef[i];
        }
        return newCoef;
    }

    poly pow(int n) const
    {
        poly ret(deque<int>(1,1)); // one
        for(int i = 0; i < n; ++i)
            ret = (ret * *this);
        return ret;
    }

    poly operator ^ (int x) const
    {
        deque<int> newCoef = coef;
        while(x > 0)
        {
            newCoef.push_front(0);
            --x;
        }
        while(x < 0 && newCoef.size())
        {
            newCoef.pop_front();
            ++x;
        }
        if(newCoef.empty()) newCoef.push_back(0);
        return newCoef;
    }

    string str() const
    {
        ostringstream outp;
        if(coef.size() == 1 && coef[0] == 0) return "0";
        for(int i = coef.size()-1; i >= 0; --i)
        {
            if(coef[i] == 0) continue;
            if(i + 1 != coef.size() && coef[i] > 0)
                outp << "+";
            if(coef[i] < 0)
                outp << "-";
            if(abs(coef[i]) != 1 || i == 0)
            {
                outp << abs(coef[i]);
            }
            if(i == 0) continue;
            outp << "x";
            if(i == 1) continue;
            outp << "^" << i;
        }
        return outp.str();
    }
};

poly operator % (poly a, poly b)
{
    if(a < b) return a;
    int addexp = a.coef.size() - b.coef.size();
    int lcm = a.coef.back() / __gcd(a.coef.back(), b.coef.back()) * b.coef.back();
    poly res = a * (lcm / a.coef.back()) - (b^addexp) * (lcm / b.coef.back());
    if(res.coef == deque<int>(1,0)) return res;
    poly ret = res % b;
    //printf("\t%s %% %s = %s\n", a.str().c_str(), b.str().c_str(), ret.str().c_str());
    return ret;

}

poly expr(const string& f);
poly term(const string& f)
{
    if(f.size() == 0)
    {
        poly one(deque<int>(1,1));
        return one;
    }
    if(f[0] == '-')
    {
        return term(f.substr(1)) * -1;
    }
    else if(f[0] == '(')
    {
        int opened = 1;
        for(int i = 1; i < f.size(); ++i)
            if(f[i] == '(')
                ++opened;
            else if(f[i] == ')')
            {
                --opened;
                if(opened == 0)
                {
                    poly head = expr(f.substr(1, i-1));
                    // check tail for exponent
                    int p = i+1;
                    if(f[p] == '^')
                    {
                        int exp = 0;
                        ++p;
                        while(isdigit(f[p]))
                        {
                            exp = exp * 10 + f[p] - '0';
                            ++p;
                        }  
                        head = head.pow(exp);
                    }
                    return head * term(f.substr(p));
                }
            }
    }
    else if(f[0] == 'x')
    {
        deque<int> coef(2);
        coef[1] = 1;
        poly head(coef);
        // check tail for exponent
        int p = 1;
        if(f[p] == '^')
        {
            int exp = 0;
            ++p;
            while(isdigit(f[p]))
            {
                exp = exp * 10 + f[p] - '0';
                ++p;
            }  
            head = head.pow(exp);
        }
        poly rest = term(f.substr(p));
        return head * rest;
    }
    int p = 0, num = 0;
    while(isdigit(f[p]))
    {
        num = num * 10 + f[p] - '0';
        ++p;
    }
    if(f[p] == '^')
    {
        int exp = 0;
        ++p;
        while(isdigit(f[p]))
        {
            exp = exp * 10 + f[p] - '0';
            ++p;
        }  
        int power = 1;
        for(int i = 0; i < exp; ++i)
            power *= num;
        num = power;
    }
    return term(f.substr(p)) * num;
}

poly expr(const string& f)
{
    int opened = 0; 
    for(int i = f.size()-1; i >= 0; --i)
        if(f[i] == ')')
            ++opened;
        else if(f[i] == '(')
            --opened;
        else if(f[i] == '+' && opened == 0)
            return expr(f.substr(0, i)) + term(f.substr(i+1));
        else if(i > 0 && f[i] == '-' && opened == 0)
            return expr(f.substr(0, i)) + term(f.substr(i+1)) * -1;
    return term(f);
}

poly gcd(const poly& a, const poly& b)
{
    poly r = a%b;
    if(r.coef == deque<int>(1,0)) return b;
    poly ret = gcd(b, a%b);
//    cout << "gcd("<<a.str()<<", "<<b.str() << ") = " << ret.str() << endl;
    return ret;
}

int main()
{
    string a, b;
    while(cin >> a >> b)
    {
        if(a == ".") break;
        poly pa = expr(a);
        poly pb = expr(b);
        //cout << "given " << pa.str() << " and " << pb.str() << endl;
        poly g = gcd(pa, pb);
        if(g.coef.back() < 0) g = g * -1;
        cout << g.normalized().str() << endl;
    }
}

J: Zigzag

The ultimate long-liner, and another geometry. Given some points on a 2D plane, find the shortest zigzag curve which passes through all these points. It might be daunting at first, however the constraints are rather severe and it greatly reduces search space. Actually, to an exponential one which is still searchable with Dijskstra.

lang:cpp
#include<iostream>
#include<algorithm>
#include<string>
#include<vector>
#include<cmath>
#include<cstdio>
#include<cassert>
#include<queue>
using namespace std;

#define FOR(i,a,b) for(int i = (a); i < (b); ++i)
#define REP(i,n) FOR(i,0,n)
#define FORE(it,x) for(typeof(x.begin()) it=x.begin();it!=x.end();++it)
#define pb push_back
#define all(x) (x).begin(),(x).end()
#define CLEAR(x,with) memset(x,with,sizeof(x))
#define sz size()
typedef long long ll;

const double eps = 1e-6;

struct Point
{
    int y, x;
    Point(int y_ = 0, int x_ = 0) : y(y_), x(x_) {}
    bool operator < (const Point& b) const { if(y != b.y) return y < b.y; return x < b.x; }
    bool operator == (const Point& b) const { return y == b.y && x == b.x; }
    Point operator + (const Point& b) const { return Point(y + b.y, x + b.x); }
    Point operator - (const Point& b) const { return Point(y - b.y, x - b.x); }
    double norm() const { return hypot(y, x); }
};

inline int sgn(int x)
{
    if(x == 0) return 0;
    if(x < 0) return -1;
    return 1;
}

int ccw(const Point& a, const Point& b, const Point& c)
{
    int by = b.y - a.y, bx = b.x - a.x;
    int cy = c.y - a.y, cx = c.x - a.x;
    return bx * cy - cx * by;
}

typedef pair<int,double> Cost;
struct State
{
    State(int v = 0, int h = 0, int l = 0) : visited(v), here(h), last(l) {}
    int visited;
    int here, last;
    bool operator < (const State& other) const
    {
        if(visited != other.visited) return visited < other.visited;
        if(here != other.here) return here < other.here;
        return last < other.last;
    }
};

bool getIntersection(const Point& p1, const Point& p2, const Point& q1, const Point& q2, double& y, double& x)
{
    int a1, b1, c1, a2, b2, c2;
    a1 = p2.y - p1.y;
    b1 = p1.x - p2.x;
    c1 = (p2.x - p1.x) * p1.y - (p2.y - p1.y) * p1.x;

    a2 = q2.y - q1.y;
    b2 = q1.x - q2.x;
    c2 = (q2.x - q1.x) * q1.y - (q2.y - q1.y) * q1.x;

    if(a2*b1 == b2*a1) return false;

    // a1*x + b1*y + c1 = 0
    // a2*x + b2*y + c2 = 0
    y = (c2*a1 - c1*a2) / double(b1*a2 - b2*a1);
    x = (c2*b1 - c1*b2) / double(a1*b2 - a2*b1);
    return true;
}

typedef pair<Cost,State> Pair;

struct Comparator
{
    bool operator () (const Pair& a, const Pair& b) const
    {
        if(a.first != b.first) return a.first > b.first;
        return a.second < b.second;
    }
};

int n;
Point p[10];
Cost c[1<<10][10][10];
double dist[10][10];

#define C(stat) c[stat.visited][stat.here][stat.last]

bool sort_by_p(int a, int b)
{
    return p[a] < p[b];
}

struct DistanceComparator
{
    double y, x;
    DistanceComparator(double y_, double x_) : y(y_), x(x_) {}
    bool operator () (int a, int b) const
    {
        return hypot(y - p[a].y, x - p[a].x) < hypot(y - p[b].y, x - p[b].x);
    }
};

int main()
{
    while(cin >> n)
    {
        if(n == 0) break;
        REP(i,n) 
        {
            int y, x;
            cin >> x >> y;
            p[i] = Point(y, x);
        }
        REP(i,n) REP(j,n) dist[i][j] = (p[i]-p[j]).norm();
        REP(i,1<<n) REP(j,n) REP(k,n) c[i][j][k] = Cost(9999, 0);
        priority_queue<pair<Cost,State>, vector<pair<Cost,State> >, Comparator > pq;

        // Find all pairwise segments..
        REP(i,n)
            FOR(j,i+1,n)
            {
                vector<int> collinear(2);
                collinear[0] = i;
                collinear[1] = j;
                REP(k,n)
                    if(i != k && j != k && ccw(p[i], p[j], p[k]) == 0)
                    {
                        if(k < i) 
                        {
                            collinear.clear();
                            break;
                        }
                        else
                            collinear.pb(k);
                    }
                if(collinear.sz == 0) continue;
                sort(all(collinear), sort_by_p);
                REP(a,collinear.sz) 
                {
                    int v = 1 << collinear[a];
                    FOR(b,a+1,collinear.sz)
                    {
                        v += 1 << collinear[b];

                        Cost cost(0, (p[collinear[b]]-p[collinear[a]]).norm());
                        State here(v, collinear[b], collinear[a]);
                        C(here) = cost;
                        pq.push(make_pair(cost, here));
                        here = State(v, collinear[a], collinear[b]);
                        C(here) = cost;
                        pq.push(make_pair(cost, here));
                    }
                }
            }
        Cost minCost(9999, 0);
        // yay, dijkstra
        while(!pq.empty())
        {
            Cost cost = pq.top().first;
            State here = pq.top().second;
            int A = here.last, B = here.here;
            if(here.visited == (1<<n)-1)                
            {
                minCost = cost;
                break;
            }
            pq.pop();
            if(C(here) < cost) continue;
            if(cost.first == 4) continue;
            // option 1: just make segment (B, i)
            REP(i,n) if((here.visited & (1<<i)) == 0)
            {
                // going straight doesn't make sense
                if(A != -1 && ccw(p[A], p[B], p[i]) == 0)
                    continue;
                State there(here.visited + (1<<i), i, B);
                Cost thereCost(cost.first + 1, cost.second + dist[i][B]);
                if(C(there) > thereCost) 
                {
                    C(there) = thereCost;                   
                    pq.push(make_pair(thereCost, there));
                }
            }
            // option 2: make a segment
            REP(i,n) if((here.visited & (1<<i)) == 0)
                FOR(j,i+1,n) if((here.visited & (1<<j)) == 0)
                {
                    // going straight doesn't make sense
                    if(ccw(p[A], p[B], p[i]) == 0 && ccw(p[A], p[B], p[j]) == 0)
                        continue;

                    // get the point of intersection
                    double iy, ix;
                    bool intersects = getIntersection(p[A], p[B], p[i], p[j], iy, ix);
                    // shit, parallel
                    if(!intersects) continue;
                    if(ccw(p[B], p[i], p[j]) != 0)
                    {
                        int minX = min(p[A].x, p[B].x), maxX = max(p[A].x, p[B].x);
                        int minY = min(p[A].y, p[B].y), maxY = max(p[A].y, p[B].y);
                        // cannot go back
                        if(minX-eps < ix && ix < maxX+eps && minY-eps < iy && iy < maxY+eps) continue;
                    }

                    // cannot be closer to tail
                    if(hypot(iy - p[A].y, ix - p[A].x) < hypot(iy - p[B].y, ix - p[B].x)) continue;

                    vector<int> collinear(2);
                    collinear[0] = i;
                    collinear[1] = j;
                    REP(k,n) if(i != k && j != k && (here.visited & (1<<k)) == 0 
                        && ccw(p[i], p[j], p[k]) == 0)
                    {
                        if(k < i)
                        {
                            collinear.clear();
                            break;
                        }
                        else
                            collinear.pb(k);
                    }
                    if(collinear.sz == 0) continue;
                    // according to sign of ccw
                    vector<int> classified[2];
                    REP(k,collinear.sz)
                    {
                        int sign = sgn(ccw(p[A], p[B], p[collinear[k]]));
                        if(sign <= 0) classified[0].pb(collinear[k]);
                        if(sign >= 0) classified[1].pb(collinear[k]);
                    }

                    // sort both classified[1] and classified[0] with the distance from the intersection point
                    DistanceComparator comp(iy, ix);
                    sort(all(classified[0]), comp);
                    sort(all(classified[1]), comp);
                    double add = hypot(p[B].y - iy, p[B].x - ix);
                    // turn left? turn right? at the intersection point
                    REP(idx,2)
                    {
                        int visited = here.visited;
                        REP(k,classified[idx].sz)
                        {
                            visited += (1<<classified[idx][k]);
                            if(k > 0)
                            {
                                State there(visited, classified[idx][k], classified[idx][k-1]);
                                double D = hypot(p[classified[idx][k]].y - iy, p[classified[idx][k]].x - ix);
                                Cost thereCost(cost.first + 1, cost.second + add + D);
                                if(C(there) > thereCost) 
                                {
                                    C(there) = thereCost;
                                    pq.push(make_pair(thereCost, there));
                                }
                            }
                        }
                    }
                }
        }
        printf("%d %.8lf\n", minCost.first, minCost.second);
    }
}
2009-05-07 16:45:38 | JM | /logs/ | 3 Comments
wook
2009-05-08 01:19:31
역시 이분은 폭풍간지군요
2009-05-09 00:49:05
아이즈 정ㅋ벅ㅋ
LIBe
2009-05-09 02:21:20
The Pretender!!!

Leave a comment

春來不似春

About

Eventstream

Pages

Guestbook

Search

Site Admin

Recent Comments