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
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. \$