diff --git a/src/gsm/Makefile.am b/src/gsm/Makefile.am
index ca850d2..cbbdcc0 100644
--- a/src/gsm/Makefile.am
+++ b/src/gsm/Makefile.am
@@ -15,6 +15,8 @@
 		 milenage/common.h milenage/crypto.h milenage/includes.h \
 		 milenage/milenage.h
 
+noinst_HEADERS += tuak/KeccakP-1600-3gpp.h tuak/tuak.h
+
 noinst_LTLIBRARIES = libgsmint.la
 lib_LTLIBRARIES = libosmogsm.la
 
@@ -30,6 +32,7 @@
 			auth_milenage.c milenage/aes-encblock.c gea.c \
 			milenage/aes-internal.c milenage/aes-internal-enc.c \
 			milenage/milenage.c gan.c ipa.c gsm0341.c apn.c \
+			tuak/KeccakP-1600-3gpp.c tuak/tuak.c auth_tuak.c \
 			gsup.c gsup_sms.c gprs_gea.c gsm0503_conv.c oap.c gsm0808_utils.c \
 			gsm23003.c gsm23236.c mncc.c bts_features.c oap_client.c \
 			gsm29118.c gsm48_rest_octets.c cbsp.c gsm48049.c \
diff --git a/src/gsm/auth_core.c b/src/gsm/auth_core.c
index ad5967d..6972cb7 100644
--- a/src/gsm/auth_core.c
+++ b/src/gsm/auth_core.c
@@ -328,6 +328,7 @@
 	{ OSMO_AUTH_ALG_XOR_3G, "XOR-3G" },
 	{ OSMO_AUTH_ALG_MILENAGE, "MILENAGE" },
 	{ OSMO_AUTH_ALG_XOR_2G, "XOR-2G" },
+	{ OSMO_AUTH_ALG_TUAK, "TUAK" },
 	{ 0, NULL }
 };
 
diff --git a/src/gsm/auth_tuak.c b/src/gsm/auth_tuak.c
new file mode 100644
index 0000000..05fbf69
--- /dev/null
+++ b/src/gsm/auth_tuak.c
@@ -0,0 +1,207 @@
+/*! \file auth_tuak.c
+ * GSM/GPRS/3G authentication core infrastructure */
+/*
+ * (C) 2023 by Harald Welte <laforge@osmocom.org>
+ *
+ * All Rights Reserved
+ *
+ * SPDX-License-Identifier: GPL-2.0+
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ */
+
+/* NOTE: TUAK offers a lot of size variability in terms of size of length of MAC_A, MAC_S,
+ * but this is not used within 3GPP.  The different sizes of Kc and RES are handled via
+ * osmo_sub_auth_data2. */
+
+#include <errno.h>
+#include <osmocom/crypt/auth.h>
+#include <osmocom/core/bits.h>
+#include "tuak/tuak.h"
+
+/*! \addtogroup auth
+ *  @{
+ */
+
+static const uint8_t *gen_opc_if_needed(const struct osmo_sub_auth_data2 *aud, uint8_t *gen_opc)
+{
+	int rc;
+
+	/* Check if we only know OP and compute OPC if required */
+	if (aud->type == OSMO_AUTH_TYPE_UMTS && aud->u.umts.opc_is_op) {
+		rc = tuak_opc_gen(gen_opc, aud->u.umts.k, aud->u.umts.k_len, aud->u.umts.opc);
+		if (rc < 0)
+			return NULL;
+		return gen_opc;
+	}
+
+	return aud->u.umts.opc;
+}
+
+static int tuak_gen_vec(struct osmo_auth_vector *vec,
+			struct osmo_sub_auth_data2 *aud,
+			const uint8_t *_rand)
+{
+	size_t res_len = vec->res_len;
+	uint64_t next_sqn;
+	uint8_t gen_opc[32];
+	const uint8_t *opc;
+	uint8_t sqn[6];
+	uint64_t ind_mask;
+	uint64_t seq_1;
+
+	switch (vec->res_len) {
+	case 4:
+	case 8:
+	case 16:
+		break;
+	default:
+		return -EINVAL;
+	}
+
+	OSMO_ASSERT(aud->algo == OSMO_AUTH_ALG_TUAK);
+
+	if (aud->u.umts.k_len != 16 && aud->u.umts.k_len != 32)
+		return -EINVAL;
+	if (aud->u.umts.opc_len != 32)
+		return -EINVAL;
+
+	opc = gen_opc_if_needed(aud, gen_opc);
+	if (!opc)
+		return -1;
+
+	/* Determine next SQN, according to 3GPP TS 33.102:
+	 * SQN consists of SEQ and a lower significant part of IND bits:
+	 *
+	 * |----------SEQ------------|
+	 * |------------------------SQN-----------|
+	 *                           |-----IND----|
+	 *
+	 * The IND part is used as "slots": e.g. a given HLR client will always
+	 * get the same IND part, called ind here, with incrementing SEQ. In
+	 * the USIM, each IND slot enforces that its SEQ are used in ascending
+	 * order -- as long as that constraint is satisfied, the SQN may jump
+	 * forwards and backwards. For example, for ind_bitlen == 5, asking the
+	 * USIM for SQN = 32, 64, 33 is allowed, because 32 and 64 are
+	 * SEQ || (ind == 0), and though 33 is below 64, it is ind == 1 and
+	 * allowed.  Not allowed would be 32, 96, 64, because 64 would go
+	 * backwards after 96, both being ind == 0.
+	 *
+	 * From the last used SQN, we want to increment SEQ + 1, and then pick
+	 * the matching IND part.
+	 *
+	 * IND size is suggested in TS 33.102 as 5 bits. SQN is 48 bits long.
+	 * If ind_bitlen is passed too large here, the algorithms will break
+	 * down. But at which point should we return an error? A sane limit
+	 * seems to be ind_bitlen == 10, but to protect against failure,
+	 * limiting ind_bitlen to 28 is enough, 28 being the number of bits
+	 * suggested for the delta in 33.102, which is discussed to still
+	 * require 2^15 > 32000 authentications to wrap the SQN back to the
+	 * start.
+	 *
+	 * Note that if a caller with ind == 1 generates N vectors, the SQN
+	 * stored after this will reflect SEQ + N. If then another caller with
+	 * ind == 2 generates another N vectors, this will then use SEQ + N
+	 * onwards and end up with SEQ + N + N. In other words, most of each
+	 * SEQ's IND slots will remain unused. When looking at SQN being 48
+	 * bits wide, after dropping ind_bitlen (say 5) from it, we will still
+	 * have a sequence range of 2^43 = 8.8e12, eight trillion sequences,
+	 * which is large enough to not bother further. With the maximum
+	 * ind_bitlen of 28 enforced below, we still get more than 1 million
+	 * sequences, which is also sufficiently large.
+	 *
+	 * An ind_bitlen of zero may be passed from legacy callers that are not
+	 * aware of the IND extension. For these, below algorithm works out as
+	 * before, simply incrementing SQN by 1.
+	 *
+	 * This is also a mechanism for tools like the osmo-auc-gen to directly
+	 * request a given SQN to be used. With ind_bitlen == 0 the caller can
+	 * be sure that this code will increment SQN by exactly one before
+	 * generating a tuple, thus a caller would simply pass
+	 * { .ind_bitlen = 0, .ind = 0, .sqn = (desired_sqn - 1) }
+	 */
+
+	if (aud->u.umts.ind_bitlen > OSMO_MILENAGE_IND_BITLEN_MAX)
+		return -2;
+
+	seq_1 = 1LL << aud->u.umts.ind_bitlen;
+	ind_mask = ~(seq_1 - 1);
+
+	/* the ind index must not affect the SEQ part */
+	if (aud->u.umts.ind >= seq_1)
+		return -3;
+
+	/* keep the incremented SQN local until gsm_milenage() succeeded. */
+	next_sqn = ((aud->u.umts.sqn + seq_1) & ind_mask) + aud->u.umts.ind;
+
+	osmo_store64be_ext(next_sqn, sqn, 6);
+
+	tuak_generate(opc, aud->u.umts.amf, aud->u.umts.k, aud->u.umts.k_len,
+		      sqn, _rand, vec->autn, vec->ik, vec->ck, vec->res, &res_len);
+
+	/* generate the GSM Kc + SRES values using C2 + C3 functions */
+	osmo_auth_c3(vec->kc, vec->ck, vec->ik);
+	osmo_auth_c2(vec->sres, vec->res, vec->res_len, 1);
+
+	vec->auth_types = OSMO_AUTH_TYPE_UMTS | OSMO_AUTH_TYPE_GSM;
+
+	/* for storage in the caller's AUC database */
+	aud->u.umts.sqn = next_sqn;
+
+	return 0;
+}
+
+static int tuak_gen_vec_auts(struct osmo_auth_vector *vec,
+			     struct osmo_sub_auth_data2 *aud,
+			     const uint8_t *auts, const uint8_t *rand_auts,
+			     const uint8_t *_rand)
+{
+	uint8_t sqn_out[6];
+	uint8_t gen_opc[32];
+	const uint8_t *opc;
+	int rc;
+
+	OSMO_ASSERT(aud->algo == OSMO_AUTH_ALG_TUAK);
+
+	if (aud->u.umts.k_len != 16 && aud->u.umts.k_len != 32)
+		return -EINVAL;
+	if (aud->u.umts.opc_len != 32)
+		return -EINVAL;
+
+	opc = gen_opc_if_needed(aud, gen_opc);
+
+	rc = tuak_auts(opc, aud->u.umts.k, sizeof(aud->u.umts.k), rand_auts, auts, sqn_out);
+	if (rc < 0)
+		return rc;
+
+	aud->u.umts.sqn_ms = osmo_load64be_ext(sqn_out, 6) >> 16;
+	/* Update our "largest used SQN" from the USIM -- milenage_gen_vec()
+	 * below will increment SQN. */
+	aud->u.umts.sqn = aud->u.umts.sqn_ms;
+
+	return tuak_gen_vec(vec, aud, _rand);
+}
+
+static struct osmo_auth_impl tuak_alg = {
+	.algo = OSMO_AUTH_ALG_TUAK,
+	.name = "TUAK (libosmogsm built-in)",
+	.priority = 1000,
+	.gen_vec = &tuak_gen_vec,
+	.gen_vec_auts = &tuak_gen_vec_auts,
+};
+
+static __attribute__((constructor)) void on_dso_load_tuak(void)
+{
+	osmo_auth_register(&tuak_alg);
+}
+
+/*! @} */
diff --git a/src/gsm/tuak/KeccakP-1600-3gpp.c b/src/gsm/tuak/KeccakP-1600-3gpp.c
new file mode 100644
index 0000000..3f5e2ad
--- /dev/null
+++ b/src/gsm/tuak/KeccakP-1600-3gpp.c
@@ -0,0 +1,176 @@
+/* -----------------------------------------------------------------------
+ * code extracted from 3GPP TS 35.231, annex E for Keccak core functions
+ * https://portal.3gpp.org/desktopmodules/Specifications/SpecificationDetails.aspx?specificationId=2402
+ *-----------------------------------------------------------------------*/
+
+/* This code may be freely used or adapted.
+*/
+
+#include "KeccakP-1600-3gpp.h"
+
+
+const uint8_t Rho[25]		= {0,1,62,28,27,36,44,6,55,20,3,10,43,25,39,41,45,
+   15,21,8,18,2,61,56,14};
+
+const uint8_t Pi[25]		= {0,6,12,18,24,3,9,10,16,22,1,7,13,19,20,4,5,11,17,
+   23,2,8,14,15,21};
+
+const uint8_t Iota[24]	= {1,146,218,112,155,33,241,89,138,136,57,42,187,203,
+   217,83,82,192,26,106,241,208,33,120};
+
+#define ROTATE64(value, n)	\
+((((uint64_t)(value))<<(n)) | (((uint64_t)(value))>>(64-(n))))
+
+/* ---------------------------------------------------------------------
+   64-bit version of Keccak_f(1600)
+   ---------------------------------------------------------------------
+*/
+void Keccak_f_64(uint64_t s[25])
+{	uint64_t t[5];
+	uint8_t i, j, round;
+
+	for(round=0; round<24; ++round)
+	{	/* Theta function */
+		for(i=0; i<5; ++i)
+			t[i] = s[i] ^ s[5+i] ^ s[10+i] ^ s[15+i] ^ s[20+i];
+		for(i=0; i<5; ++i, s+=5)
+		{	s[0] ^= t[4] ^ ROTATE64(t[1], 1);
+			s[1] ^= t[0] ^ ROTATE64(t[2], 1);
+			s[2] ^= t[1] ^ ROTATE64(t[3], 1);
+			s[3] ^= t[2] ^ ROTATE64(t[4], 1);
+			s[4] ^= t[3] ^ ROTATE64(t[0], 1);
+		}
+		s -= 25;
+
+		/* Rho function */
+		for(i=1; i<25; ++i)
+			s[i] = ROTATE64(s[i], Rho[i]);
+
+		/* Pi function */
+		for(t[1] = s[i=1]; (j=Pi[i]) > 1; s[i]=s[j], i=j);
+		s[i] = t[1];
+
+		/* Chi function */
+		for(i=0; i<5; ++i, s += 5)
+		{	t[0] = (~s[1]) & s[2];
+			t[1] = (~s[2]) & s[3];
+			t[2] = (~s[3]) & s[4];
+			t[3] = (~s[4]) & s[0];
+			t[4] = (~s[0]) & s[1];
+			for(j=0; j<5; ++j) s[j] ^= t[j];
+		}
+		s -= 25;
+
+		/* Iota function */
+		t[0] = Iota[round];
+		*s ^= (t[0] | (t[0]<<11) | (t[0]<<26) | (t[0]<<57)) 
+              & 0x800000008000808BULL; /* set & mask bits 0,1,3,7,15,31,63 */
+	}
+}
+
+
+/* ---------------------------------------------------------------------
+   8-bit version of Keccak_f(1600)
+   ---------------------------------------------------------------------
+*/
+void Keccak_f_8(uint8_t s[200])
+{	uint8_t t[40], i, j, k, round;
+
+	for(round=0; round<24; ++round)
+	{	/* Theta function */
+		for(i=0; i<40; ++i)
+			t[i]=s[i]^s[40+i]^s[80+i]^s[120+i]^s[160+i];
+		for(i=0; i<200; i+=8)
+			for(j = (i+32)%40, k=0; k<8; ++k)
+				s[i+k] ^= t[j+k];
+		for(i=0; i<40; t[i] = (t[i]<<1)|j, i+=8)
+			for(j = t[i+7]>>7, k=7; k; --k)
+				t[i+k] = (t[i+k]<<1)|(t[i+k-1]>>7);
+		for(i=0; i<200; i+=8)
+			for(j = (i+8)%40, k=0; k<8; ++k)
+				s[i+k] ^= t[j+k];
+
+		/* Rho function */
+		for(i=8; i<200; i+=8)
+		{	for(j = Rho[i>>3]>>3, k=0; k<8; ++k) 	/* j:=bytes to shift, s->t 		*/
+				t[(k+j)&7] = s[i+k];
+			for(j = Rho[i>>3]&7, k=7; k; --k) 	   /* j:=bits  to shift, t->s 	*/
+				s[i+k] = (t[k]<<j) | (t[k-1]>>(8-j));
+			s[i] = (t[0]<<j) | (t[7]>>(8-j));
+		}
+
+		/* Pi function */
+		for(k=8; k<16; ++k) t[k] = s[k];		/* =memcpy(t+8, s+8, 8) 				 	*/
+		for(i=1; (j=Pi[i])>1; i=j)
+			for(k=0; k<8; ++k)						/* =memcpy(s+(i<<3), s+(j<<3), 8)	*/
+				s[(i<<3)|k] = s[(j<<3)|k];
+		for(k=0; k<8; ++k)							/* =memcpy(s+(i<<3), t+8, 8) 		 	*/
+			s[(i<<3)|k] = t[k+8];
+
+		/* Chi function */
+		for(i=0; i<200; i+=40)
+		{	for(j=0; j<40; ++j)
+				t[j]=(~s[i+(j+8)%40]) & s[i+(j+16)%40];
+			for(j=0; j<40; ++j)	s[i+j]^=t[j];
+		}
+
+		/* Iota function */
+		k = Iota[round];
+		s[0] ^= k & 0x8B;			/* bits 0, 1, 3, 7 */
+		s[1] ^= (k<<3)&0x80;		/* bit 15 */
+		s[3] ^= (k<<2)&0x80;		/* bit 31 */
+		s[7] ^= (k<<1)&0x80;		/* bit 63 */
+
+	}
+}
+
+/* ---------------------------------------------------------------------
+   32-bit version of Keccak_f(1600)
+   ---------------------------------------------------------------------
+*/
+void Keccak_f_32(uint32_t s[50])
+{	uint32_t t[10];
+	uint8_t i, j, round, k;
+
+	for(round=0; round<24; ++round)
+	{	/* Theta function */
+		for(i=0; i<10; ++i)
+			t[i] = s[i] ^ s[10+i] ^ s[20+i] ^ s[30+i] ^ s[40+i];
+		for(i=0; i<5; ++i)
+			for(j=8, k=2; ; j%=10, k=(k+2)%10)
+			{	*s++ ^= t[j++] ^ ((t[k]<<1)|(t[k+1]>>31));
+				*s++ ^= t[j++] ^ ((t[k+1]<<1)|(t[k]>>31));
+				if(j==8) break;
+			}
+		s -= 50;
+
+		/* Rho function */
+		for(i=2; i<50; i+=2)
+		{	k = Rho[i>>1] & 0x1f;
+			t[0] = (s[i+1] << k) | (s[i] >> (32-k));
+			t[1] = (s[i] << k) | (s[i+1] >> (32-k));
+			k = Rho[i>>1] >> 5;
+			s[i] = t[1-k], s[i+1] = t[k];
+		}
+
+		/* Pi function */
+		for(i=2, t[0]=s[2], t[1]=s[3]; (j=(Pi[i>>1]<<1))>2; i=j)
+			s[i]=s[j], s[i+1]=s[j+1];
+		s[i]=t[0], s[i+1]=t[1];
+
+		/* Chi function */
+		for(i=0; i<5; ++i, s+=10)
+		{	for(j=0; j<10; ++j)
+				t[j] = (~s[(j+2)%10]) & s[(j+4)%10];
+			for(j=0; j<10; ++j)
+				s[j] ^= t[j];
+		}
+		s -= 50;
+
+		/* Iota function */
+		t[0] = Iota[round];
+		s[0] ^= (t[0] | (t[0]<<11) | (t[0]<<26)) & 0x8000808B;
+		s[1] ^= (t[0]<<25) & 0x80000000;
+	}
+}
+
diff --git a/src/gsm/tuak/KeccakP-1600-3gpp.h b/src/gsm/tuak/KeccakP-1600-3gpp.h
new file mode 100644
index 0000000..a23cc46
--- /dev/null
+++ b/src/gsm/tuak/KeccakP-1600-3gpp.h
@@ -0,0 +1,25 @@
+/* -----------------------------------------------------------------------
+ * code extracted from 3GPP TS 35.231, annex E for Keccak core functions
+ * https://portal.3gpp.org/desktopmodules/Specifications/SpecificationDetails.aspx?specificationId=2402
+ *-----------------------------------------------------------------------*/
+
+/* this is the trick to make the code cross-platform
+ * at least, Win32 / Linux */
+
+#if defined(_WIN32) || defined(__WIN32__)
+#	include <windows.h>
+#	define EXPORTIT __declspec(dllexport)
+#else
+#	define EXPORTIT
+#endif
+
+#include <stdint.h>
+
+/*------------------------------------------------------------------------
+ * KeccakP-1600-3gpp.h
+ *------------------------------------------------------------------------*/
+
+EXPORTIT void Keccak_f_8 (uint8_t s[200]);
+EXPORTIT void Keccak_f_32(uint32_t s[50]);
+EXPORTIT void Keccak_f_64(uint64_t s[25]);
+
diff --git a/src/gsm/tuak/tuak.c b/src/gsm/tuak/tuak.c
new file mode 100644
index 0000000..c044a37
--- /dev/null
+++ b/src/gsm/tuak/tuak.c
@@ -0,0 +1,372 @@
+/* (C) 2023 by Harald Welte <laforge@osmocom.org>
+ *
+ * All Rights Reserved
+ *
+ * SPDX-License-Identifier: GPL-2.0+
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ */
+
+#include <stdint.h>
+#include <string.h>
+#include <string.h>
+#include <errno.h>
+
+#include <osmocom/core/utils.h>
+
+#include "KeccakP-1600-3gpp.h"
+
+/* TUAK authentication algorithm
+ * as proposed by 3GPP as an alternative to Milenage
+ * algorithm based on SHA-3 (more exactly its KeccakP-1600 permutation)
+ * see 3GPP TS 35.231, 232 and 233 */
+
+static unsigned int g_keccak_iterations = 1;
+static const char algoname[] = "TUAK1.0";
+const uint8_t zero16[16] = { 0, };
+
+void tuak_set_keccak_iterations(unsigned int i)
+{
+	g_keccak_iterations = i;
+}
+
+/* append data from 'input' to 'buf' at 'idx', reversing byte order */
+#define PUSH_DATA(buf, idx, input, nbytes)	\
+	for (int i = nbytes-1; i >= 0; i--) {	\
+		buf[idx++] = input[i];		\
+	}
+
+/* like memcpy(), but reversing they order of bytes */
+void memcpy_reverse(uint8_t *dst, const uint8_t *src, size_t len)
+{
+	for (size_t i = 0; i < len; i++)
+		dst[i] = src[len-i-1];
+}
+
+static void tuak_core(uint8_t buf[200], const uint8_t *opc, uint8_t instance, const uint8_t *_rand,
+		      const uint8_t *amf, const uint8_t *sqn, const uint8_t *k, uint8_t k_len_bytes,
+		      unsigned int keccac_iterations)
+{
+	unsigned int idx = 0;
+
+	PUSH_DATA(buf, idx, opc, 32);
+	buf[idx++] = instance;
+	PUSH_DATA(buf, idx, algoname, strlen(algoname)); /* without trailing NUL */
+	PUSH_DATA(buf, idx, _rand, 16);
+	PUSH_DATA(buf, idx, amf, 2);
+	PUSH_DATA(buf, idx, sqn, 6);
+	PUSH_DATA(buf, idx, k, k_len_bytes);
+	memset(buf+idx, 0, 32-k_len_bytes); idx += 32-k_len_bytes;
+	buf[idx++] = 0x1f;
+	memset(buf+idx, 0, 38); idx += 38;
+	buf[idx++] = 0x80;
+	memset(buf+idx, 0, 64); idx += 64;
+	OSMO_ASSERT(idx == 200);
+
+	for (unsigned int i = 0; i < keccac_iterations; i++)
+		Keccak_f_64((uint64_t *) buf);
+}
+
+/**
+ * tuak_f1 - TUAK f1 algorithm
+ * @opc: OPc = 256-bit value derived from OP and K
+ * @k: K = 128-bit or 256-bit subscriber key
+ * @_rand: RAND = 128-bit random challenge
+ * @sqn: SQN = 48-bit sequence number
+ * @amf: AMF = 16-bit authentication management field
+ * @mac_a: Buffer for MAC-A = 64/128/256-bit network authentication code
+ * Returns: 0 on success, -1 on failure
+ */
+int tuak_f1(const uint8_t *opc, const uint8_t *k, uint8_t k_len_bytes, const uint8_t *_rand,
+	    const uint8_t *sqn, const uint8_t *amf, uint8_t *mac_a, uint8_t mac_a_len_bytes,
+	    unsigned int keccac_iterations)
+{
+	uint8_t buf[200];
+	uint8_t instance = 0x00;
+
+	switch (mac_a_len_bytes) {
+	case 8:
+		instance |= 0x08;
+		break;
+	case 16:
+		instance |= 0x10;
+		break;
+	case 32:
+		instance |= 0x20;
+		break;
+	default:
+		return -EINVAL;
+	}
+
+	switch (k_len_bytes) {
+	case 16:
+		break;
+	case 32:
+		instance |= 0x01;
+		break;
+	default:
+		return -EINVAL;
+	}
+
+	tuak_core(buf, opc, instance, _rand, amf, sqn, k, k_len_bytes, keccac_iterations);
+
+	memcpy_reverse(mac_a, buf, mac_a_len_bytes);
+
+	return 0;
+}
+
+/**
+ * tuak_f1star - TUAK f1* algorithm
+ * @opc: OPc = 256-bit value derived from OP and K
+ * @k: K = 128-bit or 256-bit subscriber key
+ * @_rand: RAND = 128-bit random challenge
+ * @sqn: SQN = 48-bit sequence number
+ * @amf: AMF = 16-bit authentication management field
+ * @mac_s: Buffer for MAC-S = 64/128/256-bit resync authentication code
+ * Returns: 0 on success, -1 on failure
+ */
+int tuak_f1star(const uint8_t *opc, const uint8_t *k, uint8_t k_len_bytes, const uint8_t *_rand,
+		const uint8_t *sqn, const uint8_t *amf, uint8_t *mac_s, uint8_t mac_s_len_bytes,
+		unsigned int keccac_iterations)
+{
+	uint8_t buf[200];
+	uint8_t instance = 0x80;
+
+	switch (mac_s_len_bytes) {
+	case 8:
+		instance |= 0x08;
+		break;
+	case 16:
+		instance |= 0x10;
+		break;
+	case 32:
+		instance |= 0x20;
+		break;
+	default:
+		return -EINVAL;
+	}
+
+	switch (k_len_bytes) {
+	case 16:
+		break;
+	case 32:
+		instance |= 0x01;
+		break;
+	default:
+		return -EINVAL;
+	}
+
+	tuak_core(buf, opc, instance, _rand, amf, sqn, k, k_len_bytes, keccac_iterations);
+
+	memcpy_reverse(mac_s, buf, mac_s_len_bytes);
+
+	return 0;
+}
+
+/**
+ * tuak_f2345 - TUAK f2, f3, f4, f5, algorithms
+ * @opc: OPc = 256-bit value derived from OP and K
+ * @k: K = 128/256-bit subscriber key
+ * @_rand: RAND = 128-bit random challenge
+ * @res: Buffer for RES = 32/64/128/256-bit signed response (f2), or %NULL
+ * @ck: Buffer for CK = 128/256-bit confidentiality key (f3), or %NULL
+ * @ik: Buffer for IK = 128/256-bit integrity key (f4), or %NULL
+ * @ak: Buffer for AK = 48-bit anonymity key (f5), or %NULL
+ * Returns: 0 on success, -1 on failure
+ */
+int tuak_f2345(const uint8_t *opc, const uint8_t *k, uint8_t k_len_bytes,
+	       const uint8_t *_rand, uint8_t *res, uint8_t res_len_bytes,
+	       uint8_t *ck, uint8_t ck_len_bytes,
+	       uint8_t *ik, uint8_t ik_len_bytes, uint8_t *ak, unsigned int keccac_iterations)
+{
+	uint8_t buf[200];
+	uint8_t instance = 0x40;
+
+	switch (res_len_bytes) {
+	case 4:
+		break;
+	case 8:
+		instance |= 0x08;
+		break;
+	case 16:
+		instance |= 0x10;
+		break;
+	case 32:
+		instance |= 0x20;
+		break;
+	default:
+		return -EINVAL;
+	}
+
+	switch (ck_len_bytes) {
+	case 16:
+		break;
+	case 32:
+		instance |= 0x04;
+		break;
+	default:
+		return -EINVAL;
+	}
+
+	switch (ik_len_bytes) {
+	case 16:
+		break;
+	case 32:
+		instance |= 0x02;
+		break;
+	default:
+		return -EINVAL;
+	}
+
+	switch (k_len_bytes) {
+	case 16:
+		break;
+	case 32:
+		instance |= 0x01;
+		break;
+	default:
+		return -EINVAL;
+	}
+
+	tuak_core(buf, opc, instance, _rand, zero16, zero16, k, k_len_bytes, keccac_iterations);
+
+	if (res)
+		memcpy_reverse(res, buf, res_len_bytes);
+
+	if (ck)
+		memcpy_reverse(ck, buf + 32, ck_len_bytes);
+
+	if (ik)
+		memcpy_reverse(ik, buf + 64, ik_len_bytes);
+
+	if (ak)
+		memcpy_reverse(ak, buf + 96, 6);
+
+	return 0;
+}
+
+/**
+ * tuak_f5star - TUAK f5* algorithm
+ * @opc: OPc = 256-bit value derived from OP and K
+ * @k: K = 128/256-bit subscriber key
+ * @_rand: RAND = 128-bit random challenge
+ * @ak: Buffer for AK = 48-bit anonymity key (f5)
+ * Returns: 0 on success, -1 on failure
+ */
+int tuak_f5star(const uint8_t *opc, const uint8_t *k, uint8_t k_len_bytes,
+		const uint8_t *_rand, uint8_t *ak, unsigned int keccac_iterations)
+{
+	uint8_t buf[200];
+	uint8_t instance = 0xc0;
+
+	switch (k_len_bytes) {
+	case 16:
+		break;
+	case 32:
+		instance += 1;
+		break;
+	default:
+		return -EINVAL;
+	}
+
+	tuak_core(buf, opc, instance, _rand, zero16, zero16, k, k_len_bytes, keccac_iterations);
+
+	memcpy_reverse(ak, buf + 96, 6);
+
+	return 0;
+}
+
+/**
+ * tuak_generate - Generate AKA AUTN,IK,CK,RES
+ * @opc: OPc = 256-bit operator variant algorithm configuration field (encr.)
+ * @amf: AMF = 16-bit authentication management field
+ * @k: K = 128/256-bit subscriber key
+ * @sqn: SQN = 48-bit sequence number
+ * @_rand: RAND = 128-bit random challenge
+ * @autn: Buffer for AUTN = 128-bit authentication token
+ * @ik: Buffer for IK = 128/256-bit integrity key (f4), or %NULL
+ * @ck: Buffer for CK = 128/256-bit confidentiality key (f3), or %NULL
+ * @res: Buffer for RES = 32/64/128-bit signed response (f2), or %NULL
+ * @res_len: Max length for res; set to used length or 0 on failure
+ */
+void tuak_generate(const uint8_t *opc, const uint8_t *amf, const uint8_t *k, uint8_t k_len_bytes,
+		   const uint8_t *sqn, const uint8_t *_rand, uint8_t *autn, uint8_t *ik,
+		   uint8_t *ck, uint8_t *res, size_t *res_len)
+{
+	int i;
+	uint8_t mac_a[8], ak[6];
+
+	if (*res_len < 4) {
+		*res_len = 0;
+		return;
+	}
+	if (tuak_f1(opc, k, k_len_bytes, _rand, sqn, amf, mac_a, sizeof(mac_a), g_keccak_iterations) ||
+	    tuak_f2345(opc, k, k_len_bytes, _rand, res, *res_len, ck, 16, ik, 16, ak, g_keccak_iterations)) {
+		*res_len = 0;
+		return;
+	}
+
+	/* AUTN = (SQN ^ AK) || AMF || MAC */
+	for (i = 0; i < 6; i++)
+		autn[i] = sqn[i] ^ ak[i];
+	memcpy(autn + 6, amf, 2);
+	memcpy(autn + 8, mac_a, 8);
+}
+
+
+/**
+ * tuak_auts - Milenage AUTS validation
+ * @opc: OPc = 256-bit operator variant algorithm configuration field (encr.)
+ * @k: K = 128/256-bit subscriber key
+ * @_rand: RAND = 128-bit random challenge
+ * @auts: AUTS = 112-bit authentication token from client
+ * @sqn: Buffer for SQN = 48-bit sequence number
+ * Returns: 0 = success (sqn filled), -1 on failure
+ */
+int tuak_auts(const uint8_t *opc, const uint8_t *k, uint8_t k_len_bytes,
+	      const uint8_t *_rand, const uint8_t *auts, uint8_t *sqn)
+{
+	uint8_t amf[2] = { 0x00, 0x00 }; /* TS 33.102 v7.0.0, 6.3.3 */
+	uint8_t ak[6], mac_s[8];
+	int i;
+
+	if (tuak_f5star(opc, k, k_len_bytes, _rand, ak, g_keccak_iterations))
+		return -1;
+	for (i = 0; i < 6; i++)
+		sqn[i] = auts[i] ^ ak[i];
+	if (tuak_f1star(opc, k, k_len_bytes, _rand, sqn, amf, mac_s, 8, g_keccak_iterations) ||
+	    memcmp(mac_s, auts + 6, 8) != 0)
+		return -1;
+	return 0;
+}
+
+int tuak_opc_gen(uint8_t *opc, const uint8_t *k, uint8_t k_len_bytes, const uint8_t *op)
+{
+	uint8_t buf[200];
+	uint8_t instance;
+
+	switch (k_len_bytes) {
+	case 16:
+		instance = 0x00;
+		break;
+	case 32:
+		instance = 0x01;
+		break;
+	default:
+		return -EINVAL;
+	}
+
+	tuak_core(buf, op, instance, zero16, zero16, zero16, k, k_len_bytes, g_keccak_iterations);
+
+	memcpy_reverse(opc, buf, 32);
+
+	return 0;
+}
diff --git a/src/gsm/tuak/tuak.h b/src/gsm/tuak/tuak.h
new file mode 100644
index 0000000..1a80822
--- /dev/null
+++ b/src/gsm/tuak/tuak.h
@@ -0,0 +1,33 @@
+#pragma once
+#include <stdint.h>
+
+/* low-level functions */
+
+int tuak_f1(const uint8_t *opc, const uint8_t *k, uint8_t k_len_bytes, const uint8_t *_rand,
+	    const uint8_t *sqn, const uint8_t *amf, uint8_t *mac_a, uint8_t mac_a_len_bytes,
+	    unsigned int keccac_iterations);
+
+int tuak_f1star(const uint8_t *opc, const uint8_t *k, uint8_t k_len_bytes, const uint8_t *_rand,
+		const uint8_t *sqn, const uint8_t *amf, uint8_t *mac_s, uint8_t mac_s_len_bytes,
+		unsigned int keccac_iterations);
+
+int tuak_f2345(const uint8_t *opc, const uint8_t *k, uint8_t k_len_bytes,
+	       const uint8_t *_rand, uint8_t *res, uint8_t res_len_bytes,
+	       uint8_t *ck, uint8_t ck_len_bytes,
+	       uint8_t *ik, uint8_t ik_len_bytes, uint8_t *ak, unsigned int keccac_iterations);
+
+int tuak_f5star(const uint8_t *opc, const uint8_t *k, uint8_t k_len_bytes,
+		const uint8_t *_rand, uint8_t *ak, unsigned int keccac_iterations);
+
+/* high-level API */
+
+void tuak_set_keccak_iterations(unsigned int i);
+
+void tuak_generate(const uint8_t *opc, const uint8_t *amf, const uint8_t *k, uint8_t k_len_bytes,
+		   const uint8_t *sqn, const uint8_t *_rand, uint8_t *autn, uint8_t *ik,
+		   uint8_t *ck, uint8_t *res, size_t *res_len);
+
+int tuak_auts(const uint8_t *opc, const uint8_t *k, uint8_t k_len_bytes,
+	      const uint8_t *_rand, const uint8_t *auts, uint8_t *sqn);
+
+int tuak_opc_gen(uint8_t *opc, const uint8_t *k, uint8_t k_len_bytes, const uint8_t *op);