Merge pull request #2 from pharpend/perf_faster

perf: add kek_fast - with unrolled uint64 implementation 800x faster

.gitignore

@@ -0,0 +1,2 @@
+*.swp
+*.beam

LICENSE

@@ -0,0 +1,15 @@
+ISC License
+
+Copyright (c) 2022 Peter Harpending
+
+Permission to use, copy, modify, and/or distribute this software for any
+purpose with or without fee is hereby granted, provided that the above
+copyright notice and this permission notice appear in all copies.
+
+THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
+REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY
+AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
+INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
+LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR
+OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
+PERFORMANCE OF THIS SOFTWARE.

README.md

@@ -0,0 +1,3 @@
+# kek
+
+Work-in-progress Erlang implementation of keccak

kek.erl

@@ -1,47 +1,157 @@
 %% @doc
 %% References
 %% 1. Helpful lecture: https://www.youtube.com/watch?v=JWskjzgiIa4
+%%    Notes: https://www.crypto-textbook.com/download/Understanding-Cryptography-Keccak.pdf
 %% 2. NIST standard: https://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.202.pdf
 %%    (btw: the double bar notation means "concatenate")
+%% 3. https://en.wikipedia.org/wiki/SHA-3
 %% @end
 -module(kek).
 
 % theta and rho steps are done
-% weird syntax errors
 
 -compile(export_all).
 
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%% TOP LEVEL API
+%%%
+%%% sha*s and shake*s
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+-spec sha3_224(Message) -> Digest
+    when Message :: bitstring(),
+         Digest  :: <<_:224>>.
+%% @doc
+%% SHA-3 with an output bit length of 224 bits.
+%% @end
+
 sha3_224(Message) ->
     sha3(224, Message).
 
+
+
+-spec sha3_256(Message) -> Digest
+    when Message :: bitstring(),
+         Digest  :: <<_:256>>.
+%% @doc
+%% SHA-3 with an output bit length of 256 bits.
+%% @end
+
 sha3_256(Message) ->
     sha3(256, Message).
 
+
+
+-spec sha3_384(Message) -> Digest
+    when Message :: bitstring(),
+         Digest  :: <<_:384>>.
+%% @doc
+%% SHA-3 with an output bit length of 384 bits.
+%% @end
+
 sha3_384(Message) ->
     sha3(384, Message).
 
+
+
+-spec sha3_512(Message) -> Digest
+    when Message :: bitstring(),
+         Digest  :: <<_:512>>.
+%% @doc
+%% SHA-3 with an output bit length of 512 bits.
+%% @end
+
 sha3_512(Message) ->
     sha3(512, Message).
 
+
+
+-spec sha3(OutputBitLength, Message) -> Digest
+    when OutputBitLength :: pos_integer(),
+         Message         :: bitstring(),
+         Digest          :: bitstring().
+%% @doc
+%% SHA-3 with an arbitrary output bit length.
+%%
+%% This means Keccak with Capacity = 2*OutputBitLength. Additionally, SHA3
+%% concatenates the bits 01 onto the end of the input, before sending the
+%% Message to keccak/3.
+%% @end
+
 sha3(OutputBitLength, Message) ->
     Capacity = 2*OutputBitLength,
     ShaMessage = <<Message/bitstring, (2#01):2>>,
     keccak(Capacity, ShaMessage, OutputBitLength).
 
 
+
+-spec shake128(Message, OutputBitLength) -> Digest
+    when Message         :: bitstring(),
+         OutputBitLength :: pos_integer(),
+         Digest          :: bitstring().
+%% @doc
+%% This is the SHAKE variable-length hash with Capacity 256 = 2*128 bits.
+%% @end
+
 shake128(Message, OutputBitLength) ->
     shake(128, Message, OutputBitLength).
 
+
+
+-spec shake256(Message, OutputBitLength) -> Digest
+    when Message         :: bitstring(),
+         OutputBitLength :: pos_integer(),
+         Digest          :: bitstring().
+%% @doc
+%% This is the SHAKE variable-length hash with Capacity 512 = 2*256 bits.
+%% @end
+
 shake256(Message, OutputBitLength) ->
     shake(256, Message, OutputBitLength).
 
+
+
+-spec shake(ShakeNumber, Message, OutputBitLength) -> Digest
+    when ShakeNumber     :: pos_integer(),
+         Message         :: bitstring(),
+         OutputBitLength :: pos_integer(),
+         Digest          :: bitstring().
+%% @doc
+%% This is the SHAKE variable-length hash with Capacity 512 = 2*ShakeNumber bits.
+%%
+%% This concatenates the bitstring 1111 onto the end of the Message before
+%% sending the message to keccak/3.
+%% @end
+
 shake(ShakeNumber, Message, OutputBitLength) ->
     Capacity = 2*ShakeNumber,
-    ShakeMessage = <<Message, (2#1111):4>>,
+    ShakeMessage = <<Message/bitstring, (2#1111):4>>,
     keccak(Capacity, ShakeMessage, OutputBitLength).
 
 
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%% OUTER KECCAK
+%%%
+%%% Keccak pads the input, absorbs it into the sponge, and squeezes the bits out
+%%% of the sponge.  The absorption and squeezing phases invoke "inner keccak",
+%%% which is the heart of the algorithm.
+%%%
+%%% - keccak/3
+%%% - pad/2
+%%% - absorb/4
+%%% - squeeze/3
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 -spec keccak(Capacity, Message, OutputBitLength) -> Digest
     when Capacity        :: pos_integer(),
          Message         :: bitstring(),
@@ -152,6 +262,7 @@ absorb(<<>>, _r, _c, FinalSponge) ->
     FinalSponge.
 
 
+
 -spec squeeze(WetSponge, OutputBitLength, BitRate) -> ResultBits
     when WetSponge       :: <<_:1600>>,
          OutputBitLength :: pos_integer(),
@@ -190,6 +301,26 @@ really_squeeze(WetSponge, OutputBitLength, BitRate, ResultAcc)->
 
 
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%% THE DREADED INNER KECCAK
+%%%
+%%% This is the "f" function that appears in all the documentation.
+%%%
+%%% The input is the 1600-bit sponge array. inner_keccak/1 sends the input
+%%% through 24 "rounds". Each round consists of the 5 Greek letter steps, each of
+%%% which is a weird transformation on the array.
+%%%
+%%% In "inner keccak", the input array is thought of as a 5x5x64 3D array. The
+%%% coordinate system is described in its own section.
+%%%
+%%% TODO: iota depends on the round index, so this code may need to be altered
+%%% slightly.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 -spec inner_keccak(Sponge) -> NewSponge
     when Sponge    :: <<_:1600>>,
          NewSponge :: <<_:1600>>.
@@ -204,33 +335,55 @@ inner_keccak(Sponge) ->
 
 -spec rounds(Sponge, NumRoundsLeft) -> ResultSponge
     when Sponge        :: <<_:1600>>,
-         NumRoundsLeft :: pos_integer(),
+         NumRoundsLeft :: non_neg_integer(),
          ResultSponge  :: <<_:1600>>.
 %% @private
 %% do however many rounds
 %% @end
-rounds(Sponge, NumRoundsLeft) when 1 =< NumRoundsLeft ->
+
-    NewSponge        = rnd(Sponge),
-    NewNumRoundsLeft = NumRoundsLeft - 1,
-    rounds(NewSponge, NewNumRoundsLeft);
 % no rounds left
 rounds(FinalSponge, 0) ->
-    FinalSponge.
+    FinalSponge;
+rounds(Sponge, NumRoundsLeft) ->
+    % NRoundsLeft = 24
+    % idx0 = 0
+    % NRoundsLeft = 1
+    % idx0 = 23
+    RoundIdx0        = 24 - NumRoundsLeft,
+    NewSponge        = rnd(RoundIdx0, Sponge),
+    NewNumRoundsLeft = NumRoundsLeft - 1,
+    rounds(NewSponge, NewNumRoundsLeft).
 
 
 
--spec rnd(Sponge) -> NewSponge
+-spec rnd(RoundIdx0, Sponge) -> NewSponge
-    when Sponge    :: <<_:1600>>,
+    when RoundIdx0 :: 0..23,
+         Sponge    :: <<_:1600>>,
          NewSponge :: <<_:1600>>.
 %% @private
 %% do a single round
 %% @private
 
-rnd(Sponge) ->
+rnd(RoundIdx0, Sponge) ->
-    iota(chi(pi(rho(theta(Sponge))))).
+    iota(RoundIdx0, chi(pi(rho(theta(Sponge))))).
 
 
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% THETA STEP
+%%
+%% For each bit in the input array,
+%% 1. take
+%%    - the bit
+%%    - the 5-bit column to the left
+%%    - the 5-bit column to the front right
+%% 2. compute the parity of their concatenation (0 if even# of 1s, 1 if odd# of
+%%    1s)
+%% 3. set the bit to that parity value
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 -spec theta(Array) -> NewArray
     when Array    :: <<_:1600>>,
          NewArray :: <<_:1600>>.
@@ -240,12 +393,12 @@ rnd(Sponge) ->
 %% @end
 
 theta(Array) ->
-    theta(Array, 0).
+    theta(Array, Array, 0).
 
 
-
+-spec theta(Array, OldArray, Idx0) -> NewArray
--spec theta(Array, Idx0) -> NewArray
     when Array    :: <<_:1600>>,
+         OldArray :: <<_:1600>>,
          Idx0     :: 0..1599,
          NewArray :: <<_:1600>>.
 %% @private
@@ -254,7 +407,7 @@ theta(Array) ->
 %% @end
 
 % done
-theta(ResultArray, 1600) ->
+theta(ResultArray, _, 1600) ->
     ResultArray;
 % do the weird permutation
 % x = left/right             -/+
@@ -262,20 +415,46 @@ theta(ResultArray, 1600) ->
 % z = outOfScreen/intoScreen -/+
 %     front/behind           -/+
 % left-handed coordinate system but what can you do
-theta(ArrayBits, ThisIdx0) ->
+theta(ArrayBits, OrigArray, ThisIdx0) ->
     <<Before:ThisIdx0, ThisBit:1, Rest/bitstring>> = ArrayBits,
     {xyz, ThisX, _ThisY, ThisZ} = idx0_to_xyz(ThisIdx0),
     XToTheLeft                  = left(ThisX),
     XToTheRight                 = right(ThisX),
     ZToTheFront                 = front(ThisZ),
-    ColumnToTheLeft             = xzth({xz, XToTheLeft, ThisZ}, ArrayBits),
+    ColumnToTheLeft             = xzth({xz, XToTheLeft, ThisZ}, OrigArray),
-    ColumnToTheFrontRight       = xzth({xz, XToTheRight, ZToTheFront}, ArrayBits),
+    ColumnToTheFrontRight       = xzth({xz, XToTheRight, ZToTheFront}, OrigArray),
     NewBit                      = parity(<<ColumnToTheLeft/bitstring, ColumnToTheFrontRight/bitstring, ThisBit:1>>),
     NewBits                     = <<Before:ThisIdx0, NewBit:1, Rest/bitstring>>,
     NewIdx0                     = ThisIdx0 + 1,
-    theta(NewBits, NewIdx0).
+    theta(NewBits, OrigArray, NewIdx0).
 
 
+-spec parity(Bits) -> Parity
+    when Bits   :: bitstring(),
+         Parity :: 0 | 1.
+%% @private
+%% Count the number of 1s in the given bitstring. Return 0 if even, 1 if odd.
+%% @end
+
+parity(Bits) ->
+    parity(Bits, 0).
+
+parity(<<0:1, Rest/bitstring>>, NOnes) -> parity(Rest, NOnes);
+parity(<<1:1, Rest/bitstring>>, NOnes) -> parity(Rest, NOnes + 1);
+parity(<<>>                   , NOnes) -> NOnes rem 2.
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% RHO STEP
+%%
+%% This step applies an affine shift to each 64-bit "lane" (fixed X,Y; Z ranges
+%% from 0 to 63).
+%%
+%% The amount of the shift is given by the "offset" table.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 -spec rho(Array) -> NewArray
     when    Array :: <<_:1600>>,
@@ -314,10 +493,10 @@ rho(Array, XY = {xy, X, Y}) ->
 
 
 
--spec rhoxy(Array, LaneXY) -> NewArray
+-spec rhoxy(Array1600, LaneXY) -> NewArray1600
-    when Array    :: <<_:1600>>,
+    when Array1600    :: <<_:1600>>,
          LaneXY       :: {xy, 0..4, 0..4},
-         NewArray :: <<_:1600>>.
+         NewArray1600 :: <<_:1600>>.
 %% @private
 %% do the rho step to a given lane
 %% @end
@@ -334,38 +513,9 @@ rhoxy(Array, ThisXY = {xy, ThisX, ThisY}) ->
     % put them on the front
     <<Foo:(64 - ThisOffset), Bar:ThisOffset>> = ThisLane,
     NewLane  = <<Bar:ThisOffset, Foo:(64 - ThisOffset)>>,
-    NewArray = replace_lane(Array, ThisXY, NewLane),
+    NewArray = xyset(ThisXY, Array, NewLane),
     NewArray.
 
-xyth(_, _) -> error(nyi).
-
--spec replace_lane(OriginalArray, LaneXY, NewLane) -> NewArray
-    when OriginalArray :: <<_:1600>>,
-         LaneXY        :: {xy, 0..4, 0..4},
-         NewLane       :: <<_:64>>,
-         NewArray      :: <<_:1600>>.
-%% @private
-%% take the original array, and swap out the lane at the given x,y coordinate
-%% with the new given lane. the lane will be represented continuously so we
-%% can do a hack
-%% @end
-
-% special case when it's the last lane
-% grab the final 64 bits off the original array and replace them with the new lane
-replace_lane(<<Pre:(1600 - 64), _:64>>, _LaneXY = {xy, 4, 4}, NewLane) ->
-    <<Pre:(1600 - 64), NewLane/bitstring>>;
-% general case, grab the shit before the lane, grab the shit after the lane
-% replace the shit in the middle
-replace_lane(OriginalArray, _LaneXY = {xy, LaneX, LaneY}, NewLane) ->
-    FirstBitOfLane_Idx0    = xyz_to_idx0({xyz, LaneX, LaneY, 0}),
-    FirstBitAfterLane_Idx0 = xyz_to_idx0({xyz, LaneX, LaneY, 63}) + 1,
-    NumberOfBitsBeforeTheLane    = FirstBitOfLane_Idx0,
-    NumberOfBitsIncludingTheLane = FirstBitAfterLane_Idx0,
-    <<PreLane:NumberOfBitsBeforeTheLane   ,         _/bitstring>> = OriginalArray,
-    <<      _:NumberOfBitsIncludingTheLane, AfterLane/bitstring>> = OriginalArray,
-    Result = <<PreLane:NumberOfBitsBeforeTheLane, NewLane/bitstring, AfterLane/bitstring>>,
-    Result.
-
 
 
 -spec offset(X, Y) -> Offset
@@ -409,17 +559,367 @@ offset(2, 3) ->  15 rem 64.
 
 
 
-pi(_Sponge) ->
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    error(nyi).
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% PI STEP
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-chi(_Sponge) ->
+-spec pi(Array1600) -> NewArray1600
-    error(nyi).
+    when Array1600    :: <<_:1600>>,
+         NewArray1600 :: <<_:1600>>.
+%% @private
+%% The effect of this step is to rearrange the lanes
+%%
+%% Result[X, Y] = Input[X + 3*Y, X]
+%%
+%% (mod 5 of course)
+%% @end
 
-iota(_Sponge) ->
+pi(Array1600) ->
-    error(nyi).
+    % what I'm going to make is a map #{{xy, X, Y} := Lane}
+    % then make a new map from the which applies the coordinate transformation
+    % then convert it back into an array
+    OriginalLaneMap = lane_map(Array1600, #{}, {xy, 0, 0}),
+    NewLaneMap      = new_lane_map(OriginalLaneMap, #{}, {xy, 0, 0}),
+    NewArray1600    = lane_map_to_arr1600(NewLaneMap, <<0:1600>>, {xy, 0, 0}),
+    NewArray1600.
 
 
 
+-spec lane_map(Array1600, MapAcc, Coord) -> LaneMap
+    when Array1600 :: <<_:1600>>,
+         MapAcc    :: #{XY := Lane},
+         Coord     :: XY,
+         LaneMap   :: #{XY := Lane},
+         XY        :: {xy, X :: 0..4, Y :: 0..4},
+         Lane      :: <<_:64>>.
+%% @private
+%% Make a map #{XY := Lane}
+%% @end
+
+% terminal case, end of array
+lane_map(Array1600, MapAcc, ThisXY = {xy, 4, 4}) ->
+    ThisLane = xyth(ThisXY, Array1600),
+    FinalMap = MapAcc#{ThisXY => ThisLane},
+    FinalMap;
+% end of Y value, set Y to 0 and increment X
+lane_map(Array1600, MapAcc, ThisXY = {xy, X, 4}) ->
+    ThisLane  = xyth(ThisXY, Array1600),
+    NewMapAcc = MapAcc#{ThisXY => ThisLane},
+    NewXY     = {xy, X + 1, 0},
+    lane_map(Array1600, NewMapAcc, NewXY);
+% general case: increment Y value
+lane_map(Array1600, MapAcc, ThisXY = {xy, X, Y}) ->
+    ThisLane  = xyth(ThisXY, Array1600),
+    NewMapAcc = MapAcc#{ThisXY => ThisLane},
+    NewXY     = {xy, X, Y + 1},
+    lane_map(Array1600, NewMapAcc, NewXY).
+
+
+
+-spec new_lane_map(LaneMap, MapAcc, Coord) -> NewLaneMap
+    when LaneMap    :: #{XY := Lane},
+         MapAcc     :: LaneMap,
+         Coord      :: XY,
+         NewLaneMap :: LaneMap,
+         XY         :: {xy, X :: 0..4, Y :: 0..4},
+         Lane       :: <<_:64>>.
+%% @private
+%% The effect of this step is to rearrange the lanes
+%%
+%% Result[X, Y] = Input[X + 3*Y, X]
+%%
+%% (mod 5 of course)
+%% @end
+
+% terminal case, end of array
+new_lane_map(OrigLaneMap, MapAcc, ThisXY = {xy, 4, 4}) ->
+    OrigXY   = xytrans(ThisXY),
+    ThisLane = maps:get(OrigXY, OrigLaneMap),
+    FinalMap = MapAcc#{ThisXY => ThisLane},
+    FinalMap;
+% end of Y value, set Y to 0 and increment X
+new_lane_map(OrigLaneMap, MapAcc, ThisXY = {xy, X, 4}) ->
+    OrigXY    = xytrans(ThisXY),
+    ThisLane  = maps:get(OrigXY, OrigLaneMap),
+    NewMapAcc = MapAcc#{ThisXY => ThisLane},
+    NewXY     = {xy, X + 1, 0},
+    new_lane_map(OrigLaneMap, NewMapAcc, NewXY);
+% general case: increment Y value
+new_lane_map(OrigLaneMap, MapAcc, ThisXY = {xy, X, Y}) ->
+    OrigXY    = xytrans(ThisXY),
+    ThisLane  = maps:get(OrigXY, OrigLaneMap),
+    NewMapAcc = MapAcc#{ThisXY => ThisLane},
+    NewXY     = {xy, X, Y + 1},
+    new_lane_map(OrigLaneMap, NewMapAcc, NewXY).
+
+
+
+-spec xytrans(ResultXY) -> InputXY
+    when ResultXY :: XY,
+         InputXY  :: XY,
+         XY       :: {xy, X :: 0..4, Y :: 0..4}.
+%% @private
+%% Result[X, Y] = Input[X + 3*Y, X]
+%%
+%% See NIST doc, pp. 14
+
+xytrans({xy, X, Y}) ->
+    {xy, (X + 3*Y) rem 5, X}.
+
+
+
+-spec lane_map_to_arr1600(LaneMap, Array1600Acc, Coord) -> Array1600
+    when LaneMap      :: #{XY := Lane},
+         Array1600Acc :: Array1600,
+         Coord        :: XY,
+         Array1600    :: <<_:1600>>,
+         XY           :: {xy, X :: 0..4, Y :: 0..4},
+         Lane         :: <<_:64>>.
+%% @private
+%% inverse of lane_map/3
+%%
+%% it would probably faster to concatenate an accumulator, but that requires
+%% iterating in the correct order, and i'm more comfortable calling xyset/3
+%% @end
+
+% terminal case, end of array
+lane_map_to_arr1600(LaneMap, Array1600Acc, ThisXY = {xy, 4, 4}) ->
+    ThisLane          = maps:get(ThisXY, LaneMap),
+    FinalArray1600Acc = xyset(ThisXY, Array1600Acc, ThisLane),
+    FinalArray1600Acc;
+% end of Y value, set Y to 0 and increment X
+lane_map_to_arr1600(LaneMap, Array1600Acc, ThisXY = {xy, X, 4}) ->
+    ThisLane        = maps:get(ThisXY, LaneMap),
+    NewArray1600Acc = xyset(ThisXY, Array1600Acc, ThisLane),
+    NewXY           = {xy, X + 1, 0},
+    lane_map_to_arr1600(LaneMap, NewArray1600Acc, NewXY);
+% general case: increment Y value
+lane_map_to_arr1600(LaneMap, Array1600Acc, ThisXY = {xy, X, Y}) ->
+    ThisLane        = maps:get(ThisXY, LaneMap),
+    NewArray1600Acc = xyset(ThisXY, Array1600Acc, ThisLane),
+    NewXY           = {xy, X, Y + 1},
+    lane_map_to_arr1600(LaneMap, NewArray1600Acc, NewXY).
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% CHI STEP
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+-spec chi(Array1600) -> NewArray1600
+    when Array1600    :: <<_:1600>>,
+         NewArray1600 :: <<_:1600>>.
+%% @private
+%% The chi step. The following transformation is applied to each bit
+%%
+%% NewBit = lxor(Bit,
+%%               land(lnot(BitToTheRight),
+%%                    Bit2ToTheRight))
+
+chi(Array1600) ->
+    chi(Array1600, Array1600, 0).
+
+
+-spec chi(Array1600, OldArray1600, Idx0) -> NewArray1600
+    when Array1600    :: <<_:1600>>,
+         OldArray1600 :: <<_:1600>>,
+         Idx0         :: non_neg_integer(),
+         NewArray1600 :: <<_:1600>>.
+%% @private
+%% The chi step. The following transformation is applied to each bit
+%%
+%% NewBit = lxor(Bit,
+%%               land(lnot(BitToTheRight),
+%%                    Bit2ToTheRight))
+%%
+%% FIXME: Could be made more efficient by operating on lanes
+
+chi(Array1600, OrigArray, ThisIdx0) when 0 =< ThisIdx0, ThisIdx0 =< 1599 ->
+    ThisXYZ      = {xyz,             ThisX  , ThisY, ThisZ} = idx0_to_xyz(ThisIdx0),
+    RightXYZ     = {xyz,       right(ThisX) , ThisY, ThisZ},
+    Right2XYZ    = {xyz, right(right(ThisX)), ThisY, ThisZ},
+    ThisBit      = xyzth(ThisXYZ  , Array1600),
+    RightBit     = xyzth(RightXYZ , OrigArray),
+    Right2Bit    = xyzth(Right2XYZ, OrigArray),
+    NewBit       = lxor(ThisBit,
+                        land(lnot(RightBit),
+                             Right2Bit)),
+    NewArray1600 = xyzset(ThisXYZ, Array1600, NewBit),
+    NewIdx0      = ThisIdx0 + 1,
+    chi(NewArray1600, OrigArray, NewIdx0);
+% terminal case
+chi(Array1600, _, 1600) ->
+    Array1600.
+
+
+lxor(0, 0) -> 0;
+lxor(0, 1) -> 1;
+lxor(1, 0) -> 1;
+lxor(1, 1) -> 0.
+
+land(0, 0) -> 0;
+land(0, 1) -> 0;
+land(1, 0) -> 0;
+land(1, 1) -> 1.
+
+lnot(0) -> 1;
+lnot(1) -> 0.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% IOTA STEP
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+-spec iota(RoundIdx0, Array1600) -> NewArray1600
+    when RoundIdx0    :: 0..23,
+         Array1600    :: <<_:1600>>,
+         NewArray1600 :: <<_:1600>>.
+%% @private
+%% iota xors the 0,0 lane by a round constant which depends on the round
+
+iota(RoundIdx0, Array1600) ->
+    <<Lane00_int:64>>     = xyth({xy, 0, 0}, Array1600),
+    ThisRoundConstant_int = round_constant_int(RoundIdx0),
+    NewLane00_bytes       = <<(Lane00_int bxor ThisRoundConstant_int):64>>,
+    xyset({xy, 0, 0}, Array1600, NewLane00_bytes).
+
+%% the algorithm is incorrect in both of these cases
+%% Well... They have a weird way of defining their "bitstrings" it is least significant bit and byte first...
+%% So just flipping them would be nice.
+
+%% textbook answer
+%% round_constant_int( 0) -> 16#0000000000000001;
+%% round_constant_int( 1) -> 16#0000000000008082;
+%% round_constant_int( 2) -> 16#800000000000808A;
+%% round_constant_int( 3) -> 16#8000000080008000;
+%% round_constant_int( 4) -> 16#000000000000808B;
+%% round_constant_int( 5) -> 16#0000000080000001;
+%% round_constant_int( 6) -> 16#8000000080008081;
+%% round_constant_int( 7) -> 16#8000000000008009;
+%% round_constant_int( 8) -> 16#000000000000008A;
+%% round_constant_int( 9) -> 16#0000000000000088;
+%% round_constant_int(10) -> 16#0000000080008009;
+%% round_constant_int(11) -> 16#000000008000000A;
+%% round_constant_int(12) -> 16#000000008000808B;
+%% round_constant_int(13) -> 16#800000000000008B;
+%% round_constant_int(14) -> 16#8000000000008089;
+%% round_constant_int(15) -> 16#8000000000008003;
+%% round_constant_int(16) -> 16#8000000000008002;
+%% round_constant_int(17) -> 16#8000000000000080;
+%% round_constant_int(18) -> 16#000000000000800A;
+%% round_constant_int(19) -> 16#800000008000000A;
+%% round_constant_int(20) -> 16#8000000080008081;
+%% round_constant_int(21) -> 16#8000000000008080;
+%% round_constant_int(22) -> 16#0000000080000001;
+%% round_constant_int(23) -> 16#8000000080008008.
+
+%% my answer
+%% TODO: In hex?
+round_constant_int( 0) -> 9223372036854775808; %% 16#8000000000000000 !!
+round_constant_int( 1) -> 4684025087442026496;
+round_constant_int( 2) -> 5836946592048873473;
+round_constant_int( 3) -> 281479271677953;
+round_constant_int( 4) -> 15060318628903649280;
+round_constant_int( 5) -> 9223372041149743104;
+round_constant_int( 6) -> 9295711110164381697;
+round_constant_int( 7) -> 10376575016438333441;
+round_constant_int( 8) -> 5836665117072162816;
+round_constant_int( 9) -> 1224979098644774912;
+round_constant_int(10) -> 10376575020733300736;
+round_constant_int(11) -> 5764607527329202176;
+round_constant_int(12) -> 15060318633198616576;
+round_constant_int(13) -> 15060037153926938625;
+round_constant_int(14) -> 10448632610476261377;
+round_constant_int(15) -> 13835339530258874369;
+round_constant_int(16) -> 4611967493404098561;
+round_constant_int(17) -> 72057594037927937;
+round_constant_int(18) -> 5764888998010945536;
+round_constant_int(19) -> 5764607527329202177;
+round_constant_int(20) -> 9295711110164381697;
+round_constant_int(21) -> 72339069014638593;
+round_constant_int(22) -> 9223372041149743104;
+round_constant_int(23) -> 1153202983878524929.
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%% INNER KECCAK COORDINATE SYSTEM
+%%%
+%%% Inner Keccak thinks of the 1600-bit input array as a 5x5x64 3D array
+%%% (5*5*64).  This section provides a variety of helper functions to talk about
+%%% the array using the X,Y,Z coordinate system.
+%%%
+%%% The coordinate system is toroidal, meaning that each coordinate is "modded
+%%% down" to be in the approprate range.  For instance, the X-coordinate "to the
+%%% right" of X=4 is X=0. And likewise, the coordinate "behind" Z=63 is Z=0. See
+%%% the section on directionality conventions.
+%%%
+%%% VOCABULARY:
+%%%
+%%% 3D state:
+%%% - The [state] is the entire 5x5x24 array
+%%%
+%%% 0D subsets of the state:
+%%% - A [bit] is a single bit in the array given by an X,Y,Z coordinate triple
+%%%   (see xyzth/3).
+%%%
+%%% 1D subsets of the state:
+%%% - a [row]
+%%%   - is a 5-bit array
+%%%   - given by a Y,Z coordinate pair in range {0..4, 0..63} (see yzth/2)
+%%%   - you should think of a row as being internally indexed with an X
+%%%     coordinate ranging in 0..4
+%%% - a [column]
+%%%   - is a 5-bit array
+%%%   - given by an X,Z coordinate pair in range {0..4, 0..63} (see xzth/2)
+%%%   - you should think of a column as being internally indexed with a Y
+%%%     coordinate ranging in 0..4
+%%% - a [lane]
+%%%   - is a 64-bit array
+%%%   - given by an X,Y coordinate pair in range {0..4, 0..4} (see xyth/2)
+%%%   - you should think of a lane as being internally indexed with a Z
+%%%     coordinate ranging in 0..63.
+%%%
+%%% 2D subsets of the state:
+%%% - a [sheet]
+%%%   - is a 5x64 array
+%%%   - given by a single X coordinate ranging in 0..4 (see xth/2)
+%%%   - you should think of a sheet as being internally indexed by a Y,Z
+%%%     coordinate pair ranging in {0..4, 0..63}.
+%%% - a [plane]
+%%%   - is a 5x64 array
+%%%   - given by a single Y coordinate ranging in 0..4 (see yth/2)
+%%%   - you should think of a plane as being internally indexed by a X,Z
+%%%     coordinate pair ranging in {0..4, 0..63}.
+%%% - a [slice]
+%%%   - is a 5x5 array
+%%%   - given by a single Z coordinate ranging in 0..63 (see zth/2)
+%%%   - you should think of a sheet as being internally indexed by a X,Y
+%%%     coordinate pair ranging in {0..4, 0..4}.
+%%%
+%%% PACKING CONVENTION:
+%%%
+%%% Each lane (64-bit long bitstring given by an {X,Y} <- {0..4, 0..4} coordinate
+%%% pair and indexed by a Z <- 0..63 coordinate.
+%%%
+%%% DIRECTIONALITY CONVENTION:
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% CONVERTING BETWEEN XYZ-INDICES AND 0-INDICES
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 -spec idx0_to_xyz(Idx0) -> XYZ
     when Idx0 :: 0..1599,
          XYZ  :: {xyz, X :: 0..4, Y :: 0..4, Z :: 0..63}.
@@ -459,6 +959,13 @@ xyz_to_idx0({xyz, X, Y, Z}) ->
 
 
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% DIRECTIONAL TRANSFORMATIONS ON SINGLE COORDINATE VALUES
+%%
+%% For instance, if you have an X-value and want to get the X-value "to the
+%% left", this section contains functions that compute such things.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 -spec left(X) -> XToTheLeft
     when X          :: 0..4,
          XToTheLeft :: 0..4.
@@ -549,12 +1056,60 @@ behind(63)                     -> 0.
 
 
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% 0D BIT ACCESSORS
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+-spec xyzth(XYZ, Array1600) -> Bit
+    when XYZ       :: {xyz, X, Y, Z},
+         Array1600 :: <<_:1600>>,
+         Bit       :: 0 | 1,
+         X         :: 0..4,
+         Y         :: 0..4,
+         Z         :: 0..63.
+%% @private
+%% Fetch the bit at the given X, Y, Z coordinate triple
+%% @end
+
+xyzth(XYZ, Array1600) ->
+    Idx0 = xyz_to_idx0(XYZ),
+    <<_Skip:Idx0, Bit:1, _Rest/bitstring>> = Array1600,
+    Bit.
+
+
+
+-spec xyzset(XYZ, Array1600, NewBit) -> NewArray1600
+    when XYZ          :: {xyz, X, Y, Z},
+         Array1600    :: <<_:1600>>,
+         NewBit       :: 0 | 1,
+         NewArray1600 :: Array1600,
+         X            :: 0..4,
+         Y            :: 0..4,
+         Z            :: 0..63.
+%% @private
+%% Replace the bit at {X, Y, Z} with the new bit
+%% @end
+
+xyzset(XYZ, Array1600, NewBit) ->
+    Idx0 = xyz_to_idx0(XYZ),
+    <<Pre:Idx0, _Bit:1, Post/bitstring>> = Array1600,
+    <<Pre:Idx0, NewBit:1, Post/bitstring>>.
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% 1D SUBSET ACCESSORS
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 -spec xzth(XZ, Bits) -> Column
     when XZ     :: {xz, X, Z},
          X      :: 0..4,
          Z      :: 0..63,
          Bits   :: <<_:1600>>,
          Column :: <<_:5>>.
+%% @private
+%% Fetch the column at the given X, Z coordinate pair
+%% @end
 
 xzth({xz, X, Z}, Bits) ->
     % just grab them one at a time
@@ -562,18 +1117,61 @@ xzth({xz, X, Z}, Bits) ->
     || Y <- lists:seq(0, 4)
     >>.
 
-xyzth(_, _) -> error(nyi).
 
 
--spec parity(bitstring()) -> 0 | 1.
+-spec xyth(XY, Array1600) -> Lane
-
+    when XY        :: {xy, X, Y},
+         Array1600 :: <<_:1600>>,
+         Lane      :: <<_:64>>,
+         X         :: 0..4,
+         Y         :: 0..4.
 %% @private
-%% count the number of 1s in the bitstring (is it even or odd?)
+%% Grab the lane at the given X, Y coordinate pair.
 %% @end
 
-parity(Bits) ->
+xyth({xy, X, Y}, Array1600) ->
-    parity(Bits, 0).
+    << <<( xyzth({xyz, X, Y, Z}, Array1600) ):1>>
+    || Z <- lists:seq(0, 63)
+    >>.
 
-parity(<<0:1, Rest/bitstring>>, NOnes) -> parity(Rest, NOnes);
+
-parity(<<1:1, Rest/bitstring>>, NOnes) -> parity(Rest, NOnes + 1);
+
-parity(<<>>                  , NOnes) -> NOnes rem 2.
+-spec xyset(LaneXY, Array1600, NewLane) -> NewArray1600
+    when Array1600    :: <<_:1600>>,
+         LaneXY       :: {xy, 0..4, 0..4},
+         NewLane      :: <<_:64>>,
+         NewArray1600 :: <<_:1600>>.
+%% @private
+%% Take the original array, and swap out the lane at the given x,y coordinate
+%% with the new given lane.
+%%
+%% The lane will be represented continuously so we can do a hack
+%% @end
+
+% special case when it's the last lane
+% grab the final 64 bits off the original array and replace them with the new lane
+xyset(_LaneXY = {xy, 4, 4}, <<Pre:(1600 - 64), _:64>>, NewLane) ->
+    <<Pre:(1600 - 64), NewLane/bitstring>>;
+% general case, grab the shit before the lane, grab the shit after the lane
+% replace the shit in the middle
+xyset(_LaneXY = {xy, LaneX, LaneY}, OriginalArray, NewLane) ->
+    FirstBitOfLane_Idx0    = xyz_to_idx0({xyz, LaneX, LaneY, 0}),
+    FirstBitAfterLane_Idx0 = xyz_to_idx0({xyz, LaneX, LaneY, 63}) + 1,
+    NumberOfBitsBeforeTheLane    = FirstBitOfLane_Idx0,
+    NumberOfBitsIncludingTheLane = FirstBitAfterLane_Idx0,
+    <<PreLane:NumberOfBitsBeforeTheLane   ,         _/bitstring>> = OriginalArray,
+    <<      _:NumberOfBitsIncludingTheLane, AfterLane/bitstring>> = OriginalArray,
+    Result = <<PreLane:NumberOfBitsBeforeTheLane, NewLane/bitstring, AfterLane/bitstring>>,
+    Result.
+
+%% Translating from Bitstring to hex table least significant bit first!!
+tabular(Bin) -> tabular(0, Bin).
+
+tabular(N, Bytes) when byte_size(Bytes) =< 16->
+    io:format("~3.10.0B: ~s\n", [N, string:join([begin <<X:8>> = <<H:1, G:1, F:1, E:1, D:1, C:1, B:1, A:1>>, io_lib:format("~2.16.0B", [X]) end
+                                                 || <<A:1, B:1, C:1, D:1, E:1, F:1, G:1, H:1>> <= Bytes], " ")]);
+tabular(N, <<Bytes:16/binary, Rest/binary>>) ->
+    tabular(N, Bytes), tabular(N + 16, Rest).
+
+h2b(Bytes) ->
+  << <<H:1, G:1, F:1, E:1, D:1, C:1, B:1, A:1>> || <<A:1, B:1, C:1, D:1, E:1, F:1, G:1, H:1>> <= Bytes >>.

kek_fast.erl

@@ -0,0 +1,593 @@
+%% @doc
+%% References
+%% 1. Helpful lecture: https://www.youtube.com/watch?v=JWskjzgiIa4
+%%    Notes: https://www.crypto-textbook.com/download/Understanding-Cryptography-Keccak.pdf
+%% 2. NIST standard: https://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.202.pdf
+%%    (btw: the double bar notation means "concatenate")
+%% 3. https://en.wikipedia.org/wiki/SHA-3
+%% @end
+-module(kek_fast).
+
+-compile([export_all, nowarn_export_all]).
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%% TOP LEVEL API
+%%%
+%%% sha*s and shake*s
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+-spec sha3_224(Message) -> Digest
+    when Message :: bitstring(),
+         Digest  :: <<_:224>>.
+%% @doc
+%% SHA-3 with an output bit length of 224 bits.
+%% @end
+
+sha3_224(Message) ->
+    sha3(224, Message).
+
+
+-spec sha3_256(Message) -> Digest
+    when Message :: bitstring(),
+         Digest  :: <<_:256>>.
+%% @doc
+%% SHA-3 with an output bit length of 256 bits.
+%% @end
+
+sha3_256(Message) ->
+    sha3(256, Message).
+
+
+-spec sha3_384(Message) -> Digest
+    when Message :: bitstring(),
+         Digest  :: <<_:384>>.
+%% @doc
+%% SHA-3 with an output bit length of 384 bits.
+%% @end
+
+sha3_384(Message) ->
+    sha3(384, Message).
+
+
+-spec sha3_512(Message) -> Digest
+    when Message :: bitstring(),
+         Digest  :: <<_:512>>.
+%% @doc
+%% SHA-3 with an output bit length of 512 bits.
+%% @end
+
+sha3_512(Message) ->
+    sha3(512, Message).
+
+
+-spec sha3(OutputBitLength, Message) -> Digest
+    when OutputBitLength :: pos_integer(),
+         Message         :: bitstring(),
+         Digest          :: bitstring().
+%% @doc
+%% SHA-3 with an arbitrary output bit length.
+%%
+%% This means Keccak with Capacity = 2*OutputBitLength. Additionally, SHA3
+%% concatenates the bits 01 onto the end of the input, before sending the
+%% Message to keccak/3.
+%% @end
+
+sha3(OutputBitLength, Message) ->
+    Capacity = 2*OutputBitLength,
+    keccak(Capacity, Message, <<2#01:2>>, OutputBitLength).
+
+
+-spec shake128(Message, OutputBitLength) -> Digest
+    when Message         :: bitstring(),
+         OutputBitLength :: pos_integer(),
+         Digest          :: bitstring().
+%% @doc
+%% This is the SHAKE variable-length hash with Capacity 256 = 2*128 bits.
+%% @end
+
+shake128(Message, OutputBitLength) ->
+    shake(128, Message, OutputBitLength).
+
+
+-spec shake256(Message, OutputBitLength) -> Digest
+    when Message         :: bitstring(),
+         OutputBitLength :: pos_integer(),
+         Digest          :: bitstring().
+%% @doc
+%% This is the SHAKE variable-length hash with Capacity 512 = 2*256 bits.
+%% @end
+
+shake256(Message, OutputBitLength) ->
+    shake(256, Message, OutputBitLength).
+
+
+-spec shake(ShakeNumber, Message, OutputBitLength) -> Digest
+    when ShakeNumber     :: pos_integer(),
+         Message         :: bitstring(),
+         OutputBitLength :: pos_integer(),
+         Digest          :: bitstring().
+%% @doc
+%% This is the SHAKE variable-length hash with Capacity 512 = 2*ShakeNumber bits.
+%%
+%% This concatenates the bitstring 1111 onto the end of the Message before
+%% sending the message to keccak/3.
+%% @end
+
+shake(ShakeNumber, Message, OutputBitLength) ->
+    Capacity = 2*ShakeNumber,
+    keccak(Capacity, Message, <<2#1111:4>>, OutputBitLength).
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%% OUTER KECCAK
+%%%
+%%% Keccak pads the input, absorbs it into the sponge, and squeezes the bits out
+%%% of the sponge.  The absorption and squeezing phases invoke "inner keccak",
+%%% which is the heart of the algorithm.
+%%%
+%%% - keccak/3
+%%% - pad/2
+%%% - absorb/4
+%%% - squeeze/3
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+-spec keccak(Capacity, Message, Delimiter, OutputBitLength) -> Digest
+    when Capacity        :: pos_integer(),
+         Message         :: bitstring(),
+         Delimiter       :: bitstring(),
+         OutputBitLength :: pos_integer(),
+         Digest          :: bitstring().
+%% @doc
+%% Note: this is Keccak 1600, the only one used in practice
+%%
+%% Capacity must be strictly less than 1600
+%% @end
+
+keccak(Capacity = _c, Message, Delimiter, OutputBitLength) ->
+    BitRate       = 1600 - Capacity,
+    PaddedMessage = pad(Message, Delimiter, BitRate),
+    InitialSponge = <<0:1600>>,
+    WetSponge     = absorb(PaddedMessage, BitRate, Capacity, InitialSponge),
+    ResultBits    = squeeze(WetSponge, OutputBitLength, BitRate),
+    ResultBits.
+
+
+-spec pad(Message, Delimiter, BitRate) -> NewMessage
+    when Message    :: bitstring(),
+         Delimiter  :: bitstring(),
+         BitRate    :: pos_integer(),
+         NewMessage :: bitstring().
+%% @private
+%% padding
+%% divide the message into r-bit blocks
+%%
+%% the message ends with 1000...0001
+%%
+%% sha3 calls this /10*1/ as in the regex
+%%
+%% Reference: https://en.wikipedia.org/wiki/SHA-3#Padding
+%% @end
+
+pad(Msg, Delimiter, BitRate) ->
+    MsgBits = bit_size(Msg),
+    DlmBits = bit_size(Delimiter),
+    <<Msg0:(MsgBits div 8)/bytes, Msg1/bitstring>> = Msg,
+    case (MsgBits + DlmBits) rem BitRate of
+        0 -> %% We add a complete RWord + flip the last chunk of the message
+            <<Msg0/binary, (rev_pad(Msg1, Delimiter, BitRate - 2))/bitstring>>;
+        N when N == BitRate - 1 -> %% Slightly retarded case
+            <<Msg0/binary, (rev_pad(Msg1, Delimiter, BitRate - 1))/bitstring>>;
+        N ->
+            <<Msg0/binary, (rev_pad(Msg1, Delimiter, BitRate - N - 2))/bitstring>>
+    end.
+
+%% Instead of reverting message bits, work with a "reversed" padding
+rev_pad(Msg, Delimiter, PadZeros) ->
+    Pad = <<Msg/bitstring, Delimiter/bitstring, 1:1, 0:PadZeros, 1:1>>,
+    << (flip_bits(X)) || <<X:8>> <= Pad >>.
+
+flip_bits(0) -> <<0:8>>;
+flip_bits(<<A:1, B:1, C:1, D:1, E:1, F:1, G:1, H:1>>) ->
+  <<H:1, G:1, F:1, E:1, D:1, C:1, B:1, A:1>>;
+flip_bits(N) -> flip_bits(<<N:8>>).
+
+-spec absorb(PaddedMessage, BitRate, Capacity, SpongeAcc) -> WetSponge
+    when PaddedMessage :: bitstring(),
+         BitRate       :: pos_integer(),
+         Capacity      :: pos_integer(),
+         SpongeAcc     :: <<_:1600>>,
+         WetSponge     :: <<_:1600>>.
+%% @private
+%% Assumptions:
+%%  1. BitRate + Capacity = 1600,
+%%  2. BitRate divides the PaddedMessage length (i.e. already have done padding)
+%% @end
+
+% can pull off r bits from the start of the message
+absorb(PaddedMessageBits, BitRate = _r, Capacity = _c, Sponge) when BitRate =< bit_size(PaddedMessageBits) ->
+    <<ThisRWord:BitRate, Rest/bitstring>> = PaddedMessageBits,
+    % we bitwise xor the sponge against the r word followed by a bunch of 0s
+    <<SpongeInt:1600>> = Sponge,
+    <<Foo:1600>>       = <<ThisRWord:BitRate, 0:Capacity>>,
+    FInputInt          = SpongeInt bxor Foo,
+    FInputBits         = <<FInputInt:1600>>,
+    NewSponge          = inner_keccak(FInputBits),
+    absorb(Rest, BitRate, Capacity, NewSponge);
+% empty string, return the sponge
+absorb(<<>>, _r, _c, FinalSponge) ->
+    FinalSponge.
+
+
+
+-spec squeeze(WetSponge, OutputBitLength, BitRate) -> ResultBits
+    when WetSponge       :: <<_:1600>>,
+         OutputBitLength :: pos_integer(),
+         BitRate         :: pos_integer(),
+         ResultBits      :: bitstring().
+%% @private
+%% squeeze the output bits out of the sponge
+%% @end
+
+%%% % simple case: bit length is less than (or equal to) the sponge size, just grab
+%%% % the first ones
+%%% % this is the case for the shas
+%%% squeeze(<<ResultBits:OutputBitLength, _Rest/bitstring>>, OutputBitLength, _BitRate) ->
+%%%     <<ResultBits:OutputBitLength>>;
+% general case: output bit length is greater than the sponge size, construct
+% accumulatively
+% this is the case for the variable-length encodings
+squeeze(WetSponge, OutputBitLength, BitRate) ->
+    InitOutputAcc = <<>>,
+    really_squeeze(WetSponge, OutputBitLength, BitRate, InitOutputAcc).
+
+% terminal case: we have enough bits in the output, return those
+really_squeeze(_WetSponge, OutputBitLength, _BitRate, FinalAccBits) when OutputBitLength =< bit_size(FinalAccBits) ->
+    <<ResultBits:OutputBitLength, _/bitstring>> = FinalAccBits,
+    <<ResultBits:OutputBitLength>>;
+% general case: need moar bits
+% in this case
+%   - we grab the first r bits of the sponge, add them to the accumulator
+%   - re-kek the sponge
+%   - try again
+really_squeeze(WetSponge, OutputBitLength, BitRate, ResultAcc)->
+    <<ThisRWord:BitRate, _/bitstring>> = WetSponge,
+    NewResultAcc                       = <<ResultAcc/bitstring, ThisRWord:BitRate>>,
+    NewWetSponge                       = inner_keccak(WetSponge),
+    really_squeeze(NewWetSponge, OutputBitLength, BitRate, NewResultAcc).
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%% THE DREADED INNER KECCAK
+%%%
+%%% This is the "f" function that appears in all the documentation.
+%%%
+%%% The input is the 1600-bit sponge array. inner_keccak/1 sends the input
+%%% through 24 "rounds". Each round consists of the 5 Greek letter steps, each of
+%%% which is a weird transformation on the array.
+%%%
+%%% Here the rounds are unrolled in terms of 64bit integers - for efficiency
+%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+inner_keccak(<<_:1600>> = State) ->
+  IntState0 = list_to_tuple([ X || <<X:64/little>> <= State ]),
+  IntState6 = inner_keccak_fast(IntState0),
+  << <<X:64/little>> || X <- tuple_to_list(IntState6) >>.
+
+inner_keccak_fast(IntState0) ->
+  IntState1 = fast_round(IntState0, {16#0000000000000001, 16#0000000000008082, 16#800000000000808A, 16#8000000080008000}),
+  IntState2 = fast_round(IntState1, {16#000000000000808B, 16#0000000080000001, 16#8000000080008081, 16#8000000000008009}),
+  IntState3 = fast_round(IntState2, {16#000000000000008A, 16#0000000000000088, 16#0000000080008009, 16#000000008000000A}),
+  IntState4 = fast_round(IntState3, {16#000000008000808B, 16#800000000000008B, 16#8000000000008089, 16#8000000000008003}),
+  IntState5 = fast_round(IntState4, {16#8000000000008002, 16#8000000000000080, 16#000000000000800A, 16#800000008000000A}),
+  fast_round(IntState5, {16#8000000080008081, 16#8000000000008080, 16#0000000080000001, 16#8000000080008008}).
+
+-define(INT64, 16#FFFFFFFFFFFFFFFF).
+-define(BSL64(X, N), ((X bsl N) band ?INT64)).
+-define(BSR64(X, N), (X bsr N)).
+
+-define(ROTL64(X, N), (?BSL64(X, N) bor ?BSR64(X, (64 - N)))).
+-define(CAN64(A, B), ((A bxor B) band A)).
+
+fast_round(As0, {RC0, RC1, RC2, RC3}) ->
+  As1 = fast_round1(As0, RC0),
+  As2 = fast_round2(As1, RC1),
+  As3 = fast_round3(As2, RC2),
+  fast_round4(As3, RC3).
+
+fast_round1({A0, A1, A2, A3, A4, A5, A6, A7, A8, A9, A10, A11, A12, A13, A14, A15, A16, A17, A18, A19, A20, A21, A22, A23, A24}, RC0) ->
+  %% Round 1
+  BC0_0 = A0 bxor A5 bxor A10 bxor A15 bxor A20,
+  BC1_0 = A1 bxor A6 bxor A11 bxor A16 bxor A21,
+  BC2_0 = A2 bxor A7 bxor A12 bxor A17 bxor A22,
+  BC3_0 = A3 bxor A8 bxor A13 bxor A18 bxor A23,
+  BC4_0 = A4 bxor A9 bxor A14 bxor A19 bxor A24,
+  D0 = BC4_0 bxor ?ROTL64(BC1_0, 1),
+  D1 = BC0_0 bxor ?ROTL64(BC2_0, 1),
+  D2 = BC1_0 bxor ?ROTL64(BC3_0, 1),
+  D3 = BC2_0 bxor ?ROTL64(BC4_0, 1),
+  D4 = BC3_0 bxor ?ROTL64(BC0_0, 1),
+
+  BC0_1 = A0 bxor D0,
+  BC1_1 = ?ROTL64(A6 bxor D1, 44),
+  BC2_1 = ?ROTL64(A12 bxor D2, 43),
+  BC3_1 = ?ROTL64(A18 bxor D3, 21),
+  BC4_1 = ?ROTL64(A24 bxor D4, 14),
+  A0_1  = BC0_1 bxor ?CAN64(BC2_1, BC1_1) bxor RC0,
+  A6_1  = BC1_1 bxor ?CAN64(BC3_1, BC2_1),
+  A12_1 = BC2_1 bxor ?CAN64(BC4_1, BC3_1),
+  A18_1 = BC3_1 bxor ?CAN64(BC0_1, BC4_1),
+  A24_1 = BC4_1 bxor ?CAN64(BC1_1, BC0_1),
+
+  BC2_2 = ?ROTL64(A10 bxor D0, 3),
+  BC3_2 = ?ROTL64(A16 bxor D1, 45),
+  BC4_2 = ?ROTL64(A22 bxor D2, 61),
+  BC0_2 = ?ROTL64(A3 bxor D3, 28),
+  BC1_2 = ?ROTL64(A9 bxor D4, 20),
+  A10_1 = BC0_2 bxor ?CAN64(BC2_2, BC1_2),
+  A16_1 = BC1_2 bxor ?CAN64(BC3_2, BC2_2),
+  A22_1 = BC2_2 bxor ?CAN64(BC4_2, BC3_2),
+  A3_1  = BC3_2 bxor ?CAN64(BC0_2, BC4_2),
+  A9_1  = BC4_2 bxor ?CAN64(BC1_2, BC0_2),
+
+  BC4_3 = ?ROTL64(A20 bxor D0, 18),
+  BC0_3 = ?ROTL64(A1 bxor D1, 1),
+  BC1_3 = ?ROTL64(A7 bxor D2, 6),
+  BC2_3 = ?ROTL64(A13 bxor D3, 25),
+  BC3_3 = ?ROTL64(A19 bxor D4, 8),
+  A20_1 = BC0_3 bxor ?CAN64(BC2_3, BC1_3),
+  A1_1  = BC1_3 bxor ?CAN64(BC3_3, BC2_3),
+  A7_1  = BC2_3 bxor ?CAN64(BC4_3, BC3_3),
+  A13_1 = BC3_3 bxor ?CAN64(BC0_3, BC4_3),
+  A19_1 = BC4_3 bxor ?CAN64(BC1_3, BC0_3),
+
+  BC1_4 = ?ROTL64(A5 bxor D0, 36),
+  BC2_4 = ?ROTL64(A11 bxor D1, 10),
+  BC3_4 = ?ROTL64(A17 bxor D2, 15),
+  BC4_4 = ?ROTL64(A23 bxor D3, 56),
+  BC0_4 = ?ROTL64(A4 bxor D4, 27),
+  A5_1  = BC0_4 bxor ?CAN64(BC2_4, BC1_4),
+  A11_1 = BC1_4 bxor ?CAN64(BC3_4, BC2_4),
+  A17_1 = BC2_4 bxor ?CAN64(BC4_4, BC3_4),
+  A23_1 = BC3_4 bxor ?CAN64(BC0_4, BC4_4),
+  A4_1  = BC4_4 bxor ?CAN64(BC1_4, BC0_4),
+
+  BC3_5 = ?ROTL64(A15 bxor D0, 41),
+  BC4_5 = ?ROTL64(A21 bxor D1, 2),
+  BC0_5 = ?ROTL64(A2 bxor D2, 62),
+  BC1_5 = ?ROTL64(A8 bxor D3, 55),
+  BC2_5 = ?ROTL64(A14 bxor D4, 39),
+  A15_1 = BC0_5 bxor ?CAN64(BC2_5, BC1_5),
+  A21_1 = BC1_5 bxor ?CAN64(BC3_5, BC2_5),
+  A2_1  = BC2_5 bxor ?CAN64(BC4_5, BC3_5),
+  A8_1  = BC3_5 bxor ?CAN64(BC0_5, BC4_5),
+  A14_1 = BC4_5 bxor ?CAN64(BC1_5, BC0_5),
+
+  {A0_1, A1_1, A2_1, A3_1, A4_1, A5_1, A6_1, A7_1, A8_1, A9_1, A10_1, A11_1, A12_1, A13_1,
+   A14_1, A15_1, A16_1, A17_1, A18_1, A19_1, A20_1, A21_1, A22_1, A23_1, A24_1}.
+
+fast_round2({A0, A1, A2, A3, A4, A5, A6, A7, A8, A9, A10, A11, A12, A13, A14, A15, A16, A17, A18, A19, A20, A21, A22, A23, A24}, RC1) ->
+  %% Round 2
+  BC0_0 = A0 bxor A5 bxor A10 bxor A15 bxor A20,
+  BC1_0 = A1 bxor A6 bxor A11 bxor A16 bxor A21,
+  BC2_0 = A2 bxor A7 bxor A12 bxor A17 bxor A22,
+  BC3_0 = A3 bxor A8 bxor A13 bxor A18 bxor A23,
+  BC4_0 = A4 bxor A9 bxor A14 bxor A19 bxor A24,
+  D0 = BC4_0 bxor ?ROTL64(BC1_0, 1),
+  D1 = BC0_0 bxor ?ROTL64(BC2_0, 1),
+  D2 = BC1_0 bxor ?ROTL64(BC3_0, 1),
+  D3 = BC2_0 bxor ?ROTL64(BC4_0, 1),
+  D4 = BC3_0 bxor ?ROTL64(BC0_0, 1),
+
+  BC0_1 = A0 bxor D0,
+  BC1_1 = ?ROTL64(A16 bxor D1, 44),
+  BC2_1 = ?ROTL64(A7 bxor D2, 43),
+  BC3_1 = ?ROTL64(A23 bxor D3, 21),
+  BC4_1 = ?ROTL64(A14 bxor D4, 14),
+  A0_1  = BC0_1 bxor ?CAN64(BC2_1, BC1_1) bxor RC1,
+  A16_1 = BC1_1 bxor ?CAN64(BC3_1, BC2_1),
+  A7_1  = BC2_1 bxor ?CAN64(BC4_1, BC3_1),
+  A23_1 = BC3_1 bxor ?CAN64(BC0_1, BC4_1),
+  A14_1 = BC4_1 bxor ?CAN64(BC1_1, BC0_1),
+
+  BC2_2 = ?ROTL64(A20 bxor D0, 3),
+  BC3_2 = ?ROTL64(A11 bxor D1, 45),
+  BC4_2 = ?ROTL64(A2 bxor D2, 61),
+  BC0_2 = ?ROTL64(A18 bxor D3, 28),
+  BC1_2 = ?ROTL64(A9 bxor D4, 20),
+  A20_1 = BC0_2 bxor ?CAN64(BC2_2, BC1_2),
+  A11_1 = BC1_2 bxor ?CAN64(BC3_2, BC2_2),
+  A2_1  = BC2_2 bxor ?CAN64(BC4_2, BC3_2),
+  A18_1 = BC3_2 bxor ?CAN64(BC0_2, BC4_2),
+  A9_1  = BC4_2 bxor ?CAN64(BC1_2, BC0_2),
+
+  BC4_3 = ?ROTL64(A15 bxor D0, 18),
+  BC0_3 = ?ROTL64(A6 bxor D1, 1),
+  BC1_3 = ?ROTL64(A22 bxor D2, 6),
+  BC2_3 = ?ROTL64(A13 bxor D3, 25),
+  BC3_3 = ?ROTL64(A4 bxor D4, 8),
+  A15_1 = BC0_3 bxor ?CAN64(BC2_3, BC1_3),
+  A6_1  = BC1_3 bxor ?CAN64(BC3_3, BC2_3),
+  A22_1 = BC2_3 bxor ?CAN64(BC4_3, BC3_3),
+  A13_1 = BC3_3 bxor ?CAN64(BC0_3, BC4_3),
+  A4_1  = BC4_3 bxor ?CAN64(BC1_3, BC0_3),
+
+  BC1_4 = ?ROTL64(A10 bxor D0, 36),
+  BC2_4 = ?ROTL64(A1 bxor D1, 10),
+  BC3_4 = ?ROTL64(A17 bxor D2, 15),
+  BC4_4 = ?ROTL64(A8 bxor D3, 56),
+  BC0_4 = ?ROTL64(A24 bxor D4, 27),
+  A10_1 = BC0_4 bxor ?CAN64(BC2_4, BC1_4),
+  A1_1  = BC1_4 bxor ?CAN64(BC3_4, BC2_4),
+  A17_1 = BC2_4 bxor ?CAN64(BC4_4, BC3_4),
+  A8_1  = BC3_4 bxor ?CAN64(BC0_4, BC4_4),
+  A24_1 = BC4_4 bxor ?CAN64(BC1_4, BC0_4),
+
+  BC3_5 = ?ROTL64(A5 bxor D0, 41),
+  BC4_5 = ?ROTL64(A21 bxor D1, 2),
+  BC0_5 = ?ROTL64(A12 bxor D2, 62),
+  BC1_5 = ?ROTL64(A3 bxor D3, 55),
+  BC2_5 = ?ROTL64(A19 bxor D4, 39),
+  A5_1  = BC0_5 bxor ?CAN64(BC2_5, BC1_5),
+  A21_1 = BC1_5 bxor ?CAN64(BC3_5, BC2_5),
+  A12_1 = BC2_5 bxor ?CAN64(BC4_5, BC3_5),
+  A3_1  = BC3_5 bxor ?CAN64(BC0_5, BC4_5),
+  A19_1 = BC4_5 bxor ?CAN64(BC1_5, BC0_5),
+
+  {A0_1, A1_1, A2_1, A3_1, A4_1, A5_1, A6_1, A7_1, A8_1, A9_1, A10_1, A11_1, A12_1, A13_1,
+   A14_1, A15_1, A16_1, A17_1, A18_1, A19_1, A20_1, A21_1, A22_1, A23_1, A24_1}.
+
+fast_round3({A0, A1, A2, A3, A4, A5, A6, A7, A8, A9, A10, A11, A12, A13, A14, A15, A16, A17, A18, A19, A20, A21, A22, A23, A24}, RC2) ->
+  %% Round 3
+  BC0_0 = A0 bxor A5 bxor A10 bxor A15 bxor A20,
+  BC1_0 = A1 bxor A6 bxor A11 bxor A16 bxor A21,
+  BC2_0 = A2 bxor A7 bxor A12 bxor A17 bxor A22,
+  BC3_0 = A3 bxor A8 bxor A13 bxor A18 bxor A23,
+  BC4_0 = A4 bxor A9 bxor A14 bxor A19 bxor A24,
+  D0 = BC4_0 bxor ?ROTL64(BC1_0, 1),
+  D1 = BC0_0 bxor ?ROTL64(BC2_0, 1),
+  D2 = BC1_0 bxor ?ROTL64(BC3_0, 1),
+  D3 = BC2_0 bxor ?ROTL64(BC4_0, 1),
+  D4 = BC3_0 bxor ?ROTL64(BC0_0, 1),
+
+  BC0_1 = A0 bxor D0,
+  BC1_1 = ?ROTL64(A11 bxor D1, 44),
+  BC2_1 = ?ROTL64(A22 bxor D2, 43),
+  BC3_1 = ?ROTL64(A8 bxor D3, 21),
+  BC4_1 = ?ROTL64(A19 bxor D4, 14),
+  A0_1  = BC0_1 bxor ?CAN64(BC2_1, BC1_1) bxor RC2,
+  A11_1 = BC1_1 bxor ?CAN64(BC3_1, BC2_1),
+  A22_1 = BC2_1 bxor ?CAN64(BC4_1, BC3_1),
+  A8_1  = BC3_1 bxor ?CAN64(BC0_1, BC4_1),
+  A19_1 = BC4_1 bxor ?CAN64(BC1_1, BC0_1),
+
+  BC2_2 = ?ROTL64(A15 bxor D0, 3),
+  BC3_2 = ?ROTL64(A1 bxor D1, 45),
+  BC4_2 = ?ROTL64(A12 bxor D2, 61),
+  BC0_2 = ?ROTL64(A23 bxor D3, 28),
+  BC1_2 = ?ROTL64(A9 bxor D4, 20),
+  A15_1 = BC0_2 bxor ?CAN64(BC2_2, BC1_2),
+  A1_1  = BC1_2 bxor ?CAN64(BC3_2, BC2_2),
+  A12_1 = BC2_2 bxor ?CAN64(BC4_2, BC3_2),
+  A23_1 = BC3_2 bxor ?CAN64(BC0_2, BC4_2),
+  A9_1  = BC4_2 bxor ?CAN64(BC1_2, BC0_2),
+
+  BC4_3 = ?ROTL64(A5 bxor D0, 18),
+  BC0_3 = ?ROTL64(A16 bxor D1, 1),
+  BC1_3 = ?ROTL64(A2 bxor D2, 6),
+  BC2_3 = ?ROTL64(A13 bxor D3, 25),
+  BC3_3 = ?ROTL64(A24 bxor D4, 8),
+  A5_1  = BC0_3 bxor ?CAN64(BC2_3, BC1_3),
+  A16_1 = BC1_3 bxor ?CAN64(BC3_3, BC2_3),
+  A2_1  = BC2_3 bxor ?CAN64(BC4_3, BC3_3),
+  A13_1 = BC3_3 bxor ?CAN64(BC0_3, BC4_3),
+  A24_1 = BC4_3 bxor ?CAN64(BC1_3, BC0_3),
+
+  BC1_4 = ?ROTL64(A20 bxor D0, 36),
+  BC2_4 = ?ROTL64(A6 bxor D1, 10),
+  BC3_4 = ?ROTL64(A17 bxor D2, 15),
+  BC4_4 = ?ROTL64(A3 bxor D3, 56),
+  BC0_4 = ?ROTL64(A14 bxor D4, 27),
+  A20_1 = BC0_4 bxor ?CAN64(BC2_4, BC1_4),
+  A6_1  = BC1_4 bxor ?CAN64(BC3_4, BC2_4),
+  A17_1 = BC2_4 bxor ?CAN64(BC4_4, BC3_4),
+  A3_1  = BC3_4 bxor ?CAN64(BC0_4, BC4_4),
+  A14_1 = BC4_4 bxor ?CAN64(BC1_4, BC0_4),
+
+  BC3_5 = ?ROTL64(A10 bxor D0, 41),
+  BC4_5 = ?ROTL64(A21 bxor D1, 2),
+  BC0_5 = ?ROTL64(A7 bxor D2, 62),
+  BC1_5 = ?ROTL64(A18 bxor D3, 55),
+  BC2_5 = ?ROTL64(A4 bxor D4, 39),
+  A10_1 = BC0_5 bxor ?CAN64(BC2_5, BC1_5),
+  A21_1 = BC1_5 bxor ?CAN64(BC3_5, BC2_5),
+  A7_1  = BC2_5 bxor ?CAN64(BC4_5, BC3_5),
+  A18_1 = BC3_5 bxor ?CAN64(BC0_5, BC4_5),
+  A4_1  = BC4_5 bxor ?CAN64(BC1_5, BC0_5),
+
+  {A0_1, A1_1, A2_1, A3_1, A4_1, A5_1, A6_1, A7_1, A8_1, A9_1, A10_1, A11_1, A12_1, A13_1,
+   A14_1, A15_1, A16_1, A17_1, A18_1, A19_1, A20_1, A21_1, A22_1, A23_1, A24_1}.
+
+fast_round4({A0, A1, A2, A3, A4, A5, A6, A7, A8, A9, A10, A11, A12, A13, A14, A15, A16, A17, A18, A19, A20, A21, A22, A23, A24}, RC3) ->
+  %% Round 4
+  BC0_0 = A0 bxor A5 bxor A10 bxor A15 bxor A20,
+  BC1_0 = A1 bxor A6 bxor A11 bxor A16 bxor A21,
+  BC2_0 = A2 bxor A7 bxor A12 bxor A17 bxor A22,
+  BC3_0 = A3 bxor A8 bxor A13 bxor A18 bxor A23,
+  BC4_0 = A4 bxor A9 bxor A14 bxor A19 bxor A24,
+  D0 = BC4_0 bxor ?ROTL64(BC1_0, 1),
+  D1 = BC0_0 bxor ?ROTL64(BC2_0, 1),
+  D2 = BC1_0 bxor ?ROTL64(BC3_0, 1),
+  D3 = BC2_0 bxor ?ROTL64(BC4_0, 1),
+  D4 = BC3_0 bxor ?ROTL64(BC0_0, 1),
+
+  BC0_1 = A0 bxor D0,
+  BC1_1 = ?ROTL64(A1 bxor D1, 44),
+  BC2_1 = ?ROTL64(A2 bxor D2, 43),
+  BC3_1 = ?ROTL64(A3 bxor D3, 21),
+  BC4_1 = ?ROTL64(A4 bxor D4, 14),
+  A0_1 = BC0_1 bxor ?CAN64(BC2_1, BC1_1) bxor RC3,
+  A1_1 = BC1_1 bxor ?CAN64(BC3_1, BC2_1),
+  A2_1 = BC2_1 bxor ?CAN64(BC4_1, BC3_1),
+  A3_1 = BC3_1 bxor ?CAN64(BC0_1, BC4_1),
+  A4_1 = BC4_1 bxor ?CAN64(BC1_1, BC0_1),
+
+  BC2_2 = ?ROTL64(A5 bxor D0, 3),
+  BC3_2 = ?ROTL64(A6 bxor D1, 45),
+  BC4_2 = ?ROTL64(A7 bxor D2, 61),
+  BC0_2 = ?ROTL64(A8 bxor D3, 28),
+  BC1_2 = ?ROTL64(A9 bxor D4, 20),
+  A5_1 = BC0_2 bxor ?CAN64(BC2_2, BC1_2),
+  A6_1 = BC1_2 bxor ?CAN64(BC3_2, BC2_2),
+  A7_1 = BC2_2 bxor ?CAN64(BC4_2, BC3_2),
+  A8_1 = BC3_2 bxor ?CAN64(BC0_2, BC4_2),
+  A9_1 = BC4_2 bxor ?CAN64(BC1_2, BC0_2),
+
+  BC4_3 = ?ROTL64(A10 bxor D0, 18),
+  BC0_3 = ?ROTL64(A11 bxor D1, 1),
+  BC1_3 = ?ROTL64(A12 bxor D2, 6),
+  BC2_3 = ?ROTL64(A13 bxor D3, 25),
+  BC3_3 = ?ROTL64(A14 bxor D4, 8),
+  A10_1 = BC0_3 bxor ?CAN64(BC2_3, BC1_3),
+  A11_1 = BC1_3 bxor ?CAN64(BC3_3, BC2_3),
+  A12_1 = BC2_3 bxor ?CAN64(BC4_3, BC3_3),
+  A13_1 = BC3_3 bxor ?CAN64(BC0_3, BC4_3),
+  A14_1 = BC4_3 bxor ?CAN64(BC1_3, BC0_3),
+
+  BC1_4 = ?ROTL64(A15 bxor D0, 36),
+  BC2_4 = ?ROTL64(A16 bxor D1, 10),
+  BC3_4 = ?ROTL64(A17 bxor D2, 15),
+  BC4_4 = ?ROTL64(A18 bxor D3, 56),
+  BC0_4 = ?ROTL64(A19 bxor D4, 27),
+  A15_1 = BC0_4 bxor ?CAN64(BC2_4, BC1_4),
+  A16_1 = BC1_4 bxor ?CAN64(BC3_4, BC2_4),
+  A17_1 = BC2_4 bxor ?CAN64(BC4_4, BC3_4),
+  A18_1 = BC3_4 bxor ?CAN64(BC0_4, BC4_4),
+  A19_1 = BC4_4 bxor ?CAN64(BC1_4, BC0_4),
+
+  BC3_5 = ?ROTL64(A20 bxor D0, 41),
+  BC4_5 = ?ROTL64(A21 bxor D1, 2),
+  BC0_5 = ?ROTL64(A22 bxor D2, 62),
+  BC1_5 = ?ROTL64(A23 bxor D3, 55),
+  BC2_5 = ?ROTL64(A24 bxor D4, 39),
+  A20_1 = BC0_5 bxor ?CAN64(BC2_5, BC1_5),
+  A21_1 = BC1_5 bxor ?CAN64(BC3_5, BC2_5),
+  A22_1 = BC2_5 bxor ?CAN64(BC4_5, BC3_5),
+  A23_1 = BC3_5 bxor ?CAN64(BC0_5, BC4_5),
+  A24_1 = BC4_5 bxor ?CAN64(BC1_5, BC0_5),
+
+  {A0_1, A1_1, A2_1, A3_1, A4_1, A5_1, A6_1, A7_1, A8_1, A9_1, A10_1, A11_1, A12_1, A13_1,
+   A14_1, A15_1, A16_1, A17_1, A18_1, A19_1, A20_1, A21_1, A22_1, A23_1, A24_1}.

rc.erl

@@ -0,0 +1,198 @@
+%% @doc
+%% this module exists to check that the table on pp. 12 of
+%% https://www.crypto-textbook.com/download/Understanding-Cryptography-Keccak.pdf
+%% is correct
+-module(rc).
+
+-compile([export_all, nowarn_export_all]).
+
+-spec little_rc(T) -> Bit
+    when T   :: non_neg_integer(),
+         Bit :: 0 | 1.
+%% copying from pp. 16 of https://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.202.pdf
+%%
+%% n.b. the || operator there i think means concatenate
+
+little_rc(T) when (T rem 255) =:= 0 ->
+    1;
+little_rc(T) ->
+    R       = <<2#1000_0000>>,
+    InitI   = 1,
+    TMod255 = T rem 255,
+    NewR = little_rc(InitI, TMod255, R),
+    <<Result:1, _/bitstring>> = NewR,
+    Result.
+
+-spec little_rc(I, MaxI, R) -> NewR
+    when I    :: pos_integer(),
+         MaxI :: pos_integer(),
+         R    :: <<_:8>>,
+         NewR :: <<_:8>>.
+
+little_rc(I, MaxI, R) when I =< MaxI ->
+    R_  = <<0:1, R/bitstring>>,
+    % need R_[0], R_[4], R_[5], R_[6], and R_[8]
+    <<R_0:1, R_123:3, R_4:1, R_5:1, R_6:1, R_7:1, R_8:1>> = R_,
+    NewR_0 = R_0 bxor R_8,
+    NewR_4 = R_4 bxor R_8,
+    NewR_5 = R_5 bxor R_8,
+    NewR_6 = R_6 bxor R_8,
+    NewR   = <<NewR_0  :1,
+                  R_123:3,
+               NewR_4  :1,
+               NewR_5  :1,
+               NewR_6  :1,
+                  R_7  :1>>,
+    NewI = I + 1,
+    little_rc(NewI, MaxI, NewR);
+little_rc(I, MaxI, R) when I > MaxI ->
+    R.
+
+
+% Table from pp. 12 of https://www.crypto-textbook.com/download/Understanding-Cryptography-Keccak.pdf
+%
+% RC[ 0] = 0x0000000000000001
+% RC[ 1] = 0x0000000000008082
+% RC[ 2] = 0x800000000000808A
+% RC[ 3] = 0x8000000080008000
+% RC[ 4] = 0x000000000000808B
+% RC[ 5] = 0x0000000080000001
+% RC[ 6] = 0x8000000080008081
+% RC[ 7] = 0x8000000000008009
+% RC[ 8] = 0x000000000000008A
+% RC[ 9] = 0x0000000000000088
+% RC[10] = 0x0000000080008009
+% RC[11] = 0x000000008000000A
+% RC[12] = 0x000000008000808B
+% RC[13] = 0x800000000000008B
+% RC[14] = 0x8000000000008089
+% RC[15] = 0x8000000000008003
+% RC[16] = 0x8000000000008002
+% RC[17] = 0x8000000000000080
+% RC[18] = 0x000000000000800A
+% RC[19] = 0x800000008000000A
+% RC[20] = 0x8000000080008081
+% RC[21] = 0x8000000000008080
+% RC[22] = 0x0000000080000001
+% RC[23] = 0x8000000080008008
+
+-spec big_rc(RoundIndex) -> BigRC
+    when RoundIndex :: 0..23,
+         BigRC      :: <<_:64>>.
+
+%% from pp. 16 of the NIST doc:
+%%  > 2. Let RC = 0w.
+%%  > 3. For j from 0 to l, let RC[2^j – 1] = rc(j + 7*ir).
+%%
+%% In this case, 0w = <<0:64>>, l = 6, ir ranges between 0..23
+%%
+%%  > 2. Let RC = 0w.
+%%  > 3. For j from 0 to 6, let RC[2^j – 1] = rc(j + 7*ir).
+
+big_rc(RoundIndex) ->
+    InitJ     = 0,
+    InitBigRC = <<0:64>>,
+    big_rc(RoundIndex, InitJ, InitBigRC).
+
+
+-spec big_rc(RoundIndex, J, BigRCAcc) -> BigRC
+    when RoundIndex :: 0..23,
+         J          :: 0..6,
+         BigRCAcc   :: BigRC,
+         BigRC      :: <<_:64>>.
+
+%%  > 3. For j from 0 to 6, let RC[2^j – 1] = rc(j + 7ir).
+big_rc(RoundIndex, J, BigRCAcc) when 0 =< J, J =< 6 ->
+    Idx0WeAreModifying                       = two_to_the(J) - 1,
+    NewBit                                   = little_rc(J + 7*RoundIndex),
+    NumSkipBits                              = Idx0WeAreModifying,
+    <<Pre:NumSkipBits, _:1, Post/bitstring>> = BigRCAcc,
+    NewJ                                     = J + 1,
+    NewBigRCAcc                              = <<Pre:NumSkipBits, NewBit:1, Post/bitstring>>,
+    big_rc(RoundIndex, NewJ, NewBigRCAcc);
+big_rc(_RoundIndex, J, BigRCAcc) when J > 6 ->
+    BigRCAcc.
+
+two_to_the(N) when 0 =< N ->
+    1 bsl N.
+
+
+%% expected values
+xrc( 0) -> <<( 16#0000000000000001 ):64>>;
+xrc( 1) -> <<( 16#0000000000008082 ):64>>;
+xrc( 2) -> <<( 16#800000000000808A ):64>>;
+xrc( 3) -> <<( 16#8000000080008000 ):64>>;
+xrc( 4) -> <<( 16#000000000000808B ):64>>;
+xrc( 5) -> <<( 16#0000000080000001 ):64>>;
+xrc( 6) -> <<( 16#8000000080008081 ):64>>;
+xrc( 7) -> <<( 16#8000000000008009 ):64>>;
+xrc( 8) -> <<( 16#000000000000008A ):64>>;
+xrc( 9) -> <<( 16#0000000000000088 ):64>>;
+xrc(10) -> <<( 16#0000000080008009 ):64>>;
+xrc(11) -> <<( 16#000000008000000A ):64>>;
+xrc(12) -> <<( 16#000000008000808B ):64>>;
+xrc(13) -> <<( 16#800000000000008B ):64>>;
+xrc(14) -> <<( 16#8000000000008089 ):64>>;
+xrc(15) -> <<( 16#8000000000008003 ):64>>;
+xrc(16) -> <<( 16#8000000000008002 ):64>>;
+xrc(17) -> <<( 16#8000000000000080 ):64>>;
+xrc(18) -> <<( 16#000000000000800A ):64>>;
+xrc(19) -> <<( 16#800000008000000A ):64>>;
+xrc(20) -> <<( 16#8000000080008081 ):64>>;
+xrc(21) -> <<( 16#8000000000008080 ):64>>;
+xrc(22) -> <<( 16#0000000080000001 ):64>>;
+xrc(23) -> <<( 16#8000000080008008 ):64>>.
+
+
+%% all false
+check() ->
+    CheckI =
+        fun(I) ->
+            io:format("I = ~p: ~p~n", [I, big_rc(I) =:= xrc(I)])
+        end,
+    lists:foreach(CheckI, lists:seq(0, 23)).
+
+%% checks if the expectation equals the outcome in the reversed case
+%% all true
+check_rev() ->
+    CheckI =
+        fun(I) ->
+            BigResult = str64(big_rc(I)),
+            ExpResult = lists:reverse(str64(xrc(I))),
+            io:format("I = ~p: ~p~n", [I, BigResult =:= ExpResult])
+        end,
+    lists:foreach(CheckI, lists:seq(0, 23)).
+
+
+
+%% format a 64 bit binary as a string (no newline)
+str64(<<N:64>>) ->
+    %% prints with 0 as the padding character
+    iolist_to_list(io_lib:format("~64.2.0B", [N])).
+
+%% prints with _ as the padding character, newline
+print64_(<<N:64>>) ->
+    io:format("~64.2._B~n", [N]).
+
+iolist_to_list(IoList) ->
+    binary_to_list(iolist_to_binary(IoList)).
+
+
+%% print out big_rc(N) -> RC(N)
+print_big_rcs() ->
+    pbr(0).
+
+% 1 digit numbers
+pbr(N) when 0 =< N, N =< 9 ->
+    <<ThisRC_int:64>> = big_rc(N),
+    io:format("round_constant_int( ~p) -> ~p;~n", [N, ThisRC_int]),
+    pbr(N + 1);
+% 2 digit numbers
+pbr(N) when 10 =< N, N =< 22 ->
+    <<ThisRC_int:64>> = big_rc(N),
+    io:format("round_constant_int(~p) -> ~p;~n", [N, ThisRC_int]),
+    pbr(N + 1);
+% terminal case
+pbr(N) when N =:= 23 ->
+    <<ThisRC_int:64>> = big_rc(N),
+    io:format("round_constant_int(~p) -> ~p.~n", [N, ThisRC_int]).