A pure-Python implementation of a subset of the Skein hash function. Written to see if I understood the spec.

skein.py

from __future__ import print_function, division
from math import ceil
from operator import xor, getitem
from functools import reduce
from struct import unpack, pack
two64 = 0x10000000000000000

def hexToBytes(hexes):
    hexes = hexes.replace(' ', '')
    return bytes(int(hexes[i:i+2], 16) for i in range(0, len(hexes), 2))

def bytesToWords(bytes):
    return list(unpack('<' + str(len(bytes) // 8) + 'Q', bytes))

def wordsToBytes(words):
    return pack('<' + str(len(words)) + 'Q', *words)

def hexPrint(data, start = '', byte = True):
    string = isinstance(data, str) or isinstance(data, bytes)
    if byte:
        print(start, ''.join('{:x}'.format(ord(byte)).zfill(2)
                             for byte in (data if string else wordsToBytes(data))))
    else:
        if string: data += '\0' * (-len(data) % 8)
        print(start, *('0x' + '{:x}'.format(ord(word)).zfill(16)
                       for word in (bytesToWords(data) if string else data)))

#first two typechecks not needed right now
#third is because v is always passed as bytes / buffer from UBI
@staticmethod
def encrypt(cls, k, t, v, append=list.append):
    if type(k) is str: k = bytesToWords(k)
    if type(t) is str: t = bytesToWords(t)
    if type(v) is str or type(v) is bytes: v = bytesToWords(v)
    R, mixes, permute, keys = cls.R, cls.mixes, cls.permute, cls.keys
    s, p, N = 0, v, cls.numwords + 1
    append(k, 0x1BD11BDAA9FC1A22 ^ reduce(xor, k))
    append(t, t[0] ^ t[1])
    for d in cls.rounds:
        if not d % 4:
            si = s + N - 1
            subkeys = [k[(s + i) % N] for i in keys]
            append(subkeys, (k[(si - 3) % N] + t[s % 3]) % two64)
            append(subkeys, (k[(si - 2) % N] + t[(s + 1) % 3]) % two64)
            append(subkeys, (k[(si - 1) % N] + s) % two64)
            v = [(vi + sk) % two64 for vi, sk in zip(v, subkeys)]
            s += 1
        Rd = R[d % 8]
        for j in mixes:
            Rdj, v1 = Rd[j // 2], v[j + 1]
            v[j] = (v[j] + v1) % two64
            v[j + 1] = (((v1 << Rdj) % two64) + (v1 >> (64 - Rdj))) ^ v[j]
        v = permute(v)
    si = s + N - 1
    subkeys = [k[(s + i) % N] for i in keys]
    append(subkeys, (k[(si - 3) % N] + t[s % 3]) % two64)
    append(subkeys, (k[(si - 2) % N] + t[(s + 1) % 3]) % two64)
    append(subkeys, (k[(si - 1) % N] + s) % two64)
    v = [(vi + sk) % two64 for vi, sk in zip(v, subkeys)]
    v = list(map(xor, v, p))
    return v

class ThreeFish(type):
    def __new__(cls, name, bases, attrs):
        attrs['mixes'] = tuple(range(0, attrs['numwords'], 2))
        attrs['rounds'] = tuple(range(attrs['numrounds']))
        attrs['keys'] = tuple(range(attrs['numwords'] - 3))
        attrs['__new__'] = encrypt
        return super(ThreeFish, cls).__new__(cls, name, bases, attrs)

class ThreeFish256(metaclass=ThreeFish):
    numwords, numrounds = 4, 72
    R = ((14, 16), (52, 57), (23, 40), (5, 37),
         (25, 33), (46, 12), (58, 22), (32, 32))
    testin = [0] * 4, [0, 0], [0] * 4
    testout = [0x94EEEA8B1F2ADA84, 0xADF103313EAE6670,
               0x952419A1F4B16D53, 0xD83F13E63C9F6B11]
    @staticmethod
    def permute(v):
        return [v[0], v[3], v[2], v[1]]

class ThreeFish512(metaclass=ThreeFish):
    __metaclass__ = ThreeFish
    numwords, numrounds = 8, 72
    R = ((46, 36, 19, 37), (33, 27, 14, 42), (17, 49, 36, 39), (44, 9, 54, 56),
        (39, 30, 34, 24), (13, 50, 10, 17), (25, 29, 39, 43), (8, 35, 56, 22))
    testin = [0] * 8, [0, 0], [0] * 8
    testout =  [0xBC2560EFC6BBA2B1, 0xE3361F162238EB40,
                0xFB8631EE0ABBD175, 0x7B9479D4C5479ED1,
                0xCFF0356E58F8C27B, 0xB1B7B08430F0E7F7,
                0xE9A380A56139ABF1, 0xBE7B6D4AA11EB47E]
    @staticmethod
    def permute(v):
        return [v[2], v[1], v[4], v[7], v[6], v[5], v[0], v[3]]

class ThreeFish1024(metaclass=ThreeFish):
    __metaclass__ = ThreeFish
    numwords, numrounds = 16, 80
    R = ((24, 13, 8, 47, 8, 17, 22, 37), (38, 19, 10, 55, 49, 18, 23, 52),
        (33, 4, 51, 13, 34, 41, 59, 17), (5, 20, 48, 41, 47, 28, 16, 25),
        (41, 9, 37, 31, 12, 47, 44, 30), (16, 34, 56, 51, 4, 53, 42, 41),
        (31, 44, 47, 46, 19, 42, 44, 25), (9, 48, 35, 52, 23, 31, 37, 20))
    testin = [0] * 16, [0, 0], [0] * 16
    testout = [0x04B3053D0A3D5CF0, 0x0136E0D1C7DD85F7,
               0x067B212F6EA78A5C, 0x0DA9C10B4C54E1C6,
               0x0F4EC27394CBACF0, 0x32437F0568EA4FD5,
               0xCFF56D1D7654B49C, 0xA2D5FB14369B2E7B,
               0x540306B460472E0B, 0x71C18254BCEA820D,
               0xC36B4068BEAF32C8, 0xFA4329597A360095,
               0xC4A36C28434A5B9A, 0xD54331444B1046CF,
               0xDF11834830B2A460, 0x1E39E8DFE1F7EE4F]
    @staticmethod
    def permute(v):
        return [v[0], v[9], v[2], v[13], v[6], v[11], v[4], v[15],
                v[10], v[7], v[12], v[3], v[14], v[5], v[8], v[1]]

class ThreeFish1024_ReducedRounds(ThreeFish1024):
    numwords, numrounds = 16, 5

threefishes = {256: ThreeFish256, 512: ThreeFish512, 1024: ThreeFish1024}

#convert M to words so can feedforward here
#also to remove need for threefish to support bytes
#this MUST take both words and bytes
#divide up M differently, don't bytepad, Nb vs words?
def ubi(H, M, Ts, encrypt, two64=two64):
    if not type(M) is bytes: M = wordsToBytes(M)
    if type(Ts) is str:
        position, positionB, treelevel, usetype = unpack('<QIxxBB', Ts)
        position += positionB * two64
    else:
        position, usetype = Ts
        position += (usetype % 0x100000000) * two64
        treelevel = (usetype >> 48) % 0x100
        usetype = (usetype >> 56) % 0x100
    assert treelevel < 128
    assert usetype < 64
    Nb, Nm, bitpad = len(H), len(M), 0
    if not type(H) is str: Nb *= 8
    if not Nm: M = b'\0' * Nb
    else: M += b'\0' * (-Nm % Nb)
    maxpos, ttr = (position + Nm, usetype * 0x100000000000000 +
                     (bitpad + treelevel) * 0x1000000000000)
    assert not maxpos > 0xffffffffffffffffffffffff
    M = tuple(M[i:i+Nb] for i in range(0, Nm, Nb))
    last_len = Nm % Nb
    last_tweak = 0xC000000000000000
    if maxpos < two64:
        if Nm > Nb:
            position += Nb
            H = encrypt(H, [position, ttr + 0x4000000000000000], M[0])
            for Mi in M[1:-1]:
                position += Nb
                H = encrypt(H, [position, ttr], Mi)
            last_tweak = 0x8000000000000000
        position += last_len or Nb
        return encrypt(H, [position, ttr + last_tweak], M[-1])
    if Nm > Nb:
        position += Nb
        H = encrypt(H, [position % two64,
                        (position >> 64) + ttr + 0x4000000000000000], M[0])
        for Mi in M[1:-1]:
            position += Nb
            H = encrypt(H, [position % two64, (position >> 64) + ttr], Mi)
        last_tweak = 0x8000000000000000
    position += last_len or Nb
    return encrypt(H, [position % two64,
                       (position >> 64) + ttr + last_tweak], M[-1])

msgtypes = {'key': 0, 'cfg': 0x0400000000000000,
                      'prs': 0x0800000000000000,
                      'PK':  0x0C00000000000000,
                      'kdf': 0x1000000000000000,
                      'non': 0x1400000000000000,
                      'msg': 0x3000000000000000,
                      'out': 0x3f00000000000000}

def simpleskein(bitsout, bits, M):
    assert not bitsout % 8
    #print('Skein', bits, 'Simple Hash')
    #hexPrint(M, 'Message:')
    encrypt = threefishes[bits]
    #This magic number is 'SHA3\01'
    G = ubi([0] * (bits >> 6), (0x0133414853, bitsout, 0, 0),
            (0, msgtypes['cfg']), encrypt)
    G = ubi(G, M, (0, msgtypes['msg']), encrypt)
    Tout = (0, msgtypes['out'])
    G = sum((ubi(G, pack('<Q', i), Tout, encrypt)
            for i in range(int(ceil(bitsout / bits)))), [])
    return wordsToBytes(G)[0:(bitsout // 8)]

plain = hexToBytes("FF FE FD FC FB FA F9 F8 F7 F6 F5 F4 F3 F2 F1 F0" +
    "EF EE ED EC EB EA E9 E8 E7 E6 E5 E4 E3 E2 E1 E0" +
    "DF DE DD DC DB DA D9 D8 D7 D6 D5 D4 D3 D2 D1 D0" +
    "CF CE CD CC CB CA C9 C8 C7 C6 C5 C4 C3 C2 C1 C0")
hash = hexToBytes("df28e916630d0b44c4a849dc9a02f07a07cb30f732318256b15d865ac4ae162f")

if __name__ == '__main__':
    debug = False
    ThreeFish1024_ReducedRounds(*ThreeFish1024.testin)
    for bits, threefish in threefishes.items():
        #hexPrint(threefish.testout, 'ThreeFish ' + str(bits) + '\nTest      :')
        output = threefish(*threefish.testin)
        assert output == threefish.testout
    assert simpleskein(256, 256, plain) == hash
    print("Success!")

I did a quick & dirty update to make this Python 3 compatible as I still get emails about this periodically. NOT UPDATED WITH NEW VERSIONS OF THE SKEIN SPEC.