Planet Scheme

July 29, 2014

Programming Praxis

Montgomery Multiplication

[ Does anybody know the proper method for writing a-bar and b-bar in plain HTML? The LaTeX symbols are \={a} and \={b}, but WordPress renders LaTeX poorly, and I prefer plain HTML. FIXED! Thanks to Jussi, who suggests writing a amp hash x304 semicolon and b amp hash x304 semicolon. ]

Montgomery multiplication is a method of computing ab mod m for positive integers a, b and m, designed for situations where there are many multiplications to be performed mod m with a small number of multipliers; in particular, it is valuable for computing xy mod m for large y. Montgomery multiplication was invented by Peter Montgomery in a 1985 article. Our exercise today will implement modular multiplication and exponentiation using Montgomery’s method.

The basic idea of Montgomery multiplication is to eliminate the division inherent in the modulo operation by translating the modulus m to a larger modulus r, with gcd(r, m) = 1, where r is a power of 2; thus, the division is rendered as a bit mask. Since r is a power of 2, the modulus m must be odd to satisfy the gcd requirement; we also require the two multipliers a and b be less than m. Then the Montgomery algorithm is as follows:

  1. Find two integers r−1 and m′ such that r r−1 m m′ = 1 by the extended Euclidean algorithm. In particular, the algorithm used is a binary version of the algorithm that performs no divisions.
  2. Translate the multipliers to “Montgomery space” by multiplying them by r (using left shifts) and reducing the product modulo m. Thus, ā = a r mod m and b̄ = b r mod m. These operations are expensive, but they are performed only once for each multiplier in a chain of multiplications.
  3. Perform the Montgomery multiplication: calculate u = a b r mod m by the sequence of operations

    t = ā * b̄
    u = (t + (t * m′ mod r) * m) / r
    if (um) then subtract m from u

    The only division is a right-shift; the final result uses subtraction instead of division because it is known that the product is less than 2m.

  4. Translate the result back from Montgomery space to an ordinary integer by ab = u r−1 mod m.

Montgomery exponentiation is similar. To compute xy mod m, select r, compute r−1 and m′, translate x to Montgomery space, use the square-and-multiply method with Montgomery multiplications to perform the exponentiation starting from 1 in Montgomery space, then translate the result back from Montgomery space. This is much faster than the calculation without Montgomery multiplication, because the cost of initializing the Montgomery space is amortized over several multiplications, each of which avoids the inherent division of the modulo operator.

Your task is to write functions that perform Montgomery multiplication and exponentiation. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

by programmingpraxis at July 29, 2014 09:00 AM

July 26, 2014

The Racket Blog

Scheme Workshop 2014

DEADLINE: 5 September 2014, (23:59 UTC-12)
LOCATION: Washington, DC (co-located with Clojure/conj)
DATE: 19 November 2014

The 2014 Scheme and Functional Programming Workshop is calling for

Submissions related to Scheme and functional programming are welcome
and encouraged. Topics of interest include but are not limited to:

  • Program-development environments, debugging, testing
  • Implementation (interpreters, compilers, tools, benchmarks, etc)
  • Syntax, macros, and hygiene
  • Distributed computing, concurrency, parallelism
  • Interoperability with other languages, FFIs
  • Continuations, modules, object systems, types
  • Theory, formal semantics, correctness
  • History, evolution and standardization of Scheme
  • Applications, experience and industrial uses of Scheme
  • Education
  • Scheme pearls (elegant, instructive uses of Scheme)

We also welcome papers related to dynamic or multiparadigmatic
languages and programming techniques.

Full papers are due 5 September 2014.
Authors will be notified by 10 October 2014.
Camera-ready versions are due 24 Oct 2014.
All deadlines are (23:59 UTC-12), "Anywhere on Earth".

For more information, please see:

See you there!

by John Clements ( at July 26, 2014 10:34 PM

July 25, 2014

Programming Praxis

Number Words

Today’s exercise is an interview question from Career Cup. When I first saw the exercise, I thought it would be easy, but I had trouble writing a solution, so maybe it will a good exercise for all you readers, as well:

Given a positive integer, return all the ways that the integer can be represented by letters using the mapping 1 -> A, 2 -> B, …, 26 -> Z. For instance, the number 1234 can be represented by the words ABCD, AWD and LCD.

Your task is to write the program to generate words from numbers. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

by programmingpraxis at July 25, 2014 09:00 AM

GNU Guix

GNU Guix 0.7 released

We are pleased to announce the next alpha release of GNU Guix, version 0.7.

This release is an important milestone for the project since it is the first to provide an image to install the GNU system from a USB stick.

Noteworthy features for the release are:

  • The GNU operating system can now be installed. Try it out!
  • To make it possible the guix system command has been augmented with new options, and support for 'operating-system' declarations has been vastly improved.
  • Programming has been simplified with the introduction of "G-expressions", which capture dependencies used by build-side expressions.
  • More than 130 packages have been added, including "big ones" like the GIMP and Maxima.

See the original announcement for details.

by Ludovic Courtès at July 25, 2014 07:28 AM

July 22, 2014

Programming Praxis

Maximum Sum Path

Sometimes here at Programming Praxis we publish exercises taken from homework in computer science classes, sometimes we publish exercises that are based on or related to problems at Project Euler, and sometimes we publish exercises that based on common technical interview questions. Today we hit the trifecta: a common interview question that is frequently assigned as homework in computer science classes and that also appears at Project Euler (twice!):

Given a binary tree with integers stored in its nodes, find the maximum sum of any path from root to leaf.

For instance, in the binary tree represented by the triangle shown below, the maximum sum path goes from the 3 at the root, to its child 7, to its grandchild 4, to the leaf node 9, a sum of 23. Since there is no other path with a greater sum, the maximum sum is 23:

  7 4
 2 4 6
8 5 9 3

Your task is to write a program to compute the maximum sum path through a binary tree; also write a variant of the program that returns the path. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

by programmingpraxis at July 22, 2014 09:00 AM

July 18, 2014

Programming Praxis

How Fermat Factored Integers

Living in the age of computers, I’ve always been fascinated by how much arithmetic it was possible to do before the advent of computers. The Egyptians knew the peasant method of multiplication when the built the pyramids. Archimedes knew the value of pi two centuries before Christ, with an approximation we still use today. Heron had an algorithm for calculating square roots two thousand years ago. And so on.

By the time of Fermat, Euler and Gauss, arithmetic was done with simple algorithms, rules of thumb, and large tables. Sadly, the tables were frequently wrong; Charles Babbage invented his Difference Engine for the purpose of automatically constructing accurate tables.

We studied Fermat’s algorithm for factoring integers in a previous exercise:

function fermat(n)
    x := floor(sqrt(n))
    r := x * x - n
    t := 2 * x + 1
    while r is not a square
        r := r + t
        t := t + 2
    x := (t - 1) / 2
    y := sqrt(r)
    return x - y, x + y

But how did Fermat actually work his algorithm? In addition to the basic arithmetic, he used some simple rules for identifying squares and a table of the squares of numbers similar to that on the next page, with the largest square at least as large as n.

Let’s consider the factorization of n = 13290059 as an example. First, Fermat removed small factors by trial division, say for factors less than 30 or 50. Then, Fermat looked up n in the table. If n appeared in the body of the table, then it is a perfect square, and the square root is its factor. But in the general case, Fermat finds x as the square root of the next smaller entry in the table. For 13290059, x = 3645, because 36452 < 13290059 < 36462. Then t = 7291 and r = -4034, which is not a perfect square (because it’s negative).

Now the iteration begins. Fermat calculated r and t in 45 steps as follows, first r, then t below it, then a sum line, then the sum r = r + t, then t incremented by 2 from the previous t, then a sum line, then a new sum r, at each step checking if r is a square:

0:  r =  -4034 = 3645 * 3645 - 13290059
    t =   7291 = 2 * 3645 + 1
1:  r =   3257 ends in 7; not a square
    t =   7293 = 7291 + 2
2:  r =  10550 ends in 50; not a square
    t =   7295 = 7293 + 2
3:  r =  17845 ends in 45; not a square
44: r = 318662 ends in 2; not a square
    t =   7381 = 7379 + 2
45: r = 326041 in table: 571 * 571

At that point, Fermat could calculate the factors of n: x = (7381 – 1) / 2 = 3690, y = 571, and the factors of n are xy = 3690 − 571 = 3119 and x + y = 3690 + 571 = 4261.

The key to Fermat’s calculation was identifying the case where r is a square. Fermat used three rules: First, the last two digits of the number had to be 00, e1, e4, 25, o6 or e9, where e is any even digit and o is any odd digit. Second, the digital root of the number had to be 1, 4, 7 or 9: sum the digits of n; if the result is less than 10, then return the sum of the digits, otherwise return the digital root of the sum of the digits. In practice, this goes much more quickly than it appears by “casting out nines” wherever they appear in the digits of n or the result. Third, if he had not already determined that the number was not a square, he looked it up in his table of squares. The whole calculation goes much more quickly than I can describe it here.

And that’s how Fermat factored integers. Of course, sometimes he failed, when it took too many iterations to find a square, or when his table of squares was too small, or when he made an error in the arithmetic. Fortunately, it was easy enough for him to confirm that any factorization he found was correct by multiplying the factors.

Your task is to write a program that factors integers the same way Fermat factored integers. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the program in the comments below.

by programmingpraxis at July 18, 2014 09:00 AM

July 15, 2014

Programming Praxis

How Many Distinct Products In A Times Table?

We all learned our times tables in elementary school. For instance, here is the standard 9 by 9 times table:

1  2  3  4  5  6  7  8  9
2  4  6  8 10 12 14 16 18
3  6  9 12 15 18 21 24 27
4  8 12 16 20 24 28 32 36
5 10 15 20 25 30 35 40 45
6 12 18 24 30 36 42 48 54
7 14 21 28 35 42 49 56 63
8 16 24 32 40 48 56 64 72
9 18 27 36 45 54 63 72 81

That table has 36 distinct products: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 24, 25, 27, 28, 30, 32, 35, 36, 40, 42, 45, 48, 49, 54, 56, 63, 64, 72 and 81.

Your task is to write a program that computes the number of distinct products in an n by n times table. There is an obvious O(n2) algorithm; can you do better? When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

by programmingpraxis at July 15, 2014 09:00 AM

July 14, 2014

Jeremy Kun

An Update on “Coloring Resilient Graphs”

A while back I announced a preprint of a paper on coloring graphs with certain resilience properties. I’m pleased to announce that it’s been accepted to the Mathematical Foundations of Computer Science 2014, which is being held in Budapest this year. Since we first published the preprint we’ve actually proved some additional results about resilience, and […]

by Jeremy Kun at July 14, 2014 03:00 PM

July 13, 2014

GNU Guix

Guix at OpenBio Codefest 2014

On Wednesday, July 9th, David Thompson gave a brief introduction to GNU Guix at the Open Bioinformatics Codefest 2014 hackathon. The objective of the Codefest is to give developers of bioinformatics software a chance to be fully focused on their projects for a few days and work in person. These developers are concerned with the reliability and reproducibility of their operating systems, and the limitations of their package management utilities.

See the slides on-line.

by David Thompson at July 13, 2014 08:26 PM

July 11, 2014

Programming Praxis

FSM Generator

A few weeks ago we had an exercise that solved the problem of removing singleton occurrences of a given character from a string while leaving multiple occurrences of the given character intact. Our solution used a finite state machine that iterated over the input string writing output as it went.

The finite state machine used in that solution was hand-crafted to solve the given problem. But it is possible to write a program that takes a definition of a finite state machine and creates a function that takes an input string and applies the finite state machine to it. Such a program isn’t hard to create and provides a useful utility for small parsing problems.

Your task is to write a program that generates a finite state machine from a simple description; the format and capabilities of the machine are up to you, but it should be at least powerful enough to solve the “remove singleton” problem. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

by programmingpraxis at July 11, 2014 09:00 AM

July 08, 2014

Programming Praxis

SQUFOF, Continued Fraction Version

Square Forms Factorization, abbreviated SQUFOF, is a special-purpose method for factoring integers. It is commonly used to split small semi-primes in the double-large-prime variants of the quadratic sieve and number field sieve factoring algorithms. We studied the “binary quadratic forms” variant of SQUFOF in a previous exercise; today we study the “continued fraction” variant of SQUFOF.

Although the mathematical basis of SQUFOF is the quadratic forms of Gauss, SQUFOF was discovered by Daniel Shanks while he was investigating the failures of the continued fraction factoring algorithm; you should read his paper, which is as good a demonstration of the process of mathematical discovery as anything I have ever read. Shanks’ algorithm expands the continued fraction of the square root of n, the number being factored, searching among the convergents for a quadratic form that is a perfect square; in this, it differs from the CFRAC algorithm, which combines multiple convergents to assemble a perfect square. We copy the algorithm from the paper by Jason Gower and Sam Wagstaff that provides the standard description of SQUFOF:

[ Note: In the description that follows, be careful to distinguish the symbols and Q; the first one has a circumflex over the Q, the second one doesn't. Some browsers render the "combining circumflex accent" readably, but some don't. On the browsers available to me, only one out of four rendered the symbol readably, sorta kinda, and only when I expanded the font to be quite large. If in doubt, look at the program on the next page, where the first symbol is named qhat and the second symbol is named q. ]

The variable N is the integer to factor, S remembers q0, q holds the current qi, P and P’ hold two consecutive values of Pi, and Q hold two consecutive values of Qi, and t is a temporary variable used in updating Q. Finally, B is an upper bound on the number of forms tested for being square forms before the algorithm gives up.

1. Initialize: Read the odd positive integer N to be factored. If N is the square of an integer, output the square root and stop. If N ≡ 1 mod 4, then set D ← 2N; otherwise, set DN. In any case, set S ← ⌊√D⌋, Q̂ ← 1, PS, QDP · P, L ← ⌊2√2√D⌋, B ← 2 · L, and i ← 0.

2. Cycle forward to find a proper square form: Steps 2a through 2e are repeated for i = 1, 2, 3, ….

2a: Set q ← ⌊(S + P) / Q⌋ and P’q · QP.

2b: If Ql, then:

If Q is even, put the pair (Q / 2, P mod (Q / 2) onto the QUEUE; otherwise, if QL / 2, then put the pair (Q, P mod Q) onto the QUEUE.

2c: Set tQ̂ + q · (PP’), Q̂ ← t, and PP’.

2d: If i is odd, go to Step 2e. If Q is not the square of an integer, go to Step 2e. Otherwise, set r ← √Q, a positive integer. If there is no pair (r, t) in the QUEUE for which r divides Pt, then go to Step 3. If r > 1 and there is a pair (r, t) in the QUEUE for which r divides Pt, then remove all pairs from the beginning of the QUEUE up to and including this pair and go to Step 2e. If r = 1 and there is a pair (1, t) in the QUEUE, then the algorithm fails because we have traversed the entire principal period of quadratic forms of discriminant 4N without finding a proper square form.

2e: Let ii + 1. If i > B, give up. Otherwise, go to Step 2a.

3. Compute an inverse square root of the square form: Set Q̂ ← r, PP + r · ⌊(SP) / r⌋, and Q ← (DP · P) / Q̂.

4. Cycle in the reverse direction to find a factor of N:

4a: Set q ← ⌊(S + P) / Q⌋ and P’q · QP.

4b: If P = P’, then go to Step 5.

4c: Set tQ̂ + q · (PP’, Q̂ ← Q, Qt, and PP’ and go to Step 4a.

5. Print the factor of N: If Q is even, set QQ / 2. Output the factor Q of N.

Later, in Section 5.2 of their paper, Gowers and Wagstaff give instructions for using a multiplier m to help find a factor of N, changing the algorithm given above as follows:

Change Step 1 to this: Read the odd positive integer N to be factored. If N is the square of an integer, output the square root and stop. If mN ≡ 1 mod 4, then set D ← 2mN; otherwise, set DmN. In any case, set S ← ⌊√D⌋, Q̂ ← 1, PS, QDP · P, L ← ⌊2√2√D⌋, B ← 2 · L, and i ← 0.

Change Step 2b to: Let gQ / gcd(Q, 2m). If gL, put (g, P mod g) onto the QUEUE. [ Note: If you're reading the linked version of the Gowers and Wagstaff paper, this formula is erroneously given as (g, B mod g); that error cost me several hours of pain. ]

In Step 5, replace Q by Q / gcd(Q, 2m) before the factor Q is output.

To find a factor, first use multiplier m = 1. If that fails, try another multiplier from the set 3, 5, 7, 11, 3 · 5 = 15, 3 · 7 = 21, 3 · 11 = 33, 5 · 7 = 35, 5 · 11 = 55, 7 · 11 = 77, 3 · 5 · 7 = 105, 3 · 5 · 11 = 165, 3 · 7 · 11 = 231, 5 · 7 · 11 = 385, 3 · 5 · 7 · 11 = 1155, or any other squarefree product of small odd primes.

Your task is to implement SQUFOF with multipliers. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

by programmingpraxis at July 08, 2014 09:00 AM

July 07, 2014

GNU Guix

GNU dmd 0.2 released

GNU dmd 0.2 has been released. It provides new features such as the ability to load new definitions for existing services, as well as bug fixes and an improved manual.

GNU dmd is a dependency-based service manager meant to be used as the init system in GNU. It is written in Guile Scheme.

by Ludovic Courtès at July 07, 2014 10:04 PM

Jeremy Kun

When Greedy Algorithms are Good Enough: Submodularity and the (1 – 1/e)-Approximation

Greedy algorithms are among the simplest and most intuitive algorithms known to humans. Their name essentially gives their description: do the thing that looks best right now, and repeat until nothing looks good anymore or you’re forced to stop. Some of the best situations in computer science are also when greedy algorithms are optimal or near-optimal. There is a […]

by Jeremy Kun at July 07, 2014 03:00 PM

July 04, 2014

Programming Praxis

Minimum And Maximum

We have today a two-part interview question:

First, write a program that takes an array of integers and returns the minimum and maximum elements of the array. Use as few comparisons as possible.

Second, write a program that takes an array of integers and returns the second minimum and maximum elements of the array; that is, the second-smallest element and the second-largest element. Again, use as few comparisons as possible.

Your task is to write the two programs described above, and state the number of comparisons each makes in the general case; be sure they are proof against malicious input. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

by programmingpraxis at July 04, 2014 09:00 AM

July 02, 2014

Jeremy Kun

AMS Network Science Mathematical Research Community

I don’t usually write promotional posts because I don’t enjoy reading them as much as I enjoy reading the technical posts. But I know that a lot of early graduate students and undergraduates read my blog, and this would be of interest to many of them. I just got back from Utah yesterday where I attended […]

by Jeremy Kun at July 02, 2014 06:47 PM