Closed3
AtCoder ARC148 F - 998244353 → 1000000007 回答例
のような演算パターンの組み合わせで
以下は自分の19行の回答例。(既にユーザ解説で18行での解法も紹介されているのでそちらも参照)
1回目のモンゴメリ変換を
2回目のモンゴメリ変換を
と2回目のモンゴメリ変換を変形して適用した場合に短縮できるかの試み。(出来てはいません)
/*
19
mul T A B
rem U T
mul U U 4915446
rem U U
mul U U 1000000007
add U U T
mul T U 8456992263029997534
rem U T
mul U U 993328907
rem U U
mul U U 18446744072709551609
add U U T
mul T U 996491785301655553
mul U T 998244353
rem U U
mul V U 696235320
rem V V
add U U V
add C T U
*/
const N: u64 = 1000000007;
const R: u64 = 998244353;
const NEG_N: u64 = 18446744072709551609; // 2 ** 64 - n;
const N_DASH: u64 = 4915446; // n * ndash % r == r - 1
const N_INV: u64 = 993328907; // n * ninv % r == 1
const R_INV: u64 = 996491785301655553; // r * rinv % (2**64) == 1
const R2: u64 = 320946142; // == r * r % n
//const RC: u64 = 8456992263029997534; // R_INV * R2 % (2**64)
const Q: u64 = 932051910; // 2**64 % R
const P: u64 = 696235320; // Q * P % R + Q == N
fn reduction1(t: u64) -> u64 {
// t := (T + (((T % R) * N_DASH % R) * N)) * R_INV (0 <= t < T/R + N)
let u = t % R;
let u = u.wrapping_mul(N_DASH);
let u = u % R;
let u = u.wrapping_mul(N);
let u = u.wrapping_add(t);
let _t = u.wrapping_mul(R_INV);
assert_eq!(_t * R % N, t % N);
_t
}
fn reduction2(t: u64) -> u64 {
// t := (T + (((T % R) * N_INV % R) * NEG_N)) * R_INV (-N < t < T/R)
let u = t % R;
let u = u.wrapping_mul(N_INV);
let u = u % R;
let u = u.wrapping_mul(NEG_N);
let u = u.wrapping_add(t);
let _t = u.wrapping_mul(R_INV);
assert_eq!(_t.wrapping_add(N) * R % N, t % N);
_t
}
fn fix(t: u64) -> u64 { // if t < 0 then return t + N else return t + 0
let u = t.wrapping_mul(R);
let u = u % R;
assert!(((t as i64) < 0 && u == Q) || ((t as i64) >= 0 && u == 0));
let v = u.wrapping_mul(P);
let v = v % R;
let u = u.wrapping_add(v);
assert!(((t as i64) < 0 && u == N) || ((t as i64) >= 0 && u == 0));
let _c = t.wrapping_add(u);
_c
}
pub fn mulmod(a: u64, b: u64) -> u64 {
let t = reduction1(a * b);
let t = reduction2(t * R2);
let c = fix(t);
assert_eq!(a * b % N, c);
c
}
18行に短縮できました。reduction2 の最後の
/*
18
mul T A B
rem U T
mul U U 4915446
rem U U
mul U U 1000000007
add U U T
mul T U 8456992263029997534
rem U T
mul U U 993328907
rem U U
mul U U 18446744072709551609
add U U T
mul T U 996491785301655553
rem U U
mul V U 696235320
rem V V
add U U V
add C T U
*/
const N: u64 = 1000000007;
const R: u64 = 998244353;
const NEG_N: u64 = 18446744072709551609; // 2 ** 64 - n;
const N_DASH: u64 = 4915446; // n * ndash % r == r - 1
const N_INV: u64 = 993328907; // n * ninv % r == 1
const R_INV: u64 = 996491785301655553; // r * rinv % (2**64) == 1
const R2: u64 = 320946142; // == r * r % n
//const RC: u64 = 8456992263029997534; // R_INV * R2 % (2**64)
const Q: u64 = 932051910; // 2**64 % R
const P: u64 = 696235320; // Q * P % R + Q == N
fn reduction1(t: u64) -> u64 {
// t := (T + (((T % R) * N_DASH % R) * N)) * R_INV (0 <= t < T/R + N)
let u = t % R;
let u = u.wrapping_mul(N_DASH);
let u = u % R;
let u = u.wrapping_mul(N);
let u = u.wrapping_add(t);
let _t = u.wrapping_mul(R_INV);
assert_eq!(_t * R % N, t % N);
_t
}
fn reduction2(t: u64) -> (u64, u64) { // t の値域が異なるモンゴメリ剰余乗算の変種を使う
// t := (T + (((T % R) * N_INV % R) * NEG_N)) * R_INV (-N < t < T/R)
let u = t % R;
let u = u.wrapping_mul(N_INV);
let u = u % R;
let u = u.wrapping_mul(NEG_N);
let u = u.wrapping_add(t);
let _t = u.wrapping_mul(R_INV);
assert_eq!(_t.wrapping_add(N) * R % N, t % N);
(_t, u)
}
fn fix(t: u64, u: u64) -> u64 { // if t < 0 then return t + N else return t + 0
//let u = t.wrapping_mul(R); // reduction2 の最後で /R しているのでその前から値を貰う
let u = u % R; // t が負なら Q=2**64%R 、 t が非負なら 0 の値を u に代入
assert!(((t as i64) < 0 && u == Q) || ((t as i64) >= 0 && u == 0));
let v = u.wrapping_mul(P);
let v = v % R;
let u = u.wrapping_add(v); // t が負なら N 、 t が非負なら 0 の値を u に代入
assert!(((t as i64) < 0 && u == N) || ((t as i64) >= 0 && u == 0));
let _c = t.wrapping_add(u); // t が負なら t+N 、 t が非負なら t の値を _c に代入
_c
}
pub fn mulmod(a: u64, b: u64) -> u64 {
let t = reduction1(a * b);
let (t, u) = reduction2(t * R2);
let c = fix(t, u);
assert_eq!(a * b % N, c);
c
}
このスクラップの内容は上記記事にまとめました。
このスクラップは2022/09/12にクローズされました
ログインするとコメントできます