起因

正规子群的群友发了一道二分, 如下, 让转为 SageMath Code

import math
from gmpy2 import next_prime, sin, cos, mpz

left = 2**63
right = 2**64
n = mpz(
    13226258685595630160447827442153267599814292160977306459359416297020102688049793771120158029953432640539514974565723831126909317859357612109049397896863392224521508922404189962829783829182814282953169
)
while left <= right:
    k = (left + right) // 2
    p = (
        2023 * k**5
        + 2022 * k**4
        - 2023 * k**3
        - 2021 * k**2
        + 2020 * k
        + 2019
        + int(k**3 * sin(k + 2023))
        - int(k**3 * cos(k**2 - 2023))
        + int(math.e ** (23) * k)
    )
    q = (
        2023 * k**5
        - 2022 * k**4
        + 2023 * k**3
        - 2021 * k**2
        + 2020 * k
        - 2019
        - int(k**2 * sin(k * 2023))
        + int(k**2 * cos(k**2 * 2023))
        + int(math.e ** (23) * k**2)
    )
    p = next_prime(p)
    q = next_prime(q)
    if p * q < n:
        left = k + 1
    elif p * q > n:
        right = k - 1
    else:
        # time < 0.1s
        print("k =", k)  # k = 17821511070084697841
        break

简单替换导入函数后发现并不符合预期

sage

import math
# from gmpy2 import next_prime, sin, cos, mpz
from sage.all import next_prime, sin, cos

left = 2**63
right = 2**64
n = (
    13226258685595630160447827442153267599814292160977306459359416297020102688049793771120158029953432640539514974565723831126909317859357612109049397896863392224521508922404189962829783829182814282953169
)

while left <= right:
    k = (left + right) // 2
    p = (
        2023 * k**5
        + 2022 * k**4
        - 2023 * k**3
        - 2021 * k**2
        + 2020 * k
        + 2019
        + int((k**3 * sin(k + 2023)))
        - int((k**3 * cos(k**2 - 2023)))
        + int((math.e ** (23) * k))
    )
    q = (
        2023 * k**5
        - 2022 * k**4
        + 2023 * k**3
        - 2021 * k**2
        + 2020 * k
        - 2019
        - int((k**2 * sin(k * 2023)))
        + int((k**2 * cos(k**2 * 2023)))
        + int((math.e ** (23) * k**2))
    )
    p = next_prime(p)
    q = next_prime(q)
    if p * q < n:
        left = k + 1
    elif p * q > n:
        right = k - 1
    else:
        print("k =", k)
        break

其并没有输出, 也就是程序并没有匹配到 else, 手动输出 $k$ 发现是对的, 但这只能说明误差不是很大.

单独输出其中一个部分发现决定符号的部分出现极大偏差

sage

from sage.all import sin

k = 13835058055282163712
sin_k = sin(k + 2023)
print(int(sin_k))  # sage.all.sin

output
0

python

from gmpy2 import sin

k = 13835058055282163712
sin_k = sin(k + 2023)
print(int(sin_k))  # gmpy2.sin()

output
-1

决定符号的值出现误差, 影响很严重

疑问

SageMath 的默认精度是多少?

在官方文档我们可以找到答案

sage

RR(sys.maxsize)

RealField()
# Real Field with 53 bits of precision

gmpy2 - python

from gmpy2 import next_prime, sin, cos, mpz, mpfr  # type: ignore
from gmpy2 import get_context  # type: ignore

context = get_context()  # https://gmpy2.readthedocs.io/en/latest/contexts.html
print(context.precision)  # 53 bits by default | IEEE 754 double precision

# from gmpy2 import context
# context().precision = int(100)  # Set precision to 100 bits

尝试

sage - 10.4 ARM64

from sage.all import sin, RealField

R64 = RealField(64)

k = 13835058055282163712
sin_k = sin(k + 2023)
print(f"sin_k (sage): {sin_k}")
print(f"sin_k (python 53-bit float): {float(sin_k)}")
print(f"sin_k (RR 53-bit): {RR(sin_k)}")
print(f"sin_k (R64 64-bit): {R64(sin_k)}")

output
sin_k (sage): sin(13835058055282165735)
sin_k (python 53-bit float): -0.9007567878871614
sin_k (RR 53-bit): -0.900756787887161
sin_k (R64 64-bit): -0.835349138937036946

结果

没有结果, 我也不知道为什么精度一样但结果不同, 也许是 sage.all.sin() 与 gmpy2.sin() 的区别

等一个看完实现的小伙伴