From 2c8277b110019c4f917cce567b8784f1f4d59e17 Mon Sep 17 00:00:00 2001 From: tk-yoshimura Date: Mon, 18 Nov 2024 03:09:49 +0900 Subject: [PATCH] fix: thread safe --- DoubleDoubleODE/BogackiShampineODESolver.cs | 11 +++++++---- DoubleDoubleODE/DormandPrinceODESolver.cs | 8 +++++--- DoubleDoubleODE/DoubleDoubleODE.csproj | 2 +- DoubleDoubleODE/EulerODESolver.cs | 5 ++++- DoubleDoubleODE/HeunODESolver.cs | 8 +++++--- DoubleDoubleODE/ODESolver.cs | 18 +++++++++--------- DoubleDoubleODE/Properties/AssemblyInfo.cs | 3 +-- DoubleDoubleODE/RungeKutta4ODESolver.cs | 10 +++++++--- DoubleDoubleODETest/DoubleDoubleODETest.csproj | 2 +- 9 files changed, 40 insertions(+), 27 deletions(-) diff --git a/DoubleDoubleODE/BogackiShampineODESolver.cs b/DoubleDoubleODE/BogackiShampineODESolver.cs index a80a965..57fec24 100644 --- a/DoubleDoubleODE/BogackiShampineODESolver.cs +++ b/DoubleDoubleODE/BogackiShampineODESolver.cs @@ -3,7 +3,8 @@ namespace DoubleDoubleODE { public class BogackiShampineODESolver : ODESolver { - private static readonly ddouble c_1d2 = (ddouble)1 / 2, c_3d4 = (ddouble)3 / 4, c_1d9 = (ddouble)1 / 9; + private const double c_3d4 = 0.75d; + private static readonly ddouble c_1d9 = (ddouble)1 / 9; public BogackiShampineODESolver(ddouble v, Func f) : base(v, f) { } @@ -21,10 +22,10 @@ public BogackiShampineODESolver(ddouble[] v, Func f) : base(v, f) { } public override void Next(ddouble h) { - ddouble h_1d2 = h * c_1d2, h_3d4 = h * c_3d4, h_1d9 = h * c_1d9; + ddouble h_1d2 = ddouble.Ldexp(h, -1), h_3d4 = h * c_3d4, h_1d9 = h * c_1d9; ddouble[] dv1 = f(v); - ddouble[] v2 = new ddouble[Params]; + ddouble[] v2 = new ddouble[Params], v_new = new ddouble[Params]; for (int i = 0; i < Params; i++) { v2[i] = v[i] + h_1d2 * dv1[i]; @@ -40,8 +41,10 @@ public override void Next(ddouble h) { ddouble[] dv3 = f(v3); for (int i = 0; i < Params; i++) { - v[i] += h_1d9 * (2 * dv1[i] + 3 * dv2[i] + 4 * dv3[i]); + v_new[i] = v[i] + h_1d9 * (ddouble.Ldexp(dv1[i] + dv2[i] + ddouble.Ldexp(dv3[i], 1), 1) + dv2[i]); } + + v = v_new; } } } diff --git a/DoubleDoubleODE/DormandPrinceODESolver.cs b/DoubleDoubleODE/DormandPrinceODESolver.cs index c3002f0..2ad49f5 100644 --- a/DoubleDoubleODE/DormandPrinceODESolver.cs +++ b/DoubleDoubleODE/DormandPrinceODESolver.cs @@ -2,7 +2,7 @@ using System; namespace DoubleDoubleODE { - public class DormandPrinceODESolver : ODESolver { + public class DormandPrinceODESolver : ODESolver { private static readonly ddouble c_1d5 = (ddouble)1 / 5, c_3d40 = (ddouble)3 / 40; private static readonly ddouble c_9d40 = (ddouble)9 / 40, c_44d45 = (ddouble)44 / 45; private static readonly ddouble c_56d15 = (ddouble)56 / 15, c_32d9 = (ddouble)32 / 9; @@ -31,7 +31,7 @@ public DormandPrinceODESolver(ddouble[] v, Func f) public override void Next(ddouble h) { ddouble[] dv1 = f(v); - ddouble[] v2 = new ddouble[Params]; + ddouble[] v2 = new ddouble[Params], v_new = new ddouble[Params]; for (int i = 0; i < Params; i++) { v2[i] = v[i] + h * c_1d5 * dv1[i]; @@ -68,8 +68,10 @@ public override void Next(ddouble h) { ddouble[] dv6 = f(v6); for (int i = 0; i < Params; i++) { - v[i] += h * (c_35d384 * dv1[i] + c_500d1113 * dv3[i] + c_125d192 * dv4[i] - c_2187d6784 * dv5[i] + c_11d84 * dv6[i]); + v_new[i] = v[i] + h * (c_35d384 * dv1[i] + c_500d1113 * dv3[i] + c_125d192 * dv4[i] - c_2187d6784 * dv5[i] + c_11d84 * dv6[i]); } + + v = v_new; } } } diff --git a/DoubleDoubleODE/DoubleDoubleODE.csproj b/DoubleDoubleODE/DoubleDoubleODE.csproj index 29bd4fc..f7c45ed 100644 --- a/DoubleDoubleODE/DoubleDoubleODE.csproj +++ b/DoubleDoubleODE/DoubleDoubleODE.csproj @@ -16,7 +16,7 @@ - + diff --git a/DoubleDoubleODE/EulerODESolver.cs b/DoubleDoubleODE/EulerODESolver.cs index 34ca5ac..5017b27 100644 --- a/DoubleDoubleODE/EulerODESolver.cs +++ b/DoubleDoubleODE/EulerODESolver.cs @@ -20,10 +20,13 @@ public EulerODESolver(ddouble[] v, Func f) public override void Next(ddouble h) { ddouble[] dv = f(v); + ddouble[] v_new = new ddouble[Params]; for (int i = 0; i < Params; i++) { - v[i] += h * dv[i]; + v_new[i] = v[i] + h * dv[i]; } + + v = v_new; } } } diff --git a/DoubleDoubleODE/HeunODESolver.cs b/DoubleDoubleODE/HeunODESolver.cs index 02ae7ac..c21f5fa 100644 --- a/DoubleDoubleODE/HeunODESolver.cs +++ b/DoubleDoubleODE/HeunODESolver.cs @@ -19,10 +19,10 @@ public HeunODESolver(ddouble[] v, Func f) : base(v, f) { } public override void Next(ddouble h) { - ddouble h_half = h / 2; + ddouble h_half = ddouble.Ldexp(h, -1); ddouble[] dv1 = f(v); - ddouble[] v2 = new ddouble[Params]; + ddouble[] v2 = new ddouble[Params], v_new = new ddouble[Params]; for (int i = 0; i < Params; i++) { v2[i] = v[i] + h * dv1[i]; @@ -31,8 +31,10 @@ public override void Next(ddouble h) { ddouble[] dv2 = f(v2); for (int i = 0; i < Params; i++) { - v[i] += h_half * (dv1[i] + dv2[i]); + v_new[i] = v[i] + h_half * (dv1[i] + dv2[i]); } + + v = v_new; } } } diff --git a/DoubleDoubleODE/ODESolver.cs b/DoubleDoubleODE/ODESolver.cs index 396f824..7269207 100644 --- a/DoubleDoubleODE/ODESolver.cs +++ b/DoubleDoubleODE/ODESolver.cs @@ -3,7 +3,7 @@ namespace DoubleDoubleODE { public abstract class ODESolver { - protected readonly ddouble[] v; + protected ddouble[] v; protected readonly Func f; public int Params => v.Length; @@ -21,23 +21,23 @@ public abstract class ODESolver { public ddouble this[int index] => v[index]; public ODESolver(ddouble v, Func f) { - this.v = new ddouble[] { v }; - this.f = (v) => new ddouble[] { f(v[0]) }; + this.v = [v]; + this.f = (v) => [f(v[0])]; } public ODESolver((ddouble x, ddouble y) v, Func f) { - this.v = new ddouble[] { v.x, v.y }; - this.f = (v) => { var (dx, dy) = f(v[0], v[1]); return new ddouble[] { dx, dy }; }; + this.v = [v.x, v.y]; + this.f = (v) => { var (dx, dy) = f(v[0], v[1]); return [dx, dy]; }; } public ODESolver((ddouble x, ddouble y, ddouble z) v, Func f) { - this.v = new ddouble[] { v.x, v.y, v.z }; - this.f = (v) => { var (dx, dy, dz) = f(v[0], v[1], v[2]); return new ddouble[] { dx, dy, dz }; }; + this.v = [v.x, v.y, v.z]; + this.f = (v) => { var (dx, dy, dz) = f(v[0], v[1], v[2]); return [dx, dy, dz]; }; } public ODESolver((ddouble x, ddouble y, ddouble z, ddouble w) v, Func f) { - this.v = new ddouble[] { v.x, v.y, v.z, v.w }; - this.f = (v) => { var (dx, dy, dz, dw) = f(v[0], v[1], v[2], v[3]); return new ddouble[] { dx, dy, dz, dw }; }; + this.v = [v.x, v.y, v.z, v.w]; + this.f = (v) => { var (dx, dy, dz, dw) = f(v[0], v[1], v[2], v[3]); return [dx, dy, dz, dw]; }; } public ODESolver(ddouble[] v, Func f) { diff --git a/DoubleDoubleODE/Properties/AssemblyInfo.cs b/DoubleDoubleODE/Properties/AssemblyInfo.cs index 62b3e61..d085f9c 100644 --- a/DoubleDoubleODE/Properties/AssemblyInfo.cs +++ b/DoubleDoubleODE/Properties/AssemblyInfo.cs @@ -1,5 +1,4 @@ using System.Reflection; -using System.Runtime.CompilerServices; using System.Runtime.InteropServices; [assembly: AssemblyTitle("DoubleDouble.ODE")] @@ -15,4 +14,4 @@ [assembly: Guid("FE7AFD95-5A6A-48AE-8519-7D37ACFDC5EB")] -[assembly: AssemblyVersion("1.4.0.*")] \ No newline at end of file +[assembly: AssemblyVersion("1.4.1.*")] \ No newline at end of file diff --git a/DoubleDoubleODE/RungeKutta4ODESolver.cs b/DoubleDoubleODE/RungeKutta4ODESolver.cs index 9ed97be..65d03d7 100644 --- a/DoubleDoubleODE/RungeKutta4ODESolver.cs +++ b/DoubleDoubleODE/RungeKutta4ODESolver.cs @@ -3,6 +3,8 @@ namespace DoubleDoubleODE { public class RungeKutta4ODESolver : ODESolver { + private static readonly ddouble c_1d6 = (ddouble)1 / 6; + public RungeKutta4ODESolver(ddouble v, Func f) : base(v, f) { } @@ -19,10 +21,10 @@ public RungeKutta4ODESolver(ddouble[] v, Func f) : base(v, f) { } public override void Next(ddouble h) { - ddouble h_1d2 = h / 2, h_1d6 = h / 6; + ddouble h_1d2 = ddouble.Ldexp(h, -1), h_1d6 = h * c_1d6; ddouble[] dv1 = f(v); - ddouble[] v2 = new ddouble[Params]; + ddouble[] v2 = new ddouble[Params], v_new = new ddouble[Params]; for (int i = 0; i < Params; i++) { v2[i] = v[i] + h_1d2 * dv1[i]; @@ -45,8 +47,10 @@ public override void Next(ddouble h) { ddouble[] dv4 = f(v4); for (int i = 0; i < Params; i++) { - v[i] += h_1d6 * (dv1[i] + 2 * (dv2[i] + dv3[i]) + dv4[i]); + v_new[i] = v[i] + h_1d6 * (dv1[i] + ddouble.Ldexp(dv2[i] + dv3[i], 1) + dv4[i]); } + + v = v_new; } } } diff --git a/DoubleDoubleODETest/DoubleDoubleODETest.csproj b/DoubleDoubleODETest/DoubleDoubleODETest.csproj index 828864f..b203b15 100644 --- a/DoubleDoubleODETest/DoubleDoubleODETest.csproj +++ b/DoubleDoubleODETest/DoubleDoubleODETest.csproj @@ -14,7 +14,7 @@ all runtime; build; native; contentfiles; analyzers; buildtransitive - +