Wednesday, September 27, 2017

Sum of multiples problem

One of the problems generalizes a Project Euler problem: what's the sum of multiples of 3 and 5 less than 1000? Seems simple enough:
sum [3, 6..999] + sum [5, 10..999]
but that's too big; it counts some values twice. With a little thought, though, you'll realize that the duplicates are the multiples of 15, so
sum [3, 6..999] + sum [5, 10..999] - sum[15, 30..999]
will do the trick. Why 15? Well, it's 3 * 5... but more on that later.

The problem makes the upper bound a parameter, and allows an array (or list, depending on the language) of factors. Looking around (which you can do on the site after you've submitted your own solution), one sees two typical approaches:
  • Go through [1..limit - 1], pick out the values that are multiples of one or more of the factors, and add up the result.
  • Create a set or hash table from [[f, 2*f..limit - 1] f <- factors] (thus getting rid of the duplicates), and add up the elements of the set/keys of the hash table.
Both work, but they have drawbacks:
  • Divisibility checking means using the relatively expensive mod, potentially repeatedly, for each value in [1..limit - 1], and as many times as possible for values that aren't summed. As limit grows, that overhead grows.
  • Sets and hash tables build up an underlying data structure that grows with each value added to it. As limit grows, so does RAM usage. (My Python version using this method, handed an upper bound much larger than the test values and set up to time with Python's timeit package, quickly sucked down immense amounts of virtual memory and had to be killed.)
There ought to be a way to generalize what we did for the initial [3, 5] case... and indeed there is. What did we do? We took the sum of the multiples of each factor and then subtracted off the sum of the duplicates. If there's just one factor, that sum is zero. If there are two, it's the sum of the multiples of... the product of the factors? Not quite!

Suppose that our factors are 6 and 21. 6 * 21 is 126, but 7 * 6 == 2 * 21 == 42. We want not the product, but the least common multiple of the two factors. 3 and 5 are relatively prime, so their least common multiple is their product (lcm m n == m * n `div` (gcd m n), and m and n are relatively prime iff their greatest common divisor is 1.) Try to fake us out, eh?

More generally, if you have factors f1, f2, ..., fn, each fi gives rise to {lcm fi fj | i < j}  as new sets of factors whose multiples are duplicates of the multiples of fi--but those in turn have duplicate issues, so it's time to recur. You won't recur forever, because the sets of factors get smaller each time.

One more thing... we don't need to evaluate sum [f, f+f, .. limit - 1] at all. There's a better way to do it, and since Young Sheldon had its premiere this week, let's return to the 19th century and Young Carl to learn why we don't have to generate that arithmetic sequence or count on GHC to do fusion. Young Carl Friedrich Gauss was assigned the problem of adding up all the integers from one to one hundred (in the first version I read, the class was assigned the problem, and only after solving it could a student go out to recess). The teacher looked, but Carl wasn't busily adding numbers. Soon he wrote down a number and handed the paper in, and to the teacher's amazement it was correct. The math book I learned it from explained Gauss's insight this way: consider the numbers from 1 to 100 written down twice: once in ascending order, and immediately under that in descending order, and consider corresponding terms in the two sequences. The first, 1 and 100, add up to 101. To get from one pair to the next, the top increases by 1 and the bottom decreases by one, so the sum remains the same: 101. There are 100 such pairs, so collectively they add up to 100 * 101 = 10100... but remember, we wrote the sequence down twice, so that sum is twice the value we want. The number young Carl wrote was 5050, and in general, the sum of the integers from 1 to n is n * (n + 1) / 2, and as an immediate consequence the sum of the multiples of f less than limit is f * n * (n + 1) `div` 2 where n = (limit - 1) `div` f.

Here's what I ended up with. Haskell has both gcd and lcm functions, I'm happy to say.

sumOfMultiples :: [Integer] -> Integer -> Integer
sumOfMultiples factors limit =
    sum (map oneFactorSum factors) -
        sum (map (`sumOfMultiples` limit) (dupFactors factors))
    where sumOneTo     n  = n * (n + 1) `div` 2
          oneFactorSum f  = f * sumOneTo ((limit - 1) `div` f)
          dupFactors   fs = map lcmHead $ take (length fs - 1) (iterate tail fs)
          lcmHead      fs = [lcm (head fs) f | f <- tail fs]

Sunday, August 06, 2017

a longest path problem has a weekly column, "The Riddler". The Riddler poses two problems in each column. The first, "Riddler Express", is intended to be an easier problem that one can solve quickly, while the second, "Riddler Classic", is more difficult. Modulo vacations, it appears each Friday, and if you submit a solution by the end of the following Sunday (Eastern Time) , you will be among those who might be given credit for the solution in the next column.

The August 28th column's Classic problem is as follows: what's the longest sequence of integers x[i] for 1 ≤  i ≤  n such that
  • for all i, 1 ≤ x[i] ≤ 100
  • for all i < n, either x[i+1] is a multiple of x[i] or x[i+1] is a factor of x[i]
  • for all i, j, x[i] = x[j] iff i = j, i.e. the x[i] are all distinct
(The last constraint avoids trivial sequences like 2, 4, 2, 4, 2, 4... which can go on forever if repeats are permitted.)

One way to characterize this problem is to look at it as a graph with a hundred vertices v[i]. Label v[i] with i, and connect v[i] and v[j], for i ≠ j, with an edge if either i is a multiple of j or j is a multiple of i. Then the problem asks for the longest path through some (or all, if it proves possible!) of the nodes of the graph following the edges with no node visited more than once. (If you solve the one-hundred node problem, the Riddler suggests you try the analogous problem with a thousand nodes.)

Alas, longest path is NP-hard. If you come up with an efficient way to solve this, your name will go down in computational complexity history. So I don't feel too bad starting out with a brute force depth-first search, and since I'm working to improve my Python skills, I start in Python.

The Riddler reports that one person ground through 30 million trials and came up with a sequence of length 55, which makes me suspect that our programs are similar. I started out with the obvious code using Python's built-in set and list types, and after over ten hours, that version had only managed to find a sequence of length 54. Then I started rewriting.
  • Step One: switch from the Python set type to using an integer (Python 3 has arbitrary precision integers, so yay). That did speed things up.
  • Step Two: rather than making lists come and go, accumulate the sequences in arrays that the recursive calls build up the sequence in. That sped things up even more.
The result: the length of known best results zips up to about the mid-40s very quickly indeed, but length 55 took at least two hours to appear. We're at four hours plus change, and 56 has yet to be seen. (UPDATE: showed up sometime between four and ten hours. Next version timestamps the output.)

Next stage will be a Go version to avoid interpreter overhead, and then dividing it up to take advantage of multiple cores, though given that the 55-number sequence is still very early on in the sequences starting with 1 (!?), I wonder how much good that will do. Here's that 55-number sequence, by the way:
1, 2, 4, 8, 16, 32, 96, 3, 6, 12, 24, 72, 36, 18, 54, 27, 81, 9, 99, 33, 66, 22, 44, 88, 11, 55, 5, 40, 80, 20, 60, 30, 90, 5, 15, 75, 25, 50, 100, 10, 70, 14, 56, 28, 84, 42, 21, 63, 7, 91, 13, 39, 78, 26, 52
The longest submitted sequences have length 77. More news as it happens. (Might experiment with how to pick the next item to try: the one with the most or least possible successors, perhaps? OTOH, if it were most, you'd always start with 1, and neither of the length 77 sequences the Riddler shows starts with 1.)

Thursday, June 15, 2017

Bikeshedding leap years

Darn it, I fell into the trap. An early exercise in Python at wants a function that in Haskell would have the type Int -> Bool and indicate whether a year is a leap year (under Gregorian rules, ignoring when various locales made the switch to keep it simple).

Some submissions are from people not yet comfortable with assigning/returning Booleans other than the constants True or False, or maybe not comfortable with the lazy nature of Python's and and or. Their versions are full of if statements. I bashed out
def is_leap_year(year): return year % (400 if year % 100 == 0 else 4) == 0
and it worked fine... but then I fell in. How about
return year % 16 == 0 if year % 25 == 0 else year % 4 == 0
Nah... no need to make the next guy, who could be me, do the prime factorization of 100 again and figure out how that ties to the leap year rules? OK, here's one way to look at it...
return (year // 100 if year % 100 == 0 else year) % 4 == 0
i.e. on century boundaries, the century has to be divisible by four, otherwise the year has to be. Well, that's one way to look at it, but one that doesn't go with the definition as clearly as one would like...

...and then the path of premature optimization that Knuth warns against beckoned.
return year % 4 == 0 and (year % 100 != 0 or year % 400 == 0)
The idea: three fourths of numbers get ruled out by a check for divisibility by four, which can be done cheaply if the compiler or interpreter notice. (That works for any power of two, not just four, which with another nod to factorization means one could write
return year % 4 == 0 and (year % 100 != 0 or year % 16 == 0)
which means, if you're lucky, only doing one divisibility check with an actual divide.)

Final decision? Leave it alone

Tuesday, June 06, 2017

Restaurants shooting themselves in the foot

It's pet peeve time, and here it is: restaurants that display graphic images of their menu on the web. Sometimes it's a PDF made up of graphics, sometimes it's just graphics.

I wish I had a nickel for each restaurant I really like whose web site suffers from this problem... though perhaps I should shoot for more than that by offering my services to fix their web sites. A few examples:
How wrong is this? Let me count the ways...
  • You say you want people to be able to search your menu? Good luck (though someday the code that lets Google Translate read text from images may get good enough to deal with the wonky display fonts people tend to use on menus or text superimposed on food photos).
  • You say you want search engines to find your restaurant's web site? Good luck again, though as above, maybe someday...
  • You say you want your web site to support travelers looking for a place to eat even when they're on the road and only able to manage a 3G or 2G connection? Then why do you make them download large graphics files to look at your menu?
  • The blind eat, too... but restaurants whose web sites don't have text menus apparently don't care. (I did an "inspect" on a menu image; its ALT attribute was, I kid you not, "Picture".

Sunday, April 02, 2017

No tutorial, I swear...

After grinding through much of the "20 Intermediate Exercises" and just about all of the Monad Challenges, Monads make much more sense to me.

Therefore, I will follow the excellent advice of Stephen Diehl, and will not write a monad-analogy tutorial. Instead, I will say: do the 20 Intermediate Exercises, and take the Monad Challenges. To paraphrase Euclid, there is no royal road to Monads; the exercises and challenges take you down that road at whose end you'll see at least part of why Monads are so useful.

If you were trying to learn group theory and people were standing around you saying "groups are like the integers",  "groups are like Rubik's Cube", "groups are like m x n matrices", "groups are like baryons and mesons", your situation would be much like the student of Haskell amidst all the monad analogy tutorials. In a sense they're backwards. All those things are groups, just as burritos et al. can at least be thought of as monads--but not all groups are any one of those things.

Yeah, I suppose this is a meta-monad tutorial. Shame on me.

Friday, March 24, 2017

Flipping out

I came across the excellent Monad Challenges, a collection of exercises designed to take you through what have become some of the standard example monads and take you through the process that motivates monads as a way to handle these seemingly diverse data structures. Both the Monad Challenges and the 20 Intermediate Exercises are well worth your time. Doing them both is helping me a lot.

That said, they don't quite match up. 20IE's banana is flip bind and apple is flip ap. (This isn't unique; the also excellent Learning Haskell from first principles has as an exercise writing bind in terms of fmap and join, but the declaration it gives has the type of flip bind. The authors do point this out.) As a result, I find myself with something that there's got to be some way to simplify, of the form

foo = flip $ (flip mumble) . (flip frotz)

I'd like to think there's some sort of distributive-like law to be found and used here... watch this space.

Friday, March 03, 2017

Trust in Schönfinkel

I'm working through the "20 Intermediate Exercises" that give Functors and Monads and such cute, non-threatening names and ask you to implement them. I've gotten most of the way through, with three or four that I'm stuck on--that's under the assumption that the inductive step from banana to banana2 will make the rest of the banana<n> obvious. (If you have sausage, it's easy to implement moppy, and vice versa, but avoiding circularity is the issue.)

So... I did a Bad Thing and looked at solutions a couple of people have posted, saying to myself, "OK, I'll look at those I've already done with an eye to more elegant expression, and look at the ones I'm having issues with and make sure I understand"... and came face to face with a question about composition.

<andy_rooney>Didja ever notice that the examples of composition you see are what we think of as functions with one argument?</andy_rooney> Well, one of the solutions involved composition that does not have that constraint, and that caused me much puzzlement until I remembered that we are, after all, using Haskell, where every function just has one argument. All it means is that the left-hand operand of (.) had better want to be passed a function...which is why the very next thing I typed at ghci just for the heck of it was stupid:

Prelude> :t ((+) . (-))
((+) . (-)) :: (Num a, Num (a -> a)) => a -> (a -> a) -> a -> a

Uh, yeah. I am pretty sure no function type is in Num, but what the heck, let's bind it to a name.

Prelude> let foo = ((+) . (-))

    Non type-variable argument in the constraint: Num (a -> a)
    (Use FlexibleContexts to permit this)
    When checking that ‘foo’ has the inferred type
      foo :: forall a. (Num a, Num (a -> a)) => a -> (a -> a) -> a -> a

Didn't expect that. Perhaps ghci isn't quite as strict when you're just asking for the type.

UPDATE: I see why the solution didn't look type correct to me. I didn't understand Hindley-Milner well enough. It must be willing to do things like decide "to make this work, type a over here must map to type a -> b over there" in such a way that all is consistent. Now I know what to learn next.

Monday, January 23, 2017

Careful with that infinite list, Eugene...

Warning: this will give away one way to solve a certain low-numbered Project Euler problem.

Any Haskell book, blog, or tutorials you come across has a good chance of including the très élégant Haskell one-liner for an infinite list containing the Fibonacci sequence:

fibs = 0 : 1 : (zipWith (+) fibs (tail fibs))

Thanks to laziness, it will only evaluate out to the last term we actually request... that was, long ago, my downfall. I wanted the terms no bigger than four million, so

filter (<= 4000000) fibs

right? Wrong. You and I know that the Fibonacci sequence is monotonically increasing, but filter doesn't, and it doesn't even notice the particular function we're filtering with to realize it could terminate thanks to monotonicity. So instead,

takeWhile (<= 4000000) fibs

is the way to go.

The particular problem has an additional constraint, because it only wants the even terms of the sequence no bigger than four million. Easy enough to do, just feed it through

filter even

and you're good. It did, though, occur to me: does the order of filtering matter? I mean, just plain

filter even fibs

is still an infinite list... but as long as we only actually evaluate finitely many values, we're safe.

filter even $ takeWhile (<= 4000000) fibs


takeWhile (<= 4000000) $ filter even fibs

both terminate and, as you'd expect, give the same result--but while the order doesn't matter to the result, it may matter for speed. Given that the first Fibonacci number is even and the next odd, the whole sequence is thus even, odd, odd, even, odd, odd, ... so

takeWhile (<= 4000000) $ filter even fibs

will only check a third as many values for being no more than four million.

A flashback and analogy

You've probably heard about how the notion of sum types (e.g. Algol 68 union s, Rust enum s, Haskell type s) and product types (e.g. tup...