Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Convolution を実装 #53

Merged
merged 10 commits into from
Sep 16, 2020
Merged

Convolution を実装 #53

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
3b095ab
Add BSF
terry-u16 Sep 15, 2020
134f5fc
Add PowMod
terry-u16 Sep 15, 2020
2fcc4b9
Add comment
terry-u16 Sep 15, 2020
ff83f8b
Add PrimitiveRoot
terry-u16 Sep 15, 2020
9b58e52
Add Butterfly
terry-u16 Sep 15, 2020
5e66626
Implement Convolution
terry-u16 Sep 15, 2020
e3fae8b
Move BSF from InternalMath to InternalBit
terry-u16 Sep 16, 2020
1edcdd5
Add cache
terry-u16 Sep 16, 2020
247a084
Make cache simple on Butterfly
terry-u16 Sep 16, 2020
b477cce
Remove InternalMath.PowMod
terry-u16 Sep 16, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 31 additions & 1 deletion AtCoderLibrary/Internal/InternalBit.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,6 @@
using System.Runtime.CompilerServices;
using System.Diagnostics;
using System.Numerics;
using System.Runtime.CompilerServices;
using System.Runtime.Intrinsics.X86;

namespace AtCoder.Internal
Expand All @@ -20,5 +22,33 @@ public static int ExtractLowestSetBit(int n)
}
return n & -n;
}

/// <summary>
/// (<paramref name="n"/> &amp; (1 &lt;&lt; x)) != 0 なる最小の非負整数 x を求めます。
/// </summary>
/// <remarks>
/// <para>BSF: Bit Scan Forward</para>
/// <para>制約: 1 ≤ <paramref name="n"/></para>
/// </remarks>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static int BSF(uint n)
{
Debug.Assert(n >= 1);
if (Bmi1.IsSupported)
{
// O(1)
return (int)Bmi1.TrailingZeroCount(n);
}
else if (Popcnt.IsSupported)
{
// O(1)
return (int)Popcnt.PopCount(~n & (n - 1));
}
else
{
// O(logn)
return BitOperations.TrailingZeroCount(n);
}
}
}
}
146 changes: 146 additions & 0 deletions AtCoderLibrary/Internal/Math/Butterfly.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,146 @@
using System;

namespace AtCoder.Internal
{
public static class Butterfly<T> where T : struct, IStaticMod
{
/// <summary>
/// sumE[i] = ies[0] * ... * ies[i - 1] * es[i]
/// </summary>
private static StaticModInt<T>[] sumE = CalcurateSumE();

/// <summary>
/// sumIE[i] = es[0] * ... * es[i - 1] * ies[i]
/// </summary>
private static StaticModInt<T>[] sumIE = CalcurateSumIE();

public static void Calculate(Span<StaticModInt<T>> a)
{
var n = a.Length;
var h = InternalMath.CeilPow2(n);

for (int ph = 1; ph <= h; ph++)
{
// ブロックサイズの半分
int w = 1 << (ph - 1);

// ブロック数
int p = 1 << (h - ph);

var now = StaticModInt<T>.Raw(1);

// 各ブロックの s 段目
for (int s = 0; s < w; s++)
{
int offset = s << (h - ph + 1);

for (int i = 0; i < p; i++)
{
var l = a[i + offset];
var r = a[i + offset + p] * now;
a[i + offset] = l + r;
a[i + offset + p] = l - r;
}
now *= sumE[InternalBit.BSF(~(uint)s)];
}
}
}

public static void CalculateInv(Span<StaticModInt<T>> a)
{
var n = a.Length;
var h = InternalMath.CeilPow2(n);

for (int ph = h; ph >= 1; ph--)
{
// ブロックサイズの半分
int w = 1 << (ph - 1);

// ブロック数
int p = 1 << (h - ph);

var iNow = StaticModInt<T>.Raw(1);

// 各ブロックの s 段目
for (int s = 0; s < w; s++)
{
int offset = s << (h - ph + 1);

for (int i = 0; i < p; i++)
{
var l = a[i + offset];
var r = a[i + offset + p];
a[i + offset] = l + r;
a[i + offset + p] = StaticModInt<T>.Raw(
unchecked((int)((ulong)(default(T).Mod + l.Value - r.Value) * (ulong)iNow.Value % default(T).Mod)));
}
iNow *= sumIE[InternalBit.BSF(~(uint)s)];
}
}
}

private static StaticModInt<T>[] CalcurateSumE()
{
int g = InternalMath.PrimitiveRoot((int)default(T).Mod);
int cnt2 = InternalBit.BSF(default(T).Mod - 1);
var e = new StaticModInt<T>(g).Pow((default(T).Mod - 1) >> cnt2);
var ie = e.Inv();

var sumE = new StaticModInt<T>[cnt2 - 2];

// es[i]^(2^(2+i)) == 1
Span<StaticModInt<T>> es = stackalloc StaticModInt<T>[cnt2 - 1];
Span<StaticModInt<T>> ies = stackalloc StaticModInt<T>[cnt2 - 1];

for (int i = es.Length - 1; i >= 0; i--)
{
// e^(2^(2+i)) == 1
es[i] = e;
ies[i] = ie;
e *= e;
ie *= ie;
}

var now = StaticModInt<T>.Raw(1);
for (int i = 0; i < sumE.Length; i++)
{
sumE[i] = es[i] * now;
now *= ies[i];
}

return sumE;
}

private static StaticModInt<T>[] CalcurateSumIE()
{
int g = InternalMath.PrimitiveRoot((int)default(T).Mod);
int cnt2 = InternalBit.BSF(default(T).Mod - 1);
var e = new StaticModInt<T>(g).Pow((default(T).Mod - 1) >> cnt2);
var ie = e.Inv();

var sumIE = new StaticModInt<T>[cnt2 - 2];

// es[i]^(2^(2+i)) == 1
Span<StaticModInt<T>> es = stackalloc StaticModInt<T>[cnt2 - 1];
Span<StaticModInt<T>> ies = stackalloc StaticModInt<T>[cnt2 - 1];

for (int i = es.Length - 1; i >= 0; i--)
{
// e^(2^(2+i)) == 1
es[i] = e;
ies[i] = ie;
e *= e;
ie *= ie;
}

var now = StaticModInt<T>.Raw(1);
for (int i = 0; i < sumIE.Length; i++)
{
sumIE[i] = ies[i] * now;
now *= es[i];
}

return sumIE;
}
}
}
84 changes: 84 additions & 0 deletions AtCoderLibrary/Internal/Math/PrimitiveRoot.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,84 @@
using System;
using System.Collections.Generic;
using System.Diagnostics;

namespace AtCoder.Internal
{
public static partial class InternalMath
{
private static readonly Dictionary<int, int> primitiveRootsCache = new Dictionary<int, int>()
{
{ 2, 1 },
{ 167772161, 3 },
{ 469762049, 3 },
{ 754974721, 11 },
{ 998244353, 3 }
};

/// <summary>
/// <paramref name="m"/> の最小の原始根を求めます。
/// </summary>
/// <remarks>
/// 制約: <paramref name="m"/> は素数
/// </remarks>
public static int PrimitiveRoot(int m)
{
Debug.Assert(m >= 2);

if (primitiveRootsCache.TryGetValue(m, out var p))
{
return p;
}

return primitiveRootsCache[m] = Calculate(m);

int Calculate(int m)
{
Span<int> divs = stackalloc int[20];
divs[0] = 2;
int cnt = 1;
int x = (m - 1) / 2;

while (x % 2 == 0)
{
x >>= 1;
}

for (int i = 3; (long)i * i <= x; i += 2)
{
if (x % i == 0)
{
divs[cnt++] = i;
while (x % i == 0)
{
x /= i;
}
}
}

if (x > 1)
{
divs[cnt++] = x;
}

for (int g = 2; ; g++)
{
bool ok = true;
for (int i = 0; i < cnt; i++)
{
if (Math.PowMod(g, (m - 1) / divs[i], m) == 1)
{
ok = false;
break;
}
}

if (ok)
{
return g;
}
}
}
}
}
}
Loading