Advertisement
Guest User

prime power sums

a guest
Jun 6th, 2020
28
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. $ cat PrimePowerSums.icn
  2. procedure main(args)
  3.  
  4. if *args == 0
  5. then args := [ 5000 ]
  6. else args[1] := integer(args[1]) |
  7. stop("Expecting an integer argument, not \"" || args[1] || "\".")
  8.  
  9. every write(prime_power_sums(args[1]))
  10.  
  11. write(&errout, "max: ", args[1], ", time: ", time_fmt(&time))
  12.  
  13. end # main
  14.  
  15.  
  16. procedure partition(lst, n)
  17.  
  18. # Return the index i into the given ordered list such that
  19. # lst[1:i] < n <= lst[i+1:].
  20.  
  21. local l, r, m
  22.  
  23. l := 1
  24. r := *lst + 1
  25.  
  26. # lst[1:l-1] < n and n <= lst[r:n]
  27.  
  28. while l < r do {
  29. m := (l + r)/2
  30. if lst[m] < n
  31. then l := m + 1
  32. else r := m
  33. }
  34.  
  35. return l - 1
  36.  
  37. end # partition
  38.  
  39.  
  40. procedure prime_power_sums(max)
  41.  
  42. # Generate all the prime power sums less then the given value. A prime power
  43. # sum is a value v with the property that p = p1^2 + p2^3 + p3^4, where p1,
  44. # p2, and p3 are primes.
  45.  
  46. local p2, p3, p4, p2s, p3c, p4q
  47. local pps, ppps
  48. local filter
  49. local max3, max2
  50.  
  51. filter := table(0)
  52.  
  53. every p4 := primes_at_most(max) do {
  54.  
  55. p4q := p4*p4*p4*p4
  56. max3 := max - p4q
  57. if max3 < 4 then break
  58.  
  59. every p3 := primes_at_most(max3) do {
  60.  
  61. p3c := p3*p3*p3
  62. ppps := p4q + p3c
  63. max2 := max - ppps
  64. if max2 < 4 then break
  65.  
  66. every p2 := primes_at_most(max2) do {
  67.  
  68. p2s := p2*p2
  69. if p2s > max2 then break
  70.  
  71. pps := ppps + p2s
  72. if filter[pps] == 0 then {
  73. filter[pps] := 1
  74. suspend pps
  75. }
  76. }
  77. }
  78. }
  79.  
  80. end # prime_power_sums
  81.  
  82.  
  83. procedure primes_at_most(max)
  84.  
  85. # Generate the sequence of primes (in increasing order) that are at most the
  86. # given value.
  87.  
  88. local i
  89.  
  90. static primes, primes_upper_bound
  91. initial {
  92. primes := []
  93. primes_upper_bound := 0
  94. }
  95.  
  96. if primes_upper_bound < max then {
  97. primes := primes_upto(max)
  98. primes_upper_bound := max
  99. }
  100.  
  101. every suspend !primes[1:partition(primes, max)]
  102.  
  103. end # primes_at_most
  104.  
  105.  
  106. procedure primes_upto(max)
  107.  
  108. # Return a list of primes (in ascending order) no greater than the given value.
  109.  
  110. local sieve, i, primes
  111.  
  112. sieve := list(max, 1)
  113.  
  114. every i := 2 to max do
  115. if sieve[i] == 1 then
  116. every sieve[i*2 to max by i] := 0
  117.  
  118. primes := []
  119. every i := 2 to max do
  120. if sieve[i] == 1 then
  121. put(primes, i)
  122.  
  123. return primes
  124.  
  125. end # primes_upto
  126.  
  127.  
  128. procedure time_fmt(t)
  129.  
  130. # Return the given time in ms as a conviently formatted string.
  131.  
  132. if t < 1000 then
  133. return string(t) || " ms"
  134.  
  135. t := integer(t/100.0 + 0.5)/10.0
  136.  
  137. if t < 60 then
  138. return string(t) || " sec"
  139.  
  140. t := integer((t/60.0)*10 + 0.5)/10.0
  141.  
  142. return string(t) || " min"
  143.  
  144. end # time_fmt
  145.  
  146.  
  147. $ icont PrimePowerSums.icn
  148. Translating:
  149. PrimePowerSums.icn:
  150. main
  151. partition
  152. prime_power_sums
  153. primes_at_most
  154. primes_upto
  155. time_fmt
  156. No errors
  157. Linking:
  158.  
  159. $ PrimePowerSums 50
  160. 28
  161. 33
  162. 49
  163. 47
  164. max: 50, time: 0 ms
  165.  
  166. $ for i in 2 3 4 5 6 ; do PrimePowerSums 5e$i | wc -l ; done
  167. max: 500, time: 0 ms
  168. 53
  169. max: 5000, time: 10 ms
  170. 395
  171. max: 50000, time: 110 ms
  172. 2579
  173. max: 500000, time: 1.2 sec
  174. 18899
  175. max: 5000000, time: 13.4 sec
  176. 138932
  177.  
  178. $ PrimePowerSums 5e7 > /dev/null
  179. max: 50000000, time: 3.1 min
  180.  
  181. $
Advertisement
Advertisement
Advertisement
RAW Paste Data Copied
Advertisement