Euler Problem 124 asks:
The radical of n, rad(n), is the product of distinct prime factors of n. For example, 504 = 23× 32 × 7, so rad(504) = 2 × 3 × 7 = 42.
If we calculate rad(n) for 1 <= n <= 10, then sort them on rad(n), and sorting on n if the radical values are equal, we get:
Let E(k) be the kth element in the sorted n column; for example, E(4) = 8 and E(6) = 9.
If rad(n) is sorted for 1 <= n <= 100000, find E(10000).
This is my 100th solution. Harumpf. All that, and I am yet but a Euler novice!
If you play with the rad() function a bit, a few things become clear: for prime p, rad(p) = p, and for all p, rad(pn) = rad(p). With those rules, you can compute a lot of rad() functions quickly. However, I couldn’t compute the remainder of the 100,000 radicals via the method shown above in Euler time (one minuteIf some one can, I’m ready to copy). So I cheated. Here is information about the radical function from “The Online Encyclopedia of Integer Sequences” or OEIS. On that page is a link to the first 100,000 computed radicals in textual format. Perfect! I downloaded that file and saved it as radn.txt. From there, bring it into E(1 to 100000, 1 to 2), sort E on column 2, and then bubble sort around E(10000,2) to get things in proper order at the point of interest. To sort E() I used blogmeister and fellow Eulerian Doug Jenkins’ SortV function. Get that here. It worked fine (Doug has since improved it). The whole routine, including the sort of 100,000 items, took a second and a half.
Here is the code that does it.
Dim E(1 To 100000, 1 To 2) As Long, Line As String
Dim TEMP(1 To 1, 1 To 2) As Long, i As Long, j As Long
Const Text As String = “D:DownloadsEuler
adn.txt”
Dim LBnd As Long, UBnd As Long
Dim Answer As Long, T As Single
T = Timer
i = 1
Open Text For Input As #1
Do While Not EOF(1)
Line Input #1, Line
j = InStr(1, Line, Chr(32)) ‘ space delimited
E(i, 1) = CLng(Left(Line, j – 1))
E(i, 2) = CLng(Right(Line, Len(Line) – j))
i = i + 1
Loop
Close #1
SortV SortRange:=E, SortBy:=2 ‘ Doug Jenkins’ function
i = 0
Do
i = i + 1 ‘ find the start of region of interest
Loop Until E(i, 2) = E(10000, 2)
LBnd = i ‘ lower bound
Do
i = i + 1 ‘ find the end of region of interest
Loop Until E(i, 2) != E(10000, 2)
UBnd = i – 1 ‘ upper bound 1 earlier
For i = LBnd To UBnd – 1 ‘ Bubble sort
For j = i + 1 To UBnd
If E(i, 1) > E(j, 1) Then
TEMP(1, 1) = E(j, 1)
E(j, 1) = E(i, 1)
E(i, 1) = TEMP(1, 1)
End If
Next j
Next i
Answer = E(10000, 1)
Debug.Print Answer; ” Time:”; Timer – T
End Sub
E(10000,2) turns out to be a pretty good year…one year after I was born. The usual angle bracket adjustments are used above.
I didn’t want to “cheat” quite that way. What I wanted to do was write code that went out on the web and pulled in the results. I’m row-challenged at 65536, and a Web-query seems to demand a range, not an array. I also couldn’t noodle out how to get past the row limit. I thought Excel might do the query via the clipboard (it does!) but even the clipboard object was truncated at 65,536 lines. I couldn’t come up with a way to populate 100,000×2 elements of an array via the internet, even though VBA doesn’t care. Is there one?
…mrt
I can’t help but think that a faster approach would be to /not/ generate the whole list, then sort; but instead to think about what the last item will be.
1) It must be greater or equal to the highest prime under 100000;
2) just check the values above that for prime factorisations.
Obviously finding that prime would be most efficient starting at 100000 and working backward. My online research shows this would be a /very/ short list of numbers to test. ;)
Hi Rick –
Not sure how that would work. There are, for instance, 74 numbers below 100,000 with rad(p) = 6. The first is 6 itself, and the last is 98,304. You’ve got to get them all next to each other (plus all the others next to their compadres) to assure E(10000,2) is the right value. Then, there are 6 numbers to sort to get to the proper E(10000,1)
Since the answer is in E(10000) I dumped the first 65536 items of the sorted array into a spreadsheet.
…mrt
Hi Michael,
I’m glad you posted another problem. I’m don’t know much about web queries, so my suggestion may not be what you’re looking for. But if you want to do the whole thing through VBA, Chip Pearson has a nice page on downloading a file from the internet at http://www.cpearson.com/Excel/DownloadFile.aspx. This way you can get the file and process the data all at once. I don’t know how to get around the row limit, other than by upgrading to 2007. The web query pasted all 100,000 values in for me, though it did take several seconds to import all of the data.
As far as computing a lot of rad() functions quickly, I have a subroutine that will take in a number and spit out an array containing the prime factorization of that number. First, I get a list of prime numbers less than the square root of 100,000 (I actually use a number slightly higher than 316, just to be on the safe side). I use the Sieve of Eratosthenes and these numbers get placed in the Global array ArrayOfPrimes. In the calling subroutine I dimension an array named something like FactorList(20,2) and a FactorCounter variable. The first dimension of the array holds the distinct prime factors, while its corresponding second dimension holds the exponent of that prime factor. I can’t imagine that any number we come actoss has more than 20 distinct prime factors, since the product of the first 20 prime numbers is ~10^26. So here is my ListPrimeFactors code, with the usual bracket adjustments.
‘puts the prime factors into the FactorList array, along with the number of occurrances
‘be sure to call GetListOfPrimes first
‘The first time you call, FactorNum should be 0 and LowPrimePosition should be 1
Dim a As Long, b As Integer, FoundFactor As Boolean
Dim SqrNumToFact As Single
SqrNumToFact = Sqr(NumberToFactor)
a = LowPrimePosition
Do Until ArrayOfPrimes(a) &GT SqrNumToFact
If NumberToFactor Mod ArrayOfPrimes(a) = 0 Then
Call AddFactor(ArrayOfPrimes(a), FactorList, FactorNum)
Call ListPrimeFactors(NumberToFactor / ArrayOfPrimes(a), FactorList, FactorNum, a)
Exit Sub
End If
a = a + 1
Loop
Call AddFactor(NumberToFactor, FactorList, FactorNum)
End Sub
It uses trial division from the list of known primes, up to the square root of the number. If it finds a factor it gets added to the FactorList array with the AddFactor (not shown) sub. The sub recurses, passing a new number to itself (new number = old number / prime factor). It also keeps track of the last prime factor, so that each time it recurses it doesn’t check the previous primes. For example if you’re factoring 33, 3 is the first prime factor you find. The sub will recurse, passing 11 as the new number to factor. Trial division will now start with 3, and not go back to 2. If you’ve tried trial division by all primes up to the square root and none divide evenly, then the number must be itself prime and it is added to the list.
This whole part, generating the primes and computing rad() for all numbers, takes under half a second. One area for improvement that I see is to use a recursive sieve to compute all of the prime factorizations. I’ve used it for other problems, but I’m not sure that the boost here is worth the effort. The idea is something like this:
Get a list of all prime numbers smaller than the maximum, say for example 10 (2,3,5,7). Start with the smallest number, 2. The PF of 2 is 2^1. Recurse, passing as arguments the prime factorization and the prime position (1 in this case). Multiply the current number by the appropriate prime and see if it’s less than the max. 2^1 * 2 = 2^2, less than 10, do it again. 2^2 * 2 = 2^3. 2^3 * 2 = 16, which is greater than 10, so here we exit. The recursion returns to 2^2, which we multiply by the next prime factor, 3. 2^2 * 3 = 12, so we exit again and it returns to 2^1. We multiply by the next prime factor, 3, and get 2^1 * 3^1 = 6. We recurse, multiply by 3, and get 18. We go back to 2^1, multiply by 5, and get 10. 2*7 = 14, so now we have the PF of all numbers containing 2. The sub then uses 3 as the first prime factor, and continues on. Sheesh, I hope that was clear. Here’s some pseudocode:
Sub GenerateAllPrimeFactors
Get a list of prime factors less than the limit
Call PrimeFactors(CurrentNumber, PFList, 1)
End Sub
Sub PrimeFactors(ByVal CurrentNumber, ByRef PFList(), ByVal PrimePosition)
For a = PrimePosition to # of Primes
If CurrentNumber * PrimeNumber(a) > Limit then exit for
Add the prime number to the PFList
Operate on the PFList if you want – compute rad(), sum of the factors, number of factors, etc.
Call PrimeFactors(CurrentNumber * PrimeNumber(a), PFList, a)
Remove the prime number from the PFList
Next a
End Sub
If you need to compute, for example, the sum of divisors for all numbers less than 10,000,000, then this is definitely the way to go. I hope that was all understandable by people other than me.
Hi Josh –
I’m happy to keep the Euler conversations going. I’m now up to 101. If there’s a particular problem you’d want to see, speak up. I’ve pretty much done the 1st 40%, except the Monopoly problem and the Sudoku problem. I’m afraid the Monopoly problem is going to be like the poker problem (big effort for a single point) and Sudoku kicked my butt some years ago before Euler, and I’m afraid it will still. There are XL solutions for Sudoku on the net. There are several problems that I solved, but not in Euler time. I just didn’t have the interest to go back and find a better way. My dice solution (#205) I think takes 6.5 minutes, and my solution rate is lagging. The real world is making some demands ;-)
Last solution was #76, which I came to see finally as a Euler Partition problem. #77 is the same, but I can’t grasp how to move my #76 solution to #77 where the steps are primes vice unitary. #77 has assumed the Sudoku position wrt my backside.
Web queries are simple, really. Record one and you pretty much have all you need to know. When I do them in a Macro, the destination range is on a TEMP worksheet I install (and remove on close). I flush the TEMP used range between queries. Before I did both, I seriously corrupted a workbook. Since, knock on formica, I haven’t.
Thank you for the steer to Chip’s function. I went and got it. Am I the only one with a subfolder inside his Excel folder labeled MVPs? That’s where those things go.
Would you post your AddFactor() routine? I’ve been wanting a routine that adds to and sorts a list on the fly, and it sounds like a good start.
…mrt
Gah, silly me – I missed a zero in the problem statement. I thought we were looking for E(100000) i.e the last number in the sequence. As 999991 is prime, the last number in the ordered list must be between 99991 and 99999.
But of course that isn’t what the problem is asking. (Hence why I feel OK now about posting actual numbers).
Rick,
If I had a nickel for every time I’ve misread a Euler problem statement…
Michael,
I’m pretty sure I’ve figured out your user name on PE, so it’s easy to see which problems you’ve solved. If you go to Statistics > BASIC, my username is jgrilly. The Monopoly problem was pretty annoying. To tell you the truth, I couldn’t get the right answer. I had to try various combinations of my top 4 squares to get the right answer. I actually enjoyed the sudoku problem. That one seemed like it revealed a useful algorithm (well, useful if you’re an Excel nerd and making your own sudoku program…). I also have a number of solutions not solved in Euler timer. 76, for example, fell to brute force after several hours. I like going to the forums afterwords and learning staggeringly simpler methods, like partition functions (which I’d never heard of before). If the forums show some solution that I’ve never seen / thought of (like 76), I usually go back and code a new solution. If the “optimal” solution is uninteresting, I just let it be. Looking at the forum for problem #77, many people adapted their solution from #76. I found that it fell to brute force in just a couple seconds, though picking the right upper bounds is key. My solution to #205 takes about 0.06 seconds. I say that not to brag, but to stimulate discussion if you like. Many of the other probability problems (213, 227, etc.) are kicking my behind. I’m currently working on #184. I have an idea, though it’s taking a while to code. I’m pretty sure I’m not going to correctly count all of the triangle reflections / rotations, and I have no idea if it’ll be anywhere close to Euler time. That’s what makes it fun though.
I guess I’ve never really needed web queries before, but they sound simple enough. I do love the macro recorder.
Here’s my AddFactor routine. It plays well with the ListPrimeFactors routine. Since the prime factors are found in order from smallest to largest, your list will always be ordered from smallest prime factor to largest prime factor. And since the routine passes the number of distinct prime factors found so far, you only need to check the last value, or add a new spot. As before, NumToAdd is the prime factor to add to the list, FactorList(20,2) is a 2D array containing the prime factors and their corresponding exponents, and NumFactors is the number of distinct prime factors found so far. Usual angle bracket adjustments made.
If NumFactors =! 0 Then
If FactorList(NumFactors, 1) = NumToAdd Then
FactorList(NumFactors, 2) = FactorList(NumFactors, 2) + 1
Exit Sub
End If
End If
NumFactors = NumFactors + 1
FactorList(NumFactors, 1) = NumToAdd
FactorList(NumFactors, 2) = 1
End Sub
If you’re looking to add to a list and then sort it on the fly, may I recommend something like Insertion Sort to start? It should work well for adding one element to an already-sorted list, and it’s pretty easy to implement.
-Josh
Hi Josh –
Problem 76 it is this weekend. If you’ve not redone it, it takes no time at all if you’re able to Google the right reference. Turns out it’s pointed to inside the forum. Same code with a tweak knocked off #78. #77 I’ll leave for you to show the way.
Then I’ll put up my brute-force #205 next week and you can really illuminate the way.
InsertionSort is exactly what I want to do. Sounds like it’s out there.
If you started with my sig, you’ve probably got me ;-)
…mrt
Michael,
Looking back through my solutions, I apparently went back and coded up the “optimal” solution for #76 and then used that for #78 as well. I’d still like to find a nice solution to #77, since I’m not really happy with my brute force approach.
A little while ago I got interested in sorting algorithms when I learned that my naive homemade sorting algorithm, which turned out to be Bubble Sort, was often given as an example of what NOT to do. My chest swelled with pride at coming up with something so inefficient on my own, but then I decided to find something better. It seems that Insertion Sort is quite good for small lists and lists that are already substantially sorted. Here’s my Insertion Sort implementation for a 1-D array of longs, with the usual angle bracket adjustments:
‘Sorts low to high
Dim a As Long, b As Long, SortTemp As Long
For a = L + 1 To R ‘a is the first unsorted key
SortTemp = SortArray(a)
For b = a To L + 1 Step -1 ‘compare it back through the sorted part if it’s bigger
If SortTemp &LT SortArray(b – 1) Then
SortArray(b) = SortArray(b – 1)
Else
Exit For
End If
Next b
SortArray(b) = SortTemp ‘it’s bigger than all keys to the left, so insert it here
Next a
End Sub
The input parameter L is the lower bound in the list you want to sort, while R is the upper bound. Since it sorts the array in place, you can also use this to only sort a small piece of a larger list. For example, I use QuickSort to do most of the heavy lifting, then switch to InsertionSort to do the small pieces. If the list is already sorted and you just want to add one element to it, you have two options. First, you can just let it go as it is. It’ll make about n comparisons, where n is the size of the list, without swapping anything. Depending on the size of the list, this may or may not matter. The other option is to change the “For a = L + 1 to R” line to say “For a = R to R”. This way you will only look at the last element.
-Josh
Hi Josh –
Following your steer to the BASIC statistics, Dick should reverse our roles. In one sense I’m only about 30 behind you, but in the other I’m about 110 behind you. That’s a lot of coding! I better get hot ;-)
There’s probably few, if any, I’ve done that you’ve not, so I’ll put #76 up anyway. I have a good reference I can point to. #77 is driving/drove me buggy. I know the answer, it’s in the OEIS. I just can’t compute the answer. The Haskell one-liners are no help and the one guy who put up some BASIC is so obfuscated that I can’t follow his juggling indices. I can often translate the Ruby, C, and Java, never the Python. Those guys say they just hacked their 76 code. Well, that’s the part I can’t get…
Thank you for the Insertion sort…looks very useful.
…mrt
You know what I am going to write, don’t you?
VBA for the cheat version? Why?
Copy the data from the web page you mentioned.
Paste the data into Excel. I put the data in cols. A and B.
Sort on: B ascending then A ascending.
Go to A10000. Get the answer.
Hi Tushar –
Yep, that’ll work if you can paste 100,000 rows. But I’m row-challenged still ;-)
I get a nasty alert that says “File not loaded completely” when I paste it. And if I can only use and sort the top 65536 rows, the correct answer is in there, but it’s not even close to A10000.
What I thought would be neat was a way to get 100000 rows off the web into an array. An XL2007 webquery probably will.
…mrt
Hi Michael,
Let’s just say that there are a few days where I don’t get any real work done…
I’m having the same problems as you – I couldn’t see how to re-use my #76 code for #77, and I couldn’t really understand the others’ solutions. There’s a link to the MathWorld article on prime partitions, but I couldn’t see how to apply it. Bah. Maybe I should go back and look at the expert solutions to #31.
-Josh
Josh –
Great thinking alike ;-)
There’s a neat recursive answer to #31 that drew Euler’s praise right at the top. It translates easily to VBA. Works very well for 8 coins. Not quite sure how to adapt it to a potful of primes.
…mrt
Michael,
Would you mind posting your translation of #31 in VBA? I keep running out of stack space, so I must be doing something wrong. Thanks.
-Josh
Hi Josh –
It runs out of stack space for me in Mac XL 2004. Runs fine in XL2002. For what it’s worth:
Sub Main()
coins(0) = 200
coins(1) = 100
coins(2) = 50
coins(3) = 20
coins(4) = 10
coins(5) = 5
coins(6) = 2
coins(7) = 1
Debug.Print findposs(200, 0)
End Sub
Function findposs(money As Long, maxcoin As Long) As Long
Dim sum As Long, i As Long
sum = 0
If maxcoin = 7 Then
findposs = 1
Exit Function
End If
For i = maxcoin To 7
If money – coins(i) = 0 Then sum = sum + 1
If money – coins(i) > 0 Then sum = sum + findposs(money – coins(i), i)
Next i
findposs = sum
End Function
…mrt
Michael, Josh,
If you are looking for a solution to 31 you may want to check
Project Euler – Problem 31
http://www.tushar-mehta.com/misc_tutorials/project_euler/euler031.html
Hi Tushar –
Thank you, but you misunderstand…we’ve both solved 31 (that code above is a different recursive solution to 31) and we’re are looking therein for hint to a solution to 77. Rather than units of currency, #77 uses prime numbers.
Trying to make a fit ;-)
…mrt
Hi Michael,
Thanks for posting the code. It works fine for me – not sure what my problem was. I managed to adapt it to solve #77. I’ll post it in your thread for #76, since it seems more appropriate there.
-Josh