Skip to content

Commit

Permalink
fix: thread safe
Browse files Browse the repository at this point in the history
  • Loading branch information
@tk-yoshimura
tk-yoshimura committed Nov 17, 2024
1 parent 62a684e commit 2c8277b
Show file tree
Hide file tree
Showing 9 changed files with 40 additions and 27 deletions.
11 changes: 7 additions & 4 deletions DoubleDoubleODE/BogackiShampineODESolver.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<ddouble, ddouble> f)
: base(v, f) { }
Expand All @@ -21,10 +22,10 @@ public BogackiShampineODESolver(ddouble[] v, Func<ddouble[], ddouble[]> 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];
Expand All @@ -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;
}
}
}
8 changes: 5 additions & 3 deletions DoubleDoubleODE/DormandPrinceODESolver.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -31,7 +31,7 @@ public DormandPrinceODESolver(ddouble[] v, Func<ddouble[], ddouble[]> 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];
Expand Down Expand Up @@ -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;
}
}
}
2 changes: 1 addition & 1 deletion DoubleDoubleODE/DoubleDoubleODE.csproj
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
</PropertyGroup>

<ItemGroup>
<PackageReference Include="TYoshimura.DoubleDouble" Version="4.0.0" />
<PackageReference Include="TYoshimura.DoubleDouble" Version="4.2.2" />
</ItemGroup>

</Project>
5 changes: 4 additions & 1 deletion DoubleDoubleODE/EulerODESolver.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,10 +20,13 @@ public EulerODESolver(ddouble[] v, Func<ddouble[], ddouble[]> 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;
}
}
}
8 changes: 5 additions & 3 deletions DoubleDoubleODE/HeunODESolver.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,10 +19,10 @@ public HeunODESolver(ddouble[] v, Func<ddouble[], ddouble[]> 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];
Expand All @@ -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;
}
}
}
18 changes: 9 additions & 9 deletions DoubleDoubleODE/ODESolver.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@

namespace DoubleDoubleODE {
public abstract class ODESolver {
protected readonly ddouble[] v;
protected ddouble[] v;
protected readonly Func<ddouble[], ddouble[]> f;

public int Params => v.Length;
Expand All @@ -21,23 +21,23 @@ public abstract class ODESolver {
public ddouble this[int index] => v[index];

public ODESolver(ddouble v, Func<ddouble, ddouble> 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<ddouble, ddouble, (ddouble dx, ddouble dy)> 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<ddouble, ddouble, ddouble, (ddouble dx, ddouble dy, ddouble dz)> 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<ddouble, ddouble, ddouble, ddouble, (ddouble dx, ddouble dy, ddouble dz, ddouble dw)> 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<ddouble[], ddouble[]> f) {
Expand Down
3 changes: 1 addition & 2 deletions DoubleDoubleODE/Properties/AssemblyInfo.cs
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
using System.Reflection;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;

[assembly: AssemblyTitle("DoubleDouble.ODE")]
Expand All @@ -15,4 +14,4 @@

[assembly: Guid("FE7AFD95-5A6A-48AE-8519-7D37ACFDC5EB")]

[assembly: AssemblyVersion("1.4.0.*")]
[assembly: AssemblyVersion("1.4.1.*")]
10 changes: 7 additions & 3 deletions DoubleDoubleODE/RungeKutta4ODESolver.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,8 @@

namespace DoubleDoubleODE {
public class RungeKutta4ODESolver : ODESolver {
private static readonly ddouble c_1d6 = (ddouble)1 / 6;

public RungeKutta4ODESolver(ddouble v, Func<ddouble, ddouble> f)
: base(v, f) { }

Expand All @@ -19,10 +21,10 @@ public RungeKutta4ODESolver(ddouble[] v, Func<ddouble[], ddouble[]> 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];
Expand All @@ -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;
}
}
}
2 changes: 1 addition & 1 deletion DoubleDoubleODETest/DoubleDoubleODETest.csproj
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
<PrivateAssets>all</PrivateAssets>
<IncludeAssets>runtime; build; native; contentfiles; analyzers; buildtransitive</IncludeAssets>
</PackageReference>
<PackageReference Include="TYoshimura.DoubleDouble" Version="4.0.0" />
<PackageReference Include="TYoshimura.DoubleDouble" Version="4.2.2" />
</ItemGroup>

<ItemGroup>
Expand Down

0 comments on commit 2c8277b

Please sign in to comment.