Achilles numbers: Difference between revisions
(Added Algol 68) |
(Added XPL0 example.) |
||
Line 730: | Line 730: | ||
10 digits: 77330 |
10 digits: 77330 |
||
11 digits: 247449 |
11 digits: 247449 |
||
</pre> |
|||
=={{header|XPL0}}== |
|||
<lang XPL0>func GCD(N, D); \Return the greatest common divisor of N and D |
|||
int N, D; \numerator and denominator |
|||
int R; |
|||
[if D > N then |
|||
[R:= D; D:= N; N:= R]; \swap D and N |
|||
while D > 0 do |
|||
[R:= rem(N/D); |
|||
N:= D; |
|||
D:= R; |
|||
]; |
|||
return N; |
|||
]; \GCD |
|||
func Totient(N); \Return the totient of N |
|||
int N, Phi, M; |
|||
[Phi:= 0; |
|||
for M:= 1 to N do |
|||
if GCD(M, N) = 1 then Phi:= Phi+1; |
|||
return Phi; |
|||
]; |
|||
func Powerful(N0); \Return 'true' if N0 is a powerful number |
|||
int N0, N, F, Q, L; |
|||
[if N0 <= 1 then return false; |
|||
N:= N0; F:= 2; |
|||
L:= sqrt(N0); |
|||
loop [Q:= N/F; |
|||
if rem(0) = 0 then \found a factor |
|||
[if rem(N0/(F*F)) then return false; |
|||
N:= Q; |
|||
if F>N then quit; |
|||
] |
|||
else [F:= F+1; |
|||
if F > L then |
|||
[if rem(N0/(N*N)) then return false; |
|||
quit; |
|||
]; |
|||
]; |
|||
]; |
|||
return true; |
|||
]; |
|||
func Achilles(N); \Return 'true' if N is an Achilles number |
|||
int N, M, A; |
|||
[if not Powerful(N) then return false; |
|||
M:= 2; |
|||
A:= M*M; |
|||
repeat loop [if A = N then return false; |
|||
if A > N then quit; |
|||
A:= A*M; |
|||
]; |
|||
M:= M+1; |
|||
A:= M*M; |
|||
until A > N; |
|||
return true; |
|||
]; |
|||
int Cnt, N, Pwr, Start; |
|||
[Cnt:= 0; |
|||
N:= 1; |
|||
loop [if Achilles(N) then |
|||
[IntOut(0, N); |
|||
Cnt:= Cnt+1; |
|||
if Cnt >= 50 then quit; |
|||
if rem(Cnt/10) then ChOut(0, 9) else CrLf(0); |
|||
]; |
|||
N:= N+1; |
|||
]; |
|||
CrLf(0); CrLf(0); |
|||
Cnt:= 0; |
|||
N:= 1; |
|||
loop [if Achilles(N) then |
|||
if Achilles(Totient(N)) then |
|||
[IntOut(0, N); |
|||
Cnt:= Cnt+1; |
|||
if Cnt >= 20 then quit; |
|||
if rem(Cnt/10) then ChOut(0, 9) else CrLf(0); |
|||
]; |
|||
N:= N+1; |
|||
]; |
|||
CrLf(0); CrLf(0); |
|||
for Pwr:= 1 to 6 do |
|||
[IntOut(0, Pwr); Text(0, ": "); |
|||
Start:= fix(Pow(10.0, float(Pwr-1))); |
|||
Cnt:= 0; |
|||
for N:= Start to Start*10-1 do |
|||
if Achilles(N) then Cnt:= Cnt+1; |
|||
IntOut(0, Cnt); CrLf(0); |
|||
]; |
|||
]</lang> |
|||
{{out}} |
|||
<pre> |
|||
72 108 200 288 392 432 500 648 675 800 |
|||
864 968 972 1125 1152 1323 1352 1372 1568 1800 |
|||
1944 2000 2312 2592 2700 2888 3087 3200 3267 3456 |
|||
3528 3872 3888 4000 4232 4500 4563 4608 5000 5292 |
|||
5324 5400 5408 5488 6075 6125 6272 6728 6912 7200 |
|||
500 864 1944 2000 2592 3456 5000 10125 10368 12348 |
|||
12500 16875 19652 19773 30375 31104 32000 33275 37044 40500 |
|||
1: 0 |
|||
2: 1 |
|||
3: 12 |
|||
4: 47 |
|||
5: 192 |
|||
6: 664 |
|||
</pre> |
</pre> |
Revision as of 15:27, 28 February 2022
This page uses content from Wikipedia. The original article was at Achilles number. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance) |
An Achilles number is a number that is powerful but imperfect. Named after Achilles, a hero of the Trojan war, who was also powerful but imperfect.
A positive integer n is a powerful number if, for every prime factor p of n, p2 is also a divisor.
In other words, every prime factor appears at least squared in the factorization.
All Achilles numbers are powerful. However, not all powerful numbers are Achilles numbers: only those that cannot be represented as mk, where m and k are positive integers greater than 1.
A strong Achilles number is an Achilles number whose Euler totient (𝜑) is also an Achilles number.
- E.G.
108 is a powerful number. Its prime factorization is 22 × 33, and thus its prime factors are 2 and 3. Both 22 = 4 and 32 = 9 are divisors of 108. However, 108 cannot be represented as mk, where m and k are positive integers greater than 1, so 108 is an Achilles number.
360 is not an Achilles number because it is not powerful. One of its prime factors is 5 but 360 is not divisible by 52 = 25.
Finally, 784 is not an Achilles number. It is a powerful number, because not only are 2 and 7 its only prime factors, but also 22 = 4 and 72 = 49 are divisors of it. Nonetheless, it is a perfect power; its square root is an even integer, so it is not an Achilles number.
500 = 22 × 53 is a strong Achilles number as its Euler totient, 𝜑(500), is 200 = 23 × 52 which is also an Achilles number.
- Task
- Find and show the first 50 Achilles numbers.
- Find and show at least the first 20 strong Achilles numbers.
- For at least 2 through 5, show the count of Achilles numbers with that many digits.
- See also
ALGOL 68
<lang algol68>BEGIN # find Achiles Numbers: numbers whose prime factors p appear at least #
# twice (i.e. if p is a prime factor, so is p^2) and cannot be # # expressed as m^k for any integer m, k > 1 # # also find strong Achiles Numbers: Achiles Numbers where the Euler's # # totient of the number is also Achiles # # returns the number of integers k where 1 <= k <= n that are mutually # # prime to n # PROC totient = ( INT n )INT: # algorithm from the second Go sample # IF n < 3 THEN 1 # in the Totient Function task # ELIF n = 3 THEN 2 ELSE INT result := n; INT v := n; INT i := 2; WHILE i * i <= v DO IF v MOD i = 0 THEN WHILE v MOD i = 0 DO v OVERAB i OD; result -:= result OVER i FI; IF i = 2 THEN i := 1 FI; i +:= 2 OD; IF v > 1 THEN result -:= result OVER v FI; result FI # totient # ; # find the numbers # INT max number = 1 000 000; # max number we will consider # PR read "primes.incl.a68" PR # include prime utilities # []BOOL prime = PRIMESIEVE max number; # construct a sieve of primes # # table of numbers, will be set to TRUE for the Achiles Numbers # [ 1 : max number ]BOOL achiles; FOR a TO UPB achiles DO achiles[ a ] := TRUE OD; # remove the numbers that don't have squared primes as factors # achiles[ 1 ] := FALSE; FOR a TO UPB achiles DO IF prime[ a ] THEN # have a prime, remove it and every multiple of it that isn't a # # multiple of a squared # INT a count := 0; FOR j FROM a BY a TO UPB achiles DO a count +:= 1; IF a count = a THEN # have a multiple of i^2, keep the number # a count := 0 ELSE # not a multiple of i^2, remove the number # achiles[ j ] := FALSE FI OD FI OD; # achiles now has TRUE for the powerful numbers, remove all m^k (m,k > 1) # FOR m FROM 2 TO UPB achiles DO INT mk := m; INT max mk = UPB achiles OVER m; # avoid overflow if INT is 32 bit # WHILE mk <= max mk DO mk *:= m; achiles[ mk ] := FALSE OD OD; # achiles now has TRUE for imperfect powerful numbers # # show the first 50 Achiles Numbers # BEGIN print( ( "First 50 Achiles Numbers:", newline ) ); INT a count := 0; FOR a WHILE a count < 50 DO IF achiles[ a ] THEN a count +:= 1; print( ( " ", whole( a, -6 ) ) ); IF a count MOD 10 = 0 THEN print( ( newline ) ) FI FI OD END; # show the first 50 Strong Achiles numbers # BEGIN print( ( "First 20 Strong Achiles Numbers:", newline ) ); INT s count := 0; FOR s WHILE s count < 20 DO IF achiles[ s ] THEN IF achiles[ totient( s ) ] THEN s count +:= 1; print( ( " ", whole( s, -6 ) ) ); IF s count MOD 10 = 0 THEN print( ( newline ) ) FI FI FI OD END; # count the number of Achiles Numbers by their digit counts # BEGIN INT a count := 0; INT power of 10 := 100; INT digit count := 2; FOR a TO UPB achiles DO IF achiles[ a ] THEN # have an Achiles Number # a count +:= 1 FI; IF a = power of 10 THEN # have reached a power of 10 # print( ( "Achiles Numbers with ", whole( digit count, 0 ) , " digits: ", whole( a count, -6 ) , newline ) ); digit count +:= 1; power of 10 *:= 10; a count := 0 FI OD END
END</lang>
- Output:
First 50 Achiles Numbers: 72 108 200 288 392 432 500 648 675 800 864 968 972 1125 1152 1323 1352 1372 1568 1800 1944 2000 2312 2592 2700 2888 3087 3200 3267 3456 3528 3872 3888 4000 4232 4500 4563 4608 5000 5292 5324 5400 5408 5488 6075 6125 6272 6728 6912 7200 First 20 Strong Achiles Numbers: 500 864 1944 2000 2592 3456 5000 10125 10368 12348 12500 16875 19652 19773 30375 31104 32000 33275 37044 40500 Achiles Numbers with 2 digits: 1 Achiles Numbers with 3 digits: 12 Achiles Numbers with 4 digits: 47 Achiles Numbers with 5 digits: 192 Achiles Numbers with 6 digits: 66
J
Implementation:
<lang J>achilles=: (*/ .>&1 * 1 = +./)@(1{__&q:)"0 strong=: achilles@(5&p:)</lang>
Task examples:
<lang J> 5 10$(#~ achilles) 1+i.10000 NB. first 50 achilles numbers
72 108 200 288 392 432 500 648 675 800 864 968 972 1125 1152 1323 1352 1372 1568 1800
1944 2000 2312 2592 2700 2888 3087 3200 3267 3456 3528 3872 3888 4000 4232 4500 4563 4608 5000 5292 5324 5400 5408 5488 6075 6125 6272 6728 6912 7200
20{.(#~ strong * achilles) 1+i.100000 NB. first twenty strong achilles numbers
500 864 1944 2000 2592 3456 5000 10125 10368 12348 12500 16875 19652 19773 30375 31104 32000 33275 37044 40500
+/achilles (+i.)/1 9*10^<:2 NB. count of two digit achilles numbers
1
+/achilles (+i.)/1 9*10^<:3
12
+/achilles (+i.)/1 9*10^<:4
47
+/achilles (+i.)/1 9*10^<:5
192
+/achilles (+i.)/1 9*10^<:6
664</lang>
Explanation of the code:
(1{__&q:) is a function which returns the non-zero powers of the prime factors of a positive integer. (__&q: returns both the primes and their factors, but here we do not care about the primes themselves.)
+./ returns the greatest common divisor of a list, and 1=+./ is true if that gcd is 1 (0 if it's false).
*/ .>&1 is true if all the values in a list are greater than 1 (0 if not).
"0 maps a function onto the individual (rank 0) items of a list or array (we use this to avoid complexities: for example if we padded our lists of prime factor powers with zeros, we could still find the gcd, but our test that the powers are greater than 1 would fail).
5&p: is euler's totient function.
(#~ predicate) list selects the elements of list where predicate is true.
Julia
<lang>using Primes
isAchilles(n) = (p = [x[2] for x in factor(n).pe]; all(>(1), p) && gcd(p) == 1)
isstrongAchilles(n) = isAchilles(n) && isAchilles(totient(n))
function teststrongachilles(nachilles = 50, nstrongachilles = 100)
# task 1 println("First $nachilles Achilles numbers:") n, found = 0, 0 while found < nachilles if isAchilles(n) found += 1 print(rpad(n, 5), found % 10 == 0 ? "\n" : "") end n += 1 end # task 2 println("\nFirst $nstrongachilles strong Achilles numbers:") n, found = 0, 0 while found < nstrongachilles if isstrongAchilles(n) found += 1 print(rpad(n, 7), found % 10 == 0 ? "\n" : "") end n += 1 end # task 3 println("\nCount of Achilles numbers for various intervals:") intervals = [10:99, 100:999, 1000:9999, 10000:99999, 100000:999999] for interval in intervals println(lpad(interval, 15), " ", count(isAchilles, interval)) end
end
teststrongachilles()
</lang>
- Output:
First 50 Achilles numbers: 72 108 200 288 392 432 500 648 675 800 864 968 972 1125 1152 1323 1352 1372 1568 1800 1944 2000 2312 2592 2700 2888 3087 3200 3267 3456 3528 3872 3888 4000 4232 4500 4563 4608 5000 5292 5324 5400 5408 5488 6075 6125 6272 6728 6912 7200 First 100 strong Achilles numbers: 500 864 1944 2000 2592 3456 5000 10125 10368 12348 12500 16875 19652 19773 30375 31104 32000 33275 37044 40500 49392 50000 52488 55296 61731 64827 67500 69984 78608 80000 81000 83349 84375 93312 108000 111132 124416 128000 135000 148176 151875 158184 162000 165888 172872 177957 197568 200000 202612 209952 219488 221184 237276 243000 246924 253125 266200 270000 273375 296352 320000 324000 333396 364500 397953 405000 432000 444528 453789 455877 493848 497664 500000 518616 533871 540000 555579 583443 605052 607500 629856 632736 648000 663552 665500 666792 675000 691488 740772 750141 790272 800000 810448 820125 831875 877952 949104 972000 987696 1000188 Count of Achilles numbers for various intervals: 10:99 1 100:999 12 1000:9999 47 10000:99999 192 100000:999999 664
Phix
You can run this online here.
with javascript_semantics requires("1.0.2") -- [prime_factors(duplicates:=2) enhancement] function achilles(integer n, maxprime=-1) sequence f = vslice(prime_factors(n,2,maxprime),2) return min(f)>=2 and gcd(f)=1 end function function totient(integer n) integer tot = n, i = 2 while i*i<=n do if rmdr(n,i)=0 then while true do n /= i if rmdr(n,i)!=0 then exit end if end while tot -= tot/i end if i += iff(i=2?1:2) end while if n>1 then tot -= tot/n end if return tot end function function strong_achilles(integer n, maxprime=-1) return achilles(n,maxprime) and achilles(totient(n),maxprime) end function atom t0 = time() integer c = iff(platform()=JS?7:8), pn = 0 -- powerful numbers must be multiples of 2 or more squared primes, -- achilles numbers must be a multiple of at least 1 cubed prime sequence p2 = repeat(0,power(10,c)) while true do pn += 1 atom pn2 = power(get_prime(pn),2) if pn2>length(p2) then exit end if for i=pn2 to length(p2) by pn2 do p2[i] += 2 end for atom pn3 = power(get_prime(pn),3) if pn3<=length(p2) then -- (avoids a for loop limit error...) for i=pn3 to length(p2) by pn3 do p2[i] = or_bits(p2[i],1) end for end if end while sequence a = {}, sa = {}, counts = repeat(0,c) bool bPrintDigits = false for p=1 to c do integer p10 = power(10,p-1), p99 = power(10,p)-1, maxprime = get_maxprime(p99) for n=p10 to p99 do integer p2n = p2[n] if p2n>4 and odd(p2n) and achilles(n,maxprime) then if length(a)<50 then a = append(a,sprintf("%4d",n)) if length(a)=50 then printf(1,"First 50 Achilles numbers:\n%s\n",{join_by(a,1,10," ")}) end if end if if length(sa)<30 and strong_achilles(n,maxprime) then sa = append(sa,sprintf("%5d",n)) if length(sa)=30 then printf(1,"First 30 strong Achilles numbers:\n%s\n",{join_by(sa,1,10," ")}) bPrintDigits = true end if end if counts[p] += 1 end if end for if bPrintDigits then for i=1 to p do if counts[i]!=0 then printf(1,"Achilles numbers with %d digits:%d\n",{i,counts[i]}) counts[i] = 0 end if end for end if end for ?elapsed(time()-t0)
- Output:
First 50 Achilles numbers: 72 108 200 288 392 432 500 648 675 800 864 968 972 1125 1152 1323 1352 1372 1568 1800 1944 2000 2312 2592 2700 2888 3087 3200 3267 3456 3528 3872 3888 4000 4232 4500 4563 4608 5000 5292 5324 5400 5408 5488 6075 6125 6272 6728 6912 7200 First 30 strong Achilles numbers: 500 864 1944 2000 2592 3456 5000 10125 10368 12348 12500 16875 19652 19773 30375 31104 32000 33275 37044 40500 49392 50000 52488 55296 61731 64827 67500 69984 78608 80000 Achilles numbers with 2 digits:1 Achilles numbers with 3 digits:12 Achilles numbers with 4 digits:47 Achilles numbers with 5 digits:192 Achilles numbers with 6 digits:664 Achilles numbers with 7 digits:2242 Achilles numbers with 8 digits:7395 "21.2s"
Unfortunately that p2 sieve is a bit of a memory hog and it won't go to 9 digits (didn't bother testing but pretty sure JavaScript would choke on 8).
Raku
Timing is going to be system / OS dependent. <lang perl6>use Prime::Factor; use Math::Root;
sub is-square-free (Int \n) {
constant @p = ^100 .map: { next unless .is-prime; .² }; for @p -> \p { return False if n %% p } True
}
sub powerful (\n, \k = 2) {
my @p; p(1, 2*k - 1); sub p (\m, \r) { @p.push(m) and return if r < k; for 1 .. (n / m).&root(r) -> \v { if r > k { next unless is-square-free(v); next unless m gcd v == 1; } p(m * v ** r, r - 1) } } @p
}
my $achilles = powerful(10**9).hyper(:500batch).grep( {
my $f = .&prime-factors.Bag; (+$f.keys > 1) && (1 == [gcd] $f.values) && (.sqrt.Int² !== $_)
} ).classify: { .chars }
my \𝜑 = 0, |(1..*).hyper.map: -> \t { t × [×] t.&prime-factors.squish.map: { 1 - 1/$_ } }
my %as = Set.new: flat $achilles.values».list;
my $strong = lazy (flat $achilles.sort».value».list».sort).grep: { ?%as{𝜑[$_]} };
put "First 50 Achilles numbers:"; put (flat $achilles.sort».value».list».sort)[^50].batch(10)».fmt("%4d").join("\n");
put "\nFirst 30 strong Achilles numbers:"; put $strong[^30].batch(10)».fmt("%5d").join("\n");
put "\nNumber of Achilles numbers with:"; say "$_ digits: " ~ +$achilles{$_} // 0 for 2..9;
printf "\n%.1f total elapsed seconds\n", now - INIT now;</lang>
- Output:
First 50 Achilles numbers: 72 108 200 288 392 432 500 648 675 800 864 968 972 1125 1152 1323 1352 1372 1568 1800 1944 2000 2312 2592 2700 2888 3087 3200 3267 3456 3528 3872 3888 4000 4232 4500 4563 4608 5000 5292 5324 5400 5408 5488 6075 6125 6272 6728 6912 7200 First 30 strong Achilles numbers: 500 864 1944 2000 2592 3456 5000 10125 10368 12348 12500 16875 19652 19773 30375 31104 32000 33275 37044 40500 49392 50000 52488 55296 61731 64827 67500 69984 78608 80000 Number of Achilles numbers with: 2 digits: 1 3 digits: 12 4 digits: 47 5 digits: 192 6 digits: 664 7 digits: 2242 8 digits: 7395 9 digits: 24008 6.1 total elapsed seconds
Could go further but slows to a crawl and starts chewing up memory in short order.
10 digits: 77330 11 digits: 247449 12 digits: 788855 410.4 total elapsed seconds
Wren
Version 1 (Brute force)
This finds the number of 6 digit Achilles numbers in 2.5 seconds, 7 digits in 51 seconds but 8 digits needs a whopping 21 minutes! <lang ecmascript>import "./math" for Int import "./seq" for Lst import "./fmt" for Fmt
var maxDigits = 8 var limit = 10.pow(maxDigits) var c = Int.primeSieve(limit-1, false)
var totient = Fn.new { |n|
var tot = n var i = 2 while (i*i <= n) { if (n%i == 0) { while(n%i == 0) n = (n/i).floor tot = tot - (tot/i).floor } if (i == 2) i = 1 i = i + 2 } if (n > 1) tot = tot - (tot/n).floor return tot
}
var isPerfectPower = Fn.new { |n|
if (n == 1) return true var x = 2 while (x * x <= n) { var y = 2 var p = x.pow(y) while (p > 0 && p <= n) { if (p == n) return true y = y + 1 p = x.pow(y) } x = x + 1 } return false
}
var isPowerful = Fn.new { |n|
while (n % 2 == 0) { var p = 0 while (n % 2 == 0) { n = (n/2).floor p = p + 1 } if (p == 1) return false } var f = 3 while (f * f <= n) { var p = 0 while (n % f == 0) { n = (n/f).floor p = p + 1 } if (p == 1) return false f = f + 2 } return n == 1
}
var isAchilles = Fn.new { |n| c[n] && isPowerful.call(n) && !isPerfectPower.call(n) }
var isStrongAchilles = Fn.new { |n|
if (!isAchilles.call(n)) return false var tot = totient.call(n) return isAchilles.call(tot)
}
System.print("First 50 Achilles numbers:") var achilles = [] var count = 0 var n = 2 while (count < 50) {
if (isAchilles.call(n)) { achilles.add(n) count = count + 1 } n = n + 1
} for (chunk in Lst.chunks(achilles, 10)) Fmt.print("$4d", chunk)
System.print("\nFirst 30 strong Achilles numbers:") var strongAchilles = [] count = 0 n = achilles[0] while (count < 30) {
if (isStrongAchilles.call(n)) { strongAchilles.add(n) count = count + 1 } n = n + 1
} for (chunk in Lst.chunks(strongAchilles, 10)) Fmt.print("$5d", chunk)
System.print("\nNumber of Achilles numbers with:") var pow = 10 for (i in 2..maxDigits) {
var count = 0 for (j in pow..pow*10-1) { if (isAchilles.call(j)) count = count + 1 } System.print("%(i) digits: %(count)") pow = pow * 10
}</lang>
- Output:
First 50 Achilles numbers: 72 108 200 288 392 432 500 648 675 800 864 968 972 1125 1152 1323 1352 1372 1568 1800 1944 2000 2312 2592 2700 2888 3087 3200 3267 3456 3528 3872 3888 4000 4232 4500 4563 4608 5000 5292 5324 5400 5408 5488 6075 6125 6272 6728 6912 7200 First 30 strong Achilles numbers: 500 864 1944 2000 2592 3456 5000 10125 10368 12348 12500 16875 19652 19773 30375 31104 32000 33275 37044 40500 49392 50000 52488 55296 61731 64827 67500 69984 78608 80000 Number of Achilles numbers with: 2 digits: 1 3 digits: 12 4 digits: 47 5 digits: 192 6 digits: 664 7 digits: 2242 8 digits: 7395
Version 2 (Much faster)
Here we make use of the fact that powerful numbers are always of the form a²b³, where a and b > 0, to generate such numbers up to a given limit. We also use an improved method for checking whether a number is a perfect power.
Ridiculously fast compared to the previous method: 8 digits now reached in 0.57 seconds, 9 digits in 4.2 seconds and 10 digits in 34 seconds. However, 11 digits takes 300 seconds and I gave up after that. <lang ecmascript>import "./set" for Set, Bag import "./math" for Int import "./seq" for Lst import "./fmt" for Fmt
var totient = Fn.new { |n|
var tot = n var i = 2 while (i*i <= n) { if (n%i == 0) { while(n%i == 0) n = (n/i).floor tot = tot - (tot/i).floor } if (i == 2) i = 1 i = i + 2 } if (n > 1) tot = tot - (tot/n).floor return tot
}
var isPerfectPower = Fn.new { |n|
if (n == 1) return true var pf = Int.primeFactors(n) var bag = Bag.new(pf) var es = bag.map { |me| me.value } if (es.count == 1) return true return Int.gcd(es.toList) != 1
}
var getAchilles = Fn.new { |minExp, maxExp|
var lower = 10.pow(minExp) var upper = 10.pow(maxExp) var achilles = Set.new() // avoids duplicates var count = 0 for (b in 1..upper.cbrt.floor) { var b3 = b * b * b for (a in 1..upper.sqrt.floor) { var p = b3 * a * a if (p >= upper) break if (p >= lower) { if (!isPerfectPower.call(p)) achilles.add(p) } } } return achilles
}
var achillesSet = getAchilles.call(1, 5) // enough for first 2 parts var achilles = achillesSet.toList achilles.sort()
System.print("First 50 Achilles numbers:") for (chunk in Lst.chunks(achilles[0..49], 10)) Fmt.print("$4d", chunk)
System.print("\nFirst 30 strong Achilles numbers:") var strongAchilles = [] var count = 0 var n = 0 while (count < 30) {
var tot = totient.call(achilles[n]) if (achillesSet.contains(tot)) { strongAchilles.add(achilles[n]) count = count + 1 } n = n + 1
} for (chunk in Lst.chunks(strongAchilles, 10)) Fmt.print("$5d", chunk)
var maxDigits = 11 System.print("\nNumber of Achilles numbers with:") for (d in 2..maxDigits) {
var ac = getAchilles.call(d-1, d).count Fmt.print("$2d digits: $d", d, ac)
}</lang>
- Output:
First 50 Achilles numbers: 72 108 200 288 392 432 500 648 675 800 864 968 972 1125 1152 1323 1352 1372 1568 1800 1944 2000 2312 2592 2700 2888 3087 3200 3267 3456 3528 3872 3888 4000 4232 4500 4563 4608 5000 5292 5324 5400 5408 5488 6075 6125 6272 6728 6912 7200 First 30 strong Achilles numbers: 500 864 1944 2000 2592 3456 5000 10125 10368 12348 12500 16875 19652 19773 30375 31104 32000 33275 37044 40500 49392 50000 52488 55296 61731 64827 67500 69984 78608 80000 Number of Achilles numbers with: 2 digits: 1 3 digits: 12 4 digits: 47 5 digits: 192 6 digits: 664 7 digits: 2242 8 digits: 7395 9 digits: 24008 10 digits: 77330 11 digits: 247449
XPL0
<lang XPL0>func GCD(N, D); \Return the greatest common divisor of N and D int N, D; \numerator and denominator int R; [if D > N then
[R:= D; D:= N; N:= R]; \swap D and N
while D > 0 do
[R:= rem(N/D); N:= D; D:= R; ];
return N; ]; \GCD
func Totient(N); \Return the totient of N int N, Phi, M; [Phi:= 0; for M:= 1 to N do
if GCD(M, N) = 1 then Phi:= Phi+1;
return Phi; ];
func Powerful(N0); \Return 'true' if N0 is a powerful number int N0, N, F, Q, L; [if N0 <= 1 then return false; N:= N0; F:= 2; L:= sqrt(N0); loop [Q:= N/F;
if rem(0) = 0 then \found a factor [if rem(N0/(F*F)) then return false; N:= Q; if F>N then quit; ] else [F:= F+1; if F > L then [if rem(N0/(N*N)) then return false; quit; ]; ]; ];
return true; ];
func Achilles(N); \Return 'true' if N is an Achilles number int N, M, A; [if not Powerful(N) then return false; M:= 2; A:= M*M; repeat loop [if A = N then return false;
if A > N then quit; A:= A*M; ]; M:= M+1; A:= M*M;
until A > N; return true; ];
int Cnt, N, Pwr, Start; [Cnt:= 0; N:= 1; loop [if Achilles(N) then
[IntOut(0, N); Cnt:= Cnt+1; if Cnt >= 50 then quit; if rem(Cnt/10) then ChOut(0, 9) else CrLf(0); ]; N:= N+1; ];
CrLf(0); CrLf(0); Cnt:= 0; N:= 1; loop [if Achilles(N) then
if Achilles(Totient(N)) then [IntOut(0, N); Cnt:= Cnt+1; if Cnt >= 20 then quit; if rem(Cnt/10) then ChOut(0, 9) else CrLf(0); ]; N:= N+1; ];
CrLf(0); CrLf(0); for Pwr:= 1 to 6 do
[IntOut(0, Pwr); Text(0, ": "); Start:= fix(Pow(10.0, float(Pwr-1))); Cnt:= 0; for N:= Start to Start*10-1 do if Achilles(N) then Cnt:= Cnt+1; IntOut(0, Cnt); CrLf(0); ];
]</lang>
- Output:
72 108 200 288 392 432 500 648 675 800 864 968 972 1125 1152 1323 1352 1372 1568 1800 1944 2000 2312 2592 2700 2888 3087 3200 3267 3456 3528 3872 3888 4000 4232 4500 4563 4608 5000 5292 5324 5400 5408 5488 6075 6125 6272 6728 6912 7200 500 864 1944 2000 2592 3456 5000 10125 10368 12348 12500 16875 19652 19773 30375 31104 32000 33275 37044 40500 1: 0 2: 1 3: 12 4: 47 5: 192 6: 664