#include <assert.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

static uint32_t f(uint32_t x, uint32_t y, uint32_t z) {
    return (x & y) | (~x & z);
}
#ifndef NDEBUG
static void TESTf(void) {
    assert(f(0, 0, 0) == 0);
    assert(f(0, 0, 1) == 1);
    assert(f(0, 1, 0) == 0);
    assert(f(0, 1, 1) == 1);
    assert(f(1, 0, 0) == 0);
    assert(f(1, 0, 1) == 0);
    assert(f(1, 1, 0) == 1);
    assert(f(1, 1, 1) == 1);
    assert(f(0x00000000, 0x00000000, 0x00000000) == 0x00000000);
    assert(f(0x00000000, 0x00000000, 0x80000000) == 0x80000000);
    assert(f(0x00000000, 0x80000000, 0x00000000) == 0x00000000);
    assert(f(0x00000000, 0x80000000, 0x80000000) == 0x80000000);
    assert(f(0x80000000, 0x00000000, 0x00000000) == 0x00000000);
    assert(f(0x80000000, 0x00000000, 0x80000000) == 0x00000000);
    assert(f(0x80000000, 0x80000000, 0x00000000) == 0x80000000);
    assert(f(0x80000000, 0x80000000, 0x80000000) == 0x80000000);
    assert(f(0xb8d95dbd, 0x917253f5, 0x3c579cfa) == 0x9456d1f7);
}
#endif

static uint32_t g(uint32_t x, uint32_t y, uint32_t z) {
    return (x & z) | (y & ~z);
}
#ifndef NDEBUG
static void TESTg(void) {
    assert(g(0, 0, 0) == 0);
    assert(g(0, 0, 1) == 0);
    assert(g(0, 1, 0) == 1);
    assert(g(0, 1, 1) == 0);
    assert(g(1, 0, 0) == 0);
    assert(g(1, 0, 1) == 1);
    assert(g(1, 1, 0) == 1);
    assert(g(1, 1, 1) == 1);
    assert(g(0x00000000, 0x00000000, 0x00000000) == 0x00000000);
    assert(g(0x00000000, 0x00000000, 0x80000000) == 0x00000000);
    assert(g(0x00000000, 0x80000000, 0x00000000) == 0x80000000);
    assert(g(0x00000000, 0x80000000, 0x80000000) == 0x00000000);
    assert(g(0x80000000, 0x00000000, 0x00000000) == 0x00000000);
    assert(g(0x80000000, 0x00000000, 0x80000000) == 0x80000000);
    assert(g(0x80000000, 0x80000000, 0x00000000) == 0x80000000);
    assert(g(0x80000000, 0x80000000, 0x80000000) == 0x80000000);
    assert(g(0xb8d95dbd, 0x917253f5, 0x3c579cfa) == 0xb9715fbd);
}
#endif

static uint32_t h(uint32_t x, uint32_t y, uint32_t z) {
    return x ^ y ^ z;
}
#ifndef NDEBUG
static void TESTh(void) {
    assert(h(0, 0, 0) == 0);
    assert(h(0, 0, 1) == 1);
    assert(h(0, 1, 0) == 1);
    assert(h(0, 1, 1) == 0);
    assert(h(1, 0, 0) == 1);
    assert(h(1, 0, 1) == 0);
    assert(h(1, 1, 0) == 0);
    assert(h(1, 1, 1) == 1);
    assert(h(0x00000000, 0x00000000, 0x00000000) == 0x00000000);
    assert(h(0x00000000, 0x00000000, 0x80000000) == 0x80000000);
    assert(h(0x00000000, 0x80000000, 0x00000000) == 0x80000000);
    assert(h(0x00000000, 0x80000000, 0x80000000) == 0x00000000);
    assert(h(0x80000000, 0x00000000, 0x00000000) == 0x80000000);
    assert(h(0x80000000, 0x00000000, 0x80000000) == 0x00000000);
    assert(h(0x80000000, 0x80000000, 0x00000000) == 0x00000000);
    assert(h(0x80000000, 0x80000000, 0x80000000) == 0x80000000);
    assert(h(0xb8d95dbd, 0x917253f5, 0x3c579cfa) == 0x15fc92b2);
}
#endif

static uint32_t i(uint32_t x, uint32_t y, uint32_t z) {
    return y ^ (x | ~z);
}
#ifndef NDEBUG
static void TESTi(void) {
    assert(i(0, 0, 0) == 0xffffffff);
    assert(i(0, 0, 1) == 0xfffffffe);
    assert(i(0, 1, 0) == 0xfffffffe);
    assert(i(0, 1, 1) == 0xffffffff);
    assert(i(1, 0, 0) == 0xffffffff);
    assert(i(1, 0, 1) == 0xffffffff);
    assert(i(1, 1, 0) == 0xfffffffe);
    assert(i(1, 1, 1) == 0xfffffffe);
    assert(i(0x00000000, 0x00000000, 0x00000000) == 0xffffffff);
    assert(i(0x00000000, 0x00000000, 0x80000000) == 0x7fffffff);
    assert(i(0x00000000, 0x80000000, 0x00000000) == 0x7fffffff);
    assert(i(0x00000000, 0x80000000, 0x80000000) == 0xffffffff);
    assert(i(0x80000000, 0x00000000, 0x00000000) == 0xffffffff);
    assert(i(0x80000000, 0x00000000, 0x80000000) == 0xffffffff);
    assert(i(0x80000000, 0x80000000, 0x00000000) == 0x7fffffff);
    assert(i(0x80000000, 0x80000000, 0x80000000) == 0x7fffffff);
    assert(i(0xb8d95dbd, 0x917253f5, 0x3c579cfa) == 0x6a8b2c48);
}
#endif

static uint32_t rotateleft(uint32_t v, uint32_t q) {
    if (q == 0) return v;
    return (uint32_t)(v << q) | (v >> (32 - q));
}
#ifndef NDEBUG
static void TESTrotate(void) {
    assert(rotateleft(0, 0) == 0);
    assert(rotateleft(0xface, 0) == 0xface);
    assert(rotateleft(0x5000, 4) == 0x50000);
    assert(rotateleft(0x5000, 5) == 0xa0000);
    assert(rotateleft(0x500000, 9) == 0xa0000000);
    assert(rotateleft(0x500000, 10) == 0x40000001);
    assert(rotateleft(0x500000, 11) == 0x80000002);
    assert(rotateleft(0x500000, 12) == 0x00000005);
    assert(rotateleft(0x500000, 13) == 0x0000000a);
    assert(rotateleft(0xf4f91816, 7) == 0x7c8c0b7a);
    assert(rotateleft(0xf4f91816, 13) == 0x2302de9f);
}
#endif

static void md5round(uint32_t *a, uint32_t b, uint32_t c, uint32_t d,
                     uint32_t k, uint32_t s, uint32_t t,
                     uint32_t (*F)(uint32_t, uint32_t, uint32_t)) {
    uint32_t tmp = *a;
    tmp += F(b, c, d);
    tmp += k;
    tmp += t;
    *a = b + rotateleft(tmp, s);
}
#ifndef NDEBUG
static void TESTround(void) {
    uint32_t a = 0x67452301;
    md5round(&a, 0xefcdab89, 0x98badcfe, 0x10325476,
             0x80000000, 7, 0xd76aa478, f);
    printf("1st round yields 0x%08x\n", a);
    assert(a == 0xefcdab89 + rotateleft((0x67452301 +
                f(0xefcdab89, 0x98badcfe, 0x10325476) +
                0x80000000 + 0xd76aa478), 7));
}
#endif

void md5mini(unsigned char *dstarr, const char *src) {
    // see https://www.ietf.org/rfc/rfc1321.txt
    static uint32_t T[64] = {0};
    if (T[0] == 0) {
        #if 0
        for (int p = 0; p < 64; p++) {
            T[p] = fabs(sin(1 + p)) * 4294967296;
            //printf("    T[%d] = 0x%08x;\n", p, T[p]);
        }
        #else
         T[0] = 0xd76aa478;  T[1] = 0xe8c7b756;  T[2] = 0x242070db;  T[3] = 0xc1bdceee;
         T[4] = 0xf57c0faf;  T[5] = 0x4787c62a;  T[6] = 0xa8304613;  T[7] = 0xfd469501;
         T[8] = 0x698098d8;  T[9] = 0x8b44f7af; T[10] = 0xffff5bb1; T[11] = 0x895cd7be;
        T[12] = 0x6b901122; T[13] = 0xfd987193; T[14] = 0xa679438e; T[15] = 0x49b40821;
        T[16] = 0xf61e2562; T[17] = 0xc040b340; T[18] = 0x265e5a51; T[19] = 0xe9b6c7aa;
        T[20] = 0xd62f105d; T[21] = 0x02441453; T[22] = 0xd8a1e681; T[23] = 0xe7d3fbc8;
        T[24] = 0x21e1cde6; T[25] = 0xc33707d6; T[26] = 0xf4d50d87; T[27] = 0x455a14ed;
        T[28] = 0xa9e3e905; T[29] = 0xfcefa3f8; T[30] = 0x676f02d9; T[31] = 0x8d2a4c8a;
        T[32] = 0xfffa3942; T[33] = 0x8771f681; T[34] = 0x6d9d6122; T[35] = 0xfde5380c;
        T[36] = 0xa4beea44; T[37] = 0x4bdecfa9; T[38] = 0xf6bb4b60; T[39] = 0xbebfbc70;
        T[40] = 0x289b7ec6; T[41] = 0xeaa127fa; T[42] = 0xd4ef3085; T[43] = 0x04881d05;
        T[44] = 0xd9d4d039; T[45] = 0xe6db99e5; T[46] = 0x1fa27cf8; T[47] = 0xc4ac5665;
        T[48] = 0xf4292244; T[49] = 0x432aff97; T[50] = 0xab9423a7; T[51] = 0xfc93a039;
        T[52] = 0x655b59c3; T[53] = 0x8f0ccc92; T[54] = 0xffeff47d; T[55] = 0x85845dd1;
        T[56] = 0x6fa87e4f; T[57] = 0xfe2ce6e0; T[58] = 0xa3014314; T[59] = 0x4e0811a1;
        T[60] = 0xf7537e82; T[61] = 0xbd3af235; T[62] = 0x2ad7d2bb; T[63] = 0xeb86d391;
        #endif
    }
    unsigned len = strlen(src);
    if (len > 55) exit(EXIT_FAILURE);
    uint8_t block[64] = {0};
    memcpy(block, src, len);
    block[len] = 0x80; // blocks [len+1], [len+2], etc are 0 from init
    block[56] = (len * 8) % 256;
    block[57] = (len * 8) / 256;
    uint32_t a = 0x67452301, b = 0xefcdab89, c = 0x98badcfe, d = 0x10325476;
    uint32_t x[16];
    for (int j = 0; j < 16; j++) {
        x[j] = (uint32_t)block[4*j] +
               (uint32_t)(block[4*j + 1] << 8) +
               (uint32_t)(block[4*j + 2] << 16) +
               (uint32_t)(block[4*j + 3] << 24);
    }
    uint32_t aa = a, bb = b, cc = c, dd = d;
    md5round(&a, b, c, d, x[ 0],  7, 0xd76aa478, f);
    md5round(&d, a, b, c, x[ 1], 12, 0xe8c7b756, f);
    md5round(&c, d, a, b, x[ 2], 17, 0x242070db, f);
    md5round(&b, c, d, a, x[ 3], 22, 0xc1bdceee, f);
    md5round(&a, b, c, d, x[ 4],  7, 0xf57c0faf, f);
    md5round(&d, a, b, c, x[ 5], 12, 0x4787c62a, f);
    md5round(&c, d, a, b, x[ 6], 17, 0xa8304613, f);
    md5round(&b, c, d, a, x[ 7], 22, 0xfd469501, f);
    md5round(&a, b, c, d, x[ 8],  7, 0x698098d8, f);
    md5round(&d, a, b, c, x[ 9], 12, 0x8b44f7af, f);
    md5round(&c, d, a, b, x[10], 17, 0xffff5bb1, f);
    md5round(&b, c, d, a, x[11], 22, 0x895cd7be, f);
    md5round(&a, b, c, d, x[12],  7, 0x6b901122, f);
    md5round(&d, a, b, c, x[13], 12, 0xfd987193, f);
    md5round(&c, d, a, b, x[14], 17, 0xa679438e, f);
    md5round(&b, c, d, a, x[15], 22, 0x49b40821, f);
    md5round(&a, b, c, d, x[ 1],  5, 0xf61e2562, g);
    md5round(&d, a, b, c, x[ 6],  9, 0xc040b340, g);
    md5round(&c, d, a, b, x[11], 14, 0x265e5a51, g);
    md5round(&b, c, d, a, x[ 0], 20, 0xe9b6c7aa, g);
    md5round(&a, b, c, d, x[ 5],  5, 0xd62f105d, g);
    md5round(&d, a, b, c, x[10],  9,  0x2441453, g);
    md5round(&c, d, a, b, x[15], 14, 0xd8a1e681, g);
    md5round(&b, c, d, a, x[ 4], 20, 0xe7d3fbc8, g);
    md5round(&a, b, c, d, x[ 9],  5, 0x21e1cde6, g);
    md5round(&d, a, b, c, x[14],  9, 0xc33707d6, g);
    md5round(&c, d, a, b, x[ 3], 14, 0xf4d50d87, g);
    md5round(&b, c, d, a, x[ 8], 20, 0x455a14ed, g);
    md5round(&a, b, c, d, x[13],  5, 0xa9e3e905, g);
    md5round(&d, a, b, c, x[ 2],  9, 0xfcefa3f8, g);
    md5round(&c, d, a, b, x[ 7], 14, 0x676f02d9, g);
    md5round(&b, c, d, a, x[12], 20, 0x8d2a4c8a, g);
    md5round(&a, b, c, d, x[ 5],  4, 0xfffa3942, h);
    md5round(&d, a, b, c, x[ 8], 11, 0x8771f681, h);
    md5round(&c, d, a, b, x[11], 16, 0x6d9d6122, h);
    md5round(&b, c, d, a, x[14], 23, 0xfde5380c, h);
    md5round(&a, b, c, d, x[ 1],  4, 0xa4beea44, h);
    md5round(&d, a, b, c, x[ 4], 11, 0x4bdecfa9, h);
    md5round(&c, d, a, b, x[ 7], 16, 0xf6bb4b60, h);
    md5round(&b, c, d, a, x[10], 23, 0xbebfbc70, h);
    md5round(&a, b, c, d, x[13],  4, 0x289b7ec6, h);
    md5round(&d, a, b, c, x[ 0], 11, 0xeaa127fa, h);
    md5round(&c, d, a, b, x[ 3], 16, 0xd4ef3085, h);
    md5round(&b, c, d, a, x[ 6], 23,  0x4881d05, h);
    md5round(&a, b, c, d, x[ 9],  4, 0xd9d4d039, h);
    md5round(&d, a, b, c, x[12], 11, 0xe6db99e5, h);
    md5round(&c, d, a, b, x[15], 16, 0x1fa27cf8, h);
    md5round(&b, c, d, a, x[ 2], 23, 0xc4ac5665, h);
    md5round(&a, b, c, d, x[ 0],  6, 0xf4292244, i);
    md5round(&d, a, b, c, x[ 7], 10, 0x432aff97, i);
    md5round(&c, d, a, b, x[14], 15, 0xab9423a7, i);
    md5round(&b, c, d, a, x[ 5], 21, 0xfc93a039, i);
    md5round(&a, b, c, d, x[12],  6, 0x655b59c3, i);
    md5round(&d, a, b, c, x[ 3], 10, 0x8f0ccc92, i);
    md5round(&c, d, a, b, x[10], 15, 0xffeff47d, i);
    md5round(&b, c, d, a, x[ 1], 21, 0x85845dd1, i);
    md5round(&a, b, c, d, x[ 8],  6, 0x6fa87e4f, i);
    md5round(&d, a, b, c, x[15], 10, 0xfe2ce6e0, i);
    md5round(&c, d, a, b, x[ 6], 15, 0xa3014314, i);
    md5round(&b, c, d, a, x[13], 21, 0x4e0811a1, i);
    md5round(&a, b, c, d, x[ 4],  6, 0xf7537e82, i);
    md5round(&d, a, b, c, x[11], 10, 0xbd3af235, i);
    md5round(&c, d, a, b, x[ 2], 15, 0x2ad7d2bb, i);
    md5round(&b, c, d, a, x[ 9], 21, 0xeb86d391, i);
    a += aa; b += bb; c += cc; d += dd;
    if (dstarr) {
        uint32_t aaa = a, bbb = b, ccc = c, ddd = d;
        *dstarr++ = aaa % 256; aaa /= 256;
        *dstarr++ = aaa % 256; aaa /= 256;
        *dstarr++ = aaa % 256; aaa /= 256;
        *dstarr++ = aaa;
        *dstarr++ = bbb % 256; bbb /= 256;
        *dstarr++ = bbb % 256; bbb /= 256;
        *dstarr++ = bbb % 256; bbb /= 256;
        *dstarr++ = bbb;
        *dstarr++ = ccc % 256; ccc /= 256;
        *dstarr++ = ccc % 256; ccc /= 256;
        *dstarr++ = ccc % 256; ccc /= 256;
        *dstarr++ = ccc;
        *dstarr++ = ddd % 256; ddd /= 256;
        *dstarr++ = ddd % 256; ddd /= 256;
        *dstarr++ = ddd % 256; ddd /= 256;
        *dstarr++ = ddd;
    }
}

#ifndef NDEBUG
void md5mini_run_all_tests(void) {
    TESTf();
    TESTg();
    TESTh();
    TESTi();
    TESTrotate();
    TESTround();
    assert(1 && "ALL TESTS PASSED");
}
#endif