从给定x值的曲线中获得一个整数值的y值

你的任务很有趣!

决定用纯c++从头开始实现我自己的解决方案。

我把你的任务分成了两个主要部分:

1. 当没有模时，即你不取方程% p为某个素数

2. 当我们有素数模量时意味着方程y^2 = x^3 + x + 24 (mod p)取% p，其中p为素数。

我选择了第二部分来实现，因为你处理的是椭圆曲线方程，它通常是取模一些素数p。

第一部分(无模)可以用两种通用算法-二进制搜索和指数搜索轻松求解。

我在循环中增加x，同时使用指数搜索和二进制搜索搜索y^2的平方根。很容易看出，如果我们增加x，那么y解决方案也会平滑地增加，因此很容易使用指数搜索，因为每个下一个"y";

众所周知，指数搜索和二分搜索都非常快，它们的复杂度为O(Log2(segment_size))，这意味着即使对于128位整数，我们在该算法中花费的基本步骤也不超过128步。

无模(x, y)点是非常稀疏的，对于给定的x我们很少能找到对应的y。它需要大量的计算能力来找到至少几个解。但是我成功地找到了你特定方程的很小的解x = 12, y = 42。更多的解决方案可能需要许多小时的8-16核机器。

我用来使程序更快的一件事是我使用了一种有趣而简单的过滤平方的技术。我取了一个数字f = 2 * 2 * 2 * 3 * 3 * 5 * 7 * 11 * 13，这个数字被用作过滤器。众所周知，如果number是平方数，那么它也是任何合数或素数的平方模数，这意味着如果N是平方数，那么N % f也可能是平方数。因此，我制作了一个大小为f的预计算表，用true标记所有可能的正方形，用false标记确定的非正方形。这个函数在我的代码中被称为FilterSquare()。

我制作了第一部分多核，这意味着具有多核的高级CPU将按比例更快地解决任务。如果你有N个核(实际上是N个硬件线程)，那么使用我的c++程序，任务的解决速度将比使用1个核快N倍。

第二部分，素数模p，可以很容易地解决，因为模一些素数几乎总是有一个解，精确地在50%的情况下。

为了提醒我们执行% p的所有操作，方程的两个部分都取% p。

这个模数任务可以分为两种情况，第一种情况是p % 4 == 3，在这种情况下，众所周知，如果n的平方根存在，那么它等于root = n^((p + 1) / 4) (mod p)，非常容易和快速的公式求出平方根。

其次，如果p % 4 == 1，那么很难找到平方根，因为这种情况下发明了Tonelli-Shanks算法，我在我的代码中完全实现了它。虽然这个算法看起来相当大，但它的计算速度仍然很快，大约在O(Log2(p))左右，这意味着128位素数p最多128步。

也很容易通过欧拉准则过滤掉没有平方根的X，这个准则说，如果我们需要找到n (mod p)的平方根，那么如果n^((p-1)/2) (mod p)是1，那么n有平方根，如果它是-1，那么它没有平方根。

为了测试目的，我随机生成任意128位X和任意128位素数p，并为这个X和p输出解。你可以看到，当你运行我的程序时，对于任何大的X和p的解都会立即输出。

您可以在我的程序中使用位大小，只需转到main()函数并将128替换为例如256，以测试256位随机数。

除了上述算法外，我还使用了两个小的额外算法-费马素数检验和模幂(使用平方或二进制幂的求幂技术)。

我的程序分为两大部分——第一部分是关于无模任务，第二部分是带模任务。你可以很容易地把我的程序切成两个文件，只要找到一段文字=== Part 2 ===，之后它就成为第二个独立的部分。

CLang和GCC下的程序编译。现在在MSVC下，它不能编译，只是因为它使用内置的128位整数__int128，这只在CLang/GCC中得到支持。但是如果你需要MSVC，我可以很容易地通过使用Boost多精度库来解决这个问题，只要告诉我。也内置__int128将比Boost的更快。

编译程序的第二部分需要Boost。它可以通过sudo apt install libboost-dev-all在Linux下轻松安装。在Windows下稍微困难一点，但仍然可以通过Chocolatey的Boost Package安装(执行choco install boost-msvc-14.3，但先安装Chocolatey)。

但是这个Boost只需要在我的程序的第二部分使用模数。你可以很容易地删除我的程序的第二部分，第一部分将没有Boost编译，只有标准库。

查看我的程序在代码后的输出示例。

上网试试!

#include <cstdint>
#include <iostream>
#include <iomanip>
#include <stdexcept>
#include <sstream>
#include <mutex>
#include <vector>
#include <future>
#include <stdexcept>
#define ASSERT_MSG(cond, msg) { if (!(cond)) { throw std::runtime_error("Assertion (" #cond ") failed at line " + std::to_string(__LINE__) + "! Msg '" + std::string(msg) + "'"); } }
#define ASSERT(cond) ASSERT_MSG(cond, "")
#define COUT(code) { std::lock_guard<std::mutex> lock(cout_mux); std::cout code; std::cout << std::flush; }
#if defined(_MSC_VER) && !defined(__clang__)
#define IS_MSVC 1
#define IS_CLANG 0
#define IS_GCC 0
#elif defined(__clang__)
#define IS_MSVC 0
#define IS_CLANG 1
#define IS_GCC 0
#else
#define IS_MSVC 0
#define IS_CLANG 0
#define IS_GCC 1
#endif
using i16 = int16_t;
using u16 = uint16_t;
using i32 = int32_t;
using u32 = uint32_t;
using i64 = int64_t;
using u64 = uint64_t;
using i128c = signed __int128;
using u128c = unsigned __int128;
std::mutex cout_mux;
template <typename T, typename F>
T BinarySearch(F const & f, T begin, T end) {
auto const end0 = end;
while (begin + 2 < end) {
auto const mid = (begin + end - 1) / 2;
if (f(mid))
end = mid + 1;
else
begin = mid + 1;
}
for (auto i = begin; i < end; ++i)
if (f(i))
return i;
return end0;
}
template <typename T, typename F>
T ExpSearch(F const & f, T begin, T end, bool back = false) {
auto const begin0 = begin, end0 = end;
T l = 0;
if (!back) {
for (l = 1;; l <<= 1) {
if (begin0 + l >= end0)
break;
if (f(begin0 + l - 1))
break;
begin = begin0 + l;
}
end = std::min(begin0 + l, end0);
} else {
for (l = 1;; l <<= 1) {
if (begin0 + l >= end0)
break;
if (!f(end0 - l))
break;
end = end0 - l + 1;
}
begin = std::max(begin0 + l, end0) - l;
}
if (begin >= end || !f(end - 1))
return end0;
return BinarySearch(f, begin, end);
}
#if defined(_MSC_VER) && !defined(__clang__)
inline u64 UDiv128(u64 h, u64 l, u64 d, u64 * r) {
return _udiv128(h, l, d, r);
}
#else
inline u64 UDiv128(u64 h, u64 l, u64 d, u64 * r) {
u64 q;
asm (R"(
.intel_syntax
mov rdx, %V[h]
mov rax, %V[l]
div %V[d]
mov %V[r], rdx
mov %V[q], rax
)"
: [q] "=r" (q), [r] "=r" (*r) : [h] "r" (h), [l] "r" (l), [d] "r" (d)
: "rax", "rdx"
);
return q;
}
#endif
template <typename YT, YT Y, typename T>
YT Mod(T x) {
static_assert(sizeof(T) <= 16);
if (x <= u64(-1))
return YT(u64(x) % Y);
u64 constexpr mod = u64(-1) / Y * Y;
u64 r, q = UDiv128(u64(x >> 64), u64(x), mod, &r);
return YT(r % Y);
}
template <typename T>
bool FilterSquare(T const & x) {
static u32 constexpr filt_primo = 2 * 2 * 2 * 3 * 3 * 5 * 7 * 11 * 13;
static auto const filt = []{
std::vector<bool> filt(filt_primo);
for (uint32_t i = 0; i < filt_primo; ++i)
filt[(uint64_t(i) * i) % filt_primo] = true;
return std::move(filt);
}();
if constexpr(IS_CLANG)
return filt[Mod<u32, filt_primo>(x)];
else
return filt[size_t(x % filt_primo)];
}
template <typename T>
std::string IntToStr(T n) {
if (n == 0)
return "0";
bool sign = !(n >= 0);
ASSERT(!sign);
//if (sign) n = -n;
u64 constexpr mod = 1'000'000'000'000'000'000ULL;
std::string r;
while (n >= mod) {
u64 const rem = u64(n % mod);
r = (std::ostringstream{} << std::setw(18) << std::setfill('0') << rem).str() + r;
n /= mod;
}
r = (std::ostringstream{} << u64(n)).str() + r;
if (sign)
r = "-" + r;
return r;
}
template <typename T>
void SolveNoMod(T x_begin, T x_end) {
using MT = i128c;
T constexpr Y_max = 1ULL << 62;
T ylast = 0;
for (T x = x_begin; x < x_end; ++x) {
MT const right = MT(x) * x * x + x + 24;
if (!FilterSquare(right)) continue;
auto const f = [&](T y){ return MT(y) * y - right; };
T const y = ExpSearch<T>([&](auto y){ return f(y) >= 0; }, ylast, Y_max);
auto const fy = f(y);
if (fy == 0)
COUT(<< "x = " << x << ", y = " << y << std::endl);
ylast = y;
}
}
template <typename T>
void SolveNoModParallel(T x_begin, T x_end) {
size_t const nthr = std::thread::hardware_concurrency();
std::vector<std::future<void>> tasks;
T constexpr block = 1 << 25;
for (T x_cur = x_begin; x_cur < x_end; x_cur += block) {
while (tasks.size() >= nthr) {
for (ptrdiff_t i = ptrdiff_t(tasks.size()) - 1; i >= 0; --i)
if (tasks[i].wait_for(std::chrono::milliseconds(1)) == std::future_status::ready) {
tasks[i].get();
tasks.erase(tasks.begin() + i);
}
std::this_thread::sleep_for(std::chrono::milliseconds(5));
std::this_thread::yield();
}
T const x_begin0 = x_cur, x_end0 = std::min<T>(x_cur + block, x_end);
if (x_begin0 >= x_end0)
continue;
tasks.push_back(std::async(std::launch::async, [&, x_begin0, x_end0]{
SolveNoMod<T>(x_begin0, x_end0);
}));
}
}
// ================ Part 2 ================
#include <random>
#include <bitset>
#include <boost/multiprecision/cpp_int.hpp>
#define bisizeof(x) (sizeof(x) * 8)
#define BiSize(T) (BiSizeS<T>::value)
using u128b = boost::multiprecision::uint128_t;
using u256 = boost::multiprecision::uint256_t;
using u512 = boost::multiprecision::uint512_t;
template <typename T> struct DWordOf;
template <> struct DWordOf<u32> { using type = u64; };
template <> struct DWordOf<u64> { using type = u128b; };
template <> struct DWordOf<u128b> { using type = u256; };
template <> struct DWordOf<u256> { using type = u512; };
template <typename T> struct BiSizeS : std::integral_constant<size_t, bisizeof(T)> {};
template <> struct BiSizeS<u32> : std::integral_constant<size_t, 32> {};
template <> struct BiSizeS<u64> : std::integral_constant<size_t, 64> {};
template <> struct BiSizeS<u128b> : std::integral_constant<size_t, 128> {};
template <> struct BiSizeS<u256> : std::integral_constant<size_t, 256> {};
template <> struct BiSizeS<u512> : std::integral_constant<size_t, 512> {};
template <typename T>
T All1() { return ~T(0); }
template <typename T>
std::string IntToStrBin(T x) {
std::ostringstream ss;
for (size_t k = 0; k < BiSize(T); k += 16) {
std::bitset<16> bs = u16(x >> k);
auto bss = bs.to_string();
std::reverse(bss.begin(), bss.end());
ss << bss << " ";
}
return ss.str();
}
template <typename T>
T PowMod(T a, T b, T c) {
using DT = typename DWordOf<T>::type;
T r = 1;
while (b) {
if (b & 1)
r = T((DT(r) * a) % c);
a = T((DT(a) * a) % c);
b >>= 1;
}
return r;
}
template <typename T>
int EulerCriterion(T n, T p) {
// https://en.wikipedia.org/wiki/Euler%27s_criterion
if (n == 0)
return 0;
auto const val = PowMod(n, (p - 1) / 2, p);
ASSERT(val == 1 || val == p - 1);
return val == 1 ? 1 : -1;
}
template <typename T>
std::vector<T> Uniquize(std::vector<T> vec) {
std::sort(vec.begin(), vec.end());
vec.erase(std::unique(vec.begin(), vec.end()), vec.end());
return std::move(vec);
}
template <typename T>
std::vector<T> SquareRoot_Mod4_3(T n, T p) {
using DT = typename DWordOf<T>::type;
ASSERT(p % 4 == 3);
n %= p;
T const r = PowMod(n, (p + 1) / 4, p);
if ((DT(r) * r) % p == n % p)
return Uniquize<T>(std::vector<T>({r, p - r}));
return {};
}
template <typename T>
std::vector<T> SquareRoot_TonelliShanks(T n, T p) {
// https://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm
// see also https://eprint.iacr.org/2020/1407.pdf
using DT = typename DWordOf<T>::type;
n %= p;
ASSERT(p % 4 == 1);
T Q = p - 1, S = 0;
while ((Q & 1) == 0) {
Q >>= 1;
++S;
}
T z = 0;
for (z = 2; z < p; ++z)
if (EulerCriterion(z, p) == -1)
break;
ASSERT_MSG(z < p, "n " + IntToStr(n) + ", p " + IntToStr(p));
T M = S, c = PowMod(z, Q, p), t = PowMod(n, Q, p), R = PowMod(n, (Q + 1) / 2, p), r = 0;
while (true) {
if (t == 0) {
r = 0;
break;
}
if (t == 1) {
r = R;
break;
}
T sqt = t;
size_t i = 0;
for (i = 1; i < M; ++i) {
sqt = T((DT(sqt) * sqt) % p);
if (sqt == 1)
break;
}
ASSERT(i < M);
ASSERT(M - i - 1 < BiSize(T));
T b = PowMod(c, T(1) << size_t(M - i - 1), p);
M = i;
c = T((DT(b) * b) % p);
t = T((DT(t) * c) % p);
R = T((DT(R) * b) % p);
}
ASSERT((DT(r) * r) % p == n);
return Uniquize<T>(std::vector<T>({r, p - r}));
}
template <typename T>
std::vector<T> SquareRoot(T n, T p) {
using DT = typename DWordOf<T>::type;
n %= p;
if (n == 0)
return {0};
if (p == 2)
return {n};
{
auto const ec = EulerCriterion(n, p);
ASSERT(ec == 1 || ec == -1);
if (ec == -1)
return {};
}
if (p % 4 == 3)
return SquareRoot_Mod4_3(n, p);
return SquareRoot_TonelliShanks(n, p);
}
auto & RNG() {
thread_local std::mt19937_64 rng{(u64(std::random_device{}()) << 32) ^ u64(std::random_device{}())};
//thread_local std::mt19937_64 rng{123};
return rng;
}
template <typename T>
T Rand() {
T r = 0;
for (size_t i = 0; i < BiSize(T); i += 64)
r ^= T(RNG()()) << i;
return r;
}
template <typename T>
size_t NumBits(T n) {
size_t cnt = 0;
while (n > 0xFF) {
n >>= 8;
cnt += 8;
}
while (n) {
n >>= 1;
++cnt;
}
return cnt;
}
template <typename T>
T RandRange(T end) {
if (end <= 1)
return 0;
size_t const bits = NumBits(end - 1);
while (true) {
auto n = Rand<T>() >> (BiSize(T) - bits);
if (n < end)
return n;
}
}
template <typename T>
struct UniformIntDistr {
UniformIntDistr(T a, T b) : a_(a), diff_(b + 1 - a) { ASSERT(a <= b); }
template <typename R>
T operator()(R & r) { return a_ + RandRange<T>(diff_); }
T a_ = 0, diff_ = 0;
};
template <typename T>
T RandUni(T a, T b) {
return UniformIntDistr<T>(a, b)(RNG());
}
template <typename T>
bool IsProbablyPrime_Fermat(T N) {
// https://en.wikipedia.org/wiki/Fermat_primality_test
size_t const ntrials = 24;
if (N <= 10)
return N == 2 || N == 3 || N == 5 || N == 7;
for (size_t trial = 0; trial < ntrials; ++trial)
if (PowMod<T>(Rand<T>() % (N - 3) + 2, N - 1, N) != 1)
return false;
return true;
}
#define IsProbablyPrime IsProbablyPrime_Fermat
template <typename T>
T GenPrime(size_t bits) {
ASSERT(bits >= 2);
if (bits == 2)
return RandUni(2, 3);
ASSERT_MSG(1 <= bits && bits <= BiSize(T), "bits " + std::to_string(bits) + ", BiSize(T) " + std::to_string(BiSize(T)));
auto const a = T(1) << (bits - 1), b = All1<T>() >> (BiSize(T) - bits);
ASSERT_MSG(a <= b, "bits " + std::to_string(bits) + ", BiSize(T) " + std::to_string(BiSize(T)) + ", an" + IntToStrBin(a) + "nbn" + IntToStrBin(b));
while (true) {
auto const p = RandUni(a, b) | 1;
if (IsProbablyPrime(p))
return p;
}
}
template <typename T>
void SolveMod(T x_begin, T mod) {
using DT = typename DWordOf<T>::type;
size_t cnt = 0;
for (T x = x_begin;; ++x) {
T const right = T((DT((DT(x) * x) % mod) * x + x + 24) % mod);
for (auto const root: SquareRoot<T>(right, mod)) {
ASSERT((DT(root) * root) % mod == right);
COUT(<< "x = " << IntToStr(x) << ", y = " << IntToStr(root) << std::endl);
++cnt;
if (cnt >= 1)
return;
}
}
}
int main() {
try {
COUT(<< "NoMod:" << std::endl);
SolveNoModParallel<i64>(0, 1ULL << 25);

for (size_t i = 0; i < 7; ++i) {
auto const p = GenPrime<u128b>(128);
COUT(<< "Mod p = " << IntToStr(p) << " (p % 4 = " << p % 4 << ")" << std::endl);
SolveMod<u128b>(Rand<u128b>(), p);
}
return 0;
} catch (std::exception const & ex) {
std::cout << "Exception: " << ex.what() << std::endl;
return -1;
}
}

输出:

NoMod:
x = 12, y = 42
Mod p = 219688491259897122736212527051092113293 (p % 4 = 1)
x = 273903578792010575565282758244483002178, y = 58545642357590275355420738836818145161
Mod p = 331730100473713163326635571798951273061 (p % 4 = 1)
x = 250565107208593852657423435334728568634, y = 61465560240598325718665905365177214482
Mod p = 252783594113540225392633573791281784563 (p % 4 = 3)
x = 288516005161218077788215822752214504957, y = 117961729460812361108079791544444521231
Mod p = 215955488059493315966305338643349539739 (p % 4 = 3)
x = 89919490689630945164185961723990553146, y = 54128779781633285229073363977029067647
Mod p = 313197233417624767187155292666519557093 (p % 4 = 1)
x = 41108118867338545127706112908846297142, y = 20447436839455724175052447276214468753
Mod p = 263976413751637208893338879428427068801 (p % 4 = 1)
x = 165225955028296917988729655751990538414, y = 62234915643079594692200530956112078114
Mod p = 291337940475282224323312279077405757331 (p % 4 = 3)
x = 160948015593193411286972051949433960251, y = 65875958798817617640182672547677071647

选择一个目标，比如y = 3;Y ^2 = 9。求一个x值x1，它的结果是y1 <= 3。求第二个x值x2，使y>= 3。假设你的方程是连续的，那么你的目标x值将在x1和x2之间。使用二进制搜索在您的机器允许的实际值范围内找到最接近的x值。

或者，尝试不同的x1和x2，直到找到它们之间至少有一个整数的一对。

考虑到在有限的计算机上表示实数的实际限制，你可能只能找到一个适用于有限小数位数的解:y = 3.0000000000120694…等等。