Эллиптическая кривая с секущей линией

Я бы гомогенизировал это x^3 + axz^2 + bz^3 - y^2z = 0 и очки P = [Px : Py : Pz] а также Q = [Qx : Qy : Qz], Тогда любая точка R = λP + µQ с (λ, µ) ≠ (0, 0) лежит на прямой, натянутой на P и Q. Если вы хотите избежать гомогенизации, вам потребуется λ + µ = 1, но это обычно приводит до дивизий, которых я бы предпочел избегать до самого конца.

Вставьте полученные координаты R в гомогенизированное уравнение эллиптической кривой, и вы получите однородное кубическое уравнение по λ и μ, то есть что-то вроде

αλ³ + βλ²μ + γλμ² + δμ³ = 0

с α, β, γ и δ в зависимости от ваших a, b, P, Q. Для μ = 0 вы получаете вектор координат, кратный P, и, поскольку однородные координаты определяют кратные, вы получаете саму точку P, которая лежит на кривой. Таким образом, μ = 0 должно удовлетворять уравнению, поэтому вы знаете α = 0 еще до того, как вычислите его. Аналогично, λ = 0 представляет собой Q, поэтому δ = 0, если эта точка лежит на кривой. Вы остались с

(βλ + γμ)λμ = 0

Конечные два фактора кодируют два известных пересечения, которые я только что упомянул. Скобка — это третье пересечение, которое вам нужно. Теперь просто выберите λ = γ и µ = −β, чтобы получить простое выражение для третьей точки пересечения.

Если вы хотите сделать гомогенизацию в конце, просто разделите первые две координаты результирующего однородного вектора координат на третью.

Если бы я не испортил мои вычисления, у вас есть

β = 3*Px^2*Qx + 2*Px*Pz*Qz*a - Py^2*Qz - 2*Py*Pz*Qy + Pz^2*Qx*a + 3*Pz^2*Qz*b
γ = 3*Px*Qx^2 + 2*Pz*Qx*Qz*a - Pz*Qy^2 - 2*Py*Qy*Qz + Px*Qz^2*a + 3*Pz*Qz^2*b

Который, как ожидается, очень симметричен в P и Q. Поэтому, по сути, вам просто нужна одна функция, и тогда вы получите β = f (P, Q) и γ = f (Q, P).

В C ++ и с полной гомогенизацией / дегомогенизацией:

inline double f(double Px, double Py, double Qx, double Qy, double a, double b) {
return 3*Px*Px*Qx + 2*Px*a - Py*Py - 2*Py*Qy + Qx*a + 3*b;
}

std::pair<double, double> third_intersection(double Px, double Py, double Qx, double Qy, double a, double b) {
double beta = f(Px, Py, Qx, Qy, a, b);
double gamma = f(Qx, Qy, Px, Py, a, b);
double denominator = gamma - beta;  // Might be zero if line PQ is an asymptote!
double x = (gamma*Px - beta*Qx) / denominator;
double y = (gamma*Py - beta*Qy) / denominator;
return std::make_pair(x, y);
}