Skip to content

Commit

Permalink
v3.0.0-atcoder3
Browse files Browse the repository at this point in the history
  • Loading branch information
kzrnm committed Jan 25, 2023
2 parents 1b2d51f + c4d21b3 commit 79218e2
Show file tree
Hide file tree
Showing 7 changed files with 197 additions and 85 deletions.
3 changes: 3 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,9 @@ All notable changes to this project will be documented in this file.
The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).

## [3.0.0-pre3] - 2023-01-26
- Multiply 128Bit in barrett reduction

## [3.0.0-pre2] - 2023-01-25
- AtCoderAnalyzer: Create operator for static abstract
- Rename DynamicModID to DynamicModIntId
Expand Down
4 changes: 2 additions & 2 deletions Directory.Build.props
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,8 @@
<RepositoryUrl>https://github.com/kzrnm/ac-library-csharp</RepositoryUrl>
<PackageReleaseNotes>https://github.com/kzrnm/ac-library-csharp/blob/main/CHANGELOG.md</PackageReleaseNotes>

<Version>3.0.0-atcoder2</Version>
<AssemblyVersion>3.0.0.2</AssemblyVersion>
<Version>3.0.0-atcoder3</Version>
<AssemblyVersion>3.0.0.3</AssemblyVersion>
<RepositoryCommit Condition="'$(GIT_COMMIT)' != ''">$(GIT_COMMIT)</RepositoryCommit>

<SignAssembly>True</SignAssembly>
Expand Down
26 changes: 9 additions & 17 deletions Source/ac-library-csharp/Math/Internal/Barrett.cs
Original file line number Diff line number Diff line change
@@ -1,8 +1,4 @@
#if NETCOREAPP3_0_OR_GREATER
using System.Runtime.Intrinsics.X86;
#endif
using System.Runtime.CompilerServices;

using System.Runtime.CompilerServices;

namespace AtCoder.Internal
{
Expand All @@ -24,19 +20,15 @@ public Barrett(uint m)
/// <paramref name="a"/> * <paramref name="b"/> mod m
/// </summary>
[MethodImpl(256)]
public uint Mul(uint a, uint b)
public uint Mul(uint a, uint b) => Reduce((ulong)a * b);

[MethodImpl(256)]
public uint Reduce(ulong z)
{
var z = (ulong)a * b;
#if NETCOREAPP3_0_OR_GREATER
if (Bmi2.X64.IsSupported)
{
var x = Bmi2.X64.MultiplyNoFlags(z, IM);
var v = unchecked((uint)(z - x * Mod));
if (Mod <= v) v += Mod;
return v;
}
#endif
return (uint)(z % Mod);
var x = InternalMath.Mul128Bit(z, IM);
var v = unchecked((uint)(z - x * Mod));
if (Mod <= v) v += Mod;
return v;
}
}
}
85 changes: 82 additions & 3 deletions Source/ac-library-csharp/Math/Internal/InternalMath.cs
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
using System.Collections.Generic;
using System.Runtime.CompilerServices;
#if NETCOREAPP3_1_OR_GREATER
using System.Runtime.Intrinsics.X86;
using System.Numerics;
#endif

Expand Down Expand Up @@ -159,14 +160,14 @@ public static bool IsPrime(int n)
return true;
}
/// <summary>
/// Same as <see cref="MathLib.PowMod(long, long, int)"/>
/// <inheritdoc cref="MathLib.PowMod(long, long, int)" />
/// </summary>
[MethodImpl(256)]
private static long PowMod(long x, long n, int m)
internal static long PowMod(long x, long n, int m)
{
Contract.Assert(0 <= n && 1 <= m, reason: $"0 <= {nameof(n)} && 1 <= {nameof(m)}");
if (m == 1) return 0;
Barrett barrett = new Barrett((uint)m);
var barrett = new Barrett((uint)m);
uint r = 1, y = (uint)SafeMod(x, m);
while (0 < n)
{
Expand Down Expand Up @@ -199,5 +200,83 @@ public static ulong FloorSumUnsigned(ulong n, ulong m, ulong a, ulong b)
(n, m, a, b) = (yMax / m, a, m, yMax % m);
}
}

/// <summary>
/// <paramref name="a"/> * <paramref name="b"/> の上位 64 ビットを返します。
/// </summary>
/// <param name="a"></param>
/// <param name="b"></param>
/// <returns></returns>
[MethodImpl(256)]
internal static ulong Mul128Bit(ulong a, ulong b)
{
#if NETCOREAPP3_1_OR_GREATER
if (Bmi2.X64.IsSupported)
return Bmi2.X64.MultiplyNoFlags(a, b);
#endif
return Mul128BitLogic(a, b);
}

/// <summary>
/// <paramref name="a"/> * <paramref name="b"/> の上位 64 ビットを返します。
/// </summary>
/// <param name="a"></param>
/// <param name="b"></param>
/// <returns></returns>
[MethodImpl(256)]
internal static ulong Mul128BitLogic(ulong a, ulong b)
{
/*
* ここでは 64 bit 整数 X の上位, 下位 32 bit をそれぞれ Xu ,Xd と表す
* A * B を計算する。
*
* 変数を下記のように定義する。
* s := 1ul << 32
* R := A * B
* L := Ad * Bd
* M1 := Ad * Bu
* M2 := Au * Bd
* H := Au * Bu
*
* A * B = s * s * H + s * (M1 + M2) + L
* であることを使って
* R = s * s * x + y
* と表せる x, y を求めたい (x, y < s * s)
*
* A * B
* = s * s * H + s * (s * M1u + M1d + s * M2u + M2d) + L
* = s * s * H + s * (s * (M1u + M2u) + M1d + M2d) + L
* = s * s * (H + M1u + M2u) + s * (M1d + M2d) + L
*
* と表せるので、繰り上がりを考慮しなければ
* x = H + M1u + M2u
* y = s * (M1d + M2d) + L = s * (M1d + M2d + Lu) + Ld
* となる
*
* 繰り上がるのは
* M1d + M2d + Lu の上位 32 bit なので
* C := M1d + M2d + Lu
* とすると
*
* x = H + M1u + M2u + Cu
* となる
*/
var au = a >> 32;
var ad = a & 0xFFFFFFFF;
var bu = b >> 32;
var bd = b & 0xFFFFFFFF;

var l = ad * bd;
var m1 = au * bd;
var m2 = ad * bu;
var h = au * bu;

var lu = l >> 32;
var m1d = m1 & 0xFFFFFFFF;
var m2d = m2 & 0xFFFFFFFF;
var c = m1d + m2d + lu;

return h + (m1 >> 32) + (m2 >> 32) + (c >> 32);
}
}
}
14 changes: 1 addition & 13 deletions Source/ac-library-csharp/Math/MathLib.PowMod.cs
Original file line number Diff line number Diff line change
Expand Up @@ -14,18 +14,6 @@ public static partial class MathLib
/// </remarks>
[MethodImpl(256)]
public static long PowMod(long x, long n, int m)
{
Contract.Assert(0 <= n && 1 <= m, reason: $"0 <= {nameof(n)} && 1 <= {nameof(m)}");
if (m == 1) return 0;
Barrett barrett = new Barrett((uint)m);
uint r = 1, y = (uint)InternalMath.SafeMod(x, m);
while (0 < n)
{
if ((n & 1) != 0) r = barrett.Mul(r, y);
y = barrett.Mul(y, y);
n >>= 1;
}
return r;
}
=> InternalMath.PowMod(x, n, m);
}
}
55 changes: 55 additions & 0 deletions Test/ac-library-csharp.Test/Internal/BarrettTest.cs
Original file line number Diff line number Diff line change
@@ -0,0 +1,55 @@
using System.Collections.Generic;
using FluentAssertions;
using Xunit;

namespace AtCoder.Internal
{
public class BarrettTest
{
[Fact]
public void Barrett()
{
for (uint m = 1; m <= 100; m++)
{
var bt = new Barrett(m);
for (uint a = 0; a < m; a++)
{
for (uint b = 0; b < m; b++)
{
bt.Mul(a, b).Should().Be((a * b) % m);
}
}
}

new Barrett(1).Mul(0, 0).Should().Be(0);
}

[Fact]
public void BarrettBorder()
{
for (uint mod = int.MaxValue; mod >= int.MaxValue - 20; mod--)
{
var bt = new Barrett(mod);
var v = new List<uint>();
for (uint i = 0; i < 10; i++)
{
v.Add(i);
v.Add(mod - i);
v.Add(mod / 2 + i);
v.Add(mod / 2 - i);
}
foreach (var a in v)
{
long a2 = a;
bt.Mul(a, bt.Mul(a, a)).Should().Be((uint)(a2 * a2 % mod * a2 % mod));
foreach (var b in v)
{
long b2 = b;
bt.Mul(a, b).Should().Be((uint)(a2 * b2 % mod));
}
}
}
}

}
}
95 changes: 45 additions & 50 deletions Test/ac-library-csharp.Test/Internal/InternalMathTest.cs
Original file line number Diff line number Diff line change
@@ -1,10 +1,16 @@
using System.Collections.Generic;
using System;
using System.Collections.Generic;
using System.Numerics;
using System.Runtime.InteropServices;
using FluentAssertions;
using Xunit;

namespace AtCoder.Internal
{
public interface IPrimitiveRootFactory : IStaticMod
{
int PrimitiveRoot();
}
public partial class InternalMathTest
{
static long Gcd(long a, long b)
Expand All @@ -25,51 +31,6 @@ private static bool IsPrimeNaive(long n)
return true;
}

[Fact]
public void Barrett()
{
for (uint m = 1; m <= 100; m++)
{
var bt = new Barrett(m);
for (uint a = 0; a < m; a++)
{
for (uint b = 0; b < m; b++)
{
bt.Mul(a, b).Should().Be((a * b) % m);
}
}
}

new Barrett(1).Mul(0, 0).Should().Be(0);
}

[Fact]
public void BarrettBorder()
{
for (uint mod = int.MaxValue; mod >= int.MaxValue - 20; mod--)
{
var bt = new Barrett(mod);
var v = new List<uint>();
for (uint i = 0; i < 10; i++)
{
v.Add(i);
v.Add(mod - i);
v.Add(mod / 2 + i);
v.Add(mod / 2 - i);
}
foreach (var a in v)
{
long a2 = a;
bt.Mul(a, bt.Mul(a, a)).Should().Be((uint)(a2 * a2 % mod * a2 % mod));
foreach (var b in v)
{
long b2 = b;
bt.Mul(a, b).Should().Be((uint)(a2 * b2 % mod));
}
}
}
}

[Fact]
public void IsPrime()
{
Expand Down Expand Up @@ -210,10 +171,44 @@ public void PrimitiveRootTest()
MathUtil.IsPrimitiveRoot(x, mods[x].PrimitiveRoot()).Should().BeTrue();
}
}
}

public interface IPrimitiveRootFactory : IStaticMod
{
int PrimitiveRoot();
#if NET7_0_OR_GREATER
[Theory]
[InlineData(3ul, 5ul)]
[InlineData((1ul << 32) - 1, 5ul)]
[InlineData(0xF0000000F0000000, 0xF0000000F0000000)]
[InlineData(0xF000000000000000, 0xF0000000F0000000)]
[InlineData(0x00000000F0000000, 0xF0000000F0000000)]
[InlineData((ulong)int.MaxValue, ulong.MaxValue)]
[InlineData((ulong)uint.MaxValue, ulong.MaxValue)]
[InlineData((ulong)int.MaxValue, (ulong)int.MaxValue)]
[InlineData((ulong)uint.MaxValue, (ulong)uint.MaxValue)]
[InlineData(ulong.MaxValue, ulong.MaxValue)]
public void Mul128Bit(ulong a, ulong b)
{
for (int i = -10; i <= 10; i++)
for (int j = -10; j <= 10; j++)
{
var x = a + (ulong)i;
var y = b + (ulong)j;
ulong expected = Math.BigMul(x, y, out _);
InternalMath.Mul128BitLogic(x, y).Should().Be(expected);
}
}
[Fact]
public void Mul128BitRandom()
{
var rnd = new Random(227);
var arr = new ulong[2];
for (int i = 0; i < 50000; i++)
{
rnd.NextBytes(MemoryMarshal.Cast<ulong, byte>(arr));
var x = arr[0];
var y = arr[1];
ulong expected = Math.BigMul(x, y, out _);
InternalMath.Mul128BitLogic(x, y).Should().Be(expected);
}
}
#endif
}
}

0 comments on commit 79218e2

Please sign in to comment.