12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788 |
- #!/usr/bin/ruby
- # Daniel "Trizen" Șuteu
- # Date: 06 September 2022
- # https://github.com/trizen
- # Generate all the squarefree Lucas pseudoprimes to the U_n(P,Q) sequence with n prime factors in a given range [a,b]. (not in sorted order)
- # See also:
- # https://en.wikipedia.org/wiki/Almost_prime
- # https://en.wikipedia.org/wiki/Lucas_sequence
- # https://trizenx.blogspot.com/2020/08/pseudoprimes-construction-methods-and.html
- func lucas_znorder(P,Q,D,n) {
- n - kronecker(D, n) -> divisors.first_by {|d| lucasUmod(P, Q, d, n) == 0 }
- }
- func squarefree_lucas_U_pseudoprimes_in_range(a, b, k, P, Q, callback) {
- a = max(k.pn_primorial, a)
- var D = (P*P - 4*Q)
- func (m, L, lo, k) {
- var hi = idiv(b,m).iroot(k)
- return nil if (lo > hi)
- if (k == 1) {
- lo = max(lo, idiv_ceil(a, m))
- lo > hi && return nil
- for j in (1, -1) {
- var t = mulmod(m.invmod(L), j, L)
- t > hi && next
- t += L*idiv_ceil(lo - t, L) if (t < lo)
- t > hi && next
- for p in (range(t, hi, L)) {
- p.is_prime || next
- with (m*p) {|n|
- with (n - kronecker(D, n)) {|w|
- if ((L `divides` w) && (lucas_znorder(P, Q, D, p) `divides` w)) {
- callback(n)
- }
- }
- }
- }
- }
- return nil
- }
- each_prime(lo, hi, {|p|
- p.divides(D) && next
- var z = lucas_znorder(P, Q, D, p)
- m.is_coprime(z) || next
- __FUNC__(m*p, lcm(L, z), p+1, k-1)
- })
- }(1, 1, 2, k)
- return callback
- }
- # Generate all the squarefree Fibonacci pseudoprimes in the range [1, 15251]
- var from = 1
- var upto = 15251
- var P = 1
- var Q = -1
- var arr = []
- for k in (2..100) {
- break if k.pn_primorial>upto
- squarefree_lucas_U_pseudoprimes_in_range(from, upto, k, P, Q, { arr << _ })
- }
- say arr.sort
- __END__
- [323, 377, 1891, 3827, 4181, 5777, 6601, 6721, 8149, 10877, 11663, 13201, 13981, 15251]
|