Implementing MD5 in AWK
It is only recently that I've discovered the Advent of Code, a coding challenge happening on December (though available all year) since 2015 and it is quite unique. I've set out to solve the 2015 puzzles in AWK (you can find my solutions here). Very soon one of them required to compute (a lot of) MD5 hashes. After some research, I found out there are basically three solutions:
 Using a pipe. This is expected to be slow, relying on an external program, but very simple to implement.
 Write a dynamic extension in C or C++. This should yield very good performances, but would be limited to GNU Awk and tedious to deploy as it would require compiling.
 Implement MD5 in AWK. This is expected to be portable and very slow, but more importantly a lot of fun.
In this post I'll go through the implementation of MD5 in AWK by first writing a GNU Awk version (because of some of its features), then porting to AWK, and finally optimizing. The main challenge is that MD5 is basically a lot of bitwise operations on 32bits integers and AWK has neither. If that sounds interesting to you, sit comfortably and read on — this is going to be a long ride.
setup()
In case you missed it I wrote all this to solve a coding challenge — for fun. Theses implementations have a lot of limitations and you should not use them for anything remotely serious (well, you should not use MD5 for anything remotely serious nowadays anyway).
Let start with a basic template that mimic the openssl md5
command output:
base.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22  { # NOTE: remember the input files inorder in the `files' array. if (nfiles == 0  files[nfiles] != FILENAME) files[++nfiles] = FILENAME; # XXX: only work with files ending with a newline, this is an OK # limitation since it is required by POSIX. content[FILENAME] = content[FILENAME] $0 "\n"; } END { # go over all the files inorder. for (i = 1; i <= nfiles; i++) { fn = files[i]; # ala `openssl md5' output. printf("MD5(%s)= %s\n", fn, md5(content[fn])); } } # our md5 implementation function md5(input) { # TODO } 
The script read every lines from the files given on the command line (or
consume the standard input) and save them into the content
variable. Because we want to display the computed checksum in the same order
as given, we have to do a little dance with files
and
nfiles
. If you are not familiar with AWK, note that
arrays are 1indexed. We've followed this practice for
files
. Also all uninitialized variables (including array values)
are 0 by default, so we don't need to explicitly initialize
nfiles
.
From now on the focus will be on the md5()
function. We'll be
following the RFC
stepbystep. It include a reference implementation in C,
Md5.c, that you should also find
here. It was really
handy to print debug the intermediate states etc. Now, this implementation
has been written a long time ago and need a trivial patch:
Md5.c.patch
1 2 3 4 5 6 7 8 9 10 11 12 13 14   Md5.c.original 20170114 21:00:21.894542600 +0100 +++ Md5.c 20170114 21:01:02.178086795 +0100 @@ 37,8 +37,10 @@ ********************************************************************** */ +#include <stdint.h> + /* typedef a 32 bit type */ typedef unsigned long int UINT4; +typedef uint32_t UINT4; /* Data structure for MD5 (Message Digest) computation */ typedef struct { 
massaging the input
For further processing we need the message to be represented as 32bits words:
Remember that arrays are 1indexed in AWK. Thus, chars
must be
walked from 1
until nbytes
included.
First, we split the input
string as an array of characters. One
byte per character is assumed (ASCII) but this could work with nonASCII text
too (you may need to LC_ALL=C
because
GNU Awk is locale aware). Then, we "convert" each
character as a byte using ord()
. Whereas many languages have an
ord()
function AWK doesn't, but fortunately we can find an
implementation from
the GNU Awk manual.
Let's hack it a bit and we get:
_ord_init()
deserve some clarification if you're not used to
AWK: the caller is not expected to actually provide a value for
i
. All variables in AWK are global except function
parameters. i
is declared as parameter only to create a
variable that is local to the function, that way _ord_init()
doesn't change the global state. While very unusual in "modern" languages
this is a very common AWK trick, so common in fact that it is documented in
the manpage:
Parameters are passed by value if scalar and by reference if array name; functions may be called recursively. Parameters are local to the function; all other variables are global. Thus local variables may be created by providing excess parameters in the function definition.
Finally, let's convert our bytes
array into 32bits words.
Remember that AWK doesn't have bitwise operators? Fortunately for us,
GNU Awk
has builtin bitwise functions.
So we'll start off with GNU Awk and then figure
out a way to implement bitwise functions later on.
Notice that we've decided to break from the 1index array convention of AWK
for the words
array here. It will make some computation less
awkward and allow us to translate "smoothly" from
Md5.c (the reference C implementation) later
on. ^{↓ Padding and Length}
How we order the input bytes into each word here is important:
[…] a sequence of bytes can be interpreted as a sequence of 32bit words, where each consecutive group of four bytes is interpreted as a word with the loworder (least significant) byte given first.
Thus we use the first two bytes to compose the "low" 16 bits part of the word, and the third and fourth bytes to compose the "high" 16 bits part.
If we wrap up we get:
step0.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57  BEGIN { _ord_init(); } { # NOTE: remember the input files inorder in the `files' array. if (nfiles == 0  files[nfiles] != FILENAME) files[++nfiles] = FILENAME; # XXX: only work with files ending with a newline, this is an OK # limitation since it is required by POSIX. content[FILENAME] = content[FILENAME] $0 "\n"; } END { # go over all the files inorder. for (i = 1; i <= nfiles; i++) { fn = files[i]; # ala `openssl md5' output. printf("MD5(%s)= %s\n", fn, md5(content[fn])); } } # our md5 implementation function md5(input, nbytes, chars, i, bytes, hi, lo, words, nwords) { # convert the input into an array of bytes using ord() on each # character. nbytes = split(input, chars, ""); for (i = 1; i <= nbytes; i++) bytes[i] = ord(chars[i]); # convert the array of bytes into an array of 32bits words. # NOTE: words is 0indexed. for (i = 1; i <= nbytes; i += 4) { hi = or(lshift(bytes[i + 3], 8), bytes[i + 2]); lo = or(lshift(bytes[i + 1], 8), bytes[i + 0]); words[nwords++] = or(lshift(hi, 16), lo); } # XXX: debug words and exit: for (i = 0; i < nwords; i++) printf "%4d: %08x\n", i, words[i]; exit; # NOTREACHED } # from https://www.gnu.org/software/gawk/manual/html_node/OrdinalFunctions.html function _ord_init( i) { for (i = 0; i < 256; i++) _ord_[sprintf("%c", i)] = i; } function ord(s) { # only first character is of interest return _ord_[substr(s, 1, 1)]; } 
diff u base.awk step0.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50   base.awk 20170220 18:24:03.945445000 +0100 +++ step0.awk 20170220 18:24:03.945454000 +0100 @@ 1,3 +1,7 @@ +BEGIN { + _ord_init(); +} + { # NOTE: remember the input files inorder in the `files' array. if (nfiles == 0  files[nfiles] != FILENAME) @@ 17,6 +21,37 @@ } # our md5 implementation function md5(input) {  # TODO +function md5(input, nbytes, chars, i, bytes, hi, lo, words, nwords) { + # convert the input into an array of bytes using ord() on each + # character. + nbytes = split(input, chars, ""); + for (i = 1; i <= nbytes; i++) + bytes[i] = ord(chars[i]); + + # convert the array of bytes into an array of 32bits words. + # NOTE: words is 0indexed. + for (i = 1; i <= nbytes; i += 4) { + hi = or(lshift(bytes[i + 3], 8), bytes[i + 2]); + lo = or(lshift(bytes[i + 1], 8), bytes[i + 0]); + words[nwords++] = or(lshift(hi, 16), lo); + } + + # XXX: debug words and exit: + for (i = 0; i < nwords; i++) + printf "%4d: %08x\n", i, words[i]; + exit; + # NOTREACHED +} + +# from https://www.gnu.org/software/gawk/manual/html_node/OrdinalFunctions.html +function _ord_init( i) +{ + for (i = 0; i < 256; i++) + _ord_[sprintf("%c", i)] = i; +} + +function ord(s) +{ + # only first character is of interest + return _ord_[substr(s, 1, 1)]; } 
We've added the variables used in the md5()
function body to the
signature at line 27 and a bit of debug code (lines 41 — 45). Let's make a
quick test with 1234abcd\n
:
This should be familiar (see ascii(7)), In
particular 0a
is the ASCII code for \n
.
Let's run our script:
Alright! The byte order in each words is as expected: loworder (least significant) byte given first.
Padding Bits and Length
We'll go through "3.1 Step 1. Append Padding Bits" and "3.2 Step 2. Append Length" of the RFC in one go, as I feel like it's easier to visualize and understand both at once. Basically, we need to concatenate
 the input message,
 some padding,
 the input message length.
So that the result looks like this:
The padding is at least a '1' bit, and then as many as '0' bits needed so that the full result can be split into chunks of 512 bits. The trick is to take into account that 64 bits will be "used" by the length.
The loop at the end may be confusing if you don't remember that in AWK all uninitialized variables (including array values) are 0 by default. As a result, we "only" need to increment the array size to "append zeros" to the array.
The length is represented on 64 bits as two 32bits words (lower word given first).
We truncate the length if it is bigger than 2^64 as the RFC mandate, just in case you feel like feeding more than 2097152 TiB to AWK.
Let's wrap up again:
step1.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77  BEGIN { _ord_init(); } { # NOTE: remember the input files inorder in the `files' array. if (nfiles == 0  files[nfiles] != FILENAME) files[++nfiles] = FILENAME; # XXX: only work with files ending with a newline, this is an OK # limitation since it is required by POSIX. content[FILENAME] = content[FILENAME] $0 "\n"; } END { # go over all the files inorder. for (i = 1; i <= nfiles; i++) { fn = files[i]; # ala `openssl md5' output. printf("MD5(%s)= %s\n", fn, md5(content[fn])); } } # our md5 implementation function md5(input, nbytes, chars, i, bytes, hi, lo, words, nwords) { # convert the input into an array of bytes using ord() on each # character. nbytes = split(input, chars, ""); for (i = 1; i <= nbytes; i++) bytes[i] = ord(chars[i]); # convert the array of bytes into an array of 32bits words. # NOTE: words is 0indexed. for (i = 1; i <= nbytes; i += 4) { hi = or(lshift(bytes[i + 3], 8), bytes[i + 2]); lo = or(lshift(bytes[i + 1], 8), bytes[i + 0]); words[nwords++] = or(lshift(hi, 16), lo); } # Step 1. Append Padding Bits if (nbytes % 4 == 0) { # the input size is congruent modulo 32, we need a new word to # store the first '1' padding bit. words[nwords++] = 0x80; } else { # append a '1' bit in the byte just after the last input byte. words[nwords  1] = or(words[nwords  1], lshift(0x80, (nbytes % 4) * 8)); } # "fill" the remaining bytes with 0 until we're just shy two words of # having 16Word Blocks. while ((nwords % 16) != 14) nwords++; # Step 2. Append Length hi = rshift(nbytes * 8, 32); lo = (nbytes * 8)  lshift(hi, 32); words[nwords++] = lo; words[nwords++] = and(hi, 0xffffffff); # truncate to 32 bits # XXX: debug words and exit: for (i = 0; i < nwords; i++) printf "%4d: %08x\n", i, words[i]; exit; # NOTREACHED } # from https://www.gnu.org/software/gawk/manual/html_node/OrdinalFunctions.html function _ord_init( i) { for (i = 0; i < 256; i++) _ord_[sprintf("%c", i)] = i; } function ord(s) { # only first character is of interest return _ord_[substr(s, 1, 1)]; } 
diff u step0.awk step1.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29   step0.awk 20170220 18:24:03.945454000 +0100 +++ step1.awk 20170220 18:24:03.945458000 +0100 @@ 36,6 +36,26 @@ words[nwords++] = or(lshift(hi, 16), lo); } + # Step 1. Append Padding Bits + if (nbytes % 4 == 0) { + # the input size is congruent modulo 32, we need a new word to + # store the first '1' padding bit. + words[nwords++] = 0x80; + } else { + # append a '1' bit in the byte just after the last input byte. + words[nwords  1] = or(words[nwords  1], lshift(0x80, (nbytes % 4) * 8)); + } + # "fill" the remaining bytes with 0 until we're just shy two words of + # having 16Word Blocks. + while ((nwords % 16) != 14) + nwords++; + + # Step 2. Append Length + hi = rshift(nbytes * 8, 32); + lo = (nbytes * 8)  lshift(hi, 32); + words[nwords++] = lo; + words[nwords++] = and(hi, 0xffffffff); # truncate to 32 bits + # XXX: debug words and exit: for (i = 0; i < nwords; i++) printf "%4d: %08x\n", i, words[i]; 
Using our previous 1234abcd\n
example, the result should look
like:
Let's test our script:
We get as expected 16 32bits words (i.e. 512 bits) in total,
0x80
right after \n
and some zeros to complete the
padding, and finally the length in the two last words.
MD Buffer initialization, Message processing and Output
I will pass quickly through the MD5 rounds and setup, highlighting the challenging parts.
Step 3. is a trivial translation from Md5.c (lines 150 — 155).
Step 4. was easy too. The main loop and x
setup was
adapted directly from the RFC, while the rounds were inspired from
Md5.c. Here having words
,
state
and x
0indexed (despite the AWK 1indexed
arrays tradition) was key to translate almost verbatim the processing loop.
^{↑ massaging the input}
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85  # Step 4. Process Message in 16Word Blocks # Process each 16word block. for (i = 0; i < nwords; i += 16) { # Copy block i into x. for (j = 0; j < 16; j++) x[j] = words[i + j]; a = state[0]; b = state[1]; c = state[2]; d = state[3]; # Round 1 a = FF(a, b, c, d, x[ 0], S11, 0xd76aa478); d = FF(d, a, b, c, x[ 1], S12, 0xe8c7b756); c = FF(c, d, a, b, x[ 2], S13, 0x242070db); b = FF(b, c, d, a, x[ 3], S14, 0xc1bdceee); a = FF(a, b, c, d, x[ 4], S11, 0xf57c0faf); d = FF(d, a, b, c, x[ 5], S12, 0x4787c62a); c = FF(c, d, a, b, x[ 6], S13, 0xa8304613); b = FF(b, c, d, a, x[ 7], S14, 0xfd469501); a = FF(a, b, c, d, x[ 8], S11, 0x698098d8); d = FF(d, a, b, c, x[ 9], S12, 0x8b44f7af); c = FF(c, d, a, b, x[10], S13, 0xffff5bb1); b = FF(b, c, d, a, x[11], S14, 0x895cd7be); a = FF(a, b, c, d, x[12], S11, 0x6b901122); d = FF(d, a, b, c, x[13], S12, 0xfd987193); c = FF(c, d, a, b, x[14], S13, 0xa679438e); b = FF(b, c, d, a, x[15], S14, 0x49b40821); # Round 2 a = GG(a, b, c, d, x[ 1], S21, 0xf61e2562); d = GG(d, a, b, c, x[ 6], S22, 0xc040b340); c = GG(c, d, a, b, x[11], S23, 0x265e5a51); b = GG(b, c, d, a, x[ 0], S24, 0xe9b6c7aa); a = GG(a, b, c, d, x[ 5], S21, 0xd62f105d); d = GG(d, a, b, c, x[10], S22, 0x2441453); c = GG(c, d, a, b, x[15], S23, 0xd8a1e681); b = GG(b, c, d, a, x[ 4], S24, 0xe7d3fbc8); a = GG(a, b, c, d, x[ 9], S21, 0x21e1cde6); d = GG(d, a, b, c, x[14], S22, 0xc33707d6); c = GG(c, d, a, b, x[ 3], S23, 0xf4d50d87); b = GG(b, c, d, a, x[ 8], S24, 0x455a14ed); a = GG(a, b, c, d, x[13], S21, 0xa9e3e905); d = GG(d, a, b, c, x[ 2], S22, 0xfcefa3f8); c = GG(c, d, a, b, x[ 7], S23, 0x676f02d9); b = GG(b, c, d, a, x[12], S24, 0x8d2a4c8a); # Round 3 a = HH(a, b, c, d, x[ 5], S31, 0xfffa3942); d = HH(d, a, b, c, x[ 8], S32, 0x8771f681); c = HH(c, d, a, b, x[11], S33, 0x6d9d6122); b = HH(b, c, d, a, x[14], S34, 0xfde5380c); a = HH(a, b, c, d, x[ 1], S31, 0xa4beea44); d = HH(d, a, b, c, x[ 4], S32, 0x4bdecfa9); c = HH(c, d, a, b, x[ 7], S33, 0xf6bb4b60); b = HH(b, c, d, a, x[10], S34, 0xbebfbc70); a = HH(a, b, c, d, x[13], S31, 0x289b7ec6); d = HH(d, a, b, c, x[ 0], S32, 0xeaa127fa); c = HH(c, d, a, b, x[ 3], S33, 0xd4ef3085); b = HH(b, c, d, a, x[ 6], S34, 0x4881d05); a = HH(a, b, c, d, x[ 9], S31, 0xd9d4d039); d = HH(d, a, b, c, x[12], S32, 0xe6db99e5); c = HH(c, d, a, b, x[15], S33, 0x1fa27cf8); b = HH(b, c, d, a, x[ 2], S34, 0xc4ac5665); # Round 4 a = II(a, b, c, d, x[ 0], S41, 0xf4292244); d = II(d, a, b, c, x[ 7], S42, 0x432aff97); c = II(c, d, a, b, x[14], S43, 0xab9423a7); b = II(b, c, d, a, x[ 5], S44, 0xfc93a039); a = II(a, b, c, d, x[12], S41, 0x655b59c3); d = II(d, a, b, c, x[ 3], S42, 0x8f0ccc92); c = II(c, d, a, b, x[10], S43, 0xffeff47d); b = II(b, c, d, a, x[ 1], S44, 0x85845dd1); a = II(a, b, c, d, x[ 8], S41, 0x6fa87e4f); d = II(d, a, b, c, x[15], S42, 0xfe2ce6e0); c = II(c, d, a, b, x[ 6], S43, 0xa3014314); b = II(b, c, d, a, x[13], S44, 0x4e0811a1); a = II(a, b, c, d, x[ 4], S41, 0xf7537e82); d = II(d, a, b, c, x[11], S42, 0xbd3af235); c = II(c, d, a, b, x[ 2], S43, 0x2ad7d2bb); b = II(b, c, d, a, x[ 9], S44, 0xeb86d391); state[0] = mod32bits(state[0] + a); state[1] = mod32bits(state[1] + b); state[2] = mod32bits(state[2] + c); state[3] = mod32bits(state[3] + d); } 
We'll need to setup the shift constants (S11
, S12
,
S13
and so on) and also implement the round functions (namely
FF
, GG
, HH
, and II
).
An interesting point is that when state
is updated (lines 81 —
84) we have to use "modulo2^32 addition". As it will be the case
for all additions, we'll write a mod32bits()
helper function. If
you remember we already had to implement "32bits truncation" for the "high
32bits" part of length:
This will do fine for now:
Output
Similarly to the last part of MD5Final()
, the 32bits
state
words is "converted" into a digest
of 16
bytes. Then, the hexadecimal string representation of the hash is built and
returned.
It could be made shorter (and faster) but I've left it in two steps as I find it more readable. The performance penalty should not be significant as this part is O(1), in other words it run in constant time regardless of the input message length.
Round functions
Again, overall easy translation from the reference implementation:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51  function F(x, y, z) { return or(and(x, y), and(not(x), z)); } function G(x, y, z) { return or(and(x, z), and(y, not(z))); } function H(x, y, z) { return xor(x, xor(y, z)); } function I(x, y, z) { return xor(y, or(x, not(z))); } function FF(a, b, c, d, x, s, ac) { a = mod32bits(a + F(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function GG(a, b, c, d, x, s, ac) { a = mod32bits(a + G(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function HH(a, b, c, d, x, s, ac) { a = mod32bits(a + H(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function II(a, b, c, d, x, s, ac) { a = mod32bits(a + I(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function ROTATE_LEFT(x, n) { return or(mod32bits(lshift(x, n)), rshift(x, 32  n)); } function not(x) { return mod32bits(compl(x)); } 
We need to be careful with both lshift()
(line 46) and
compl()
(line 50) as their return value may be greater than
2^32:
For all of these functions, first the doubleprecision floatingpoint value is converted to the widest C unsigned integer type, then the bitwise operation is performed. If the result cannot be represented exactly as a C double, leading nonzero bits are removed one by one until it can be represented exactly. The result is then converted back into a C double. (If you don’t understand this paragraph, don’t worry about it.)
"the widest C unsigned integer type" is likely to be represented on more than
32bits, and so we're ensuring in both ROTATE_LEFT()
and
not()
to truncate the result to 32bits.
Shift constants
Simple and inspired by _ord_init()
that we previously adapted:
Finally it should work! Here is the final GNU Awk code:
step2.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239  BEGIN { _ord_init(); _md5_init(); } { # NOTE: remember the input files inorder in the `files' array. if (nfiles == 0  files[nfiles] != FILENAME) files[++nfiles] = FILENAME; # XXX: only work with files ending with a newline, this is an OK # limitation since it is required by POSIX. content[FILENAME] = content[FILENAME] $0 "\n"; } END { # go over all the files inorder. for (i = 1; i <= nfiles; i++) { fn = files[i]; # ala `openssl md5' output. printf("MD5(%s)= %s\n", fn, md5(content[fn])); } } # our md5 implementation function md5(input, nbytes, chars, i, bytes, hi, lo, words, nwords, state, a, b, c, d, j, x, digest, ret) { # convert the input into an array of bytes using ord() on each # character. nbytes = split(input, chars, ""); for (i = 1; i <= nbytes; i++) bytes[i] = ord(chars[i]); # convert the array of bytes into an array of 32bits words. # NOTE: words is 0indexed. for (i = 1; i <= nbytes; i += 4) { hi = or(lshift(bytes[i + 3], 8), bytes[i + 2]); lo = or(lshift(bytes[i + 1], 8), bytes[i + 0]); words[nwords++] = or(lshift(hi, 16), lo); } # Step 1. Append Padding Bits if (nbytes % 4 == 0) { # the input size is congruent modulo 32, we need a new word to # store the first '1' padding bit. words[nwords++] = 0x80; } else { # append a '1' bit in the byte just after the last input byte. words[nwords  1] = or(words[nwords  1], lshift(0x80, (nbytes % 4) * 8)); } # "fill" the remaining bytes with 0 until we're just shy two words of # having 16Word Blocks. while ((nwords % 16) != 14) nwords++; # Step 2. Append Length hi = rshift(nbytes * 8, 32); lo = (nbytes * 8)  lshift(hi, 32); words[nwords++] = lo; words[nwords++] = mod32bits(hi); # truncate to 32 bits # Step 3. Initialize MD Buffer state[0] = 0x67452301; state[1] = 0xefcdab89; state[2] = 0x98badcfe; state[3] = 0x10325476; # Step 4. Process Message in 16Word Blocks # Process each 16word block. for (i = 0; i < nwords; i += 16) { # Copy block i into x. for (j = 0; j < 16; j++) x[j] = words[i + j]; a = state[0]; b = state[1]; c = state[2]; d = state[3]; # Round 1 a = FF(a, b, c, d, x[ 0], S11, 0xd76aa478); d = FF(d, a, b, c, x[ 1], S12, 0xe8c7b756); c = FF(c, d, a, b, x[ 2], S13, 0x242070db); b = FF(b, c, d, a, x[ 3], S14, 0xc1bdceee); a = FF(a, b, c, d, x[ 4], S11, 0xf57c0faf); d = FF(d, a, b, c, x[ 5], S12, 0x4787c62a); c = FF(c, d, a, b, x[ 6], S13, 0xa8304613); b = FF(b, c, d, a, x[ 7], S14, 0xfd469501); a = FF(a, b, c, d, x[ 8], S11, 0x698098d8); d = FF(d, a, b, c, x[ 9], S12, 0x8b44f7af); c = FF(c, d, a, b, x[10], S13, 0xffff5bb1); b = FF(b, c, d, a, x[11], S14, 0x895cd7be); a = FF(a, b, c, d, x[12], S11, 0x6b901122); d = FF(d, a, b, c, x[13], S12, 0xfd987193); c = FF(c, d, a, b, x[14], S13, 0xa679438e); b = FF(b, c, d, a, x[15], S14, 0x49b40821); # Round 2 a = GG(a, b, c, d, x[ 1], S21, 0xf61e2562); d = GG(d, a, b, c, x[ 6], S22, 0xc040b340); c = GG(c, d, a, b, x[11], S23, 0x265e5a51); b = GG(b, c, d, a, x[ 0], S24, 0xe9b6c7aa); a = GG(a, b, c, d, x[ 5], S21, 0xd62f105d); d = GG(d, a, b, c, x[10], S22, 0x2441453); c = GG(c, d, a, b, x[15], S23, 0xd8a1e681); b = GG(b, c, d, a, x[ 4], S24, 0xe7d3fbc8); a = GG(a, b, c, d, x[ 9], S21, 0x21e1cde6); d = GG(d, a, b, c, x[14], S22, 0xc33707d6); c = GG(c, d, a, b, x[ 3], S23, 0xf4d50d87); b = GG(b, c, d, a, x[ 8], S24, 0x455a14ed); a = GG(a, b, c, d, x[13], S21, 0xa9e3e905); d = GG(d, a, b, c, x[ 2], S22, 0xfcefa3f8); c = GG(c, d, a, b, x[ 7], S23, 0x676f02d9); b = GG(b, c, d, a, x[12], S24, 0x8d2a4c8a); # Round 3 a = HH(a, b, c, d, x[ 5], S31, 0xfffa3942); d = HH(d, a, b, c, x[ 8], S32, 0x8771f681); c = HH(c, d, a, b, x[11], S33, 0x6d9d6122); b = HH(b, c, d, a, x[14], S34, 0xfde5380c); a = HH(a, b, c, d, x[ 1], S31, 0xa4beea44); d = HH(d, a, b, c, x[ 4], S32, 0x4bdecfa9); c = HH(c, d, a, b, x[ 7], S33, 0xf6bb4b60); b = HH(b, c, d, a, x[10], S34, 0xbebfbc70); a = HH(a, b, c, d, x[13], S31, 0x289b7ec6); d = HH(d, a, b, c, x[ 0], S32, 0xeaa127fa); c = HH(c, d, a, b, x[ 3], S33, 0xd4ef3085); b = HH(b, c, d, a, x[ 6], S34, 0x4881d05); a = HH(a, b, c, d, x[ 9], S31, 0xd9d4d039); d = HH(d, a, b, c, x[12], S32, 0xe6db99e5); c = HH(c, d, a, b, x[15], S33, 0x1fa27cf8); b = HH(b, c, d, a, x[ 2], S34, 0xc4ac5665); # Round 4 a = II(a, b, c, d, x[ 0], S41, 0xf4292244); d = II(d, a, b, c, x[ 7], S42, 0x432aff97); c = II(c, d, a, b, x[14], S43, 0xab9423a7); b = II(b, c, d, a, x[ 5], S44, 0xfc93a039); a = II(a, b, c, d, x[12], S41, 0x655b59c3); d = II(d, a, b, c, x[ 3], S42, 0x8f0ccc92); c = II(c, d, a, b, x[10], S43, 0xffeff47d); b = II(b, c, d, a, x[ 1], S44, 0x85845dd1); a = II(a, b, c, d, x[ 8], S41, 0x6fa87e4f); d = II(d, a, b, c, x[15], S42, 0xfe2ce6e0); c = II(c, d, a, b, x[ 6], S43, 0xa3014314); b = II(b, c, d, a, x[13], S44, 0x4e0811a1); a = II(a, b, c, d, x[ 4], S41, 0xf7537e82); d = II(d, a, b, c, x[11], S42, 0xbd3af235); c = II(c, d, a, b, x[ 2], S43, 0x2ad7d2bb); b = II(b, c, d, a, x[ 9], S44, 0xeb86d391); state[0] = mod32bits(state[0] + a); state[1] = mod32bits(state[1] + b); state[2] = mod32bits(state[2] + c); state[3] = mod32bits(state[3] + d); } for (i = j = 0; j < 16; j += 4) { digest[j + 0] = and(state[i], 0xff); digest[j + 1] = and(rshift(state[i], 8), 0xff); digest[j + 2] = and(rshift(state[i], 16), 0xff); digest[j + 3] = and(rshift(state[i++], 24), 0xff); } for (i = 0; i < 16; i++) ret = sprintf("%s%02x", ret, digest[i]); return ret; } function F(x, y, z) { return or(and(x, y), and(not(x), z)); } function G(x, y, z) { return or(and(x, z), and(y, not(z))); } function H(x, y, z) { return xor(x, xor(y, z)); } function I(x, y, z) { return xor(y, or(x, not(z))); } function FF(a, b, c, d, x, s, ac) { a = mod32bits(a + F(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function GG(a, b, c, d, x, s, ac) { a = mod32bits(a + G(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function HH(a, b, c, d, x, s, ac) { a = mod32bits(a + H(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function II(a, b, c, d, x, s, ac) { a = mod32bits(a + I(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function ROTATE_LEFT(x, n) { return or(mod32bits(lshift(x, n)), rshift(x, 32  n)); } function mod32bits(x) { return and(x, 0xffffffff); } function not(x) { return mod32bits(compl(x)); } # from https://www.gnu.org/software/gawk/manual/html_node/OrdinalFunctions.html function _ord_init( i) { for (i = 0; i < 256; i++) _ord_[sprintf("%c", i)] = i; } function ord(s) { # only first character is of interest return _ord_[substr(s, 1, 1)]; } function _md5_init() { # MD5 shift constants setup. S11 = 7; S12 = 12; S13 = 17; S14 = 22; S21 = 5; S22 = 9; S23 = 14; S24 = 20; S31 = 4; S32 = 11; S33 = 16; S34 = 23; S41 = 6; S42 = 10; S43 = 15; S44 = 21; } 
diff u step1.awk step2.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206   step1.awk 20170621 09:43:06.590286000 +0200 +++ step2.awk 20191024 09:35:57.779452000 +0200 @@ 1,5 +1,6 @@ BEGIN { _ord_init(); + _md5_init(); } { @@ 21,7 +22,8 @@ } # our md5 implementation function md5(input, nbytes, chars, i, bytes, hi, lo, words, nwords) { +function md5(input, nbytes, chars, i, bytes, hi, lo, words, nwords, state, + a, b, c, d, j, x, digest, ret) { # convert the input into an array of bytes using ord() on each # character. nbytes = split(input, chars, ""); @@ 54,15 +56,167 @@ hi = rshift(nbytes * 8, 32); lo = (nbytes * 8)  lshift(hi, 32); words[nwords++] = lo;  words[nwords++] = and(hi, 0xffffffff); # truncate to 32 bits + words[nwords++] = mod32bits(hi); # truncate to 32 bits  # XXX: debug words and exit:  for (i = 0; i < nwords; i++)  printf "%4d: %08x\n", i, words[i];  exit;  # NOTREACHED + # Step 3. Initialize MD Buffer + state[0] = 0x67452301; + state[1] = 0xefcdab89; + state[2] = 0x98badcfe; + state[3] = 0x10325476; + + # Step 4. Process Message in 16Word Blocks + # Process each 16word block. + for (i = 0; i < nwords; i += 16) { + # Copy block i into x. + for (j = 0; j < 16; j++) + x[j] = words[i + j]; + a = state[0]; b = state[1]; c = state[2]; d = state[3]; + + # Round 1 + a = FF(a, b, c, d, x[ 0], S11, 0xd76aa478); + d = FF(d, a, b, c, x[ 1], S12, 0xe8c7b756); + c = FF(c, d, a, b, x[ 2], S13, 0x242070db); + b = FF(b, c, d, a, x[ 3], S14, 0xc1bdceee); + a = FF(a, b, c, d, x[ 4], S11, 0xf57c0faf); + d = FF(d, a, b, c, x[ 5], S12, 0x4787c62a); + c = FF(c, d, a, b, x[ 6], S13, 0xa8304613); + b = FF(b, c, d, a, x[ 7], S14, 0xfd469501); + a = FF(a, b, c, d, x[ 8], S11, 0x698098d8); + d = FF(d, a, b, c, x[ 9], S12, 0x8b44f7af); + c = FF(c, d, a, b, x[10], S13, 0xffff5bb1); + b = FF(b, c, d, a, x[11], S14, 0x895cd7be); + a = FF(a, b, c, d, x[12], S11, 0x6b901122); + d = FF(d, a, b, c, x[13], S12, 0xfd987193); + c = FF(c, d, a, b, x[14], S13, 0xa679438e); + b = FF(b, c, d, a, x[15], S14, 0x49b40821); + + # Round 2 + a = GG(a, b, c, d, x[ 1], S21, 0xf61e2562); + d = GG(d, a, b, c, x[ 6], S22, 0xc040b340); + c = GG(c, d, a, b, x[11], S23, 0x265e5a51); + b = GG(b, c, d, a, x[ 0], S24, 0xe9b6c7aa); + a = GG(a, b, c, d, x[ 5], S21, 0xd62f105d); + d = GG(d, a, b, c, x[10], S22, 0x2441453); + c = GG(c, d, a, b, x[15], S23, 0xd8a1e681); + b = GG(b, c, d, a, x[ 4], S24, 0xe7d3fbc8); + a = GG(a, b, c, d, x[ 9], S21, 0x21e1cde6); + d = GG(d, a, b, c, x[14], S22, 0xc33707d6); + c = GG(c, d, a, b, x[ 3], S23, 0xf4d50d87); + b = GG(b, c, d, a, x[ 8], S24, 0x455a14ed); + a = GG(a, b, c, d, x[13], S21, 0xa9e3e905); + d = GG(d, a, b, c, x[ 2], S22, 0xfcefa3f8); + c = GG(c, d, a, b, x[ 7], S23, 0x676f02d9); + b = GG(b, c, d, a, x[12], S24, 0x8d2a4c8a); + + # Round 3 + a = HH(a, b, c, d, x[ 5], S31, 0xfffa3942); + d = HH(d, a, b, c, x[ 8], S32, 0x8771f681); + c = HH(c, d, a, b, x[11], S33, 0x6d9d6122); + b = HH(b, c, d, a, x[14], S34, 0xfde5380c); + a = HH(a, b, c, d, x[ 1], S31, 0xa4beea44); + d = HH(d, a, b, c, x[ 4], S32, 0x4bdecfa9); + c = HH(c, d, a, b, x[ 7], S33, 0xf6bb4b60); + b = HH(b, c, d, a, x[10], S34, 0xbebfbc70); + a = HH(a, b, c, d, x[13], S31, 0x289b7ec6); + d = HH(d, a, b, c, x[ 0], S32, 0xeaa127fa); + c = HH(c, d, a, b, x[ 3], S33, 0xd4ef3085); + b = HH(b, c, d, a, x[ 6], S34, 0x4881d05); + a = HH(a, b, c, d, x[ 9], S31, 0xd9d4d039); + d = HH(d, a, b, c, x[12], S32, 0xe6db99e5); + c = HH(c, d, a, b, x[15], S33, 0x1fa27cf8); + b = HH(b, c, d, a, x[ 2], S34, 0xc4ac5665); + + # Round 4 + a = II(a, b, c, d, x[ 0], S41, 0xf4292244); + d = II(d, a, b, c, x[ 7], S42, 0x432aff97); + c = II(c, d, a, b, x[14], S43, 0xab9423a7); + b = II(b, c, d, a, x[ 5], S44, 0xfc93a039); + a = II(a, b, c, d, x[12], S41, 0x655b59c3); + d = II(d, a, b, c, x[ 3], S42, 0x8f0ccc92); + c = II(c, d, a, b, x[10], S43, 0xffeff47d); + b = II(b, c, d, a, x[ 1], S44, 0x85845dd1); + a = II(a, b, c, d, x[ 8], S41, 0x6fa87e4f); + d = II(d, a, b, c, x[15], S42, 0xfe2ce6e0); + c = II(c, d, a, b, x[ 6], S43, 0xa3014314); + b = II(b, c, d, a, x[13], S44, 0x4e0811a1); + a = II(a, b, c, d, x[ 4], S41, 0xf7537e82); + d = II(d, a, b, c, x[11], S42, 0xbd3af235); + c = II(c, d, a, b, x[ 2], S43, 0x2ad7d2bb); + b = II(b, c, d, a, x[ 9], S44, 0xeb86d391); + + state[0] = mod32bits(state[0] + a); + state[1] = mod32bits(state[1] + b); + state[2] = mod32bits(state[2] + c); + state[3] = mod32bits(state[3] + d); + } + + for (i = j = 0; j < 16; j += 4) { + digest[j + 0] = and(state[i], 0xff); + digest[j + 1] = and(rshift(state[i], 8), 0xff); + digest[j + 2] = and(rshift(state[i], 16), 0xff); + digest[j + 3] = and(rshift(state[i++], 24), 0xff); + } + for (i = 0; i < 16; i++) + ret = sprintf("%s%02x", ret, digest[i]); + return ret; } +function F(x, y, z) { + return or(and(x, y), and(not(x), z)); +} + +function G(x, y, z) { + return or(and(x, z), and(y, not(z))); +} + +function H(x, y, z) { + return xor(x, xor(y, z)); +} + +function I(x, y, z) { + return xor(y, or(x, not(z))); +} + +function FF(a, b, c, d, x, s, ac) { + a = mod32bits(a + F(b, c, d) + x + ac); + a = ROTATE_LEFT(a, s); + a = mod32bits(a + b); + return a; +} + +function GG(a, b, c, d, x, s, ac) { + a = mod32bits(a + G(b, c, d) + x + ac); + a = ROTATE_LEFT(a, s); + a = mod32bits(a + b); + return a; +} + +function HH(a, b, c, d, x, s, ac) { + a = mod32bits(a + H(b, c, d) + x + ac); + a = ROTATE_LEFT(a, s); + a = mod32bits(a + b); + return a; +} + +function II(a, b, c, d, x, s, ac) { + a = mod32bits(a + I(b, c, d) + x + ac); + a = ROTATE_LEFT(a, s); + a = mod32bits(a + b); + return a; +} + +function ROTATE_LEFT(x, n) { + return or(mod32bits(lshift(x, n)), rshift(x, 32  n)); +} + +function mod32bits(x) { + return and(x, 0xffffffff); +} + +function not(x) { + return mod32bits(compl(x)); +} + # from https://www.gnu.org/software/gawk/manual/html_node/OrdinalFunctions.html function _ord_init( i) { @@ 74,4 +228,12 @@ { # only first character is of interest return _ord_[substr(s, 1, 1)]; +} + +function _md5_init() { + # MD5 shift constants setup. + S11 = 7; S12 = 12; S13 = 17; S14 = 22; + S21 = 5; S22 = 9; S23 = 14; S24 = 20; + S31 = 4; S32 = 11; S33 = 16; S34 = 23; + S41 = 6; S42 = 10; S43 = 15; S44 = 21; } 
Now, we need a bunch of ASCII files for testing. Arbitrarily, I'll use the SQLite sources files (i.e. files from the src/ directory) because it's enough data without being too much (about 6.5MB for 200k sloc).
We start by computing the MD5 sum for each file using openssl md5
(because we mimic its output):
Let's do the same using our script, and finally diff
the two
outputs to ensure that our implementation compute the hash correctly:
Victory! Our implementation seems correct and take about 16s (on my machine) to compute all the MD5 hashes. The timings show that our implementation is roughly three order of magnitude slower than openssl. It is actually even worst as we can see that the openssl command was not using the CPU about one third of the time (probably waiting on I/O).
Now we're not trying to race the openssl command of course, but it gives a good idea of how impractical our implementation is for use beyond a coding challenge. Furthermore, our implementation is going to be only slower from now on.
From GNU Awk to AWK
In order to focus on our initial problem of implementing the MD5 algorithm, we've used two GNU Awk features unavailable in AWK:

hexadecimal numbers literal (e.g.,
0xff
) 
bitwise functions:
or
,and
,xor
,compl
,lshift
, andrshift
.
Now that we've successfully wrote the MD5 part, let's "port" our implementation to AWK.
hexadecimal numbers and wrappers
Replacing hexadecimal numbers is going to be trivial, just a pinch of Perl oneliner and some handmade alignment fixes:
We'll also "wrap" the bitwise functions that we intent to reimplement, so that their names don't clash with GNU Awk builtin functions. Now the script looks like this:
step3.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259  BEGIN { _ord_init(); _md5_init(); } { # NOTE: remember the input files inorder in the `files' array. if (nfiles == 0  files[nfiles] != FILENAME) files[++nfiles] = FILENAME; # XXX: only work with files ending with a newline, this is an OK # limitation since it is required by POSIX. content[FILENAME] = content[FILENAME] $0 "\n"; } END { # go over all the files inorder. for (i = 1; i <= nfiles; i++) { fn = files[i]; # ala `openssl md5' output. printf("MD5(%s)= %s\n", fn, md5(content[fn])); } } # our md5 implementation function md5(input, nbytes, chars, i, bytes, hi, lo, words, nwords, state, a, b, c, d, j, x, digest, ret) { # convert the input into an array of bytes using ord() on each # character. nbytes = split(input, chars, ""); for (i = 1; i <= nbytes; i++) bytes[i] = ord(chars[i]); # convert the array of bytes into an array of 32bits words. # NOTE: words is 0indexed. for (i = 1; i <= nbytes; i += 4) { hi = bw_or(bw_lshift(bytes[i + 3], 8), bytes[i + 2]); lo = bw_or(bw_lshift(bytes[i + 1], 8), bytes[i + 0]); words[nwords++] = bw_or(bw_lshift(hi, 16), lo); } # Step 1. Append Padding Bits if (nbytes % 4 == 0) { # the input size is congruent modulo 32, we need a new word to # store the first '1' padding bit. words[nwords++] = 128; # 0x80 } else { # append a '1' bit in the byte just after the last input byte. words[nwords  1] = bw_or(words[nwords  1], bw_lshift(128, (nbytes % 4) * 8)); # 0x80 } # "fill" the remaining bytes with 0 until we're just shy two words of # having 16Word Blocks. while ((nwords % 16) != 14) nwords++; # Step 2. Append Length hi = bw_rshift(nbytes * 8, 32); lo = (nbytes * 8)  bw_lshift(hi, 32); words[nwords++] = lo; words[nwords++] = mod32bits(hi); # truncate to 32 bits # Step 3. Initialize MD Buffer state[0] = 1732584193; # 0x67452301 state[1] = 4023233417; # 0xefcdab89 state[2] = 2562383102; # 0x98badcfe state[3] = 271733878; # 0x10325476 # Step 4. Process Message in 16Word Blocks # Process each 16word block. for (i = 0; i < nwords; i += 16) { # Copy block i into x. for (j = 0; j < 16; j++) x[j] = words[i + j]; a = state[0]; b = state[1]; c = state[2]; d = state[3]; # Round 1 a = FF(a, b, c, d, x[ 0], S11, 3614090360); # 0xd76aa478 d = FF(d, a, b, c, x[ 1], S12, 3905402710); # 0xe8c7b756 c = FF(c, d, a, b, x[ 2], S13, 606105819); # 0x242070db b = FF(b, c, d, a, x[ 3], S14, 3250441966); # 0xc1bdceee a = FF(a, b, c, d, x[ 4], S11, 4118548399); # 0xf57c0faf d = FF(d, a, b, c, x[ 5], S12, 1200080426); # 0x4787c62a c = FF(c, d, a, b, x[ 6], S13, 2821735955); # 0xa8304613 b = FF(b, c, d, a, x[ 7], S14, 4249261313); # 0xfd469501 a = FF(a, b, c, d, x[ 8], S11, 1770035416); # 0x698098d8 d = FF(d, a, b, c, x[ 9], S12, 2336552879); # 0x8b44f7af c = FF(c, d, a, b, x[10], S13, 4294925233); # 0xffff5bb1 b = FF(b, c, d, a, x[11], S14, 2304563134); # 0x895cd7be a = FF(a, b, c, d, x[12], S11, 1804603682); # 0x6b901122 d = FF(d, a, b, c, x[13], S12, 4254626195); # 0xfd987193 c = FF(c, d, a, b, x[14], S13, 2792965006); # 0xa679438e b = FF(b, c, d, a, x[15], S14, 1236535329); # 0x49b40821 # Round 2 a = GG(a, b, c, d, x[ 1], S21, 4129170786); # 0xf61e2562 d = GG(d, a, b, c, x[ 6], S22, 3225465664); # 0xc040b340 c = GG(c, d, a, b, x[11], S23, 643717713); # 0x265e5a51 b = GG(b, c, d, a, x[ 0], S24, 3921069994); # 0xe9b6c7aa a = GG(a, b, c, d, x[ 5], S21, 3593408605); # 0xd62f105d d = GG(d, a, b, c, x[10], S22, 38016083); # 0x2441453 c = GG(c, d, a, b, x[15], S23, 3634488961); # 0xd8a1e681 b = GG(b, c, d, a, x[ 4], S24, 3889429448); # 0xe7d3fbc8 a = GG(a, b, c, d, x[ 9], S21, 568446438); # 0x21e1cde6 d = GG(d, a, b, c, x[14], S22, 3275163606); # 0xc33707d6 c = GG(c, d, a, b, x[ 3], S23, 4107603335); # 0xf4d50d87 b = GG(b, c, d, a, x[ 8], S24, 1163531501); # 0x455a14ed a = GG(a, b, c, d, x[13], S21, 2850285829); # 0xa9e3e905 d = GG(d, a, b, c, x[ 2], S22, 4243563512); # 0xfcefa3f8 c = GG(c, d, a, b, x[ 7], S23, 1735328473); # 0x676f02d9 b = GG(b, c, d, a, x[12], S24, 2368359562); # 0x8d2a4c8a # Round 3 a = HH(a, b, c, d, x[ 5], S31, 4294588738); # 0xfffa3942 d = HH(d, a, b, c, x[ 8], S32, 2272392833); # 0x8771f681 c = HH(c, d, a, b, x[11], S33, 1839030562); # 0x6d9d6122 b = HH(b, c, d, a, x[14], S34, 4259657740); # 0xfde5380c a = HH(a, b, c, d, x[ 1], S31, 2763975236); # 0xa4beea44 d = HH(d, a, b, c, x[ 4], S32, 1272893353); # 0x4bdecfa9 c = HH(c, d, a, b, x[ 7], S33, 4139469664); # 0xf6bb4b60 b = HH(b, c, d, a, x[10], S34, 3200236656); # 0xbebfbc70 a = HH(a, b, c, d, x[13], S31, 681279174); # 0x289b7ec6 d = HH(d, a, b, c, x[ 0], S32, 3936430074); # 0xeaa127fa c = HH(c, d, a, b, x[ 3], S33, 3572445317); # 0xd4ef3085 b = HH(b, c, d, a, x[ 6], S34, 76029189); # 0x4881d05 a = HH(a, b, c, d, x[ 9], S31, 3654602809); # 0xd9d4d039 d = HH(d, a, b, c, x[12], S32, 3873151461); # 0xe6db99e5 c = HH(c, d, a, b, x[15], S33, 530742520); # 0x1fa27cf8 b = HH(b, c, d, a, x[ 2], S34, 3299628645); # 0xc4ac5665 # Round 4 a = II(a, b, c, d, x[ 0], S41, 4096336452); # 0xf4292244 d = II(d, a, b, c, x[ 7], S42, 1126891415); # 0x432aff97 c = II(c, d, a, b, x[14], S43, 2878612391); # 0xab9423a7 b = II(b, c, d, a, x[ 5], S44, 4237533241); # 0xfc93a039 a = II(a, b, c, d, x[12], S41, 1700485571); # 0x655b59c3 d = II(d, a, b, c, x[ 3], S42, 2399980690); # 0x8f0ccc92 c = II(c, d, a, b, x[10], S43, 4293915773); # 0xffeff47d b = II(b, c, d, a, x[ 1], S44, 2240044497); # 0x85845dd1 a = II(a, b, c, d, x[ 8], S41, 1873313359); # 0x6fa87e4f d = II(d, a, b, c, x[15], S42, 4264355552); # 0xfe2ce6e0 c = II(c, d, a, b, x[ 6], S43, 2734768916); # 0xa3014314 b = II(b, c, d, a, x[13], S44, 1309151649); # 0x4e0811a1 a = II(a, b, c, d, x[ 4], S41, 4149444226); # 0xf7537e82 d = II(d, a, b, c, x[11], S42, 3174756917); # 0xbd3af235 c = II(c, d, a, b, x[ 2], S43, 718787259); # 0x2ad7d2bb b = II(b, c, d, a, x[ 9], S44, 3951481745); # 0xeb86d391 state[0] = mod32bits(state[0] + a); state[1] = mod32bits(state[1] + b); state[2] = mod32bits(state[2] + c); state[3] = mod32bits(state[3] + d); } for (i = j = 0; j < 16; j += 4) { digest[j + 0] = bw_and(state[i], 255); # 0xff digest[j + 1] = bw_and(bw_rshift(state[i], 8), 255); # 0xff digest[j + 2] = bw_and(bw_rshift(state[i], 16), 255); # 0xff digest[j + 3] = bw_and(bw_rshift(state[i++], 24), 255); # 0xff } for (i = 0; i < 16; i++) ret = sprintf("%s%02x", ret, digest[i]); return ret; } function F(x, y, z) { return bw_or(bw_and(x, y), bw_and(bw_not(x), z)); } function G(x, y, z) { return bw_or(bw_and(x, z), bw_and(y, bw_not(z))); } function H(x, y, z) { return bw_xor(x, bw_xor(y, z)); } function I(x, y, z) { return bw_xor(y, bw_or(x, bw_not(z))); } function FF(a, b, c, d, x, s, ac) { a = mod32bits(a + F(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function GG(a, b, c, d, x, s, ac) { a = mod32bits(a + G(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function HH(a, b, c, d, x, s, ac) { a = mod32bits(a + H(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function II(a, b, c, d, x, s, ac) { a = mod32bits(a + I(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function ROTATE_LEFT(x, n) { return bw_or(mod32bits(bw_lshift(x, n)), bw_rshift(x, 32  n)); } function mod32bits(x) { return bw_and(x, 4294967295); # 0xffffffff } function bw_not(x) { return mod32bits(compl(x)); } function bw_lshift(x, n) { return lshift(x, n); } function bw_rshift(x, n) { return rshift(x, n); } function bw_and(x, y) { return and(x, y); } function bw_or(x, y) { return or(x, y); } function bw_xor(x, y) { return xor(x, y); } # from https://www.gnu.org/software/gawk/manual/html_node/OrdinalFunctions.html function _ord_init( i) { for (i = 0; i < 256; i++) _ord_[sprintf("%c", i)] = i; } function ord(s) { # only first character is of interest return _ord_[substr(s, 1, 1)]; } function _md5_init() { # MD5 shift constants setup. S11 = 7; S12 = 12; S13 = 17; S14 = 22; S21 = 5; S22 = 9; S23 = 14; S24 = 20; S31 = 4; S32 = 11; S33 = 16; S34 = 23; S41 = 6; S42 = 10; S43 = 15; S44 = 21; } 
diff u step2.awk step3.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269   step2.awk 20191024 09:35:57.779452000 +0200 +++ step3.awk 20170621 09:43:06.599989000 +0200 @@ 33,19 +33,19 @@ # convert the array of bytes into an array of 32bits words. # NOTE: words is 0indexed. for (i = 1; i <= nbytes; i += 4) {  hi = or(lshift(bytes[i + 3], 8), bytes[i + 2]);  lo = or(lshift(bytes[i + 1], 8), bytes[i + 0]);  words[nwords++] = or(lshift(hi, 16), lo); + hi = bw_or(bw_lshift(bytes[i + 3], 8), bytes[i + 2]); + lo = bw_or(bw_lshift(bytes[i + 1], 8), bytes[i + 0]); + words[nwords++] = bw_or(bw_lshift(hi, 16), lo); } # Step 1. Append Padding Bits if (nbytes % 4 == 0) { # the input size is congruent modulo 32, we need a new word to # store the first '1' padding bit.  words[nwords++] = 0x80; + words[nwords++] = 128; # 0x80 } else { # append a '1' bit in the byte just after the last input byte.  words[nwords  1] = or(words[nwords  1], lshift(0x80, (nbytes % 4) * 8)); + words[nwords  1] = bw_or(words[nwords  1], bw_lshift(128, (nbytes % 4) * 8)); # 0x80 } # "fill" the remaining bytes with 0 until we're just shy two words of # having 16Word Blocks. @@ 53,16 +53,16 @@ nwords++; # Step 2. Append Length  hi = rshift(nbytes * 8, 32);  lo = (nbytes * 8)  lshift(hi, 32); + hi = bw_rshift(nbytes * 8, 32); + lo = (nbytes * 8)  bw_lshift(hi, 32); words[nwords++] = lo; words[nwords++] = mod32bits(hi); # truncate to 32 bits # Step 3. Initialize MD Buffer  state[0] = 0x67452301;  state[1] = 0xefcdab89;  state[2] = 0x98badcfe;  state[3] = 0x10325476; + state[0] = 1732584193; # 0x67452301 + state[1] = 4023233417; # 0xefcdab89 + state[2] = 2562383102; # 0x98badcfe + state[3] = 271733878; # 0x10325476 # Step 4. Process Message in 16Word Blocks # Process each 16word block. @@ 73,76 +73,76 @@ a = state[0]; b = state[1]; c = state[2]; d = state[3]; # Round 1  a = FF(a, b, c, d, x[ 0], S11, 0xd76aa478);  d = FF(d, a, b, c, x[ 1], S12, 0xe8c7b756);  c = FF(c, d, a, b, x[ 2], S13, 0x242070db);  b = FF(b, c, d, a, x[ 3], S14, 0xc1bdceee);  a = FF(a, b, c, d, x[ 4], S11, 0xf57c0faf);  d = FF(d, a, b, c, x[ 5], S12, 0x4787c62a);  c = FF(c, d, a, b, x[ 6], S13, 0xa8304613);  b = FF(b, c, d, a, x[ 7], S14, 0xfd469501);  a = FF(a, b, c, d, x[ 8], S11, 0x698098d8);  d = FF(d, a, b, c, x[ 9], S12, 0x8b44f7af);  c = FF(c, d, a, b, x[10], S13, 0xffff5bb1);  b = FF(b, c, d, a, x[11], S14, 0x895cd7be);  a = FF(a, b, c, d, x[12], S11, 0x6b901122);  d = FF(d, a, b, c, x[13], S12, 0xfd987193);  c = FF(c, d, a, b, x[14], S13, 0xa679438e);  b = FF(b, c, d, a, x[15], S14, 0x49b40821); + a = FF(a, b, c, d, x[ 0], S11, 3614090360); # 0xd76aa478 + d = FF(d, a, b, c, x[ 1], S12, 3905402710); # 0xe8c7b756 + c = FF(c, d, a, b, x[ 2], S13, 606105819); # 0x242070db + b = FF(b, c, d, a, x[ 3], S14, 3250441966); # 0xc1bdceee + a = FF(a, b, c, d, x[ 4], S11, 4118548399); # 0xf57c0faf + d = FF(d, a, b, c, x[ 5], S12, 1200080426); # 0x4787c62a + c = FF(c, d, a, b, x[ 6], S13, 2821735955); # 0xa8304613 + b = FF(b, c, d, a, x[ 7], S14, 4249261313); # 0xfd469501 + a = FF(a, b, c, d, x[ 8], S11, 1770035416); # 0x698098d8 + d = FF(d, a, b, c, x[ 9], S12, 2336552879); # 0x8b44f7af + c = FF(c, d, a, b, x[10], S13, 4294925233); # 0xffff5bb1 + b = FF(b, c, d, a, x[11], S14, 2304563134); # 0x895cd7be + a = FF(a, b, c, d, x[12], S11, 1804603682); # 0x6b901122 + d = FF(d, a, b, c, x[13], S12, 4254626195); # 0xfd987193 + c = FF(c, d, a, b, x[14], S13, 2792965006); # 0xa679438e + b = FF(b, c, d, a, x[15], S14, 1236535329); # 0x49b40821 # Round 2  a = GG(a, b, c, d, x[ 1], S21, 0xf61e2562);  d = GG(d, a, b, c, x[ 6], S22, 0xc040b340);  c = GG(c, d, a, b, x[11], S23, 0x265e5a51);  b = GG(b, c, d, a, x[ 0], S24, 0xe9b6c7aa);  a = GG(a, b, c, d, x[ 5], S21, 0xd62f105d);  d = GG(d, a, b, c, x[10], S22, 0x2441453);  c = GG(c, d, a, b, x[15], S23, 0xd8a1e681);  b = GG(b, c, d, a, x[ 4], S24, 0xe7d3fbc8);  a = GG(a, b, c, d, x[ 9], S21, 0x21e1cde6);  d = GG(d, a, b, c, x[14], S22, 0xc33707d6);  c = GG(c, d, a, b, x[ 3], S23, 0xf4d50d87);  b = GG(b, c, d, a, x[ 8], S24, 0x455a14ed);  a = GG(a, b, c, d, x[13], S21, 0xa9e3e905);  d = GG(d, a, b, c, x[ 2], S22, 0xfcefa3f8);  c = GG(c, d, a, b, x[ 7], S23, 0x676f02d9);  b = GG(b, c, d, a, x[12], S24, 0x8d2a4c8a); + a = GG(a, b, c, d, x[ 1], S21, 4129170786); # 0xf61e2562 + d = GG(d, a, b, c, x[ 6], S22, 3225465664); # 0xc040b340 + c = GG(c, d, a, b, x[11], S23, 643717713); # 0x265e5a51 + b = GG(b, c, d, a, x[ 0], S24, 3921069994); # 0xe9b6c7aa + a = GG(a, b, c, d, x[ 5], S21, 3593408605); # 0xd62f105d + d = GG(d, a, b, c, x[10], S22, 38016083); # 0x2441453 + c = GG(c, d, a, b, x[15], S23, 3634488961); # 0xd8a1e681 + b = GG(b, c, d, a, x[ 4], S24, 3889429448); # 0xe7d3fbc8 + a = GG(a, b, c, d, x[ 9], S21, 568446438); # 0x21e1cde6 + d = GG(d, a, b, c, x[14], S22, 3275163606); # 0xc33707d6 + c = GG(c, d, a, b, x[ 3], S23, 4107603335); # 0xf4d50d87 + b = GG(b, c, d, a, x[ 8], S24, 1163531501); # 0x455a14ed + a = GG(a, b, c, d, x[13], S21, 2850285829); # 0xa9e3e905 + d = GG(d, a, b, c, x[ 2], S22, 4243563512); # 0xfcefa3f8 + c = GG(c, d, a, b, x[ 7], S23, 1735328473); # 0x676f02d9 + b = GG(b, c, d, a, x[12], S24, 2368359562); # 0x8d2a4c8a # Round 3  a = HH(a, b, c, d, x[ 5], S31, 0xfffa3942);  d = HH(d, a, b, c, x[ 8], S32, 0x8771f681);  c = HH(c, d, a, b, x[11], S33, 0x6d9d6122);  b = HH(b, c, d, a, x[14], S34, 0xfde5380c);  a = HH(a, b, c, d, x[ 1], S31, 0xa4beea44);  d = HH(d, a, b, c, x[ 4], S32, 0x4bdecfa9);  c = HH(c, d, a, b, x[ 7], S33, 0xf6bb4b60);  b = HH(b, c, d, a, x[10], S34, 0xbebfbc70);  a = HH(a, b, c, d, x[13], S31, 0x289b7ec6);  d = HH(d, a, b, c, x[ 0], S32, 0xeaa127fa);  c = HH(c, d, a, b, x[ 3], S33, 0xd4ef3085);  b = HH(b, c, d, a, x[ 6], S34, 0x4881d05);  a = HH(a, b, c, d, x[ 9], S31, 0xd9d4d039);  d = HH(d, a, b, c, x[12], S32, 0xe6db99e5);  c = HH(c, d, a, b, x[15], S33, 0x1fa27cf8);  b = HH(b, c, d, a, x[ 2], S34, 0xc4ac5665); + a = HH(a, b, c, d, x[ 5], S31, 4294588738); # 0xfffa3942 + d = HH(d, a, b, c, x[ 8], S32, 2272392833); # 0x8771f681 + c = HH(c, d, a, b, x[11], S33, 1839030562); # 0x6d9d6122 + b = HH(b, c, d, a, x[14], S34, 4259657740); # 0xfde5380c + a = HH(a, b, c, d, x[ 1], S31, 2763975236); # 0xa4beea44 + d = HH(d, a, b, c, x[ 4], S32, 1272893353); # 0x4bdecfa9 + c = HH(c, d, a, b, x[ 7], S33, 4139469664); # 0xf6bb4b60 + b = HH(b, c, d, a, x[10], S34, 3200236656); # 0xbebfbc70 + a = HH(a, b, c, d, x[13], S31, 681279174); # 0x289b7ec6 + d = HH(d, a, b, c, x[ 0], S32, 3936430074); # 0xeaa127fa + c = HH(c, d, a, b, x[ 3], S33, 3572445317); # 0xd4ef3085 + b = HH(b, c, d, a, x[ 6], S34, 76029189); # 0x4881d05 + a = HH(a, b, c, d, x[ 9], S31, 3654602809); # 0xd9d4d039 + d = HH(d, a, b, c, x[12], S32, 3873151461); # 0xe6db99e5 + c = HH(c, d, a, b, x[15], S33, 530742520); # 0x1fa27cf8 + b = HH(b, c, d, a, x[ 2], S34, 3299628645); # 0xc4ac5665 # Round 4  a = II(a, b, c, d, x[ 0], S41, 0xf4292244);  d = II(d, a, b, c, x[ 7], S42, 0x432aff97);  c = II(c, d, a, b, x[14], S43, 0xab9423a7);  b = II(b, c, d, a, x[ 5], S44, 0xfc93a039);  a = II(a, b, c, d, x[12], S41, 0x655b59c3);  d = II(d, a, b, c, x[ 3], S42, 0x8f0ccc92);  c = II(c, d, a, b, x[10], S43, 0xffeff47d);  b = II(b, c, d, a, x[ 1], S44, 0x85845dd1);  a = II(a, b, c, d, x[ 8], S41, 0x6fa87e4f);  d = II(d, a, b, c, x[15], S42, 0xfe2ce6e0);  c = II(c, d, a, b, x[ 6], S43, 0xa3014314);  b = II(b, c, d, a, x[13], S44, 0x4e0811a1);  a = II(a, b, c, d, x[ 4], S41, 0xf7537e82);  d = II(d, a, b, c, x[11], S42, 0xbd3af235);  c = II(c, d, a, b, x[ 2], S43, 0x2ad7d2bb);  b = II(b, c, d, a, x[ 9], S44, 0xeb86d391); + a = II(a, b, c, d, x[ 0], S41, 4096336452); # 0xf4292244 + d = II(d, a, b, c, x[ 7], S42, 1126891415); # 0x432aff97 + c = II(c, d, a, b, x[14], S43, 2878612391); # 0xab9423a7 + b = II(b, c, d, a, x[ 5], S44, 4237533241); # 0xfc93a039 + a = II(a, b, c, d, x[12], S41, 1700485571); # 0x655b59c3 + d = II(d, a, b, c, x[ 3], S42, 2399980690); # 0x8f0ccc92 + c = II(c, d, a, b, x[10], S43, 4293915773); # 0xffeff47d + b = II(b, c, d, a, x[ 1], S44, 2240044497); # 0x85845dd1 + a = II(a, b, c, d, x[ 8], S41, 1873313359); # 0x6fa87e4f + d = II(d, a, b, c, x[15], S42, 4264355552); # 0xfe2ce6e0 + c = II(c, d, a, b, x[ 6], S43, 2734768916); # 0xa3014314 + b = II(b, c, d, a, x[13], S44, 1309151649); # 0x4e0811a1 + a = II(a, b, c, d, x[ 4], S41, 4149444226); # 0xf7537e82 + d = II(d, a, b, c, x[11], S42, 3174756917); # 0xbd3af235 + c = II(c, d, a, b, x[ 2], S43, 718787259); # 0x2ad7d2bb + b = II(b, c, d, a, x[ 9], S44, 3951481745); # 0xeb86d391 state[0] = mod32bits(state[0] + a); state[1] = mod32bits(state[1] + b); @@ 151,10 +151,10 @@ } for (i = j = 0; j < 16; j += 4) {  digest[j + 0] = and(state[i], 0xff);  digest[j + 1] = and(rshift(state[i], 8), 0xff);  digest[j + 2] = and(rshift(state[i], 16), 0xff);  digest[j + 3] = and(rshift(state[i++], 24), 0xff); + digest[j + 0] = bw_and(state[i], 255); # 0xff + digest[j + 1] = bw_and(bw_rshift(state[i], 8), 255); # 0xff + digest[j + 2] = bw_and(bw_rshift(state[i], 16), 255); # 0xff + digest[j + 3] = bw_and(bw_rshift(state[i++], 24), 255); # 0xff } for (i = 0; i < 16; i++) ret = sprintf("%s%02x", ret, digest[i]); @@ 162,19 +162,19 @@ } function F(x, y, z) {  return or(and(x, y), and(not(x), z)); + return bw_or(bw_and(x, y), bw_and(bw_not(x), z)); } function G(x, y, z) {  return or(and(x, z), and(y, not(z))); + return bw_or(bw_and(x, z), bw_and(y, bw_not(z))); } function H(x, y, z) {  return xor(x, xor(y, z)); + return bw_xor(x, bw_xor(y, z)); } function I(x, y, z) {  return xor(y, or(x, not(z))); + return bw_xor(y, bw_or(x, bw_not(z))); } function FF(a, b, c, d, x, s, ac) { @@ 206,15 +206,35 @@ } function ROTATE_LEFT(x, n) {  return or(mod32bits(lshift(x, n)), rshift(x, 32  n)); + return bw_or(mod32bits(bw_lshift(x, n)), bw_rshift(x, 32  n)); } function mod32bits(x) {  return and(x, 0xffffffff); + return bw_and(x, 4294967295); # 0xffffffff } function not(x) { +function bw_not(x) { return mod32bits(compl(x)); +} + +function bw_lshift(x, n) { + return lshift(x, n); +} + +function bw_rshift(x, n) { + return rshift(x, n); +} + +function bw_and(x, y) { + return and(x, y); +} + +function bw_or(x, y) { + return or(x, y); +} + +function bw_xor(x, y) { + return xor(x, y); } # from https://www.gnu.org/software/gawk/manual/html_node/OrdinalFunctions.html 
We'll assume that the hexadecimal change did not impact performances, but "wrapping" the bitwise functions is expected to incur an overhead as they are intensively used by the MD5 rounds:
If we compare to the step2.awk timing (about 16s), this version is ~56% slower (!) Hang on though, as I'm sure we can do worse ;)
lshift & rshift
We need to implement the bitwise functions based on what AWK has to offer.
Because we feel like we have a hint, let's start with lshift
and
rshift
.
A classic "bit trick" in C (and others languages) to perform multiplication
and division on unsigned integer by a power of two is to use logical
left shift and logical right shift (respectively). Here we can do the
opposite, i.e. using multiplication and division to implement our
bw_lshift
and bw_rshift
functions:
step4.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259  BEGIN { _ord_init(); _md5_init(); } { # NOTE: remember the input files inorder in the `files' array. if (nfiles == 0  files[nfiles] != FILENAME) files[++nfiles] = FILENAME; # XXX: only work with files ending with a newline, this is an OK # limitation since it is required by POSIX. content[FILENAME] = content[FILENAME] $0 "\n"; } END { # go over all the files inorder. for (i = 1; i <= nfiles; i++) { fn = files[i]; # ala `openssl md5' output. printf("MD5(%s)= %s\n", fn, md5(content[fn])); } } # our md5 implementation function md5(input, nbytes, chars, i, bytes, hi, lo, words, nwords, state, a, b, c, d, j, x, digest, ret) { # convert the input into an array of bytes using ord() on each # character. nbytes = split(input, chars, ""); for (i = 1; i <= nbytes; i++) bytes[i] = ord(chars[i]); # convert the array of bytes into an array of 32bits words. # NOTE: words is 0indexed. for (i = 1; i <= nbytes; i += 4) { hi = bw_or(bw_lshift(bytes[i + 3], 8), bytes[i + 2]); lo = bw_or(bw_lshift(bytes[i + 1], 8), bytes[i + 0]); words[nwords++] = bw_or(bw_lshift(hi, 16), lo); } # Step 1. Append Padding Bits if (nbytes % 4 == 0) { # the input size is congruent modulo 32, we need a new word to # store the first '1' padding bit. words[nwords++] = 128; # 0x80 } else { # append a '1' bit in the byte just after the last input byte. words[nwords  1] = bw_or(words[nwords  1], bw_lshift(128, (nbytes % 4) * 8)); # 0x80 } # "fill" the remaining bytes with 0 until we're just shy two words of # having 16Word Blocks. while ((nwords % 16) != 14) nwords++; # Step 2. Append Length hi = bw_rshift(nbytes * 8, 32); lo = (nbytes * 8)  bw_lshift(hi, 32); words[nwords++] = lo; words[nwords++] = mod32bits(hi); # truncate to 32 bits # Step 3. Initialize MD Buffer state[0] = 1732584193; # 0x67452301 state[1] = 4023233417; # 0xefcdab89 state[2] = 2562383102; # 0x98badcfe state[3] = 271733878; # 0x10325476 # Step 4. Process Message in 16Word Blocks # Process each 16word block. for (i = 0; i < nwords; i += 16) { # Copy block i into x. for (j = 0; j < 16; j++) x[j] = words[i + j]; a = state[0]; b = state[1]; c = state[2]; d = state[3]; # Round 1 a = FF(a, b, c, d, x[ 0], S11, 3614090360); # 0xd76aa478 d = FF(d, a, b, c, x[ 1], S12, 3905402710); # 0xe8c7b756 c = FF(c, d, a, b, x[ 2], S13, 606105819); # 0x242070db b = FF(b, c, d, a, x[ 3], S14, 3250441966); # 0xc1bdceee a = FF(a, b, c, d, x[ 4], S11, 4118548399); # 0xf57c0faf d = FF(d, a, b, c, x[ 5], S12, 1200080426); # 0x4787c62a c = FF(c, d, a, b, x[ 6], S13, 2821735955); # 0xa8304613 b = FF(b, c, d, a, x[ 7], S14, 4249261313); # 0xfd469501 a = FF(a, b, c, d, x[ 8], S11, 1770035416); # 0x698098d8 d = FF(d, a, b, c, x[ 9], S12, 2336552879); # 0x8b44f7af c = FF(c, d, a, b, x[10], S13, 4294925233); # 0xffff5bb1 b = FF(b, c, d, a, x[11], S14, 2304563134); # 0x895cd7be a = FF(a, b, c, d, x[12], S11, 1804603682); # 0x6b901122 d = FF(d, a, b, c, x[13], S12, 4254626195); # 0xfd987193 c = FF(c, d, a, b, x[14], S13, 2792965006); # 0xa679438e b = FF(b, c, d, a, x[15], S14, 1236535329); # 0x49b40821 # Round 2 a = GG(a, b, c, d, x[ 1], S21, 4129170786); # 0xf61e2562 d = GG(d, a, b, c, x[ 6], S22, 3225465664); # 0xc040b340 c = GG(c, d, a, b, x[11], S23, 643717713); # 0x265e5a51 b = GG(b, c, d, a, x[ 0], S24, 3921069994); # 0xe9b6c7aa a = GG(a, b, c, d, x[ 5], S21, 3593408605); # 0xd62f105d d = GG(d, a, b, c, x[10], S22, 38016083); # 0x2441453 c = GG(c, d, a, b, x[15], S23, 3634488961); # 0xd8a1e681 b = GG(b, c, d, a, x[ 4], S24, 3889429448); # 0xe7d3fbc8 a = GG(a, b, c, d, x[ 9], S21, 568446438); # 0x21e1cde6 d = GG(d, a, b, c, x[14], S22, 3275163606); # 0xc33707d6 c = GG(c, d, a, b, x[ 3], S23, 4107603335); # 0xf4d50d87 b = GG(b, c, d, a, x[ 8], S24, 1163531501); # 0x455a14ed a = GG(a, b, c, d, x[13], S21, 2850285829); # 0xa9e3e905 d = GG(d, a, b, c, x[ 2], S22, 4243563512); # 0xfcefa3f8 c = GG(c, d, a, b, x[ 7], S23, 1735328473); # 0x676f02d9 b = GG(b, c, d, a, x[12], S24, 2368359562); # 0x8d2a4c8a # Round 3 a = HH(a, b, c, d, x[ 5], S31, 4294588738); # 0xfffa3942 d = HH(d, a, b, c, x[ 8], S32, 2272392833); # 0x8771f681 c = HH(c, d, a, b, x[11], S33, 1839030562); # 0x6d9d6122 b = HH(b, c, d, a, x[14], S34, 4259657740); # 0xfde5380c a = HH(a, b, c, d, x[ 1], S31, 2763975236); # 0xa4beea44 d = HH(d, a, b, c, x[ 4], S32, 1272893353); # 0x4bdecfa9 c = HH(c, d, a, b, x[ 7], S33, 4139469664); # 0xf6bb4b60 b = HH(b, c, d, a, x[10], S34, 3200236656); # 0xbebfbc70 a = HH(a, b, c, d, x[13], S31, 681279174); # 0x289b7ec6 d = HH(d, a, b, c, x[ 0], S32, 3936430074); # 0xeaa127fa c = HH(c, d, a, b, x[ 3], S33, 3572445317); # 0xd4ef3085 b = HH(b, c, d, a, x[ 6], S34, 76029189); # 0x4881d05 a = HH(a, b, c, d, x[ 9], S31, 3654602809); # 0xd9d4d039 d = HH(d, a, b, c, x[12], S32, 3873151461); # 0xe6db99e5 c = HH(c, d, a, b, x[15], S33, 530742520); # 0x1fa27cf8 b = HH(b, c, d, a, x[ 2], S34, 3299628645); # 0xc4ac5665 # Round 4 a = II(a, b, c, d, x[ 0], S41, 4096336452); # 0xf4292244 d = II(d, a, b, c, x[ 7], S42, 1126891415); # 0x432aff97 c = II(c, d, a, b, x[14], S43, 2878612391); # 0xab9423a7 b = II(b, c, d, a, x[ 5], S44, 4237533241); # 0xfc93a039 a = II(a, b, c, d, x[12], S41, 1700485571); # 0x655b59c3 d = II(d, a, b, c, x[ 3], S42, 2399980690); # 0x8f0ccc92 c = II(c, d, a, b, x[10], S43, 4293915773); # 0xffeff47d b = II(b, c, d, a, x[ 1], S44, 2240044497); # 0x85845dd1 a = II(a, b, c, d, x[ 8], S41, 1873313359); # 0x6fa87e4f d = II(d, a, b, c, x[15], S42, 4264355552); # 0xfe2ce6e0 c = II(c, d, a, b, x[ 6], S43, 2734768916); # 0xa3014314 b = II(b, c, d, a, x[13], S44, 1309151649); # 0x4e0811a1 a = II(a, b, c, d, x[ 4], S41, 4149444226); # 0xf7537e82 d = II(d, a, b, c, x[11], S42, 3174756917); # 0xbd3af235 c = II(c, d, a, b, x[ 2], S43, 718787259); # 0x2ad7d2bb b = II(b, c, d, a, x[ 9], S44, 3951481745); # 0xeb86d391 state[0] = mod32bits(state[0] + a); state[1] = mod32bits(state[1] + b); state[2] = mod32bits(state[2] + c); state[3] = mod32bits(state[3] + d); } for (i = j = 0; j < 16; j += 4) { digest[j + 0] = bw_and(state[i], 255); # 0xff digest[j + 1] = bw_and(bw_rshift(state[i], 8), 255); # 0xff digest[j + 2] = bw_and(bw_rshift(state[i], 16), 255); # 0xff digest[j + 3] = bw_and(bw_rshift(state[i++], 24), 255); # 0xff } for (i = 0; i < 16; i++) ret = sprintf("%s%02x", ret, digest[i]); return ret; } function F(x, y, z) { return bw_or(bw_and(x, y), bw_and(bw_not(x), z)); } function G(x, y, z) { return bw_or(bw_and(x, z), bw_and(y, bw_not(z))); } function H(x, y, z) { return bw_xor(x, bw_xor(y, z)); } function I(x, y, z) { return bw_xor(y, bw_or(x, bw_not(z))); } function FF(a, b, c, d, x, s, ac) { a = mod32bits(a + F(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function GG(a, b, c, d, x, s, ac) { a = mod32bits(a + G(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function HH(a, b, c, d, x, s, ac) { a = mod32bits(a + H(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function II(a, b, c, d, x, s, ac) { a = mod32bits(a + I(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function ROTATE_LEFT(x, n) { return bw_or(mod32bits(bw_lshift(x, n)), bw_rshift(x, 32  n)); } function mod32bits(x) { return bw_and(x, 4294967295); # 0xffffffff } function bw_not(x) { return mod32bits(compl(x)); } function bw_lshift(x, n) { return x * (2 ^ n); } function bw_rshift(x, n) { return int(x / (2 ^ n)); } function bw_and(x, y) { return and(x, y); } function bw_or(x, y) { return or(x, y); } function bw_xor(x, y) { return xor(x, y); } # from https://www.gnu.org/software/gawk/manual/html_node/OrdinalFunctions.html function _ord_init( i) { for (i = 0; i < 256; i++) _ord_[sprintf("%c", i)] = i; } function ord(s) { # only first character is of interest return _ord_[substr(s, 1, 1)]; } function _md5_init() { # MD5 shift constants setup. S11 = 7; S12 = 12; S13 = 17; S14 = 22; S21 = 5; S22 = 9; S23 = 14; S24 = 20; S31 = 4; S32 = 11; S33 = 16; S34 = 23; S41 = 6; S42 = 10; S43 = 15; S44 = 21; } 
diff u step3.awk step4.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16   step3.awk 20170220 21:30:13.502377000 +0100 +++ step4.awk 20170220 21:30:23.464046000 +0100 @@ 218,11 +218,11 @@ } function bw_lshift(x, n) {  return lshift(x, n); + return x * (2 ^ n); } function bw_rshift(x, n) {  return rshift(x, n); + return int(x / (2 ^ n)); } function bw_and(x, y) { 
Surprisingly (at least to me), AWK has an exponentiation operator
^
making computing power of two trivial.
Note that since GNU Awk lshift
may yield a number that is greater than 2^32 we already ensure to
mod32bits
its result where needed (in ROTATE_LEFT
).
Thus, our bw_lshift
doesn't have to hold this invariant. Because
AWK division is a "real division" (i.e. not an integer division) we need to
drop the fractional part from bw_rshift
's result, hence the
int()
truncation.
The running time increased by "only" about 3%, which is less that I
personally expected. My best guess is that lshift
and
rshift
are actually slower than what "we think" (we usually see
them as fast), because of the number representation dance that
GNU Awk has to do (from double
to
uintmax_t
and then back to double
).
and, or, xor & not
Intuitively I thought about iterating through each bit using %
and shifting. Although this should "do the trick" it is awfully slow:
example of bw_and() bit by bit
1 2 3 4 5 6 7 8  function bw_and(x, y, i, r) { for (i = 0; i < 32; i++) { r += ((x % 2) && (y % 2)) * (2 ^ i); x = int(x / 2); y = int(y / 2); } return r; } 
Line 3 can easily be adapted to implement bw_or
,
bw_xor
and bw_not
.
Seeking a better solution, I found this helpful StackOverflow answer. The AND solution described looks promising: by using a lookup table for each combination of 4 x 4 bits of input we only have to loop 8 times (instead of 32). Additionally, XOR and OR are elegantly implemented on top of AND.
Instead of using a precomputed lookup table, I've opted for bootstrapping it
using the slow bitbybit method. This will be more consistent with
_ord_init
, yield cleaner and smaller code, and just feels more
in tune with the "pure AWK" vision that is the raison d'être of this
implementation.
Here is our version, translated in AWK (using multiple array subscripts) and tweaked a bit (because the StackOverflow answer was for 16bits numbers and we use 32bits).
step5.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280  BEGIN { _ord_init(); _bitwise_init(); _md5_init(); } { # NOTE: remember the input files inorder in the `files' array. if (nfiles == 0  files[nfiles] != FILENAME) files[++nfiles] = FILENAME; # XXX: only work with files ending with a newline, this is an OK # limitation since it is required by POSIX. content[FILENAME] = content[FILENAME] $0 "\n"; } END { # go over all the files inorder. for (i = 1; i <= nfiles; i++) { fn = files[i]; # ala `openssl md5' output. printf("MD5(%s)= %s\n", fn, md5(content[fn])); } } # our md5 implementation function md5(input, nbytes, chars, i, bytes, hi, lo, words, nwords, state, a, b, c, d, j, x, digest, ret) { # convert the input into an array of bytes using ord() on each # character. nbytes = split(input, chars, ""); for (i = 1; i <= nbytes; i++) bytes[i] = ord(chars[i]); # convert the array of bytes into an array of 32bits words. # NOTE: words is 0indexed. for (i = 1; i <= nbytes; i += 4) { hi = bw_or(bw_lshift(bytes[i + 3], 8), bytes[i + 2]); lo = bw_or(bw_lshift(bytes[i + 1], 8), bytes[i + 0]); words[nwords++] = bw_or(bw_lshift(hi, 16), lo); } # Step 1. Append Padding Bits if (nbytes % 4 == 0) { # the input size is congruent modulo 32, we need a new word to # store the first '1' padding bit. words[nwords++] = 128; # 0x80 } else { # append a '1' bit in the byte just after the last input byte. words[nwords  1] = bw_or(words[nwords  1], bw_lshift(128, (nbytes % 4) * 8)); # 0x80 } # "fill" the remaining bytes with 0 until we're just shy two words of # having 16Word Blocks. while ((nwords % 16) != 14) nwords++; # Step 2. Append Length hi = bw_rshift(nbytes * 8, 32); lo = (nbytes * 8)  bw_lshift(hi, 32); words[nwords++] = lo; words[nwords++] = mod32bits(hi); # truncate to 32 bits # Step 3. Initialize MD Buffer state[0] = 1732584193; # 0x67452301 state[1] = 4023233417; # 0xefcdab89 state[2] = 2562383102; # 0x98badcfe state[3] = 271733878; # 0x10325476 # Step 4. Process Message in 16Word Blocks # Process each 16word block. for (i = 0; i < nwords; i += 16) { # Copy block i into x. for (j = 0; j < 16; j++) x[j] = words[i + j]; a = state[0]; b = state[1]; c = state[2]; d = state[3]; # Round 1 a = FF(a, b, c, d, x[ 0], S11, 3614090360); # 0xd76aa478 d = FF(d, a, b, c, x[ 1], S12, 3905402710); # 0xe8c7b756 c = FF(c, d, a, b, x[ 2], S13, 606105819); # 0x242070db b = FF(b, c, d, a, x[ 3], S14, 3250441966); # 0xc1bdceee a = FF(a, b, c, d, x[ 4], S11, 4118548399); # 0xf57c0faf d = FF(d, a, b, c, x[ 5], S12, 1200080426); # 0x4787c62a c = FF(c, d, a, b, x[ 6], S13, 2821735955); # 0xa8304613 b = FF(b, c, d, a, x[ 7], S14, 4249261313); # 0xfd469501 a = FF(a, b, c, d, x[ 8], S11, 1770035416); # 0x698098d8 d = FF(d, a, b, c, x[ 9], S12, 2336552879); # 0x8b44f7af c = FF(c, d, a, b, x[10], S13, 4294925233); # 0xffff5bb1 b = FF(b, c, d, a, x[11], S14, 2304563134); # 0x895cd7be a = FF(a, b, c, d, x[12], S11, 1804603682); # 0x6b901122 d = FF(d, a, b, c, x[13], S12, 4254626195); # 0xfd987193 c = FF(c, d, a, b, x[14], S13, 2792965006); # 0xa679438e b = FF(b, c, d, a, x[15], S14, 1236535329); # 0x49b40821 # Round 2 a = GG(a, b, c, d, x[ 1], S21, 4129170786); # 0xf61e2562 d = GG(d, a, b, c, x[ 6], S22, 3225465664); # 0xc040b340 c = GG(c, d, a, b, x[11], S23, 643717713); # 0x265e5a51 b = GG(b, c, d, a, x[ 0], S24, 3921069994); # 0xe9b6c7aa a = GG(a, b, c, d, x[ 5], S21, 3593408605); # 0xd62f105d d = GG(d, a, b, c, x[10], S22, 38016083); # 0x2441453 c = GG(c, d, a, b, x[15], S23, 3634488961); # 0xd8a1e681 b = GG(b, c, d, a, x[ 4], S24, 3889429448); # 0xe7d3fbc8 a = GG(a, b, c, d, x[ 9], S21, 568446438); # 0x21e1cde6 d = GG(d, a, b, c, x[14], S22, 3275163606); # 0xc33707d6 c = GG(c, d, a, b, x[ 3], S23, 4107603335); # 0xf4d50d87 b = GG(b, c, d, a, x[ 8], S24, 1163531501); # 0x455a14ed a = GG(a, b, c, d, x[13], S21, 2850285829); # 0xa9e3e905 d = GG(d, a, b, c, x[ 2], S22, 4243563512); # 0xfcefa3f8 c = GG(c, d, a, b, x[ 7], S23, 1735328473); # 0x676f02d9 b = GG(b, c, d, a, x[12], S24, 2368359562); # 0x8d2a4c8a # Round 3 a = HH(a, b, c, d, x[ 5], S31, 4294588738); # 0xfffa3942 d = HH(d, a, b, c, x[ 8], S32, 2272392833); # 0x8771f681 c = HH(c, d, a, b, x[11], S33, 1839030562); # 0x6d9d6122 b = HH(b, c, d, a, x[14], S34, 4259657740); # 0xfde5380c a = HH(a, b, c, d, x[ 1], S31, 2763975236); # 0xa4beea44 d = HH(d, a, b, c, x[ 4], S32, 1272893353); # 0x4bdecfa9 c = HH(c, d, a, b, x[ 7], S33, 4139469664); # 0xf6bb4b60 b = HH(b, c, d, a, x[10], S34, 3200236656); # 0xbebfbc70 a = HH(a, b, c, d, x[13], S31, 681279174); # 0x289b7ec6 d = HH(d, a, b, c, x[ 0], S32, 3936430074); # 0xeaa127fa c = HH(c, d, a, b, x[ 3], S33, 3572445317); # 0xd4ef3085 b = HH(b, c, d, a, x[ 6], S34, 76029189); # 0x4881d05 a = HH(a, b, c, d, x[ 9], S31, 3654602809); # 0xd9d4d039 d = HH(d, a, b, c, x[12], S32, 3873151461); # 0xe6db99e5 c = HH(c, d, a, b, x[15], S33, 530742520); # 0x1fa27cf8 b = HH(b, c, d, a, x[ 2], S34, 3299628645); # 0xc4ac5665 # Round 4 a = II(a, b, c, d, x[ 0], S41, 4096336452); # 0xf4292244 d = II(d, a, b, c, x[ 7], S42, 1126891415); # 0x432aff97 c = II(c, d, a, b, x[14], S43, 2878612391); # 0xab9423a7 b = II(b, c, d, a, x[ 5], S44, 4237533241); # 0xfc93a039 a = II(a, b, c, d, x[12], S41, 1700485571); # 0x655b59c3 d = II(d, a, b, c, x[ 3], S42, 2399980690); # 0x8f0ccc92 c = II(c, d, a, b, x[10], S43, 4293915773); # 0xffeff47d b = II(b, c, d, a, x[ 1], S44, 2240044497); # 0x85845dd1 a = II(a, b, c, d, x[ 8], S41, 1873313359); # 0x6fa87e4f d = II(d, a, b, c, x[15], S42, 4264355552); # 0xfe2ce6e0 c = II(c, d, a, b, x[ 6], S43, 2734768916); # 0xa3014314 b = II(b, c, d, a, x[13], S44, 1309151649); # 0x4e0811a1 a = II(a, b, c, d, x[ 4], S41, 4149444226); # 0xf7537e82 d = II(d, a, b, c, x[11], S42, 3174756917); # 0xbd3af235 c = II(c, d, a, b, x[ 2], S43, 718787259); # 0x2ad7d2bb b = II(b, c, d, a, x[ 9], S44, 3951481745); # 0xeb86d391 state[0] = mod32bits(state[0] + a); state[1] = mod32bits(state[1] + b); state[2] = mod32bits(state[2] + c); state[3] = mod32bits(state[3] + d); } for (i = j = 0; j < 16; j += 4) { digest[j + 0] = bw_and(state[i], 255); # 0xff digest[j + 1] = bw_and(bw_rshift(state[i], 8), 255); # 0xff digest[j + 2] = bw_and(bw_rshift(state[i], 16), 255); # 0xff digest[j + 3] = bw_and(bw_rshift(state[i++], 24), 255); # 0xff } for (i = 0; i < 16; i++) ret = sprintf("%s%02x", ret, digest[i]); return ret; } function F(x, y, z) { return bw_or(bw_and(x, y), bw_and(bw_not(x), z)); } function G(x, y, z) { return bw_or(bw_and(x, z), bw_and(y, bw_not(z))); } function H(x, y, z) { return bw_xor(x, bw_xor(y, z)); } function I(x, y, z) { return bw_xor(y, bw_or(x, bw_not(z))); } function FF(a, b, c, d, x, s, ac) { a = mod32bits(a + F(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function GG(a, b, c, d, x, s, ac) { a = mod32bits(a + G(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function HH(a, b, c, d, x, s, ac) { a = mod32bits(a + H(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function II(a, b, c, d, x, s, ac) { a = mod32bits(a + I(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function ROTATE_LEFT(x, n) { return bw_or(mod32bits(bw_lshift(x, n)), bw_rshift(x, 32  n)); } function mod32bits(x) { return bw_and(x, 4294967295); # 0xffffffff } function bw_not(x) { return bw_xor(x, 4294967295); # 0xffffffff } function bw_lshift(x, n) { return x * (2 ^ n); } function bw_rshift(x, n) { return int(x / (2 ^ n)); } function bw_and(x, y, i, r) { for (i = 0; i < 32; i += 4) { r = r / (2 ^ 4) + andlookup[x % 16, y % 16] * (2 ^ 28); x = int(x / (2 ^ 4)); y = int(y / (2 ^ 4)); } return r; } function bw_or(x, y) { return bw_xor(bw_xor(x, y), bw_and(x, y)); } function bw_xor(x, y) { return (x + y  2 * bw_and(x, y)); } # from https://www.gnu.org/software/gawk/manual/html_node/OrdinalFunctions.html function _ord_init( i) { for (i = 0; i < 256; i++) _ord_[sprintf("%c", i)] = i; } function ord(s) { # only first character is of interest return _ord_[substr(s, 1, 1)]; } function _bitwise_init( a, b, x, y, i) { # generate the lookup table for bw_and(). for (a = 0; a < 16; a++) { for (b = 0; b < 16; b++) { x = a; y = b; for (i = 0; i < 4; i++) { andlookup[a, b] += ((x % 2) && (y % 2)) * (2 ^ i); x = int(x / 2); y = int(y / 2); } } } } function _md5_init() { # MD5 shift constants setup. S11 = 7; S12 = 12; S13 = 17; S14 = 22; S21 = 5; S22 = 9; S23 = 14; S24 = 20; S31 = 4; S32 = 11; S33 = 16; S34 = 23; S41 = 6; S42 = 10; S43 = 15; S44 = 21; } 
diff u step4.awk step5.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66   step4.awk 20170220 21:30:23.464046000 +0100 +++ step5.awk 20170221 14:58:31.499103000 +0100 @@ 1,5 +1,6 @@ BEGIN { _ord_init(); + _bitwise_init(); _md5_init(); } @@ 214,7 +215,7 @@ } function bw_not(x) {  return mod32bits(compl(x)); + return bw_xor(x, 4294967295); # 0xffffffff } function bw_lshift(x, n) { @@ 225,16 +226,21 @@ return int(x / (2 ^ n)); } function bw_and(x, y) {  return and(x, y); +function bw_and(x, y, i, r) { + for (i = 0; i < 32; i += 4) { + r = r / (2 ^ 4) + andlookup[x % 16, y % 16] * (2 ^ 28); + x = int(x / (2 ^ 4)); + y = int(y / (2 ^ 4)); + } + return r; } function bw_or(x, y) {  return or(x, y); + return bw_xor(bw_xor(x, y), bw_and(x, y)); } function bw_xor(x, y) {  return xor(x, y); + return (x + y  2 * bw_and(x, y)); } # from https://www.gnu.org/software/gawk/manual/html_node/OrdinalFunctions.html @@ 250,6 +256,21 @@ return _ord_[substr(s, 1, 1)]; } +function _bitwise_init( a, b, x, y, i) { + # generate the lookup table for bw_and(). + for (a = 0; a < 16; a++) { + for (b = 0; b < 16; b++) { + x = a; + y = b; + for (i = 0; i < 4; i++) { + andlookup[a, b] += ((x % 2) && (y % 2)) * (2 ^ i); + x = int(x / 2); + y = int(y / 2); + } + } + } +} + function _md5_init() { # MD5 shift constants setup. S11 = 7; S12 = 12; S13 = 17; S14 = 22; 
Let's test our changes:
Now the good news is that our script doesn't use GNU Awk bitwise function anymore, and consequently we can also test with onetrueawk too:
Yay! We've succeeded in our attempt to implement MD5 in AWK! That's a big step as we can test with several AWK implementations from now on. At this point though, the performances are horrible even for my coding challenge context: the script is roughly 20 times slower than the previous version. Also worth noting, GNU Awk is about two times faster than onetrueawk for this test.
Improvements
While we didn't worry too much about performances in our porting steps, we've kept an eye on the impact of each change. Thus, we should have a good idea of what should be improved, but let's run our script through GNU Awk's profiler anyway (yes, it has one):
For something more visual, here is a (handmade by yours truly) callgraph:
step5.dot
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16  digraph G { ROTATE_LEFT > bw_lshift; ROTATE_LEFT > bw_rshift; ROTATE_LEFT > mod32bits; ROTATE_LEFT > bw_or; mod32bits > bw_and; bw_not > bw_xor; bw_or > bw_xor; bw_or > bw_xor; bw_or > bw_and; bw_xor > bw_and; bw_and [color="#dd4444"]; bw_lshift [color="#268bd2"]; bw_rshift [color="#268bd2"]; } 
We see that there is a ton of pressure on bw_and
. Sadly, the
profiler is only able to report how many time each function was called. The
running time of each function would have been very useful, but we'll have to
"play it by ear". From the timings we gathered so far we think that:
 From step3.awk we learned that calling functions is costly.

step4.awk showed that both
bw_lshift
andbw_rshift
performances are close to before porting, so they're not on our radar to optimize (hence the blue color in the callgraph). 
Since step5.awk
bw_and
has become really slow (hence the red color in the callgraph), but we think the memory vs runtime tradeoff is acceptable.
bw_or
Obviously our bw_or
function can be improved: it calls
bw_xor
twice and bw_and
once, for a total of "3
times bw_and
cost". The most straightforward refactoring is to
mimic bw_and
implementation and should make bw_or
three time faster:
step6.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286  BEGIN { _ord_init(); _bitwise_init(); _md5_init(); } { # NOTE: remember the input files inorder in the `files' array. if (nfiles == 0  files[nfiles] != FILENAME) files[++nfiles] = FILENAME; # XXX: only work with files ending with a newline, this is an OK # limitation since it is required by POSIX. content[FILENAME] = content[FILENAME] $0 "\n"; } END { # go over all the files inorder. for (i = 1; i <= nfiles; i++) { fn = files[i]; # ala `openssl md5' output. printf("MD5(%s)= %s\n", fn, md5(content[fn])); } } # our md5 implementation function md5(input, nbytes, chars, i, bytes, hi, lo, words, nwords, state, a, b, c, d, j, x, digest, ret) { # convert the input into an array of bytes using ord() on each # character. nbytes = split(input, chars, ""); for (i = 1; i <= nbytes; i++) bytes[i] = ord(chars[i]); # convert the array of bytes into an array of 32bits words. # NOTE: words is 0indexed. for (i = 1; i <= nbytes; i += 4) { hi = bw_or(bw_lshift(bytes[i + 3], 8), bytes[i + 2]); lo = bw_or(bw_lshift(bytes[i + 1], 8), bytes[i + 0]); words[nwords++] = bw_or(bw_lshift(hi, 16), lo); } # Step 1. Append Padding Bits if (nbytes % 4 == 0) { # the input size is congruent modulo 32, we need a new word to # store the first '1' padding bit. words[nwords++] = 128; # 0x80 } else { # append a '1' bit in the byte just after the last input byte. words[nwords  1] = bw_or(words[nwords  1], bw_lshift(128, (nbytes % 4) * 8)); # 0x80 } # "fill" the remaining bytes with 0 until we're just shy two words of # having 16Word Blocks. while ((nwords % 16) != 14) nwords++; # Step 2. Append Length hi = bw_rshift(nbytes * 8, 32); lo = (nbytes * 8)  bw_lshift(hi, 32); words[nwords++] = lo; words[nwords++] = mod32bits(hi); # truncate to 32 bits # Step 3. Initialize MD Buffer state[0] = 1732584193; # 0x67452301 state[1] = 4023233417; # 0xefcdab89 state[2] = 2562383102; # 0x98badcfe state[3] = 271733878; # 0x10325476 # Step 4. Process Message in 16Word Blocks # Process each 16word block. for (i = 0; i < nwords; i += 16) { # Copy block i into x. for (j = 0; j < 16; j++) x[j] = words[i + j]; a = state[0]; b = state[1]; c = state[2]; d = state[3]; # Round 1 a = FF(a, b, c, d, x[ 0], S11, 3614090360); # 0xd76aa478 d = FF(d, a, b, c, x[ 1], S12, 3905402710); # 0xe8c7b756 c = FF(c, d, a, b, x[ 2], S13, 606105819); # 0x242070db b = FF(b, c, d, a, x[ 3], S14, 3250441966); # 0xc1bdceee a = FF(a, b, c, d, x[ 4], S11, 4118548399); # 0xf57c0faf d = FF(d, a, b, c, x[ 5], S12, 1200080426); # 0x4787c62a c = FF(c, d, a, b, x[ 6], S13, 2821735955); # 0xa8304613 b = FF(b, c, d, a, x[ 7], S14, 4249261313); # 0xfd469501 a = FF(a, b, c, d, x[ 8], S11, 1770035416); # 0x698098d8 d = FF(d, a, b, c, x[ 9], S12, 2336552879); # 0x8b44f7af c = FF(c, d, a, b, x[10], S13, 4294925233); # 0xffff5bb1 b = FF(b, c, d, a, x[11], S14, 2304563134); # 0x895cd7be a = FF(a, b, c, d, x[12], S11, 1804603682); # 0x6b901122 d = FF(d, a, b, c, x[13], S12, 4254626195); # 0xfd987193 c = FF(c, d, a, b, x[14], S13, 2792965006); # 0xa679438e b = FF(b, c, d, a, x[15], S14, 1236535329); # 0x49b40821 # Round 2 a = GG(a, b, c, d, x[ 1], S21, 4129170786); # 0xf61e2562 d = GG(d, a, b, c, x[ 6], S22, 3225465664); # 0xc040b340 c = GG(c, d, a, b, x[11], S23, 643717713); # 0x265e5a51 b = GG(b, c, d, a, x[ 0], S24, 3921069994); # 0xe9b6c7aa a = GG(a, b, c, d, x[ 5], S21, 3593408605); # 0xd62f105d d = GG(d, a, b, c, x[10], S22, 38016083); # 0x2441453 c = GG(c, d, a, b, x[15], S23, 3634488961); # 0xd8a1e681 b = GG(b, c, d, a, x[ 4], S24, 3889429448); # 0xe7d3fbc8 a = GG(a, b, c, d, x[ 9], S21, 568446438); # 0x21e1cde6 d = GG(d, a, b, c, x[14], S22, 3275163606); # 0xc33707d6 c = GG(c, d, a, b, x[ 3], S23, 4107603335); # 0xf4d50d87 b = GG(b, c, d, a, x[ 8], S24, 1163531501); # 0x455a14ed a = GG(a, b, c, d, x[13], S21, 2850285829); # 0xa9e3e905 d = GG(d, a, b, c, x[ 2], S22, 4243563512); # 0xfcefa3f8 c = GG(c, d, a, b, x[ 7], S23, 1735328473); # 0x676f02d9 b = GG(b, c, d, a, x[12], S24, 2368359562); # 0x8d2a4c8a # Round 3 a = HH(a, b, c, d, x[ 5], S31, 4294588738); # 0xfffa3942 d = HH(d, a, b, c, x[ 8], S32, 2272392833); # 0x8771f681 c = HH(c, d, a, b, x[11], S33, 1839030562); # 0x6d9d6122 b = HH(b, c, d, a, x[14], S34, 4259657740); # 0xfde5380c a = HH(a, b, c, d, x[ 1], S31, 2763975236); # 0xa4beea44 d = HH(d, a, b, c, x[ 4], S32, 1272893353); # 0x4bdecfa9 c = HH(c, d, a, b, x[ 7], S33, 4139469664); # 0xf6bb4b60 b = HH(b, c, d, a, x[10], S34, 3200236656); # 0xbebfbc70 a = HH(a, b, c, d, x[13], S31, 681279174); # 0x289b7ec6 d = HH(d, a, b, c, x[ 0], S32, 3936430074); # 0xeaa127fa c = HH(c, d, a, b, x[ 3], S33, 3572445317); # 0xd4ef3085 b = HH(b, c, d, a, x[ 6], S34, 76029189); # 0x4881d05 a = HH(a, b, c, d, x[ 9], S31, 3654602809); # 0xd9d4d039 d = HH(d, a, b, c, x[12], S32, 3873151461); # 0xe6db99e5 c = HH(c, d, a, b, x[15], S33, 530742520); # 0x1fa27cf8 b = HH(b, c, d, a, x[ 2], S34, 3299628645); # 0xc4ac5665 # Round 4 a = II(a, b, c, d, x[ 0], S41, 4096336452); # 0xf4292244 d = II(d, a, b, c, x[ 7], S42, 1126891415); # 0x432aff97 c = II(c, d, a, b, x[14], S43, 2878612391); # 0xab9423a7 b = II(b, c, d, a, x[ 5], S44, 4237533241); # 0xfc93a039 a = II(a, b, c, d, x[12], S41, 1700485571); # 0x655b59c3 d = II(d, a, b, c, x[ 3], S42, 2399980690); # 0x8f0ccc92 c = II(c, d, a, b, x[10], S43, 4293915773); # 0xffeff47d b = II(b, c, d, a, x[ 1], S44, 2240044497); # 0x85845dd1 a = II(a, b, c, d, x[ 8], S41, 1873313359); # 0x6fa87e4f d = II(d, a, b, c, x[15], S42, 4264355552); # 0xfe2ce6e0 c = II(c, d, a, b, x[ 6], S43, 2734768916); # 0xa3014314 b = II(b, c, d, a, x[13], S44, 1309151649); # 0x4e0811a1 a = II(a, b, c, d, x[ 4], S41, 4149444226); # 0xf7537e82 d = II(d, a, b, c, x[11], S42, 3174756917); # 0xbd3af235 c = II(c, d, a, b, x[ 2], S43, 718787259); # 0x2ad7d2bb b = II(b, c, d, a, x[ 9], S44, 3951481745); # 0xeb86d391 state[0] = mod32bits(state[0] + a); state[1] = mod32bits(state[1] + b); state[2] = mod32bits(state[2] + c); state[3] = mod32bits(state[3] + d); } for (i = j = 0; j < 16; j += 4) { digest[j + 0] = bw_and(state[i], 255); # 0xff digest[j + 1] = bw_and(bw_rshift(state[i], 8), 255); # 0xff digest[j + 2] = bw_and(bw_rshift(state[i], 16), 255); # 0xff digest[j + 3] = bw_and(bw_rshift(state[i++], 24), 255); # 0xff } for (i = 0; i < 16; i++) ret = sprintf("%s%02x", ret, digest[i]); return ret; } function F(x, y, z) { return bw_or(bw_and(x, y), bw_and(bw_not(x), z)); } function G(x, y, z) { return bw_or(bw_and(x, z), bw_and(y, bw_not(z))); } function H(x, y, z) { return bw_xor(x, bw_xor(y, z)); } function I(x, y, z) { return bw_xor(y, bw_or(x, bw_not(z))); } function FF(a, b, c, d, x, s, ac) { a = mod32bits(a + F(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function GG(a, b, c, d, x, s, ac) { a = mod32bits(a + G(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function HH(a, b, c, d, x, s, ac) { a = mod32bits(a + H(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function II(a, b, c, d, x, s, ac) { a = mod32bits(a + I(b, c, d) + x + ac); a = ROTATE_LEFT(a, s); a = mod32bits(a + b); return a; } function ROTATE_LEFT(x, n) { return bw_or(mod32bits(bw_lshift(x, n)), bw_rshift(x, 32  n)); } function mod32bits(x) { return bw_and(x, 4294967295); # 0xffffffff } function bw_not(x) { return bw_xor(x, 4294967295); # 0xffffffff } function bw_lshift(x, n) { return x * (2 ^ n); } function bw_rshift(x, n) { return int(x / (2 ^ n)); } function bw_and(x, y, i, r) { for (i = 0; i < 32; i += 4) { r = r / (2 ^ 4) + bw_lookup["and", x % 16, y % 16] * (2 ^ 28); x = int(x / (2 ^ 4)); y = int(y / (2 ^ 4)); } return r; } function bw_or(x, y, i, r) { for (i = 0; i < 32; i += 4) { r = r / (2 ^ 4) + bw_lookup["or", x % 16, y % 16] * (2 ^ 28); x = int(x / (2 ^ 4)); y = int(y / (2 ^ 4)); } return r; } function bw_xor(x, y) { return (x + y  2 * bw_and(x, y)); } # from https://www.gnu.org/software/gawk/manual/html_node/OrdinalFunctions.html function _ord_init( i) { for (i = 0; i < 256; i++) _ord_[sprintf("%c", i)] = i; } function ord(s) { # only first character is of interest return _ord_[substr(s, 1, 1)]; } function _bitwise_init( a, b, x, y, i) { # generate the bw_lookup table used by bw_and() and bw_or(). for (a = 0; a < 16; a++) { for (b = 0; b < 16; b++) { x = a; y = b; for (i = 0; i < 4; i++) { bw_lookup["and", a, b] += ((x % 2) && (y % 2)) * (2 ^ i); bw_lookup["or", a, b] += ((x % 2)  (y % 2)) * (2 ^ i); x = int(x / 2); y = int(y / 2); } } } } function _md5_init() { # MD5 shift constants setup. S11 = 7; S12 = 12; S13 = 17; S14 = 22; S21 = 5; S22 = 9; S23 = 14; S24 = 20; S31 = 4; S32 = 11; S33 = 16; S34 = 23; S41 = 6; S42 = 10; S43 = 15; S44 = 21; } 
diff u step5.awk step6.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43   step5.awk 20170221 14:58:31.499103000 +0100 +++ step6.awk 20170304 12:57:02.936165000 +0100 @@ 228,15 +228,20 @@ function bw_and(x, y, i, r) { for (i = 0; i < 32; i += 4) {  r = r / (2 ^ 4) + andlookup[x % 16, y % 16] * (2 ^ 28); + r = r / (2 ^ 4) + bw_lookup["and", x % 16, y % 16] * (2 ^ 28); x = int(x / (2 ^ 4)); y = int(y / (2 ^ 4)); } return r; } function bw_or(x, y) {  return bw_xor(bw_xor(x, y), bw_and(x, y)); +function bw_or(x, y, i, r) { + for (i = 0; i < 32; i += 4) { + r = r / (2 ^ 4) + bw_lookup["or", x % 16, y % 16] * (2 ^ 28); + x = int(x / (2 ^ 4)); + y = int(y / (2 ^ 4)); + } + return r; } function bw_xor(x, y) { @@ 257,13 +262,14 @@ } function _bitwise_init( a, b, x, y, i) {  # generate the lookup table for bw_and(). + # generate the bw_lookup table used by bw_and() and bw_or(). for (a = 0; a < 16; a++) { for (b = 0; b < 16; b++) { x = a; y = b; for (i = 0; i < 4; i++) {  andlookup[a, b] += ((x % 2) && (y % 2)) * (2 ^ i); + bw_lookup["and", a, b] += ((x % 2) && (y % 2)) * (2 ^ i); + bw_lookup["or", a, b] += ((x % 2)  (y % 2)) * (2 ^ i); x = int(x / 2); y = int(y / 2); } 
Although we could easily do the same for bw_not
and
bw_xor
, their "perceived bw_and
cost" is 1, and so
there is not much to be gain. Let's see the impact of our change:
That's a 35% improvement, well worth it. Let's focus now on the next "obvious" target.
mod32bits
Originally, we wrote mod32bits
using and
because it
was "fast". Now the modulo %
operator is very likely to be
faster than bw_and
and we can easily just get rid of the
mod32bits
function:
step7.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282  BEGIN { _ord_init(); _bitwise_init(); _md5_init(); } { # NOTE: remember the input files inorder in the `files' array. if (nfiles == 0  files[nfiles] != FILENAME) files[++nfiles] = FILENAME; # XXX: only work with files ending with a newline, this is an OK # limitation since it is required by POSIX. content[FILENAME] = content[FILENAME] $0 "\n"; } END { # go over all the files inorder. for (i = 1; i <= nfiles; i++) { fn = files[i]; # ala `openssl md5' output. printf("MD5(%s)= %s\n", fn, md5(content[fn])); } } # our md5 implementation function md5(input, nbytes, chars, i, bytes, hi, lo, words, nwords, state, a, b, c, d, j, x, digest, ret) { # convert the input into an array of bytes using ord() on each # character. nbytes = split(input, chars, ""); for (i = 1; i <= nbytes; i++) bytes[i] = ord(chars[i]); # convert the array of bytes into an array of 32bits words. # NOTE: words is 0indexed. for (i = 1; i <= nbytes; i += 4) { hi = bw_or(bw_lshift(bytes[i + 3], 8), bytes[i + 2]); lo = bw_or(bw_lshift(bytes[i + 1], 8), bytes[i + 0]); words[nwords++] = bw_or(bw_lshift(hi, 16), lo); } # Step 1. Append Padding Bits if (nbytes % 4 == 0) { # the input size is congruent modulo 32, we need a new word to # store the first '1' padding bit. words[nwords++] = 128; # 0x80 } else { # append a '1' bit in the byte just after the last input byte. words[nwords  1] = bw_or(words[nwords  1], bw_lshift(128, (nbytes % 4) * 8)); # 0x80 } # "fill" the remaining bytes with 0 until we're just shy two words of # having 16Word Blocks. while ((nwords % 16) != 14) nwords++; # Step 2. Append Length hi = bw_rshift(nbytes * 8, 32); lo = (nbytes * 8)  bw_lshift(hi, 32); words[nwords++] = lo; words[nwords++] = hi % (2 ^ 32); # truncate to 32 bits # Step 3. Initialize MD Buffer state[0] = 1732584193; # 0x67452301 state[1] = 4023233417; # 0xefcdab89 state[2] = 2562383102; # 0x98badcfe state[3] = 271733878; # 0x10325476 # Step 4. Process Message in 16Word Blocks # Process each 16word block. for (i = 0; i < nwords; i += 16) { # Copy block i into x. for (j = 0; j < 16; j++) x[j] = words[i + j]; a = state[0]; b = state[1]; c = state[2]; d = state[3]; # Round 1 a = FF(a, b, c, d, x[ 0], S11, 3614090360); # 0xd76aa478 d = FF(d, a, b, c, x[ 1], S12, 3905402710); # 0xe8c7b756 c = FF(c, d, a, b, x[ 2], S13, 606105819); # 0x242070db b = FF(b, c, d, a, x[ 3], S14, 3250441966); # 0xc1bdceee a = FF(a, b, c, d, x[ 4], S11, 4118548399); # 0xf57c0faf d = FF(d, a, b, c, x[ 5], S12, 1200080426); # 0x4787c62a c = FF(c, d, a, b, x[ 6], S13, 2821735955); # 0xa8304613 b = FF(b, c, d, a, x[ 7], S14, 4249261313); # 0xfd469501 a = FF(a, b, c, d, x[ 8], S11, 1770035416); # 0x698098d8 d = FF(d, a, b, c, x[ 9], S12, 2336552879); # 0x8b44f7af c = FF(c, d, a, b, x[10], S13, 4294925233); # 0xffff5bb1 b = FF(b, c, d, a, x[11], S14, 2304563134); # 0x895cd7be a = FF(a, b, c, d, x[12], S11, 1804603682); # 0x6b901122 d = FF(d, a, b, c, x[13], S12, 4254626195); # 0xfd987193 c = FF(c, d, a, b, x[14], S13, 2792965006); # 0xa679438e b = FF(b, c, d, a, x[15], S14, 1236535329); # 0x49b40821 # Round 2 a = GG(a, b, c, d, x[ 1], S21, 4129170786); # 0xf61e2562 d = GG(d, a, b, c, x[ 6], S22, 3225465664); # 0xc040b340 c = GG(c, d, a, b, x[11], S23, 643717713); # 0x265e5a51 b = GG(b, c, d, a, x[ 0], S24, 3921069994); # 0xe9b6c7aa a = GG(a, b, c, d, x[ 5], S21, 3593408605); # 0xd62f105d d = GG(d, a, b, c, x[10], S22, 38016083); # 0x2441453 c = GG(c, d, a, b, x[15], S23, 3634488961); # 0xd8a1e681 b = GG(b, c, d, a, x[ 4], S24, 3889429448); # 0xe7d3fbc8 a = GG(a, b, c, d, x[ 9], S21, 568446438); # 0x21e1cde6 d = GG(d, a, b, c, x[14], S22, 3275163606); # 0xc33707d6 c = GG(c, d, a, b, x[ 3], S23, 4107603335); # 0xf4d50d87 b = GG(b, c, d, a, x[ 8], S24, 1163531501); # 0x455a14ed a = GG(a, b, c, d, x[13], S21, 2850285829); # 0xa9e3e905 d = GG(d, a, b, c, x[ 2], S22, 4243563512); # 0xfcefa3f8 c = GG(c, d, a, b, x[ 7], S23, 1735328473); # 0x676f02d9 b = GG(b, c, d, a, x[12], S24, 2368359562); # 0x8d2a4c8a # Round 3 a = HH(a, b, c, d, x[ 5], S31, 4294588738); # 0xfffa3942 d = HH(d, a, b, c, x[ 8], S32, 2272392833); # 0x8771f681 c = HH(c, d, a, b, x[11], S33, 1839030562); # 0x6d9d6122 b = HH(b, c, d, a, x[14], S34, 4259657740); # 0xfde5380c a = HH(a, b, c, d, x[ 1], S31, 2763975236); # 0xa4beea44 d = HH(d, a, b, c, x[ 4], S32, 1272893353); # 0x4bdecfa9 c = HH(c, d, a, b, x[ 7], S33, 4139469664); # 0xf6bb4b60 b = HH(b, c, d, a, x[10], S34, 3200236656); # 0xbebfbc70 a = HH(a, b, c, d, x[13], S31, 681279174); # 0x289b7ec6 d = HH(d, a, b, c, x[ 0], S32, 3936430074); # 0xeaa127fa c = HH(c, d, a, b, x[ 3], S33, 3572445317); # 0xd4ef3085 b = HH(b, c, d, a, x[ 6], S34, 76029189); # 0x4881d05 a = HH(a, b, c, d, x[ 9], S31, 3654602809); # 0xd9d4d039 d = HH(d, a, b, c, x[12], S32, 3873151461); # 0xe6db99e5 c = HH(c, d, a, b, x[15], S33, 530742520); # 0x1fa27cf8 b = HH(b, c, d, a, x[ 2], S34, 3299628645); # 0xc4ac5665 # Round 4 a = II(a, b, c, d, x[ 0], S41, 4096336452); # 0xf4292244 d = II(d, a, b, c, x[ 7], S42, 1126891415); # 0x432aff97 c = II(c, d, a, b, x[14], S43, 2878612391); # 0xab9423a7 b = II(b, c, d, a, x[ 5], S44, 4237533241); # 0xfc93a039 a = II(a, b, c, d, x[12], S41, 1700485571); # 0x655b59c3 d = II(d, a, b, c, x[ 3], S42, 2399980690); # 0x8f0ccc92 c = II(c, d, a, b, x[10], S43, 4293915773); # 0xffeff47d b = II(b, c, d, a, x[ 1], S44, 2240044497); # 0x85845dd1 a = II(a, b, c, d, x[ 8], S41, 1873313359); # 0x6fa87e4f d = II(d, a, b, c, x[15], S42, 4264355552); # 0xfe2ce6e0 c = II(c, d, a, b, x[ 6], S43, 2734768916); # 0xa3014314 b = II(b, c, d, a, x[13], S44, 1309151649); # 0x4e0811a1 a = II(a, b, c, d, x[ 4], S41, 4149444226); # 0xf7537e82 d = II(d, a, b, c, x[11], S42, 3174756917); # 0xbd3af235 c = II(c, d, a, b, x[ 2], S43, 718787259); # 0x2ad7d2bb b = II(b, c, d, a, x[ 9], S44, 3951481745); # 0xeb86d391 state[0] = (state[0] + a) % (2 ^ 32); state[1] = (state[1] + b) % (2 ^ 32); state[2] = (state[2] + c) % (2 ^ 32); state[3] = (state[3] + d) % (2 ^ 32); } for (i = j = 0; j < 16; j += 4) { digest[j + 0] = bw_and(state[i], 255); # 0xff digest[j + 1] = bw_and(bw_rshift(state[i], 8), 255); # 0xff digest[j + 2] = bw_and(bw_rshift(state[i], 16), 255); # 0xff digest[j + 3] = bw_and(bw_rshift(state[i++], 24), 255); # 0xff } for (i = 0; i < 16; i++) ret = sprintf("%s%02x", ret, digest[i]); return ret; } function F(x, y, z) { return bw_or(bw_and(x, y), bw_and(bw_not(x), z)); } function G(x, y, z) { return bw_or(bw_and(x, z), bw_and(y, bw_not(z))); } function H(x, y, z) { return bw_xor(x, bw_xor(y, z)); } function I(x, y, z) { return bw_xor(y, bw_or(x, bw_not(z))); } function FF(a, b, c, d, x, s, ac) { a = (a + F(b, c, d) + x + ac) % (2 ^ 32); a = ROTATE_LEFT(a, s); a = (a + b) % (2 ^ 32); return a; } function GG(a, b, c, d, x, s, ac) { a = (a + G(b, c, d) + x + ac) % (2 ^ 32); a = ROTATE_LEFT(a, s); a = (a + b) % (2 ^ 32); return a; } function HH(a, b, c, d, x, s, ac) { a = (a + H(b, c, d) + x + ac) % (2 ^ 32); a = ROTATE_LEFT(a, s); a = (a + b) % (2 ^ 32); return a; } function II(a, b, c, d, x, s, ac) { a = (a + I(b, c, d) + x + ac) % (2 ^ 32); a = ROTATE_LEFT(a, s); a = (a + b) % (2 ^ 32); return a; } function ROTATE_LEFT(x, n) { return bw_or((bw_lshift(x, n)) % (2 ^ 32), bw_rshift(x, 32  n)); } function bw_not(x) { return bw_xor(x, 4294967295); # 0xffffffff } function bw_lshift(x, n) { return x * (2 ^ n); } function bw_rshift(x, n) { return int(x / (2 ^ n)); } function bw_and(x, y, i, r) { for (i = 0; i < 32; i += 4) { r = r / (2 ^ 4) + bw_lookup["and", x % 16, y % 16] * (2 ^ 28); x = int(x / (2 ^ 4)); y = int(y / (2 ^ 4)); } return r; } function bw_or(x, y, i, r) { for (i = 0; i < 32; i += 4) { r = r / (2 ^ 4) + bw_lookup["or", x % 16, y % 16] * (2 ^ 28); x = int(x / (2 ^ 4)); y = int(y / (2 ^ 4)); } return r; } function bw_xor(x, y) { return (x + y  2 * bw_and(x, y)); } # from https://www.gnu.org/software/gawk/manual/html_node/OrdinalFunctions.html function _ord_init( i) { for (i = 0; i < 256; i++) _ord_[sprintf("%c", i)] = i; } function ord(s) { # only first character is of interest return _ord_[substr(s, 1, 1)]; } function _bitwise_init( a, b, x, y, i) { # generate the bw_lookup table used by bw_and() and bw_or(). for (a = 0; a < 16; a++) { for (b = 0; b < 16; b++) { x = a; y = b; for (i = 0; i < 4; i++) { bw_lookup["and", a, b] += ((x % 2) && (y % 2)) * (2 ^ i); bw_lookup["or", a, b] += ((x % 2)  (y % 2)) * (2 ^ i); x = int(x / 2); y = int(y / 2); } } } } function _md5_init() { # MD5 shift constants setup. S11 = 7; S12 = 12; S13 = 17; S14 = 22; S21 = 5; S22 = 9; S23 = 14; S24 = 20; S31 = 4; S32 = 11; S33 = 16; S34 = 23; S41 = 6; S42 = 10; S43 = 15; S44 = 21; } 
diff u step6.awk step7.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75   step6.awk 20170304 12:57:02.936165000 +0100 +++ step7.awk 20170307 12:08:08.666000000 +0100 @@ 57,7 +57,7 @@ hi = bw_rshift(nbytes * 8, 32); lo = (nbytes * 8)  bw_lshift(hi, 32); words[nwords++] = lo;  words[nwords++] = mod32bits(hi); # truncate to 32 bits + words[nwords++] = hi % (2 ^ 32); # truncate to 32 bits # Step 3. Initialize MD Buffer state[0] = 1732584193; # 0x67452301 @@ 145,10 +145,10 @@ c = II(c, d, a, b, x[ 2], S43, 718787259); # 0x2ad7d2bb b = II(b, c, d, a, x[ 9], S44, 3951481745); # 0xeb86d391  state[0] = mod32bits(state[0] + a);  state[1] = mod32bits(state[1] + b);  state[2] = mod32bits(state[2] + c);  state[3] = mod32bits(state[3] + d); + state[0] = (state[0] + a) % (2 ^ 32); + state[1] = (state[1] + b) % (2 ^ 32); + state[2] = (state[2] + c) % (2 ^ 32); + state[3] = (state[3] + d) % (2 ^ 32); } for (i = j = 0; j < 16; j += 4) { @@ 179,39 +179,35 @@ } function FF(a, b, c, d, x, s, ac) {  a = mod32bits(a + F(b, c, d) + x + ac); + a = (a + F(b, c, d) + x + ac) % (2 ^ 32); a = ROTATE_LEFT(a, s);  a = mod32bits(a + b); + a = (a + b) % (2 ^ 32); return a; } function GG(a, b, c, d, x, s, ac) {  a = mod32bits(a + G(b, c, d) + x + ac); + a = (a + G(b, c, d) + x + ac) % (2 ^ 32); a = ROTATE_LEFT(a, s);  a = mod32bits(a + b); + a = (a + b) % (2 ^ 32); return a; } function HH(a, b, c, d, x, s, ac) {  a = mod32bits(a + H(b, c, d) + x + ac); + a = (a + H(b, c, d) + x + ac) % (2 ^ 32); a = ROTATE_LEFT(a, s);  a = mod32bits(a + b); + a = (a + b) % (2 ^ 32); return a; } function II(a, b, c, d, x, s, ac) {  a = mod32bits(a + I(b, c, d) + x + ac); + a = (a + I(b, c, d) + x + ac) % (2 ^ 32); a = ROTATE_LEFT(a, s);  a = mod32bits(a + b); + a = (a + b) % (2 ^ 32); return a; } function ROTATE_LEFT(x, n) {  return bw_or(mod32bits(bw_lshift(x, n)), bw_rshift(x, 32  n)); }  function mod32bits(x) {  return bw_and(x, 4294967295); # 0xffffffff + return bw_or((bw_lshift(x, n)) % (2 ^ 32), bw_rshift(x, 32  n)); } function bw_not(x) { 
Again, about 35% faster and arguably a bit more readable too.
ROTATE_LEFT
By taking a look at ROTATE_LEFT
, we see that we can use addition
instead of OR'ing between the high and low part of the result, because there
is no "overlap" between the operands "interesting bits" (i.e. no chance of
carry). In other words, addition and OR in this particular case will yield
the same result (but the addition will be faster than our slow
bw_or
implementation).
We can apply this logic in a few other places and also to bw_and
by using modulo instead when we only care about the last few bits:
step8.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284  BEGIN { _ord_init(); _bitwise_init(); _md5_init(); } { # NOTE: remember the input files inorder in the `files' array. if (nfiles == 0  files[nfiles] != FILENAME) files[++nfiles] = FILENAME; # XXX: only work with files ending with a newline, this is an OK # limitation since it is required by POSIX. content[FILENAME] = content[FILENAME] $0 "\n"; } END { # go over all the files inorder. for (i = 1; i <= nfiles; i++) { fn = files[i]; # ala `openssl md5' output. printf("MD5(%s)= %s\n", fn, md5(content[fn])); } } # our md5 implementation function md5(input, nbytes, chars, i, bytes, hi, lo, words, nwords, state, a, b, c, d, j, x, digest, ret) { # convert the input into an array of bytes using ord() on each # character. nbytes = split(input, chars, ""); for (i = 1; i <= nbytes; i++) bytes[i] = ord(chars[i]); # convert the array of bytes into an array of 32bits words. # NOTE: words is 0indexed. for (i = 1; i <= nbytes; i += 4) { hi = bw_lshift(bytes[i + 3], 8) + bytes[i + 2]; lo = bw_lshift(bytes[i + 1], 8) + bytes[i + 0]; words[nwords++] = bw_lshift(hi, 16) + lo; } # Step 1. Append Padding Bits if (nbytes % 4 == 0) { # the input size is congruent modulo 32, we need a new word to # store the first '1' padding bit. words[nwords++] = 128; # 0x80 } else { # append a '1' bit in the byte just after the last input byte. words[nwords  1] = words[nwords  1] + bw_lshift(128, (nbytes % 4) * 8); # 0x80 } # "fill" the remaining bytes with 0 until we're just shy two words of # having 16Word Blocks. while ((nwords % 16) != 14) nwords++; # Step 2. Append Length hi = bw_rshift(nbytes * 8, 32); lo = (nbytes * 8)  bw_lshift(hi, 32); words[nwords++] = lo; words[nwords++] = hi % (2 ^ 32); # truncate to 32 bits # Step 3. Initialize MD Buffer state[0] = 1732584193; # 0x67452301 state[1] = 4023233417; # 0xefcdab89 state[2] = 2562383102; # 0x98badcfe state[3] = 271733878; # 0x10325476 # Step 4. Process Message in 16Word Blocks # Process each 16word block. for (i = 0; i < nwords; i += 16) { # Copy block i into x. for (j = 0; j < 16; j++) x[j] = words[i + j]; a = state[0]; b = state[1]; c = state[2]; d = state[3]; # Round 1 a = FF(a, b, c, d, x[ 0], S11, 3614090360); # 0xd76aa478 d = FF(d, a, b, c, x[ 1], S12, 3905402710); # 0xe8c7b756 c = FF(c, d, a, b, x[ 2], S13, 606105819); # 0x242070db b = FF(b, c, d, a, x[ 3], S14, 3250441966); # 0xc1bdceee a = FF(a, b, c, d, x[ 4], S11, 4118548399); # 0xf57c0faf d = FF(d, a, b, c, x[ 5], S12, 1200080426); # 0x4787c62a c = FF(c, d, a, b, x[ 6], S13, 2821735955); # 0xa8304613 b = FF(b, c, d, a, x[ 7], S14, 4249261313); # 0xfd469501 a = FF(a, b, c, d, x[ 8], S11, 1770035416); # 0x698098d8 d = FF(d, a, b, c, x[ 9], S12, 2336552879); # 0x8b44f7af c = FF(c, d, a, b, x[10], S13, 4294925233); # 0xffff5bb1 b = FF(b, c, d, a, x[11], S14, 2304563134); # 0x895cd7be a = FF(a, b, c, d, x[12], S11, 1804603682); # 0x6b901122 d = FF(d, a, b, c, x[13], S12, 4254626195); # 0xfd987193 c = FF(c, d, a, b, x[14], S13, 2792965006); # 0xa679438e b = FF(b, c, d, a, x[15], S14, 1236535329); # 0x49b40821 # Round 2 a = GG(a, b, c, d, x[ 1], S21, 4129170786); # 0xf61e2562 d = GG(d, a, b, c, x[ 6], S22, 3225465664); # 0xc040b340 c = GG(c, d, a, b, x[11], S23, 643717713); # 0x265e5a51 b = GG(b, c, d, a, x[ 0], S24, 3921069994); # 0xe9b6c7aa a = GG(a, b, c, d, x[ 5], S21, 3593408605); # 0xd62f105d d = GG(d, a, b, c, x[10], S22, 38016083); # 0x2441453 c = GG(c, d, a, b, x[15], S23, 3634488961); # 0xd8a1e681 b = GG(b, c, d, a, x[ 4], S24, 3889429448); # 0xe7d3fbc8 a = GG(a, b, c, d, x[ 9], S21, 568446438); # 0x21e1cde6 d = GG(d, a, b, c, x[14], S22, 3275163606); # 0xc33707d6 c = GG(c, d, a, b, x[ 3], S23, 4107603335); # 0xf4d50d87 b = GG(b, c, d, a, x[ 8], S24, 1163531501); # 0x455a14ed a = GG(a, b, c, d, x[13], S21, 2850285829); # 0xa9e3e905 d = GG(d, a, b, c, x[ 2], S22, 4243563512); # 0xfcefa3f8 c = GG(c, d, a, b, x[ 7], S23, 1735328473); # 0x676f02d9 b = GG(b, c, d, a, x[12], S24, 2368359562); # 0x8d2a4c8a # Round 3 a = HH(a, b, c, d, x[ 5], S31, 4294588738); # 0xfffa3942 d = HH(d, a, b, c, x[ 8], S32, 2272392833); # 0x8771f681 c = HH(c, d, a, b, x[11], S33, 1839030562); # 0x6d9d6122 b = HH(b, c, d, a, x[14], S34, 4259657740); # 0xfde5380c a = HH(a, b, c, d, x[ 1], S31, 2763975236); # 0xa4beea44 d = HH(d, a, b, c, x[ 4], S32, 1272893353); # 0x4bdecfa9 c = HH(c, d, a, b, x[ 7], S33, 4139469664); # 0xf6bb4b60 b = HH(b, c, d, a, x[10], S34, 3200236656); # 0xbebfbc70 a = HH(a, b, c, d, x[13], S31, 681279174); # 0x289b7ec6 d = HH(d, a, b, c, x[ 0], S32, 3936430074); # 0xeaa127fa c = HH(c, d, a, b, x[ 3], S33, 3572445317); # 0xd4ef3085 b = HH(b, c, d, a, x[ 6], S34, 76029189); # 0x4881d05 a = HH(a, b, c, d, x[ 9], S31, 3654602809); # 0xd9d4d039 d = HH(d, a, b, c, x[12], S32, 3873151461); # 0xe6db99e5 c = HH(c, d, a, b, x[15], S33, 530742520); # 0x1fa27cf8 b = HH(b, c, d, a, x[ 2], S34, 3299628645); # 0xc4ac5665 # Round 4 a = II(a, b, c, d, x[ 0], S41, 4096336452); # 0xf4292244 d = II(d, a, b, c, x[ 7], S42, 1126891415); # 0x432aff97 c = II(c, d, a, b, x[14], S43, 2878612391); # 0xab9423a7 b = II(b, c, d, a, x[ 5], S44, 4237533241); # 0xfc93a039 a = II(a, b, c, d, x[12], S41, 1700485571); # 0x655b59c3 d = II(d, a, b, c, x[ 3], S42, 2399980690); # 0x8f0ccc92 c = II(c, d, a, b, x[10], S43, 4293915773); # 0xffeff47d b = II(b, c, d, a, x[ 1], S44, 2240044497); # 0x85845dd1 a = II(a, b, c, d, x[ 8], S41, 1873313359); # 0x6fa87e4f d = II(d, a, b, c, x[15], S42, 4264355552); # 0xfe2ce6e0 c = II(c, d, a, b, x[ 6], S43, 2734768916); # 0xa3014314 b = II(b, c, d, a, x[13], S44, 1309151649); # 0x4e0811a1 a = II(a, b, c, d, x[ 4], S41, 4149444226); # 0xf7537e82 d = II(d, a, b, c, x[11], S42, 3174756917); # 0xbd3af235 c = II(c, d, a, b, x[ 2], S43, 718787259); # 0x2ad7d2bb b = II(b, c, d, a, x[ 9], S44, 3951481745); # 0xeb86d391 state[0] = (state[0] + a) % (2 ^ 32); state[1] = (state[1] + b) % (2 ^ 32); state[2] = (state[2] + c) % (2 ^ 32); state[3] = (state[3] + d) % (2 ^ 32); } for (i = j = 0; j < 16; j += 4) { digest[j + 0] = state[i] % (2 ^ 8); digest[j + 1] = bw_rshift(state[i], 8) % (2 ^ 8); digest[j + 2] = bw_rshift(state[i], 16) % (2 ^ 8); digest[j + 3] = bw_rshift(state[i++], 24) % (2 ^ 8); } for (i = 0; i < 16; i++) ret = sprintf("%s%02x", ret, digest[i]); return ret; } function F(x, y, z) { return bw_or(bw_and(x, y), bw_and(bw_not(x), z)); } function G(x, y, z) { return bw_or(bw_and(x, z), bw_and(y, bw_not(z))); } function H(x, y, z) { return bw_xor(x, bw_xor(y, z)); } function I(x, y, z) { return bw_xor(y, bw_or(x, bw_not(z))); } function FF(a, b, c, d, x, s, ac) { a = (a + F(b, c, d) + x + ac) % (2 ^ 32); a = ROTATE_LEFT(a, s); a = (a + b) % (2 ^ 32); return a; } function GG(a, b, c, d, x, s, ac) { a = (a + G(b, c, d) + x + ac) % (2 ^ 32); a = ROTATE_LEFT(a, s); a = (a + b) % (2 ^ 32); return a; } function HH(a, b, c, d, x, s, ac) { a = (a + H(b, c, d) + x + ac) % (2 ^ 32); a = ROTATE_LEFT(a, s); a = (a + b) % (2 ^ 32); return a; } function II(a, b, c, d, x, s, ac) { a = (a + I(b, c, d) + x + ac) % (2 ^ 32); a = ROTATE_LEFT(a, s); a = (a + b) % (2 ^ 32); return a; } function ROTATE_LEFT(x, n, l, r) { l = bw_lshift(x, n) % (2 ^ 32); r = bw_rshift(x, 32  n); return (r + l); } function bw_not(x) { return bw_xor(x, 4294967295); # 0xffffffff } function bw_lshift(x, n) { return x * (2 ^ n); } function bw_rshift(x, n) { return int(x / (2 ^ n)); } function bw_and(x, y, i, r) { for (i = 0; i < 32; i += 4) { r = r / (2 ^ 4) + bw_lookup["and", x % 16, y % 16] * (2 ^ 28); x = int(x / (2 ^ 4)); y = int(y / (2 ^ 4)); } return r; } function bw_or(x, y, i, r) { for (i = 0; i < 32; i += 4) { r = r / (2 ^ 4) + bw_lookup["or", x % 16, y % 16] * (2 ^ 28); x = int(x / (2 ^ 4)); y = int(y / (2 ^ 4)); } return r; } function bw_xor(x, y) { return (x + y  2 * bw_and(x, y)); } # from https://www.gnu.org/software/gawk/manual/html_node/OrdinalFunctions.html function _ord_init( i) { for (i = 0; i < 256; i++) _ord_[sprintf("%c", i)] = i; } function ord(s) { # only first character is of interest return _ord_[substr(s, 1, 1)]; } function _bitwise_init( a, b, x, y, i) { # generate the bw_lookup table used by bw_and() and bw_or(). for (a = 0; a < 16; a++) { for (b = 0; b < 16; b++) { x = a; y = b; for (i = 0; i < 4; i++) { bw_lookup["and", a, b] += ((x % 2) && (y % 2)) * (2 ^ i); bw_lookup["or", a, b] += ((x % 2)  (y % 2)) * (2 ^ i); x = int(x / 2); y = int(y / 2); } } } } function _md5_init() { # MD5 shift constants setup. S11 = 7; S12 = 12; S13 = 17; S14 = 22; S21 = 5; S22 = 9; S23 = 14; S24 = 20; S31 = 4; S32 = 11; S33 = 16; S34 = 23; S41 = 6; S42 = 10; S43 = 15; S44 = 21; } 
diff u step7.awk step8.awk
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52   step7.awk 20170307 12:08:08.666000000 +0100 +++ step8.awk 20170307 15:16:41.668779000 +0100 @@ 34,9 +34,9 @@ # convert the array of bytes into an array of 32bits words. # NOTE: words is 0indexed. for (i = 1; i <= nbytes; i += 4) {  hi = bw_or(bw_lshift(bytes[i + 3], 8), bytes[i + 2]);  lo = bw_or(bw_lshift(bytes[i + 1], 8), bytes[i + 0]);  words[nwords++] = bw_or(bw_lshift(hi, 16), lo); + hi = bw_lshift(bytes[i + 3], 8) + bytes[i + 2]; + lo = bw_lshift(bytes[i + 1], 8) + bytes[i + 0]; + words[nwords++] = bw_lshift(hi, 16) + lo; } # Step 1. Append Padding Bits @@ 46,7 +46,7 @@ words[nwords++] = 128; # 0x80 } else { # append a '1' bit in the byte just after the last input byte.  words[nwords  1] = bw_or(words[nwords  1], bw_lshift(128, (nbytes % 4) * 8)); # 0x80 + words[nwords  1] = words[nwords  1] + bw_lshift(128, (nbytes % 4) * 8); # 0x80 } # "fill" the remaining bytes with 0 until we're just shy two words of # having 16Word Blocks. @@ 152,10 +152,10 @@ } for (i = j = 0; j < 16; j += 4) {  digest[j + 0] = bw_and(state[i], 255); # 0xff  digest[j + 1] = bw_and(bw_rshift(state[i], 8), 255); # 0xff  digest[j + 2] = bw_and(bw_rshift(state[i], 16), 255); # 0xff  digest[j + 3] = bw_and(bw_rshift(state[i++], 24), 255); # 0xff + digest[j + 0] = state[i] % (2 ^ 8); + digest[j + 1] = bw_rshift(state[i], 8) % (2 ^ 8); + digest[j + 2] = bw_rshift(state[i], 16) % (2 ^ 8); + digest[j + 3] = bw_rshift(state[i++], 24) % (2 ^ 8); } for (i = 0; i < 16; i++) ret = sprintf("%s%02x", ret, digest[i]); @@ 206,8 +206,10 @@ return a; } function ROTATE_LEFT(x, n) {  return bw_or((bw_lshift(x, n)) % (2 ^ 32), bw_rshift(x, 32  n)); +function ROTATE_LEFT(x, n, l, r) { + l = bw_lshift(x, n) % (2 ^ 32); + r = bw_rshift(x, 32  n); + return (r + l); } function bw_not(x) { 
30% faster although this was a tricky patch. Overall we've made very good progress on performances, going from about nine minutes to two and a half.
Wrapup
At that point I am happy with the result. I guess it could be even more optimized but I'm not sure where to start (we could inline more functions as we've seen that function call are costly but that will come with a heavy readability penalty).
Through this post we've implemented MD5 from the RFC in GNU Awk to simplify the problem at first, then refactored the code to port it to "standard" AWK, and finally made optimizations. You can also find the full code on GitHub. Thank you for reading and I hope you had as much fun as I did :)