%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Copyright: 2010 Integrated Sytems Laboratory, ETH Zurich
%%            http://www.iis.ee.ethz.ch/~sha3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = NTT(M, f)

y = zeros(1,128);
y(1:64) = M;

FFT64_8_8_Twiddle = [
    1,    1,    1,    1,    1,    1,    1,    1;
    1,    2,    4,    8,   16,   32,   64,  128;
    1,   60,    2,  120,    4,  -17,    8,  -34;
    1,  120,    8,  -68,   64,  -30,   -2,   17; 
    1,   46,   60,  -67,    2,   92,  120,  123; 
    1,   92,  -17,  -22,   32,  117,  -30,   67; 
    1,  -67,  120,  -73,    8,  -22,  -68,  -70; 
    1,  123,  -34,  -70,  128,   67,   17,   35];

FFT128_2_64_Twiddle =  [
    1, -118,   46,  -31,   60,  116,  -67,  -61, ...
    2,   21,   92,  -62,  120,  -25,  123, -122, ...
    4,   42,  -73, -124,  -17,  -50,  -11,   13, ...
    8,   84,  111,    9,  -34, -100,  -22,   26, ...
   16,  -89,  -35,   18,  -68,   57,  -44,   52, ...
   32,   79,  -70,   36,  121,  114,  -88,  104, ...
   64,  -99,  117,   72,  -15,  -29,   81,  -49, ...
  128,   59,  -23, -113,  -30,  -58,  -95,  -98];





P = [1
 65
 33
 97
 17
 81
 49
113
  9
 73
 41
105
 25
 89
 57
121
  5
 69
 37
101
 21
 85
 53
117
 13
 77
 45
109
 29
 93
 61
125
  3
 67
 35
 99
 19
 83
 51
115
 11
 75
 43
107
 27
 91
 59
123
  7
 71
 39
103
 23
 87
 55
119
 15
 79
 47
111
 31
 95
 63
127
  2
 66
 34
 98
 18
 82
 50
114
 10
 74
 42
106
 26
 90
 58
122
  6
 70
 38
102
 22
 86
 54
118
 14
 78
 46
110
 30
 94
 62
126
  4
 68
 36
100
 20
 84
 52
116
 12
 76
 44
108
 28
 92
 60
124
  8
 72
 40
104
 24
 88
 56
120
 16
 80
 48
112
 32
 96
 64
128];

Pf = [
  1
 65
 33
 97
 17
 81
 49
113
  9
 73
 41
105
 25
 89
 57
121
  5
 69
 37
101
 21
 85
 53
117
 13
 77
 45
109
 29
 93
 61
125
  3
 67
 35
 99
 19
 83
 51
115
 11
 75
 43
107
 27
 91
 59
123
  7
 71
 39
103
 23
 87
 55
119
 15
 79
 47
111
 31
 95
 63
127
  2
 66
 34
 98
 18
 82
 50
114
 10
 74
 42
106
 26
 90
 58
122
  6
 70
 38
102
 22
 86
 54
118
 14
 78
 46
110
 30
 94
 62
126
  4
 68
 36
100
 20
 84
 52
116
 12
 76
 44
108
 28
 92
 60
124
  8
 72
 40
104
 24
 88
 56
120
 16
 80
 48
112
 32
 96
 64
128];

yo = y;

tmp = y(64);

for i=0:62,
    y(64+i+1) = reduce(y(i+1) * FFT128_2_64_Twiddle(i+1));
end

if f == 1,
  tmp2 = y(62);
  y(62)  = reduce(tmp2 + 1);
  y(126) = reduce((tmp2 - 1) * FFT128_2_64_Twiddle(62));
end

y(63+1)  = reduce(tmp + 1);
y(127+1) = reduce((tmp - 1) * FFT128_2_64_Twiddle(63+1));

ym = zeros(16,8);
for i = 1:16,
    ym(i,1:8) = int16(y(i*8-8+[1:8]));
end

% first -------------------------------------------------------------------

for i=0:7,
    ym(:,i+1) = FFT_8(ym(:,i+1));
end
    
for i=1:8
    for j = 1:8,
        ym(i,j) = reduce(ym(i,j)*FFT64_8_8_Twiddle(i,j));
    end
end

for i=0:7,
    ym(i+1,:) = FFT_8(ym(i+1,:));
end

% Second ------------------------------------------------------------------

for i=0:7,
    ym(9:end,i+1) = FFT_8(ym(9:end,i+1));
end
    
for i=1:8
    for j = 1:8,
        ym(8+i,j) = reduce(ym(8+i,j)*FFT64_8_8_Twiddle(i,j));
    end
end

for i=0:7,
    ym(i+9,:) = FFT_8(ym(i+9,:));
end


y2 = [ym(1,:) ym(2,:) ym(3,:) ym(4,:) ...
     ym(5,:) ym(6,:) ym(7,:) ym(8,:) ...
     ym(9,:) ym(10,:) ym(11,:) ym(12,:) ...
     ym(13,:) ym(14,:) ym(15,:) ym(16,:)];
 
% 
    y = y2(P(1:128));
%  else
%      y = y2(Pf(1:128));
%  end