数学研发论坛

 找回密码
 欢迎注册
查看: 216|回复: 19

[讨论] 素性测试中的平方数判定问题

[复制链接]
发表于 2019-10-8 09:01:13 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?欢迎注册

x
不论是BSPW中的Lucas测试,还是我发的Quadratic Frobenius测试,都要求测试整数n不能是完全平方数,否则用Jacobi symbol找不到参数

我想到个办法,就是用n mod m余数的方法,当选择特定m时候,余数是很少的
比如m = 120, 当n不含有2,3,5素因子时候,如果n是完全平方数,那么n模m的剩余只有1,49两个
这是暴力搜索m到65536得到的结果
其中length是满足与number互素的整数的平方剩余的个数,list里是具体的值

(length = 2, number = 120, list = Any[1, 49])
(length = 6, number = 420, list = Any[1, 109, 121, 169, 289, 361])
(length = 6, number = 840, list = Any[1, 121, 169, 289, 361, 529])
(length = 30, number = 4620, list = Any[1, 169, 289, 361, 421, 529, 709, 841, 949, 961, 1369, 1549, 1621, 1681, 1849, 2209, 2269, 2641, 2689, 2809, 2941, 3061, 3301, 3469, 3481, 3529, 3721, 4141, 4321, 4489])
(length = 30, number = 9240, list = Any[1, 169, 289, 361, 529, 841, 961, 1369, 1681, 1849, 2209, 2641, 2689, 2809, 3481, 3529, 3721, 4321, 4489, 5041, 5329, 5569, 6169, 6241, 6889, 7561, 7681, 7921, 8089, 8761])
(length = 180, number = 60060, list = Any[1, 289, 361, 529, 841, 961, 1369, 1621, 1681, 1849, 2209, 2809, 2941, 3301, 3481, 3721, 4489, 5041, 5149, 5329, 5461, 5581, 5989, 6241, 6301, 6829, 6889, 7309, 7921, 8089, 8761, 8941, 9109, 9409, 9949, 10189, 10201, 10609, 10789, 10921, 11449, 11509, 11881, 12301, 12541, 12769, 13381, 13729, 13861, 14221, 14569, 14821, 15409, 16069, 16129, 16501, 16669, 17161, 17341,
18001, 18769, 18841, 18901, 19009, 19189, 19321, 20029, 20101, 20329, 20749, 21121, 21421, 21961, 22201, 22621, 22801, 23269, 23461, 23521, 24049, 24469, 24649, 24781, 25741, 25789, 26569, 26581, 27421, 27589, 27889, 28081, 28141, 28669, 28681, 29929, 30781, 31021, 31201, 31249, 32041, 32341, 32509, 32629, 32761, 33049, 33289, 34189, 34609, 35149, 35281, 36061, 36481, 36661, 37129, 37249, 37489, 37801, 37909, 38509, 38581, 38809, 39601, 39649, 39901, 40429, 40681, 41869, 41941, 42949, 43261, 43429, 44269, 44521, 44641, 45049, 45061, 45109, 45301, 46069, 46201, 46489, 46621, 47161, 48049, 48409, 48889, 49009, 49141, 49261, 49501, 49669, 49729, 49921, 50521, 50821, 50989, 51349, 51529, 51661, 51769, 52441, 53089, 53509, 53629, 53881, 54289, 54349, 54961, 55441, 55969, 56281, 56809, 56989, 57061, 57121, 58081, 58249, 58501, 58969, 59929])
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2019-10-8 09:05:32 | 显示全部楼层
上面的数据里,有价值的是
(length = 2, number = 120, list = Any[1, 49]) 120=2^3*3*5
(length = 6, number = 840, list = Any[1, 121, 169, 289, 361, 529]) 840=2^3*3*5*7
最后一组,候选太多了,虽然数据比较好,待测试数除非很大,否则不建议选择
(length = 180, number = 60060, list = Any[1, 289, 361, 529, 841, 961, 1369, 1621, 1681, 1849, 2209, 2809, 2941, 3301, 3481, 3721, 4489, 5041, 5149, 5329, 5461, 5581, 5989, 6241, 6301, 6829, 6889, 7309, 7921, 8089, 8761, 8941, 9109, 9409, 9949, 10189, 10201, 10609, 10789, 10921, 11449, 11509, 11881, 12301, 12541, 12769, 13381, 13729, 13861, 14221, 14569, 14821, 15409, 16069, 16129, 16501, 16669, 17161, 17341,
18001, 18769, 18841, 18901, 19009, 19189, 19321, 20029, 20101, 20329, 20749, 21121, 21421, 21961, 22201, 22621, 22801, 23269, 23461, 23521, 24049, 24469, 24649, 24781, 25741, 25789, 26569, 26581, 27421, 27589, 27889, 28081, 28141, 28669, 28681, 29929, 30781, 31021, 31201, 31249, 32041, 32341, 32509, 32629, 32761, 33049, 33289, 34189, 34609, 35149, 35281, 36061, 36481, 36661, 37129, 37249, 37489, 37801, 37909, 38509, 38581, 38809, 39601, 39649, 39901, 40429, 40681, 41869, 41941, 42949, 43261, 43429, 44269, 44521, 44641, 45049, 45061, 45109, 45301, 46069, 46201, 46489, 46621, 47161, 48049, 48409, 48889, 49009, 49141, 49261, 49501, 49669, 49729, 49921, 50521, 50821, 50989, 51349, 51529, 51661, 51769, 52441, 53089, 53509, 53629, 53881, 54289, 54349, 54961, 55441, 55969, 56281, 56809, 56989, 57061, 57121, 58081, 58249, 58501, 58969, 59929])
60060=2^2*3*5*7*11*13
这组可以用
(length = 30, number = 9240, list = Any[1, 169, 289, 361, 529, 841, 961, 1369, 1681, 1849, 2209, 2641, 2689, 2809, 3481, 3529, 3721, 4321, 4489, 5041, 5329, 5569, 6169, 6241, 6889, 7561, 7681, 7921, 8089, 8761])  9240=2^3*3*5*7*11
代替,180/60060=0.002997,30/9240=0.0003247,差距很小
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2019-10-8 09:16:19 | 显示全部楼层
  1. global cnt = (length = 65536, number = 65537, list = [])
  2. println("Begin")
  3. for m in range(120, step = 2, stop = 65536)
  4.     global s = []
  5.     for i in range(1, stop = m-1)
  6.         if gcd(i, m) == 1
  7.             t = i * i % m
  8.             push!(s, t)
  9.         end
  10.     end
  11.     sort!(s)
  12.     unique!(s)
  13.     if cnt.length//cnt.number > length(s)//m
  14.         global cnt = (length = length(s), number = m, list = s)
  15.         println(cnt)
  16.     end
  17. end
  18. println("End")
  19.         
复制代码

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2019-10-8 10:05:22 | 显示全部楼层
求平方根是很快的
https://gmplib.org/manual/Integer-Roots.html
函数mpz_sqrt
得到平方根以后,在自相乘,再判断结果是否等于原整数即可
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2019-10-8 10:16:30 | 显示全部楼层
  1. function isSquareOrDivBy2357(n)
  2.     if gcd(n, 840) != 1
  3.         return true
  4.     end
  5.     t = n % 840
  6.     if t in [1, 121, 169, 289, 361, 529]
  7.         sq = isqrt(n)
  8.         return sq * sq == n
  9.     else
  10.         return false
  11.     end
  12. end

  13. for i in range(1, stop = 65535)
  14.     t = i * i
  15.     if !isSquareOrDivBy2357(t)
  16.         print("error at $i * $i = $t")
  17.     end
  18. end

  19. t = big(10)^67 - 1
  20. println("$t: ", isSquareOrDivBy2357(t))

  21. t = t * t
  22. println("$t: ", isSquareOrDivBy2357(t))
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2019-10-8 10:42:10 | 显示全部楼层
mathe 发表于 2019-10-8 10:05
求平方根是很快的
https://gmplib.org/manual/Integer-Roots.html
函数mpz_sqrt

我知道,但是我想先筛选掉绝大部分不是平方的数
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2019-10-8 10:51:07 | 显示全部楼层
  1. function isSquareOrDivBy2357(n)
  2.     if gcd(n, 840) != 1
  3.         return true
  4.     end
  5.     t = n % 840
  6.     if t in [1, 121, 169, 289, 361, 529]
  7.         sq = isqrt(n)
  8.         return sq * sq == n
  9.     else
  10.         return false
  11.     end
  12. end

  13. function isSquare1(n)
  14.     sq = isqrt(n)
  15.     return sq * sq == n
  16. end

  17. @time for i in range(1, stop = 65535)
  18.     t = i * i
  19.     if !isSquareOrDivBy2357(t)
  20.         print("error at $i * $i = $t")
  21.     end
  22. end

  23. @time for i in range(1, stop = 65535)
  24.     t = i * i
  25.     if !isSquare1(t)
  26.         print("error at $i * $i = $t")
  27.     end
  28. end

  29. t = t * t
  30. @time println("$t: ", isSquareOrDivBy2357(t))
  31. @time println("$t: ", isSquare1(t))

复制代码


好吧,结果果然不如直接开平方~
0.004386 seconds (14.98 k allocations: 1.828 MiB)
  0.000710 seconds
9999999999999999999999999999999999999999999999999999999999999999996000000000000000000000000000000000000000000000000000000000000000000599999999999999999999999999999999999999999999999999999999999999999960000000000000000000000000000000000000000000000000000000000000000001: true
  0.072349 seconds (197.96 k allocations: 10.581 MiB, 15.72% gc time)
9999999999999999999999999999999999999999999999999999999999999999996000000000000000000000000000000000000000000000000000000000000000000599999999999999999999999999999999999999999999999999999999999999999960000000000000000000000000000000000000000000000000000000000000000001: true
  0.018368 seconds (2.21 k allocations: 106.926 KiB)

点评

后来我想了想,这种测试执行时间的测试方法不对  发表于 2019-10-9 08:23
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2019-10-8 17:12:04 | 显示全部楼层
借题,问一下,可有此类算法:可快速判定一个大整数,有/无 含非 1 的平方因子?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2019-10-8 17:48:16 | 显示全部楼层
验证n是完全平方数
sqrt(n)求一下看看是不是整数,不就出结果了吗?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2019-10-8 18:44:44 | 显示全部楼层
gxqcn 发表于 2019-10-8 17:12
借题,问一下,可有此类算法:可快速判定一个大整数,有/无 含非 1 的平方因子?

3年前我问过,也再问一下。

https://bbs.emath.ac.cn/thread-9242-1-1.html
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2019-10-19 17:52 , Processed in 0.056931 second(s), 17 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表