纯lsp用Miller-Rabin算法求千万亿级别（16位数字）的质数，速度为毫秒级

tryhi

注：该算法核心源于网络
(try-isPrime 4330632126147343 8)
4330632126147343这是用这个函数求出来目前最大的质数，你们可以试试传统算法，lsp验证一个这样级别的质数需要多少时间，反正我的电脑估计一个星期都可能跑不出来，几乎达到双精度的极限，感觉再往上就要用字符串来模拟数值计算了
;;快速乘取模
;;快速求出（a*b)mod m 的值
(defun quickMulMod(a b m / ret)
        (setq ret 0)
        (while (not(zerop b))
                (if (not(zerop(rem b 2)))
                        (setq ret(rem(+ a ret)m))
                )
                (setq
                        b(fix(/ b 2))
                        a(rem(* a 2.)m)
                )
        )
        (fix ret)
)

;;快速幂取模
;;快速求出（a^b)mod m 的值
(defun quickPowMod (a b m / ret)
        (setq ret 1)
        (while (not(zerop b))
                (if (not(zerop(rem b 2)))
                        (setq ret(quickMulMod ret a m))
                )
                (setq
                        b(fix(/ b 2))
                        a(quickMulMod a a m)
                )
        )
        ret
)

;;判断一个整数是否为质数
;;lsp中最大的整数为2147483647，超过则自动转为实数
;;本函数求出来的最大质数为4330632126147343，几乎达到双精度的极限
;;参数为一个整数，第二个参数为猜测次数，建议用8
(defun try-isPrime (n tt / a d i j k k2 r ret x)
        (setq tt (min (- n 3) tt))
        (if (or(> n 4330632126147343)(< n 2))(princ"\必须输入3-4330632126147343的整数"))
        (if (= n 2)
                T
                (progn
                        (setq d(1- n)
                                r 0
                        )
                        (while (zerop(rem d 2))
                                (setq r(1+ r)d(/ d 2))
                        )
                        (setq i 0 k t ret t)
                        (while (and k(<= (setq i(1+ i)) tt))
                                (setq a(try_Ran-fix 2(- n 2))
                                        x(quickPowMod a d n)
                                )
                                (if(not(or(= x 1)(= x (1- n))))
                                        (progn
                                                (setq j 0 k2 t)
                                                (while (and k2(<= (setq j(1+ j)) r))
                                                        (setq x(quickMulMod x x n))
                                                        (if (= x(1- n))
                                                                (setq k2 nil)
                                                        )
                                                )
                                                (if k2
                                                        (setq k nil ret nil )
                                                )
                                        )
                                       
                                )
                        )
                        ret
                )
        )
)

(try-isPrime 4330632126147343 8)

tigcat

用cad2010，加载函数，控制台输入后提示“必须输入3-4330632126147343的整数nil”，海哥给的那个数用这个函数是不是不能验证啊？这个程序只能验证，我觉得加个功能，改成a0-b0，然后princ “符合的质数”
把海哥函数套了下，写了这个代码
(defun c:qzs (/ a0 b0 lst pp)
  (setq a0 (getreal "\n请输入起始正整数(不能<2):"))
  (setq b0 (getreal "\n请输入结束正整数(有限值，但相信你不会这么输这么大:"))
  (setq lst nil)
  (while (<= a0 b0)
    (if        (try-isprime a0 8)
      (setq lst (cons a0 lst))
    )
    (setq a0 (1+ a0))
  )
  (setq lst (reverse lst))
  (if  lst
    (progn
      (princ "\n符合要求的数如下(正确率99.99%,程序是靠猜测的):")
      (foreach pp lst
        (princ (strcat "\n" (rtos pp 2 0)))
      )
    )
    (princ "\n未发现符合要求的数(正确率99.99%,程序是靠猜测的)")
  )
  (prin1)
)
(prompt "\n***********<c:qzs>****求解质数")
(prin1)

tryhi

mahuan1279 发表于 2021-6-24 09:56
参数tt是起啥作用？  (setq i 0 k t ret t)似乎不对，应该是 (setq i 0 k tt ret tt),，不然k和ret变量值为 ...

怎么可能不对，本来就是要是逻辑值为T，k跟ret设计就是一个布尔值，你赋值为tt想干嘛（虽然你改了也能跑，但属于瞎改）？？
其实这个是Miller-Rabin算法，他验证一个质数并不是100%的，当tt取值为8时，正确率为，99.99999%，即1-4^(-8)

highflybird

Miller-Rabin素数测试算法介绍：
https://blog.csdn.net/forever_dreams/article/details/82314237
下面的文章也有参考价值：
https://www.cnblogs.com/zwfymqz/p/8150969.html
Miller Rabin质数判定算法_wang_yue_xi的博客-CSDN博客
https://blog.csdn.net/weixin_45321190/article/details/105100930

tryhi

(defun try_Ran-fix(n m)(fix(+ n (rem (getvar "CPUTICKS")(- m n -1)))))

补一个函数

cghdy

海哥优秀

mahuan1279

参数tt是起啥作用？  (setq i 0 k t ret t)似乎不对，应该是 (setq i 0 k tt ret tt),，不然k和ret变量值为逻辑值T了。

mahuan1279

tryhi 发表于 2021-6-24 10:14
怎么可能不对，本来就是要是逻辑值为T，k跟ret设计就是一个布尔值，你赋值为tt想干嘛（虽然你改了也能跑 ...

看到你前面自定义函数里设置ret为数值变量，还以为是同一个变量。

mahuan1279

已经验证100万以内没有反例。

tigcat

tigcat 发表于 2021-6-26 16:11
用cad2010，加载函数，控制台输入后提示“必须输入3-4330632126147343的整数nil”，海哥给的那个数用这个函 ...

是的，没看清，不过还是没明白的是，超出范围的整数是要怎么用这个函数验证是不是质数呢？