6 Derived from public domain code by D. J. Bernstein.
13   unsigned int j;  in add()  local
16   for (j = 0;j < 31;++j) { u += a[j] + b[j]; out[j] = u & 255; u >>= 8; }  in add()
22   unsigned int j;  in sub()  local
25   for (j = 0;j < 31;++j) {  in sub()
26     u += a[j] + 65280 - b[j];  in sub()
27     out[j] = u & 255;  in sub()
30   u += a[31] - b[31];  in sub()
36   unsigned int j;  in squeeze()  local
39   for (j = 0;j < 31;++j) { u += a[j]; a[j] = u & 255; u >>= 8; }  in squeeze()
42   for (j = 0;j < 31;++j) { u += a[j]; a[j] = u & 255; u >>= 8; }  in squeeze()
53   unsigned int j;  in freeze()  local
56   for (j = 0;j < 32;++j) aorig[j] = a[j];  in freeze()
58   negative = -((a[31] >> 7) & 1);  in freeze()
59   for (j = 0;j < 32;++j) a[j] ^= negative & (aorig[j] ^ a[j]);  in freeze()
65   unsigned int j;  in mult()  local
70     for (j = 0;j <= i;++j) u += a[j] * b[i - j];  in mult()
71     for (j = i + 1;j < 32;++j) u += 38 * a[j] * b[i + 32 - j];  in mult()
79   unsigned int j;  in mult121665()  local
83   for (j = 0;j < 31;++j) { u += 121665 * a[j]; out[j] = u & 255; u >>= 8; }  in mult121665()
86   for (j = 0;j < 31;++j) { u += out[j]; out[j] = u & 255; u >>= 8; }  in mult121665()
87   u += out[j]; out[j] = u;  in mult121665()
93   unsigned int j;  in square()  local
98     for (j = 0;j < i - j;++j) u += a[j] * a[i - j];  in square()
99     for (j = i + 1;j < i + 32 - j;++j) u += 38 * a[j] * a[i + 32 - j];  in square()
112   unsigned int j;  in select()  local
116   bminus1 = b - 1;  in select()
117   for (j = 0;j < 64;++j) {  in select()
118     t = bminus1 & (r[j] ^ s[j]);  in select()
119     p[j] = s[j] ^ t;  in select()
120     q[j] = r[j] ^ t;  in select()
141   unsigned int j;  in mainloop()  local
145   for (j = 0;j < 32;++j) xzm1[j] = work[j];  in mainloop()
147   for (j = 33;j < 64;++j) xzm1[j] = 0;  in mainloop()
150   for (j = 1;j < 64;++j) xzm[j] = 0;  in mainloop()
152   for (pos = 254;pos >= 0;--pos) {  in mainloop()
177   for (j = 0;j < 64;++j) work[j] = xzm[j];  in mainloop()
200   /* 2^5 - 2^0 = 31 */ mult(z2_5_0,t0,z9);  in recip()
202   /* 2^6 - 2^1 */ square(t0,z2_5_0);  in recip()
203   /* 2^7 - 2^2 */ square(t1,t0);  in recip()
204   /* 2^8 - 2^3 */ square(t0,t1);  in recip()
205   /* 2^9 - 2^4 */ square(t1,t0);  in recip()
206   /* 2^10 - 2^5 */ square(t0,t1);  in recip()
207   /* 2^10 - 2^0 */ mult(z2_10_0,t0,z2_5_0);  in recip()
209   /* 2^11 - 2^1 */ square(t0,z2_10_0);  in recip()
210   /* 2^12 - 2^2 */ square(t1,t0);  in recip()
211   /* 2^20 - 2^10 */ for (i = 2;i < 10;i += 2) { square(t0,t1); square(t1,t0); }  in recip()
212   /* 2^20 - 2^0 */ mult(z2_20_0,t1,z2_10_0);  in recip()
214   /* 2^21 - 2^1 */ square(t0,z2_20_0);  in recip()
215   /* 2^22 - 2^2 */ square(t1,t0);  in recip()
216   /* 2^40 - 2^20 */ for (i = 2;i < 20;i += 2) { square(t0,t1); square(t1,t0); }  in recip()
217   /* 2^40 - 2^0 */ mult(t0,t1,z2_20_0);  in recip()
219   /* 2^41 - 2^1 */ square(t1,t0);  in recip()
220   /* 2^42 - 2^2 */ square(t0,t1);  in recip()
221   /* 2^50 - 2^10 */ for (i = 2;i < 10;i += 2) { square(t1,t0); square(t0,t1); }  in recip()
222   /* 2^50 - 2^0 */ mult(z2_50_0,t0,z2_10_0);  in recip()
224   /* 2^51 - 2^1 */ square(t0,z2_50_0);  in recip()
225   /* 2^52 - 2^2 */ square(t1,t0);  in recip()
226   /* 2^100 - 2^50 */ for (i = 2;i < 50;i += 2) { square(t0,t1); square(t1,t0); }  in recip()
227   /* 2^100 - 2^0 */ mult(z2_100_0,t1,z2_50_0);  in recip()
229   /* 2^101 - 2^1 */ square(t1,z2_100_0);  in recip()
230   /* 2^102 - 2^2 */ square(t0,t1);  in recip()
231   /* 2^200 - 2^100 */ for (i = 2;i < 100;i += 2) { square(t1,t0); square(t0,t1); }  in recip()
232   /* 2^200 - 2^0 */ mult(t1,t0,z2_100_0);  in recip()
234   /* 2^201 - 2^1 */ square(t0,t1);  in recip()
235   /* 2^202 - 2^2 */ square(t1,t0);  in recip()
236   /* 2^250 - 2^50 */ for (i = 2;i < 50;i += 2) { square(t0,t1); square(t1,t0); }  in recip()
237   /* 2^250 - 2^0 */ mult(t0,t1,z2_50_0);  in recip()
239   /* 2^251 - 2^1 */ square(t1,t0);  in recip()
240   /* 2^252 - 2^2 */ square(t0,t1);  in recip()
241   /* 2^253 - 2^3 */ square(t1,t0);  in recip()
242   /* 2^254 - 2^4 */ square(t0,t1);  in recip()
243   /* 2^255 - 2^5 */ square(t1,t0);  in recip()
244   /* 2^255 - 21 */ mult(out,t1,z11);  in recip()