中国剰余定理 (CRT) 使いました記念

// 約数一覧を返します
vector<ll> calc_devisors(ll a){
    vector<ll> ret;
    for(ll i=1; i*i<=a; i++){
        if(a%i!=0) continue;
        ret.push_back(i);
        ll another = a/i;
        if(i!=another) ret.push_back(another);
    }
    sort(ALL(ret), greater<ll>()); // 使いやすいように大きい順にソートしておく
    return ret;
}

// drken's lib
// https://qiita.com/drken/items/ae02240cd1f8edfc86fd
// 負の数にも対応した mod
// 例えば -17 を 5 で割った余りは本当は 3 (-17 ≡ 3 (mod. 5))
// しかし単に -17 % 5 では -2 になってしまう
inline long long super_mod(long long a, long long m) {
    return (a % m + m) % m;
}

// 拡張 Euclid の互除法
// ap + bq = gcd(a, b) となる (p, q) を求め、d = gcd(a, b) をリターンします
long long extGcd(long long a, long long b, long long &p, long long &q) {  
    if (b == 0) { p = 1; q = 0; return a; }  
    long long d = extGcd(b, a%b, q, p);  
    q -= a/b * p;  
    return d;  
}

// 中国剰余定理
// リターン値を (r, m) とすると解は x ≡ r (mod. m)
// 解なしの場合は (0, -1) をリターン
// 入力:
// m1で割った余りがb1
// m2で割った余りがb2
PII ChineseRem(long long b1, long long m1, long long b2, long long m2) {
  long long p, q;
  long long d = extGcd(m1, m2, p, q); // p is inv of m1/d (mod. m2/d)
  if ((b2 - b1) % d != 0) return make_pair(0, -1);
  long long m = m1 * (m2/d); // lcm of (m1, m2)
  long long tmp = (b2 - b1) / d * p % (m2/d);
  long long r = super_mod(b1 + m1 * tmp, m);
  return make_pair(r, m);
}

int main(){
    cin.tie(0);
    ios::sync_with_stdio(false);

    // input
    ll N; 
    cin>>N;
    
    auto Y = calc_devisors(2*N);
    ll mi=inf;

    for(ll a : Y){
      ll b = 2*N/a;
      // aで割って0
      // bで割って-1であるものは?
      auto pa = ChineseRem(0,a,-1,b);
      if(pa.first==0)continue;
      ll k = pa.first;
      ll m = pa.second;
      chmin(mi,k);
    }
    p(mi);
    
    return 0;
}

python3で整数どうしの割り算に注意。片方が負のときの挙動がC++と違う

>>> -1/3
-0.3333333333333333

>>> -1//3
-1

>>> 1//3
0
  • スラッシュ1つだと小数点になることがある。これは基本
  • スラッシュ2つだと整数の割り算で答えも整数になるが、-1//3が-1になっており、これはC++だと0が返ってくる。挙動の違いに注意

除算値を超えない最大の整数

実際に問題を解きたい人は

点Pから円に引いた2本の接線が求まるが、その2つの接点を求めよ (polar) "Tangent to a Circle" AOJ

f:id:peroon:20200925053205p:plain

typedef complex<ld> C;

// ベクトルの長さ
ld length(C v){
  return sqrt(v.real()*v.real() + v.imag()*v.imag());
}

// 2点間の距離
ld distance_between_two_points(C a, C b){
  return length(a-b);
}

// 単位ベクトル
C unit_vector(C v){
  return v / length(v);
}

int main(){
    cin.tie(0);
    ios::sync_with_stdio(false);

    // input
    ld px,py;cin>>px>>py;
    C p = {px,py};

    ld cx,cy,r;cin>>cx>>cy>>r;
    C c = {cx,cy};

    // p-c間距離
    ld d = distance_between_two_points(p,c);
    ld sin_theta = r / d;
    ld theta = asin(sin_theta);

    C v = c-p;
    v *= polar(1.0, (double)theta);
    auto u = unit_vector(v);
    auto ans0 = p + u * d * cos(theta);

    theta *= -1;
    v = c-p;
    v *= polar(1.0, (double)theta);
    u = unit_vector(v);
    auto ans1 = p + u * d * cos(theta);

    vector<pair<ld,ld>> V;
    V.push_back({ans0.real(), ans0.imag()});
    V.push_back({ans1.real(), ans1.imag()});
    sort(ALL(V));

    rep(i,2){
      ld x = V[i].first;
      ld y = V[i].second;
      printf("%.10Lf %.10Lf\n", x, y);
    }
    
    return 0;
}

外積で凸性判定

f:id:peroon:20200925035135p:plain

typedef complex<ld> C;

// 外積
ld cross(C a, C b){
  return a.real()*b.imag() - a.imag()*b.real();
}

// 凸性判定
// https://imagingsolution.net/math/calc_n_point_area/
bool is_convex(vector<C> V){
  ll N = V.size();
  V.push_back(V[0]);
  V.push_back(V[1]);
  rep(i,N){
    C a = V[i+1]-V[i];
    C b = V[i+2]-V[i+1];
    if(cross(a,b)<0) return false;
  }
  return true;
}

int main(){
    cin.tie(0);
    ios::sync_with_stdio(false);

    // input
    ll N; cin>>N;
    vector<C> V;
    rep(i,N){
      ld x,y;cin>>x>>y;
      V.push_back({x,y});
    }
    
    if(is_convex(V)){
      p(1);
    }else{
      p(0);
    }
    
    return 0;
}

多角形の面積(凸じゃなくても可)は外積で求まる

typedef complex<ld> C;

// 外積
ld cross(C a, C b){
  return a.real()*b.imag() - a.imag()*b.real();
}

// 多角形の面積
// 凸に限らない
// https://imagingsolution.net/math/calc_n_point_area/
ld polygon_area(vector<C> V){
  ll N = V.size();
  // 原点,0,1
  // 原点,1,2
  // ...
  // 原点,N-1,N(0)
  V.push_back(V[0]);
  ld sum=0;
  C g = {0,0}; // 原点
  rep(i,N){
    C a = V[i]-g;
    C b = V[i+1]-g;
    sum += cross(a,b);
  }
  sum = abs(sum);
  return sum/2;
}

int main(){
    cin.tie(0);
    ios::sync_with_stdio(false);

    // input
    ll N; cin>>N;
    vector<C> V;
    rep(i,N){
      ld x,y;cin>>x>>y;
      V.push_back({x,y});
    }
    
    ld area = polygon_area(V);
    printf("%.1Lf\n",area);
    
    return 0;
}