I will also add the scalar values that cause the anomaly on the secp256k1 field prime p, i.e., the cases where p equals n. That will probably be a first as well.
from math import isqrt
def extract_small_square_factor(n: int, bound: int = 20000):
sf = 1
m = n
sieve = [True] * (bound + 1)
primes = []
for i in range(2, bound + 1):
if sieve[i]:
primes.append(i)
ii = i * i
if ii <= bound:
step = i
sieve[ii:bound + 1:step] = [False] * (((bound - ii) // step) + 1)
for q in primes:
q2 = q * q
if q2 > m:
break
while m % q2 == 0:
m //= q2
sf *= q
return sf, m
def report_anomalous_difficulty(p: int, bound: int = 20000):
Da = 4 * p - 1 # = -(1 - 4p)
sf, rem = extract_small_square_factor(Da, bound=bound)
print(f"p = {p}")
print(f"Delta = 1 - 4p = -({Da})")
print(f"small square factor: {sf * sf} (sf={sf})")
print(f"remaining cofactor bit length (Da / sf^2): {rem.bit_length()}")
sp = isqrt(p)
print(f"expected count of order==p cases ~ sqrt(p)/4 ≈ {sp // 4}")
print(f"hit probability for random k ~ 1/(4*sqrt(p)) ≈ 1/{4 * sp}")
print()
for P in (
79,
967,
1303,
115792089237316195423570985008687907853269984665640564039457584007908834671663,
): report_anomalous_difficulty(P)
p = 79
Delta = 1 - 4p = -(315)
small square factor: 9 (sf=3)
remaining cofactor bit length (Da / sf^2): 6
expected count of order==p cases ~ sqrt(p)/4 ≈ 2
hit probability for random k ~ 1/(4*sqrt(p)) ≈ 1/32
p = 967
Delta = 1 - 4p = -(3867)
small square factor: 1 (sf=1)
remaining cofactor bit length (Da / sf^2): 12
expected count of order==p cases ~ sqrt(p)/4 ≈ 7
hit probability for random k ~ 1/(4*sqrt(p)) ≈ 1/124
p = 1303
Delta = 1 - 4p = -(5211)
small square factor: 9 (sf=3)
remaining cofactor bit length (Da / sf^2): 10
expected count of order==p cases ~ sqrt(p)/4 ≈ 9
hit probability for random k ~ 1/(4*sqrt(p)) ≈ 1/144
p = 115792089237316195423570985008687907853269984665640564039457584007908834671663
Delta = 1 - 4p = -(463168356949264781694283940034751631413079938662562256157830336031635338686651)
small square factor: 9 (sf=3)
remaining cofactor bit length (Da / sf^2): 255
expected count of order==p cases ~ sqrt(p)/4 ≈ 85070591730234615865843651857942052863
hit probability for random k ~ 1/(4*sqrt(p)) ≈ 1/1361129467683753853853498429727072845820
This script gives a quick, practical summary of how rare and difficult the “order equals p” (anomalous) target is for each prime p you test.
In the output,The “square factor” line shows whether there is any small square factor hiding inside the relevant value tied to p (some primes have a small square factor, some do not).
The “bit length” line indicates how large the remaining cofactor is; the larger it is, the more complex the structure typically is.
The “expected count” line is a rough estimate of how many such special cases you might expect for that p.
The “hit probability for random k” line shows how unlikely it is to hit this special case by choosing k at random.
Overall; for small p it can be feasible to search, but for huge p (like the secp field prime) random searching becomes practically hopeless.
I have never come across a curve with order 115792089237316195423570985008687907853269984665640564039457584007908834671663!
I think it has never been shared.
Curves of type p=n are vulnerable to various attack methods, both known and not widely known, such as Smart's Attack (p-adic math). These curves are not used because they are known to be weak. However, this will continue until a bridging method is discovered between the weak curve and the original curve.
import random
from math import isqrt
def legendre(a, p):
a %= p
if a == 0:
return 0
return 1 if pow(a, (p - 1) // 2, p) == 1 else -1
def tonelli_shanks(n, p):
n %= p
if n == 0:
return 0
if pow(n, (p - 1) // 2, p) != 1:
return None
if p % 4 == 3:
return pow(n, (p + 1) // 4, p)
q = p - 1
s = 0
while q % 2 == 0:
s += 1
q //= 2
z = 2
while pow(z, (p - 1) // 2, p) != p - 1:
z += 1
m = s
c = pow(z, q, p)
t = pow(n, q, p)
r = pow(n, (q + 1) // 2, p)
while t != 1:
i = 1
t2i = (t * t) % p
while i < m and t2i != 1:
t2i = (t2i * t2i) % p
i += 1
b = pow(c, 1 << (m - i - 1), p)
r = (r * b) % p
t = (t * b * b) % p
c = (b * b) % p
m = i
return r
def singular_ks(p, a=1):
# 4a^3 + 27k^2 ≡ 0 (mod p) => k^2 ≡ -4a^3 * inv(27)
inv27 = pow(27, -1, p)
rhs = (-4 * (a**3) * inv27) % p
r = tonelli_shanks(rhs, p)
if r is None:
return []
return sorted({r % p, (-r) % p})
def anomalous_ks_small_p(p, a=1):
base_count = [0] * p
for x in range(p):
v = (x * x % p * x + a * x) % p
base_count[v] += 1
chi = [0] * p
for r in range(p):
chi[r] = legendre(r, p)
hits = []
singular = []
for k in range(p):
if (4 * (a ** 3) + 27 * (k * k)) % p == 0:
singular.append(k)
continue
s = 0
for v, c in enumerate(base_count):
if c:
s += c * chi[(v + k) % p]
if s == -1: # order == p
hits.append(k)
return hits, singular
def anomalous_ks_for_y2_eq_x3_plus_ax_plus_k(p, a=1, small_p_limit=2_000_000):
if p <= small_p_limit:
return anomalous_ks_small_p(p, a=a)
sing = singular_ks(p, a=a)
return [], sing
Ps = (79, 967, 1303, 115792089237316195423570985008687907853269984665640564039457584007908834671663)
for p in Ps:
hits, sing = anomalous_ks_for_y2_eq_x3_plus_ax_plus_k(p, a=1)
print(f"p={p}\n order==p k'ları: {hits}\n singular: {sing}\n")
p=79
order==p k'ları: [14, 43, 53]
singular: [28, 51]
p=967
order==p k'ları: [451, 549, 587, 703, 705, 888]
singular: [259, 708]
p=1303
order==p k'ları: [35, 47, 142, 154, 240, 395, 503, 579, 754, 764, 841, 941, 1095, 1124, 1175, 1202, 1215, 1223]
singular: [332, 971]
p=115792089237316195423570985008687907853269984665640564039457584007908834671663
order==p k'ları: []
singular: [11842912595111486228085625214058158124919313799234703649739607102117954628595, 103949176642204709195485359794629749728350670866405860389717976905790880043068]
This code cannot be used to determine the b coefficient of large curves. I have developed a different and faster technique for this purpose; I am keeping the method and the b coefficients I found confidential at this stage, but I will share these values at the end of the process.
