#803557 (Rust) No.144 エラトステネスのざる

提出ソース

結果

問題	No.144 エラトステネスのざる
ユーザー	Ryota Bannai
提出日時	2022-10-02 00:33:59
言語	Rust (1.77.0 + proconio)
結果	AC
実行時間	64 ms / 2,000 ms
コード長	28,033 bytes
コンパイル時間	12,355 ms
コンパイル使用メモリ	402,076 KB
実行使用メモリ	6,944 KB
最終ジャッジ日時	2024-06-06 14:46:08
合計ジャッジ時間	12,595 ms
ジャッジサーバーID （参考情報）	judge4 / judge5

このコードへのチャレンジ
（要ログイン）

テストケース

テストケース表示

入力	結果	実行時間実行使用メモリ
testcase_00	AC	1 ms 6,816 KB
testcase_01	AC	1 ms 6,944 KB
testcase_02	AC	1 ms 6,944 KB
testcase_03	AC	1 ms 6,940 KB
testcase_04	AC	1 ms 6,940 KB
testcase_05	AC	1 ms 6,944 KB
testcase_06	AC	1 ms 6,944 KB
testcase_07	AC	1 ms 6,940 KB
testcase_08	AC	1 ms 6,944 KB
testcase_09	AC	0 ms 6,940 KB
testcase_10	AC	1 ms 6,944 KB
testcase_11	AC	1 ms 6,944 KB
testcase_12	AC	1 ms 6,944 KB
testcase_13	AC	55 ms 6,940 KB
testcase_14	AC	64 ms 6,944 KB
testcase_15	AC	63 ms 6,940 KB
testcase_16	AC	61 ms 6,940 KB
testcase_17	AC	62 ms 6,940 KB
testcase_18	AC	62 ms 6,940 KB
testcase_19	AC	49 ms 6,940 KB

権限があれば一括ダウンロードができます

コンパイルメッセージ

warning: unused import: `marker::Chars`
 --> src/main.rs:8:32
  |
8 | use proconio::{fastout, input, marker::Chars};
  |                                ^^^^^^^^^^^^^
  |
  = note: `#[warn(unused_imports)]` on by default

ソースコード

raw source code

pub use __cargo_equip::prelude::*;

/**
 * @cpg_dirspec eratosthenes_zaru
 *
 * cpg run -p src/bin/dp/probability/eratosthenes_zaru.rs
 */
use proconio::{fastout, input, marker::Chars};
// use std::cmp::{max, min};
// use superslice::Ext;
// use ac_library_rs::modint::ModInt998244353 as Mint;
// use superslice::{self, Ext};
// use derive_new::new;
// use indexmap::indexmap;
// use std::collections::{BTreeMap, BTreeSet};
// type Map = BTreeMap<String, usize>;
// type Set = BTreeSet<(usize, usize)>;
// use easy_ext::ext;
// use std::collections::{BinaryHeap, VecDeque};
// use collection::{geo_lib::*, utils::read};

/**
 * 確率DP HSI
 *
 * https://yukicoder.me/problems/no/144
 *
 * tags: #DP #確率DP
 *
 * 2つリストがある. これらをA、Bをする.
 * Aは自然数x {x|2<=x<=N} であり、このリストから、最小mi を選んでBに移すことを考える.
 * この時、確率p でmi の倍数をAから削除する.
 * Aが空{} になるまで、移動すること操作を行うとき、Bの要素数の期待値を求めよ.
 *
 *
 * 「Bの要素数の期待値」について何を求めるか？
 * -> 各要素がB に移される期待値の総和
 * -> 各要素は独立に計算できて、それぞれカウントは１だから、各要素がB に移される時の確率を求めれば良い！（ei=1*pi=pi {pi: i∈N}）
 * -> 倍数が削除されるというのは、削除される p の余事象（1-p）
 * -> 約数全てが削除されないなら、自身が最小となる時にB へ移される
 *
 * 6 が B に移される
 * =2,3 の倍数が削除されない (1-p)^2 2:= １以上約数 - 1(=自分)
 *
*/

#[fastout]
fn main() {
    input! {
      n: usize,
      p: f64
    }

    let mut div = vec![0; n + 1];

    let mut e = 0.;
    for k in 2..=n {
        for m in (k..=n).step_by(k) {
            div[m] += 1;
        }
        e += (1. - p).powf((div[k] - 1) as f64); // 約数が移動されない確率で良い
    }

    println!("{:.5}", e);
}

// The following code was expanded by `cargo-equip`.

///  # Bundled libraries
/// 
///  - `collection 0.1.0 (path+████████████████████████████████████████████████████████████████████████)` published in **missing** licensed under **missing** as `crate::__cargo_equip::crates::collection`
#[cfg_attr(any(), rustfmt::skip)]
#[allow(unused)]
mod __cargo_equip {
    pub(crate) mod crates {
        pub mod collection {pub mod geo_lib { use std :: { cmp :: { Ordering , Ordering :: { Equal , Greater , Less } , } , f64 :: { consts :: PI , NAN } , ops :: { Add , Div , Mul , Sub } , } ; # [derive (Clone , Copy , Debug)] pub struct Vector { pub x : f64 , pub y : f64 , } impl PartialEq for Vector { fn eq (& self , other : & Self) -> bool { let eps = 1e-10 ; (self . x - other . x) . abs () < eps && (self . y - other . y) . abs () < eps } } impl Ord for Vector { fn cmp (& self , other : & Self) -> Ordering { let eps = 1e-10 ; if (self . x - other . x) . abs () < eps { if (self . y - other . y) . abs () < eps { Equal } else if self . y < other . y { Less } else { Greater } } else if self . x < other . x { Less } else { Greater } } } impl PartialOrd for Vector { fn partial_cmp (& self , other : & Self) -> Option < Ordering > { Some (self . cmp (other)) } } impl Eq for Vector { } impl Sub for Vector { type Output = Self ; fn sub (self , other : Vector) -> Self { Self :: new (self . x - other . x , self . y - other . y) } } impl Add for Vector { type Output = Self ; fn add (self , other : Self) -> Self { Self :: new (self . x + other . x , self . y + other . y) } } impl Mul < f64 > for Vector { type Output = Self ; fn mul (self , k : f64) -> Self { Self :: new (self . x * k , self . y * k) } } impl Div < f64 > for Vector { type Output = Self ; fn div (self , k : f64) -> Self { Self :: new (self . x / k , self . y / k) } } impl Vector { pub fn new (x : f64 , y : f64) -> Self { Self { x , y } } pub fn from ((x , y) : (f64 , f64)) -> Self { Self { x , y } } pub fn dot (self , other : Vector) -> f64 { VectorFns :: dot (self , other) } pub fn cross (self , other : Vector) -> f64 { VectorFns :: cross (self , other) } pub fn norm (self) -> f64 { VectorFns :: norm (self) } pub fn cmp_y (self , other : Vector) -> Ordering { VectorFns :: cmp_y (self , other) } pub fn unit (self) -> Self { VectorFns :: unit (Vector :: new (0.0 , 0.0) , self) } } pub struct VectorFns { } impl VectorFns { pub fn unit (v1 : Vector , v2 : Vector) -> Vector { let nv = v2 . sub (v1) ; nv . div (Self :: norm (nv)) } pub fn norm (v : Vector) -> f64 { Self :: dot (v , v) . sqrt () } pub fn abs (v1 : Vector , v2 : Vector) -> f64 { Self :: dot (v1 , v2) . sqrt () } pub fn dot (v1 : Vector , v2 : Vector) -> f64 { v1 . x * v2 . x + v1 . y * v2 . y } pub fn rad (v1 : Vector , v2 : Vector) -> f64 { (Self :: dot (v1 , v2) / (Self :: norm (v1) * Self :: norm (v2))) . acos () } pub fn cross (v1 : Vector , v2 : Vector) -> f64 { v1 . x * v2 . y - v1 . y * v2 . x } pub fn cmp_y (v1 : Vector , v2 : Vector) -> Ordering { let eps = 1e-10 ; if (v1 . y - v2 . y) . abs () < eps { if (v1 . x - v2 . x) . abs () < eps { Equal } else if v1 . x < v2 . x { Less } else { Greater } } else if v1 . y < v2 . y { Less } else { Greater } } pub fn is_orthogonal (v1 : Vector , v2 : Vector , u1 : Vector , u2 : Vector) -> bool { let eps = 1e-10 ; let nv = v1 . sub (v2) ; let nu = u1 . sub (u2) ; Self :: dot (nv , nu) . abs () < eps } pub fn is_parallel (v1 : Vector , v2 : Vector , u1 : Vector , u2 : Vector) -> bool { let eps = 1e-10 ; let nv = v1 . sub (v2) ; let nu = u1 . sub (u2) ; Self :: cross (nv , nu) . abs () < eps } pub fn projection (v : Vector , v1 : Vector , v2 : Vector) -> Vector { let base = v2 . sub (v1) ; let hypo = v . sub (v1) ; let t = Self :: dot (hypo , base) / Self :: norm (base) ; let r = t / Self :: norm (base) ; v1 . add (base . mul (r)) } pub fn reflection (v : Vector , v1 : Vector , v2 : Vector) -> Vector { v . add ((Self :: projection (v , v1 , v2) . sub (v)) . mul (2.0)) } pub fn distance_vv (v1 : Vector , v2 : Vector) -> f64 { Self :: abs (v1 , v2) } pub fn distance_lv (v : Vector , v1 : Vector , v2 : Vector) -> f64 { let nv = v2 . sub (v1) ; let nu = v . sub (v1) ; (Self :: cross (nv , nu) / Self :: norm (nv)) . abs () } pub fn distance_sv (v : Vector , v1 : Vector , v2 : Vector) -> f64 { if Self :: dot (v2 . sub (v1) , v . sub (v1)) < 0.0 { Self :: norm (v . sub (v1)) } else if Self :: dot (v1 . sub (v2) , v . sub (v2)) < 0.0 { Self :: norm (v . sub (v2)) } else { Self :: distance_lv (v , v1 , v2) } } pub fn distance_ss (v1 : Vector , v2 : Vector , u1 : Vector , u2 : Vector) -> f64 { if Self :: intersect (v1 , v2 , u1 , u2) { 0.0 } else { let mut mi = std :: f64 :: MAX ; for & d in & [Self :: distance_sv (u1 , v1 , v2) , Self :: distance_sv (u2 , v1 , v2) , Self :: distance_sv (v1 , u1 , u2) , Self :: distance_sv (v2 , u1 , u2) ,] { if d < mi { mi = d } } mi } } pub fn intersect (v1 : Vector , v2 : Vector , u1 : Vector , u2 : Vector) -> bool { let place1 = Self :: place (u1 , v1 , v2) ; let place2 = Self :: place (u2 , v1 , v2) ; let place3 = Self :: place (v1 , u1 , u2) ; let place4 = Self :: place (v2 , u1 , u2) ; if vec ! [place1 , place2 , place3 , place4] . iter () . any (| & p | p == 7 || p == 11 || p == 17 || p == 19 || p == 23) { true } else { place1 * place2 == 3 && place3 * place4 == 3 } } pub fn place (v : Vector , v1 : Vector , v2 : Vector) -> usize { let eps = 1e-10 ; let cross = Self :: cross (v2 . sub (v1) , v . sub (v1)) ; let dot = Self :: dot (v . sub (v1) , v2 . sub (v1)) ; let bnorm = v2 . sub (v1) . norm () ; let anorm = v . sub (v1) . norm () ; if cross > 1e-10 { 1 } else if cross < - 1e-10 { 3 } else if dot < - 1e-10 { 5 } else if bnorm < 1e-10 && anorm < 1e-10 { 17 } else if bnorm < 1e-10 { 19 } else if anorm < 1e-10 { 23 } else if (anorm - bnorm) . abs () < eps { 11 } else if anorm < bnorm { 7 } else { 13 } } pub fn point_at_intersection_on_ss (v1 : Vector , v2 : Vector , u1 : Vector , u2 : Vector) -> Vector { if ! Self :: intersect (v1 , v2 , u1 , u2) { return Vector :: new (NAN , NAN) ; } Self :: point_at_intersection_on_ll (v1 , v2 , u1 , u2) } pub fn point_at_intersection_on_ll (v1 : Vector , v2 : Vector , u1 : Vector , u2 : Vector) -> Vector { let base = v2 . sub (v1) ; let h1 = (Self :: cross (base , u1 . sub (v1)) / Self :: norm (base)) . abs () ; let h2 = (Self :: cross (base , u2 . sub (v1)) / Self :: norm (base)) . abs () ; let t = h1 / (h1 + h2) ; let nv = u2 . sub (u1) . mul (t) ; u1 . add (nv) } pub fn polar_on_r (r : f64 , rad : f64) -> Vector { Vector :: new (r * rad . cos () , r * rad . sin ()) } pub fn polar_on_v (v : Vector , rad : f64) -> Vector { Vector :: new (v . x * rad . cos () - v . y * rad . sin () , v . y * rad . cos () + v . x * rad . sin () ,) } } # [derive (Clone , Copy , Debug , PartialEq)] pub struct Circle { pub c : Vector , pub r : f64 , } impl Circle { pub fn new (c : Vector , r : f64) -> Self { Self { c , r } } pub fn area (& self) -> f64 { self . r * self . r * PI } } pub struct CircleFns { } impl CircleFns { pub fn is_intersect_line (c : Circle , v1 : Vector , v2 : Vector) -> bool { let d = VectorFns :: distance_lv (c . c , v1 , v2) ; d <= c . r } pub fn is_intersect_circles (c1 : Circle , c2 : Circle) -> usize { let eps = 1e-10 ; let nv = c2 . c . sub (c1 . c) ; let d = VectorFns :: norm (nv) ; if (d - (c1 . r + c2 . r)) . abs () < eps { 3 } else if ((d + c1 . r) - c2 . r) . abs () < eps || ((d + c2 . r) - c1 . r) . abs () < eps { 1 } else if d + c1 . r < c2 . r || d + c2 . r < c1 . r { 0 } else if d > c1 . r + c2 . r { 4 } else { 2 } } pub fn points_at_intersection_segment_from_two_vectors (c : Circle , v1 : Vector , v2 : Vector ,) -> Vec < Vector > { let cp = Self :: points_at_intersection_line_from_two_vectors (c , v1 , v2) ; if cp [0] . x . is_nan () || cp [0] . y . is_nan () { return cp ; } let mut res = vec ! [] ; let p1 = VectorFns :: place (cp [0] , v1 , v2) ; let p2 = VectorFns :: place (cp [1] , v1 , v2) ; if vec ! [7 , 11 , 17 , 23] . iter () . any (| & x | x == p1) { res . push (cp [0]) ; } if vec ! [7 , 11 , 17 , 23] . iter () . any (| & x | x == p2) { res . push (cp [1]) ; } res } pub fn points_at_intersection_line_from_two_vectors (c : Circle , v1 : Vector , v2 : Vector ,) -> Vec < Vector > { if ! Self :: is_intersect_line (c , v1 , v2) { return vec ! [Vector :: new (NAN , NAN)] ; } let pr = VectorFns :: projection (c . c , v1 , v2) ; let e = VectorFns :: unit (v1 , v2) ; let nv = pr . sub (c . c) ; let base = (c . r * c . r - VectorFns :: dot (nv , nv)) . sqrt () ; let nu = e . mul (base) ; vec ! [pr . sub (nu) , pr . add (nu)] } pub fn points_at_intersection_line_from_le (c : Circle , mut le : LinearEquation) -> Vec < Vector > { le = le . normalize () . unwrap () ; let (a , b , k , x0 , y0 , r) = (le . a , le . b , le . k , c . c . x , c . c . y , c . r) ; let d = (a * x0 + b * y0 + k) . abs () ; if d > r { return vec ! [Vector :: new (NAN , NAN)] ; } let cmn = (c . r * c . r - d * d) . sqrt () ; vec ! [Vector :: new (a * d - b * cmn + x0 , b * d + a * cmn + y0) , Vector :: new (a * d + b * cmn + x0 , b * d - a * cmn + y0) ,] } fn points_at_intersection_line_from_le_2 (c : Circle , le : LinearEquation) -> Vec < Vector > { let mut res = Vec :: with_capacity (2) ; let r = le . e . dot (c . c - le . v) ; let p_mid = le . v + le . e * r ; let d = (c . r * c . r - VectorFns :: dot (p_mid - c . c , p_mid - c . c)) . sqrt () ; if d > 0. { res . push (p_mid + le . e * d) ; res . push (p_mid - le . e * d) ; } else { res . push (Vector :: new (NAN , NAN)) ; } res } pub fn tangent_point_from_polar (c : Circle , v : Vector) -> Vec < Vector > { let (x0 , y0 , r , x1 , y1) = (c . c . x , c . c . y , c . r , v . x , v . y) ; let a = x1 - x0 ; let b = y1 - y0 ; let k = x0 * x0 + y0 * y0 - x1 * x0 - y1 * y0 - r * r ; let le = LinearEquation :: new (a , b , k) ; Self :: points_at_intersection_line_from_le (c , le) } pub fn points_at_intersection_circles (c1 : Circle , c2 : Circle) -> Vec < Vector > { let nv = c2 . c . sub (c1 . c) ; let d = VectorFns :: norm (nv) ; if d > c1 . r + c2 . r { return vec ! [Vector :: new (NAN , NAN)] ; } let a = ((c1 . r * c1 . r + d * d - c2 . r * c2 . r) / (2.0 * c1 . r * d)) . acos () ; let t = nv . y . atan2 (nv . x) ; vec ! [c1 . c . add (VectorFns :: polar_on_r (c1 . r , t + a)) , c1 . c . add (VectorFns :: polar_on_r (c1 . r , t - a)) ,] } pub fn incircle (v1 : Vector , v2 : Vector , v3 : Vector) -> Circle { let nv1 = v1 . sub (v2) ; let nv2 = v2 . sub (v3) ; let nv3 = v3 . sub (v1) ; let r = VectorFns :: cross (v1 . sub (v3) , v2 . sub (v3)) . abs () / (nv1 . norm () + nv2 . norm () + nv3 . norm ()) ; let k = nv3 . norm () / (nv1 . norm () + nv3 . norm ()) ; let j = k * nv2 . norm () ; let v4 = v3 . add (nv2 . mul (j / nv2 . norm ())) ; let m = nv3 . norm () / (j + nv3 . norm ()) ; let center = v1 . add (v4 . sub (v1) . mul (m)) ; Circle :: new (center , r) } pub fn outcircle (v1 : Vector , v2 : Vector , v3 : Vector) -> Circle { let (x1 , x2 , x3 , y1 , y2 , y3) = (v1 . x , v2 . x , v3 . x , v1 . y , v2 . y , v3 . y) ; let cmn1 = x2 * x2 + y2 * y2 - x1 * x1 - y1 * y1 ; let cmn2 = x2 * x2 + y2 * y2 - x3 * x3 - y3 * y3 ; let q = (cmn1 * (x2 - x3) - cmn2 * (x2 - x1)) / (2. * ((y2 - y1) * (x2 - x3) - (y2 - y3) * (x2 - x1))) ; let p = (cmn1 * (y2 - y3) - cmn2 * (y2 - y1)) / (2. * ((y2 - y3) * (x2 - x1) - (y2 - y1) * (x2 - x3))) ; Circle :: new (Vector :: new (p , q) , v1 . sub (Vector :: new (p , q)) . norm ()) } pub fn common_line (c1 : Circle , c2 : Circle) -> LinearEquation { let a = - 2. * c1 . c . x ; let b = - 2. * c1 . c . y ; let c = c1 . c . x * c1 . c . x + c1 . c . y * c1 . c . y - c1 . r * c1 . r ; let d = - 2. * c2 . c . x ; let e = - 2. * c2 . c . y ; let f = c2 . c . x * c2 . c . x + c2 . c . y * c2 . c . y - c2 . r * c2 . r ; LinearEquation :: new (a - d , b - e , c - f) } pub fn tangent_point_from_polar_2 (c : Circle , v : Vector) -> Vec < Vector > { let center = v - c . c ; let r = (center . x * center . x + center . y * center . y - c . r * c . r) . sqrt () ; let new_c = Circle :: new (v , r) ; Self :: points_at_intersection_circles (c , new_c) } pub fn tangent_circle (c1 : Circle , c2 : Circle) -> Vec < Vector > { let (r1 , r2 , v) = (c1 . r , c2 . r , c2 . c - c1 . c) ; let d = v . norm () ; let mut res = vec ! [] ; if d > r1 + r2 { let polar = c1 . c + v * (r1 / (r1 + r2)) ; let pt = Self :: tangent_point_from_polar_2 (c1 , polar) ; res . extend (pt) ; } if d == r1 + r2 || d == (r1 - r2) . abs () { let pt = Self :: points_at_intersection_circles (c1 , c2) ; res . push (pt [0]) ; } if d > (r1 - r2) . abs () { if r1 == r2 { let n = Vector :: new (v . y / d , - v . x / d) ; res . extend (vec ! [c1 . c + n * r1 , c1 . c - n * r1]) ; } else { let polar = c1 . c + v * (r1 / (r1 - r2)) ; let pt = Self :: tangent_point_from_polar_2 (c1 , polar) ; res . extend (pt) ; } } res } pub fn area_of_two_cricles (c1 : Circle , c2 : Circle) -> f64 { let pos = Self :: is_intersect_circles (c1 , c2) ; if pos == 3 || pos == 4 { return 0. ; } else if pos == 0 || pos == 1 { if c1 . r < c2 . r { return c1 . area () ; } else { return c2 . area () ; } } let cps = Self :: points_at_intersection_circles (c1 , c2) ; let (cp1 , cp2) = (cps [0] , cps [1]) ; let place1 = VectorFns :: place (c1 . c , cp1 , cp2) ; let place2 = VectorFns :: place (c2 . c , cp1 , cp2) ; let cross = VectorFns :: cross ; let dot = VectorFns :: dot ; let s1 = cross (cp2 - c1 . c , cp1 - c1 . c) . abs () * 0.5 ; let ang1 = cross (cp1 - c1 . c , cp2 - c1 . c) . atan2 (dot (cp1 - c1 . c , cp2 - c1 . c)) . abs () ; let cs1 = c1 . r * c1 . r * ang1 * 0.5 ; let s2 = cross (cp2 - c2 . c , cp1 - c2 . c) . abs () * 0.5 ; let ang2 = cross (cp1 - c2 . c , cp2 - c2 . c) . atan2 (dot (cp1 - c2 . c , cp2 - c2 . c)) . abs () ; let cs2 = c2 . r * c2 . r * ang2 * 0.5 ; if place1 * place2 == 3 { cs1 + cs2 - s1 - s2 } else if c1 . r < c2 . r { let cs1 = c1 . r * c1 . r * (2. * PI - ang1) * 0.5 ; cs1 + cs2 + s1 - s2 } else { let cs2 = c2 . r * c2 . r * (2. * PI - ang2) * 0.5 ; cs1 + cs2 - s1 + s2 } } pub fn area_of_circle_polygon (c : Circle , p : Polygon) -> f64 { # [derive (Clone , Copy , Debug , PartialEq , PartialOrd , Eq , Ord)] enum TYPE { Cross , Edge , } let mut ps = vec ! [] ; let mut area = 0.0 ; let n = p . len () ; for i in 0 .. n { let ni = (i + 1) % n ; Self :: points_at_intersection_segment_from_two_vectors (c , p [i] , p [ni]) . iter () . filter (| cp | ! cp . x . is_nan () && ! cp . y . is_nan ()) . for_each (| & x | ps . push ((x , TYPE :: Cross , true))) ; let is_in = p [ni] . sub (c . c) . norm () <= c . r ; ps . push ((p [ni] , TYPE :: Edge , is_in)) ; } let cross = VectorFns :: cross ; let dot = VectorFns :: dot ; let m = ps . len () ; for i in 0 .. m { let ni = (i + 1) % m ; let (p , _ , is_in) = ps [i] ; let (np , _ , nis_in) = ps [ni] ; if is_in && nis_in { area += cross (p - c . c , np - c . c) * 0.5 ; } else { let ang = cross (p - c . c , np - c . c) . atan2 (dot (p - c . c , np - c . c)) ; area += c . r * c . r * ang * 0.5 ; } } area } } pub type Polygon = Vec < Vector > ; pub struct PolygonFns { } impl PolygonFns { pub fn area (p : Polygon) -> f64 { let n = p . len () ; let mut area = 0.0 ; for i in 0 .. n { let ni = (i + 1) % n ; area += VectorFns :: cross (p [i] , p [ni]) / 2.0 ; } area } pub fn contain_point_opt (p : Polygon , v : Vector) -> usize { let eps = 1e-10 ; let n = p . len () ; let mut cn = 0 ; for i in 0 .. n { let ni = (i + 1) % n ; if (p [i] . y <= v . y && v . y < p [ni] . y) || (p [ni] . y <= v . y && v . y < p [i] . y) { let rat = (v . y - p [i] . y) / (p [ni] . y - p [i] . y) ; let qx = p [i] . x + rat * (p [ni] . x - p [i] . x) ; if (qx - v . x) . abs () > eps && v . x < qx { cn += 1 ; } } } if cn % 2 == 1 { 1 } else { 3 } } pub fn contain_point (p : Polygon , v1 : Vector) -> usize { let eps = 1e-10 ; let n = p . len () ; let mut cn = 0 ; for i in 0 .. n { let mut v2 = p [i] ; let mut v = p [(i + 1) % n] ; if [5 , 17 , 19 , 23] . contains (& VectorFns :: place (v , v1 , v2)) { return 5 ; } if v2 . y > v . y { std :: mem :: swap (& mut v , & mut v2) ; } if VectorFns :: place (v , v1 , v2) == 1 && (v2 . sub (v1) . y < eps && eps < v . sub (v1) . y) { cn += 1 ; } } if cn % 2 == 1 { 1 } else { 3 } } pub fn contain_circle (p : Polygon , c : Circle) -> bool { let n = p . len () ; let mut ans = Self :: contain_point (p . clone () , c . c) == 1 ; if ! ans { return ans ; } for i in 0 .. n { let ni = (i + 1) % n ; let d = VectorFns :: distance_lv (c . c , p [i] , p [ni]) ; ans &= d >= c . r ; } ans } pub fn convex_hull (mut p : Polygon) -> (Vec < Vector > , Vec < Vector >) { p . sort () ; let n = p . len () ; let mut up = vec ! [p [0] , p [1]] ; let mut low = vec ! [p [n - 1] , p [n - 2]] ; for & v in p [2 .. n] . iter () { let mut k = up . len () ; while k >= 2 && VectorFns :: place (v , up [k - 2] , up [k - 1]) == 1 { up . pop () ; k -= 1 ; } up . push (v) ; } for & v in p [0 .. n - 2] . iter () . rev () { let mut k = low . len () ; while k >= 2 && VectorFns :: place (v , low [k - 2] , low [k - 1]) == 1 { low . pop () ; k -= 1 ; } low . push (v) ; } (up , low) } pub fn diameter (p : Polygon) -> f64 { let (up , mut low) = Self :: convex_hull (p) ; let mut vs = up ; let wl = low . len () ; vs . append (& mut low [1 .. wl - 1] . to_vec ()) ; let n = vs . len () ; let mut ma = 0. ; let (mut pos1 , mut pos2) = (0 , 1) ; let (mut ma_pos1 , mut ma_pos2) = (0 , 1) ; loop { let d = vs [pos2] . sub (vs [pos1]) . norm () ; if ma < d { ma = d ; ma_pos1 = pos1 ; ma_pos2 = pos2 ; } ; let s1 = vs [(pos1 + 1) % n] . sub (vs [pos1]) ; let s2 = vs [(pos2 + 1) % n] . sub (vs [pos2]) ; if VectorFns :: cross (s1 , s2) > 0. { pos1 = (pos1 + 1) % n ; } else { pos2 = (pos2 + 1) % n ; } if pos1 == ma_pos1 && pos2 == ma_pos2 { break ; } } ma } pub fn convex_cut (p : Polygon , v1 : Vector , v2 : Vector) -> f64 { let nv1 = v1 . add (v1 . sub (v2) . mul (10000.)) ; let nv2 = v2 . add (v2 . sub (v1) . mul (10000.)) ; let n = p . len () ; let mut cc = 0 ; let mut vs = vec ! [p [0]] ; let mut edge = 0 ; for i in 0 .. n { let cur = p [i] ; let next = p [(i + 1) % n] ; let dist1 = VectorFns :: distance_lv (cur , nv1 , nv2) ; let dist2 = VectorFns :: distance_lv (next , nv1 , nv2) ; if dist1 == 0. && dist2 == 0. { let cos = VectorFns :: dot (next . sub (cur) , nv2 . sub (nv1)) ; if cos < 0. { return 0. ; } else { return Self :: area (p) ; } } if VectorFns :: intersect (nv1 , nv2 , cur , next) { let cross_point = VectorFns :: point_at_intersection_on_ll (nv1 , nv2 , cur , next) ; if cross_point != p [i] { vs . push (cross_point) ; cc += 1 ; } if cross_point == p [i] && cc % 2 == 1 { edge = i ; } } if cc % 2 == 0 { vs . push (p [(i + 1) % n]) ; } } if cc % 2 == 1 { let cur = p [edge] ; let next = p [(edge + 1) % n] ; let sin = VectorFns :: cross (nv2 . sub (nv1) , next . sub (cur)) ; if sin < 0. { return 0. ; } else { return Self :: area (p) ; } } let place = VectorFns :: place (p [0] , nv1 , nv2) ; let area = Self :: area (p) ; let part = Self :: area (vs) ; if place == 1 { part } else { area - part } } } pub mod manhattan_geo { use crate :: __cargo_equip :: crates :: collection :: geo_lib :: Vector ; use std :: cmp :: { Ordering , Ordering :: { Equal , Greater , Less } , } ; use std :: collections :: { BTreeMap , BTreeSet } ; type Set < T > = BTreeSet < T > ; type Point = Vector ; # [derive (Clone , Copy , Debug)] struct Num (f64) ; impl Num { pub fn new (num : f64) -> Self { Self (num) } } impl Eq for Num { } impl PartialEq for Num { fn eq (& self , other : & Self) -> bool { let eps = 1e-10 ; (self . 0 - other . 0) . abs () < eps } } impl Ord for Num { fn cmp (& self , other : & Self) -> Ordering { if self == other { Equal } else if self . 0 < other . 0 { Less } else { Greater } } } impl PartialOrd for Num { fn partial_cmp (& self , other : & Self) -> Option < Ordering > { Some (self . cmp (other)) } } # [derive (Clone , Copy , Debug , PartialEq , PartialOrd , Eq , Ord)] enum TYPE { BOTTOM , LEFT , RIGHT , TOP , } # [derive (Clone , Copy , Debug , PartialEq)] struct EndPoint { p : Point , pos : usize , t : TYPE , } impl EndPoint { pub fn new (p : Point , pos : usize , t : TYPE) -> Self { Self { p , pos , t } } } pub fn plane_sweep (mut segs : Vec < (Point , Point) >) -> usize { use crate :: __cargo_equip :: crates :: collection :: geo_lib :: manhattan_geo :: TYPE :: * ; let mut eps = vec ! [] ; for (i , (p1 , p2)) in segs . iter_mut () . enumerate () { if p1 > p2 { std :: mem :: swap (p1 , p2) ; } if Num :: new (p1 . y) == Num :: new (p2 . y) { eps . push (EndPoint :: new (* p1 , i , LEFT)) ; eps . push (EndPoint :: new (* p2 , i , RIGHT)) ; } else { eps . push (EndPoint :: new (* p1 , i , BOTTOM)) ; eps . push (EndPoint :: new (* p2 , i , TOP)) ; } } eps . sort_by (| p1 , p2 | { let y1 = Num :: new (p1 . p . y) ; let y2 = Num :: new (p2 . p . y) ; if y1 == y2 { p1 . t . cmp (& p2 . t) } else { y1 . cmp (& y2) } }) ; let mut bt = Set :: new () ; let mut cnt = 0 ; for ep in eps { let (t , pos , p) = (ep . t , ep . pos , ep . p) ; if t == TOP { bt . remove (& Num :: new (p . x)) ; } else if t == BOTTOM { bt . insert (Num :: new (p . x)) ; } else if t == LEFT { let begin = segs [pos] . 0 . x ; let end = segs [pos] . 1 . x ; let a = bt . range (Num :: new (begin) ..= Num :: new (end)) ; cnt += a . count () ; } } cnt } } # [derive (Clone , Copy , Debug , PartialEq)] pub struct LinearEquation { a : f64 , b : f64 , k : f64 , v : Vector , e : Vector , } impl LinearEquation { pub fn new (a : f64 , b : f64 , k : f64) -> Self { let a_point_on_line = Vector :: new (b , - a) ; let e = a_point_on_line / a_point_on_line . norm () ; let t = if a != 0. { (- k / a , 0.) } else if b != 0. { (0. , - k / b) } else { panic ! ("no line") } ; let v = Vector :: from (t) ; Self { a , b , k , v , e } } pub fn one_y (& self) -> Option < Self > { let eps = 1e-10 ; if self . b . abs () < eps { return None ; } Some (Self :: new (self . a / self . b , 1. , self . k / self . b)) } pub fn one_x (& self) -> Option < Self > { let eps = 1e-10 ; if self . a . abs () < eps { return None ; } Some (Self :: new (1. , self . b / self . a , self . k / self . a)) } pub fn normalize (& self) -> Option < Self > { let eps = 1e-10 ; let n = (self . a * self . a + self . b * self . b) . sqrt () ; if n < eps { return None ; } Some (Self :: new (self . a / n , self . b / n , self . k / n)) } } pub struct LinearEquationFns { } impl LinearEquationFns { pub fn slope (v1 : Vector , v2 : Vector) -> f64 { let eps = 1e-10 ; let dx = v2 . x - v1 . x ; if dx . abs () < eps { 0. } else { (v2 . y - v1 . y) / dx } } pub fn solve (le1 : LinearEquation , le2 : LinearEquation) -> Vector { let eps = 1e-10 ; if (le1 . b . abs () < eps && le2 . b . abs () < eps) || (le1 . a . abs () < eps && le2 . a . abs () < eps) { Vector :: new (NAN , NAN) } else if le1 . b . abs () < eps && le2 . a . abs () < eps { Vector :: new (- le1 . k / le1 . a , - le2 . k / le2 . b) } else if le1 . a . abs () < eps && le2 . b . abs () < eps { Vector :: new (- le2 . k / le2 . a , - le1 . k / le1 . b) } else if le1 . a . abs () < eps { let y = - le1 . k / le1 . b ; let nle = le2 . one_x () . unwrap () ; let x = - (nle . k + nle . b * y) ; Vector :: new (x , y) } else if le1 . b . abs () < eps { let x = - le1 . k / le1 . a ; let nle = le2 . one_y () . unwrap () ; let y = - (nle . k + nle . a * x) ; Vector :: new (x , y) } else if le2 . a . abs () < eps { let y = - le2 . k / le2 . b ; let nle = le1 . one_x () . unwrap () ; let x = - (nle . k + nle . b * y) ; Vector :: new (x , y) } else if le2 . b . abs () < eps { let x = - le2 . k / le2 . a ; let nle = le1 . one_y () . unwrap () ; let y = - (nle . k + nle . a * x) ; Vector :: new (x , y) } else { let nle1 = le1 . one_y () . unwrap () ; let nle2 = le2 . one_y () . unwrap () ; let x = - (nle1 . k - nle2 . k) / (nle1 . a - nle2 . a) ; let y = - (nle1 . a * x + nle1 . k) ; Vector :: new (x , y) } } } } pub mod nt_lib { use std :: collections :: { BTreeMap , BTreeSet } ; pub fn prime (ma : usize) -> Vec < bool > { let mut ps = vec ! [2] ; for i in (3 ..= ma) . step_by (2) { let mut is_prime = true ; for p in & ps { if p * p > i { break ; } if i % p == 0 { is_prime = false ; break ; } } if is_prime { ps . push (i) ; } } let mut table = vec ! [false ; ma + 1] ; for x in ps { table [x] = true ; } table } pub fn factorize (n : usize) -> BTreeMap < usize , usize > { let mut factors = BTreeMap :: < usize , usize > :: new () ; let mut rest = n ; for i in 2 .. n { if i > rest { break ; } if i * i > n && factors . is_empty () { factors . insert (n , 1) ; break ; } while rest != 0 && rest % i == 0 { rest /= i ; if let Some (x) = factors . get_mut (& i) { * x += 1 ; } else { factors . insert (i , 1) ; } } } factors } pub fn euler_phi (n : usize) -> usize { let mut ans = n as f64 ; let factors = factorize (n) ; for (k , _) in factors . iter () { ans *= 1. - 1. / * k as f64 ; } ans as usize } pub fn gcd (xs : Vec < usize >) -> usize { if xs . is_empty () { return 0 ; } let mut ret = xs [0] ; for mut a in xs { let mut b = ret ; while a % b != 0 { let tmp = b ; b = a % b ; a = tmp ; } ret = b ; } ret } pub fn lcm (xs : Vec < usize >) -> usize { if xs . is_empty () { return 0 ; } let mut a = xs [0] ; for b in xs { a = a * b / gcd (vec ! [a , b]) ; } a } pub fn ext_gcd (a : isize , b : isize) -> (isize , isize) { if b == 0 { return (0 , 1) ; } let (x , y) = ext_gcd (b , a % b) ; (y - a / b * x , x) } } pub mod utils { use std :: io ; pub fn read < T : std :: str :: FromStr > () -> Vec < T > { let mut buf = String :: new () ; io :: stdin () . read_line (& mut buf) . unwrap () ; buf . trim () . split (' ') . flat_map (str :: parse) . collect () } }}
    }

    pub(crate) mod macros {
        pub mod collection {}
    }

    pub(crate) mod prelude {pub use crate::__cargo_equip::crates::*;}

    mod preludes {
        pub mod collection {}
    }
}