問題はこちら。
問題の主眼
計算量の削減パートは公式解説に譲るとして、実質的に解けばよい問題は以下になります。
高橋くんは座標 $(r_t, c_t)$ から方向 $d_t$ に、青木くんは座標 $(r_a, c_a)$ から方向 $d_a$ にそれぞれ同じ速さで長さ $l$ だけ移動します。
方向は上下左右の4通りです。
2人が「同時に」同じ座標にいることがあるかどうか、lの大きさに依存しない方法で求めてください。
そこまで難しくないけど面倒です。面倒だなあと思いながらモニタを眺めていたら時間内に実装できませんでした。
心を込めて16通りの場合分けを書けば解けるのですが、やりたくないので、ほかの方法を考えます。
問題の言い換え
高橋くんと青木くんを点と捉えれば、2点が2次元座標内を動く問題とみなすことができるので、思いつくのは2人の軌跡が交差するかどうかです。
ただし、線分が交差するかどうかだけでは、タイミングの問題があり、二人が「同時に」同じ座標にいるかどうかはわかりません。
そこで、問題を3次元に拡張します。
時間tを変数として導入し、高橋くんの移動の軌跡を $(r_t, c_t, t_t)$ から $(r_t + dr_t, c_t + dc_t, t_t + l)$ への線分として表現します。
青木くんの軌跡も同様に表すことで、もとの問題は、3次元空間上でこの2つの線分が交わるか？という問題に帰着されました。
2つの線分に共有点があれば、このとき2人は同じ時間で、同じ座標にいることがわかります。
実は、この判定方法は完全ではなく、共有点の座標が整数でない場合があるのですが、この問題の場合はパリティで判断ができるため、前処理で弾いておくことで適用できます。
3次元の線分が交差するか判定する方法
2次元の場合の求め方は、たとえばこちらの記事にあります。
https://qiita.com/zu_rin/items/e04fdec4e3dec6072104
これを拡張することを考えます。
2つの線分が同一平面上にあるか判定する
3次元の場合、「ねじれの位置」の場合があるので、まず2つの線分が同一平面上にあることが必要です。
実は、$\vec{u} = \vec{AB}, \vec{v} = \vec{CD}$ が同一平面上にあるか判定する方法は、
$$
(\vec{c} - \vec{a}) \cdot (\vec{u} \times \vec{v}) = 0
$$
です。私は知らなかったのですがChatGPTが教えてくれました。
線分 $AC$ と、ベクトル $\vec{u}, \vec{v}$ のなす平面が直交しているか、ということですね。なるほど～。
2つの線分が直交するか判定する
同一平面上にあることがわかったので、上の方法が使えそうです。
平面の方程式を用いて適当な2次元の座標に帰着させて上の方法を適用してもよさそうですが、ここでは外積を用いる方法に着目します。
上の記事で紹介されている考え方は、「線分のどちら側にあるかは外積の符号で判断できる」というもので、条件は以下です。
$$
AB \times AC \cdot AB \times AD < 0 かつ
CD \times CA \cdot CD \times CB < 0
$$
2次元のベクトルの場合、外積はスカラー値なので、ドットはスカラー積になることに注意してください。
3次元のベクトルの外積はベクトルですが、考え方は応用できそうです。これらは平面と垂直なので、方向ベクトルが等しいことがわかります。つまり、2つのベクトルの関係は、「同じ向き」か「真逆の向き」で、これらの向きが違う場合に「違う側にある」といえそうです。
したがって、2つの外積の内積を取り、符号を見ることで、同じ判定方法が使えそうです。
注意すべき点として、どれかの端点が他の線分の上にある場合、外積は0となりこの条件を満たしません。そのようなケースは別途考慮する必要があります。
実装
実装は結構面倒で、オーバーフローの問題があります。
3次元のベクトルの外積は、各項がint程度であっても外積はlong程度まで数値が大きくなります。内積に至ってはlongですらオーバーフローするケースがあり、Int128が必要です。
intでこれですから、Vector<long>に対する内積はBigIntegerじゃないと表現できないケースがあり、きれいな実装が作れていません。（一応upsolveしたのですが、かなり汚いコードになってしまった……)
doubleならばオーバーフローの心配はほぼありませんが、今度は浮動小数点による誤差の扱いがあり、doubleで作ればうまくいくよ、とは筆者には言い切れません。
一応、筆者の実装を掲載しますが、高速化のための黒魔術がひどく可読性を下げており、INumberで作るのはちょっとつらみが深いかも。
using System.Numerics;
public readonly record struct Point3D<T>(T X, T Y, T Z) where T : INumber<T>
{
public static Vec3D<T> operator -(Point3D<T> left, Point3D<T> right)
{
return Vec3D<T>.FromPoints(right, left);
}
}
public readonly record struct Vec3D<T>(T X, T Y, T Z) where T : INumber<T>
{
public static Vec3D<T> FromPoints(Point3D<T> start, Point3D<T> end)
{
return new Vec3D<T>(end.X - start.X, end.Y - start.Y, end.Z - start.Z);
}
public Vec3D<TResult> As<TResult>() where TResult : INumber<TResult>
{
return new Vec3D<TResult>(
TResult.CreateChecked(this.X),
TResult.CreateChecked(this.Y),
TResult.CreateChecked(this.Z)
);
}
/// <summary>
/// 2つのベクトルの内積を計算します。
/// </summary>
/// <param name="v"></param>
/// <returns></returns>
public T InnerProduct(Vec3D<T> v)
{
return this.X * v.X
+ this.Y * v.Y
+ this.Z * v.Z;
}
/// <summary>
/// 2つのベクトルの内積を計算します。
/// </summary>
/// <param name="v"></param>
/// <returns></returns>
public TResult InnerProduct<TResult>(Vec3D<T> v) where TResult : INumber<TResult>
{
//JIT最適化で消えてくれてちょっと速くなることを期待……
if (typeof(T) == typeof(long) && typeof(TResult) == typeof(BigInteger))
{
return (TResult)(object)(
(BigInteger)(long)(object)this.X * (BigInteger)(long)(object)v.X
+ (BigInteger)(long)(object)this.Y * (BigInteger)(long)(object)v.Y
+ (BigInteger)(long)(object)this.Z * (BigInteger)(long)(object)v.Z
);
}
if (typeof(T) == typeof(Int128) && typeof(TResult) == typeof(BigInteger))
{
return (TResult)(object)(
(BigInteger)(Int128)(object)this.X * (BigInteger)(Int128)(object)v.X
+ (BigInteger)(Int128)(object)this.Y * (BigInteger)(Int128)(object)v.Y
+ (BigInteger)(Int128)(object)this.Z * (BigInteger)(Int128)(object)v.Z
);
}
return TResult.CreateChecked(this.X) * TResult.CreateChecked(v.X)
+ TResult.CreateChecked(this.Y) * TResult.CreateChecked(v.Y)
+ TResult.CreateChecked(this.Z) * TResult.CreateChecked(v.Z);
}
/// <summary>
/// 2つのベクトルの内積を計算し、符号のみを返します。
/// オーバーフローしないように頑張ります。
/// </summary>
/// <param name="v"></param>
/// <returns></returns>
public int InnerProductSign(Vec3D<T> v)
{
//このifは多分JIT最適化で消える たぶん
if (typeof(T) == typeof(Int128))
{
return InnerProduct<BigInteger>(v).Sign;
}
else if (typeof(T) == typeof(long))
{
return InnerProduct<BigInteger>(v).Sign;
}
else if (typeof(T) == typeof(int))
{
return Int128.Sign(InnerProduct<Int128>(v));
}
return T.Sign(this.InnerProduct(v));
}
/// <summary>
/// 2つのベクトルの外積を計算します。
/// </summary>
/// <param name="v"></param>
/// <returns></returns>
public Vec3D<T> OuterProduct(Vec3D<T> v)
{
return new Vec3D<T>(
this.Y * v.Z - this.Z * v.Y,
this.Z * v.X - this.X * v.Z,
this.X * v.Y - this.Y * v.X
);
}
/// <summary>
/// 2つのベクトルの外積を計算します。
/// </summary>
/// <param name="v"></param>
/// <returns></returns>
public Vec3D<TResult> OuterProduct<TResult>(Vec3D<T> v) where TResult : INumber<TResult>
{
//JIT最適化で消えてくれてちょっと速くなることを期待……
if (typeof(T) == typeof(long) && typeof(TResult) == typeof(Int128))
{
return new Vec3D<TResult>(
(TResult)(object)((Int128)(long)(object)(this.Y) * (Int128)(long)(object)(v.Z) - (Int128)(long)(object)(this.Z) * (Int128)(long)(object)(v.Y)),
(TResult)(object)((Int128)(long)(object)(this.Z) * (Int128)(long)(object)(v.X) - (Int128)(long)(object)(this.X) * (Int128)(long)(object)(v.Z)),
(TResult)(object)((Int128)(long)(object)(this.X) * (Int128)(long)(object)(v.Y) - (Int128)(long)(object)(this.Y) * (Int128)(long)(object)(v.X))
);
}
return new Vec3D<TResult>(
TResult.CreateChecked(this.Y) * TResult.CreateChecked(v.Z) - TResult.CreateChecked(this.Z) * TResult.CreateChecked(v.Y),
TResult.CreateChecked(this.Z) * TResult.CreateChecked(v.X) - TResult.CreateChecked(this.X) * TResult.CreateChecked(v.Z),
TResult.CreateChecked(this.X) * TResult.CreateChecked(v.Y) - TResult.CreateChecked(this.Y) * TResult.CreateChecked(v.X)
);
}
public static Vec3D<T> operator -(Vec3D<T> left, Vec3D<T> right)
{
return new Vec3D<T>(left.X - right.X, left.Y - right.Y, left.Z - right.Z);
}
}
public static class Geometry
{
/// <summary>
/// 2つの線分がクロスするかどうかを判定(端点が重なる場合はFalse)
/// </summary>
/// <param name="A"></param>
/// <param name="B"></param>
/// <param name="C"></param>
/// <param name="D"></param>
public static bool IsCross<T>(
Point3D<T> A,
Point3D<T> B,
Point3D<T> C,
Point3D<T> D
) where T : INumber<T>
{
if (!IsSamePlane(A, B, C, D))
{
return false;
}
var AB = B - A;
var CD = D - C;
var AC = C - A;
var AD = D - A;
var CA = A - C;
var CB = B - C;
var ABxAC = AB.OuterProduct<Int128>(AC);
var ABxAD = AB.OuterProduct<Int128>(AD);
var CDxCA = CD.OuterProduct<Int128>(CA);
var CDxCB = CD.OuterProduct<Int128>(CB);
return ABxAC.InnerProductSign(ABxAD) < 0
&& CDxCA.InnerProductSign(CDxCB) < 0;
}
/// <summary>
/// 2つの線分が同一平面上にあるかを判定
/// </summary>
/// <param name="A"></param>
/// <param name="B"></param>
/// <param name="C"></param>
/// <param name="D"></param>
public static bool IsSamePlane<T>(
Point3D<T> A,
Point3D<T> B,
Point3D<T> C,
Point3D<T> D
) where T : INumber<T>
{
var u = B - A;
var v = D - C;
//(b-a)・(u × v) == 0 が条件
var ca = (C - A).As<Int128>();
return ca.InnerProductSign(u.OuterProduct<Int128>(v)) == 0;
}
}
まとめ
- 3次元での線分の交差判定は実装できるよ
- でも数字が大きくなりすぎてオーバーフローしない保証がちょっとつらいよ
- この程度の問題のために持ち出す大道具ではなかったかも