naoya_t@hatenablog

いわゆるチラシノウラであります

SRM734 Div1 Easy(300): TheRoundCityDiv1

最近の(出てない)SRMの300点問題。agwたんが話してたので見てみた。

円の中に含まれる(原点は含まないが円周上は含む)グリッドの点で、原点から「見える」点の数を求める問題。
x,y座標(ともに整数)が互いに素であるペアの数。r≦10^6で。
(r=1なら4を返して終わり。)
対称な8つの区画(各象限を45度でさらに半分に切ったもの)に切って、区画の線上の原点から見える点(8つ)は自明として、区画内(オンラインなものは除く)にある点を数える。(それをあとで8倍して足す)。
x>0,y>0x^2+y^2\le r^2y\lt xな整数の組み合わせなんだけど、全部見ていたらO(r^2)で時間が足りない。
あるyに対し、それより大きくて \sqrt{r^2-y^2}以内のxのうち、x,yが互いに素なもの、すなわちxがyの約数のどれかで割り切れたりしない数の個数を数える。
例えばy=6なら 2の倍数、3の倍数 を排除していけばいいのだけれど、それだと 2x3=6 の倍数は2回排除されることになる。そこで包除原理の出番。
rが100万までだと素因数の数は高々7つ(重複除く)だし、(大抵は2つか3つの素数で出来てる)、全パターン見ても高々127通りだしなんとかなるでしょう…
素因数分解そのものはsieveするときに最小(ないし最大)の素因数を記録しておくことでするすると分解できる。


system testで落ちてて、何でだろうと思ったら\sqrt{r^2-y^2}を計算するときにrとyがintのままだった。そこだけlong longに直して再投稿してAC。(いつもながら最後の詰めが甘い系のケアレスミス…)

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef vector<int> vi;

#define sz(a)  int((a).size())
#define pb  push_back
#define rep(var,n)  for(int var=0;var<(n);++var)
#define ALL(c)  (c).begin(),(c).end()

#include <vector>
using namespace std;

vector<int> primes;
vector<int> smallest_prime_factor;

int sieve(int nmax){
    primes.clear();
    smallest_prime_factor.assign(nmax+1, 0);

    for (int n=2; n<=nmax; ++n) {
        if (!smallest_prime_factor[n]) {
            primes.push_back(n);
            for (int kn=n; kn<=nmax; kn+=n) {
                if (!smallest_prime_factor[kn]) {  // ←このチェックをせずに上書きを続けると largest_prime_factor[] が出来上がる
                    smallest_prime_factor[kn] = n;
                }
            }
        }
    }
    return primes.size();
}

// xを素因数分解
vector<int> factorize(int x, bool uniq=false) {
    vector<int> prime_factors;
    while (x > 1) {
        int s = smallest_prime_factor[x];
        if (!s) break;
        if (uniq) {
            if (prime_factors.empty() || prime_factors.back() != s)
                prime_factors.push_back(s);
        } else {
            prime_factors.push_back(s);
        }
        x /= s;
    }
    return prime_factors;
}

inline int bai(int from, int to, int x){
    return to/x - (from-1)/x;
}

inline int bitpart(vi& fy, int pat){
    int L = fy.size();
    int x = 1;
    for (int i=0,m=1; i<L; ++i,m<<=1) {
        if (m&pat) x *= fy[i];
    }
    return x;
}

class TheRoundCityDiv1 {
    public:
    long long calc(int y, int xmax, vi& fy) {
        // |fy| <= 7
        ll ans = xmax - y; //(y-1);

        int fyL = fy.size();
        int P = 1 << fyL;
        ll a = 0;
        for(int pat=1; pat<P; ++pat){
            int b = __builtin_popcount(pat);
            int pm = (b % 2) ? -1 : 1;
            int ff = bitpart(fy, pat);
            a += pm * bai(y+1, xmax, ff);
        }
        return ans + a;
    }
    long long find(int r) {
        if (r == 1) return 4LL;

        // r >= 2
        ll a = 0;

        a += r-2; // for x=1
        int rSq = r * sqrt(0.5);
        int P = sieve(rSq+1);
        for (int y=2; y<=rSq; ++y) {
            vector<int> fy = factorize(y, true);
            // x > y
            int xmax = sqrt((ll)r*r - (ll)y*y);  /// ←←ここ(ll)するの忘れててfailed
            a += calc(y, xmax, fy);
        }
        return 8LL + a*8;
    }
};