prog.sf 2.9 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162
  1. #!/usr/bin/ruby
  2. # Square array A(n, k) read by antidiagonals downwards: smallest base-n strong Fermat pseudoprime with k distinct prime factors for k, n >= 2.
  3. # See also:
  4. # https://oeis.org/A271873
  5. #`(
  6. # PARI/GP program:
  7. strong_check(p, base, e, r) = my(tv=valuation(p-1, 2)); tv > e && Mod(base, p)^((p-1)>>(tv-e)) == r;
  8. strong_fermat_psp(A, B, k, base) = A=max(A, vecprod(primes(k))); (f(m, l, lo, k, e, r) = my(list=List()); my(hi=sqrtnint(B\m, k)); if(lo > hi, return(list)); if(k==1, forstep(p=lift(1/Mod(m, l)), hi, l, if(isprimepower(p) && gcd(m*base, p) == 1 && strong_check(p, base, e, r), my(n=m*p); if(n >= A && (n-1) % znorder(Mod(base, p)) == 0, listput(list, n)))), forprime(p=lo, hi, base%p == 0 && next; strong_check(p, base, e, r) || next; my(z=znorder(Mod(base, p))); gcd(m,z) == 1 || next; my(q=p, v=m*p); while(v <= B, list=concat(list, f(v, lcm(l, z), p+1, k-1, e, r)); q *= p; Mod(base, q)^z == 1 || break; v *= p))); list); my(res=f(1, 1, 2, k, 0, 1)); for(v=0, logint(B, 2), res=concat(res, f(1, 1, 2, k, v, -1))); vecsort(Set(res));
  9. T(n, k) = if(n < 2, return()); my(x=vecprod(primes(k)), y=2*x); while(1, my(v=strong_fermat_psp(x, y, k, n)); if(#v >= 1, return(v[1])); x=y+1; y=2*x);
  10. print_table(n, k) = for(x=2, n, for(y=2, k, print1(T(x, y), ", ")); print(""));
  11. for(k=2, 9, for(n=2, k, print1(T(n, k-n+2)", "))); \\ ~~~~
  12. )
  13. func T(n,k) { # OEIS: A271873
  14. if (n < 2) {
  15. return []
  16. }
  17. var x = 1
  18. var y = 2*x
  19. loop {
  20. var v = k.strong_fermat_psp(n,x,y)
  21. v.len >= 1 && return v[0]
  22. x = y+1
  23. y = 2*x
  24. }
  25. }
  26. say "\nTable:\n"
  27. {|x| {|y| T(x,y) }.map(2..6) }.map(2..6).each{.say}
  28. say "\nTerms:\n"
  29. {|k|
  30. {|n|
  31. print(T(n, k - n + 2), ", ")
  32. }.map(2..k)
  33. }.map(2..14).flat,
  34. say ''
  35. __END__
  36. Table:
  37. [2047, 15841, 800605, 293609485, 10761055201]
  38. [703, 8911, 152551, 41341321, 12283706701]
  39. [341, 4371, 129921, 9224391, 2592053871]
  40. [781, 24211, 4382191, 381347461, 9075517561]
  41. [217, 29341, 3405961, 557795161, 333515107081]
  42. Terms:
  43. 2047, 15841, 703, 800605, 8911, 341, 293609485, 152551, 4371, 781, 10761055201, 41341321, 129921, 24211, 217, 5478598723585, 12283706701, 9224391, 4382191, 29341, 325, 713808066913201, 1064404682551, 2592053871, 381347461, 3405961, 58825, 65, 90614118359482705, 19142275066201, 201068525791, 9075517561, 557795161, 1611805, 15841, 91, 5993318051893040401, 31475449738947061, 15804698567581, 2459465259031, 333515107081, 299048101, 205465, 1729, 91, 24325630440506854886701, 5405254780334022901, 8211666556476901, 1009614375233101, 17075314684351, 451698887431, 62338081, 63973, 1729, 133, 27146803388402594456683201, 649053318057640906411, 285574887575981641, 97383398181018301, 2805619369385521, 37388680793101, 10761055201, 20821801, 251251, 50997, 91, 4365221464536367089854499301, 512362775244582336338551, 350010251220384686701, 11290262679772872301, 364410024106892401, 5290810955703001, 1387708021141, 130027051, 67976623, 111055, 1729, 85,