diff options
Diffstat (limited to 'kara.c')
-rw-r--r-- | kara.c | 301 |
1 files changed, 301 insertions, 0 deletions
@@ -0,0 +1,301 @@ | |||
1 | #include <stdio.h> | ||
2 | #include <stdint.h> | ||
3 | #include <string.h> | ||
4 | |||
5 | typedef uint8_t leg_t; | ||
6 | typedef uint16_t dleg_t; | ||
7 | |||
8 | enum { g_legs = 256 }; | ||
9 | |||
10 | #define dump_int(A,B,C) dump_int2(A,B,C) | ||
11 | static void dump_int2( uint8_t const * p, int l, int nl ) | ||
12 | { | ||
13 | while( l-- ) | ||
14 | printf( "%02X", *(p++) ); | ||
15 | if( nl ) putchar( 10 ); | ||
16 | } | ||
17 | |||
18 | //#define printf(...) | ||
19 | //#define putchar(A) | ||
20 | |||
21 | static int mp_cmp( leg_t const *a, leg_t const *b, int legs ) | ||
22 | { | ||
23 | int leg; | ||
24 | for( leg=0; leg<legs; ++leg ) | ||
25 | { | ||
26 | if( a[leg] < b[leg] ) return -1; | ||
27 | if( a[leg] > b[leg] ) return 1; | ||
28 | } | ||
29 | return 0; | ||
30 | } | ||
31 | |||
32 | static void add( leg_t *oper, leg_t const *oper2, int legs ) | ||
33 | { | ||
34 | dleg_t acc = 0; | ||
35 | while( legs-- > 0 ) | ||
36 | { | ||
37 | acc += oper[legs] + oper2[legs]; | ||
38 | oper[legs] = (leg_t)acc; | ||
39 | acc >>= 8*sizeof(leg_t); | ||
40 | } | ||
41 | while( acc ) | ||
42 | { | ||
43 | acc += oper[legs]; | ||
44 | oper[legs--] = (leg_t)acc; | ||
45 | acc >>= 8*sizeof(leg_t); | ||
46 | } | ||
47 | } | ||
48 | |||
49 | static void adm( leg_t *oper, leg_t const *oper2, int legs, int flegs ) | ||
50 | { | ||
51 | dleg_t acc = 0; | ||
52 | while( legs-- > 0 ) | ||
53 | { | ||
54 | acc += oper[legs] + oper2[legs]; | ||
55 | oper[legs] = (leg_t)acc; | ||
56 | acc >>= 8*sizeof(leg_t); | ||
57 | } | ||
58 | while( acc && flegs-- ) | ||
59 | { | ||
60 | acc += oper[legs]; | ||
61 | oper[legs--] = (leg_t)acc; | ||
62 | acc >>= 8*sizeof(leg_t); | ||
63 | } | ||
64 | } | ||
65 | |||
66 | static void sub( leg_t *oper, leg_t const *oper2, int legs ) | ||
67 | { | ||
68 | int borrow = 0, borrow_temp; | ||
69 | while( legs-- ) | ||
70 | { | ||
71 | leg_t temp = oper[legs] - oper2[legs]; | ||
72 | borrow_temp = temp > oper[legs]; | ||
73 | oper[legs] = temp - borrow; | ||
74 | borrow = borrow_temp | ( oper[legs] > temp ); | ||
75 | } | ||
76 | while( borrow ) | ||
77 | { | ||
78 | leg_t temp = oper[legs] - borrow; | ||
79 | borrow = temp > oper[legs]; | ||
80 | oper[legs--] = temp; | ||
81 | } | ||
82 | printf( "sub: %d\n", legs ); | ||
83 | } | ||
84 | |||
85 | /* result is guaranteed to be initialized with enough legs prepended to take the carry */ | ||
86 | static void mp_mul_uint_add( leg_t *result, leg_t const *oper, leg_t fac, int legs ) | ||
87 | { | ||
88 | dleg_t acc = 0; | ||
89 | while( legs-- ) | ||
90 | { | ||
91 | acc += (dleg_t)result[legs] + (dleg_t)oper[legs] * (dleg_t)fac; | ||
92 | result[legs] = (leg_t)acc; | ||
93 | acc >>= 8*sizeof(leg_t); | ||
94 | } | ||
95 | while( acc ) | ||
96 | { | ||
97 | acc += result[legs]; | ||
98 | result[legs--] = (leg_t)acc; | ||
99 | acc >>= 8*sizeof(leg_t); | ||
100 | } | ||
101 | } | ||
102 | |||
103 | static void mul( leg_t * result, leg_t const *oper, leg_t const *oper2, int legs ) | ||
104 | { | ||
105 | int leg = legs; | ||
106 | while( leg-- ) | ||
107 | mp_mul_uint_add( result + 1 + leg, oper, oper2[leg], legs ); | ||
108 | } | ||
109 | |||
110 | static int sub_mod( leg_t * result, leg_t const *a, leg_t const * b, int legs ) | ||
111 | { | ||
112 | int borrow = 0, borrow_temp; | ||
113 | while( legs-- ) | ||
114 | { | ||
115 | leg_t temp = a[legs] - b[legs]; | ||
116 | borrow_temp = temp > a[legs]; | ||
117 | result[legs] = temp - borrow; | ||
118 | borrow = borrow_temp | (result[legs] > temp ); | ||
119 | } | ||
120 | return borrow; | ||
121 | } | ||
122 | |||
123 | static void negate( leg_t * p, int legs ) | ||
124 | { | ||
125 | while( legs && !p[legs-1] ) p[legs--]^=(leg_t)-1; | ||
126 | if( !legs ) return; | ||
127 | p[legs-1]--; | ||
128 | while( legs--) p[legs]^=(leg_t)-1; | ||
129 | } | ||
130 | |||
131 | void kara2( leg_t* p, leg_t *a, leg_t *b, int len ) | ||
132 | { | ||
133 | printf( "L%d karaloop: a*b ", len ); dump_int( a, len, 0 ); putchar('*'); dump_int( b, len, 1 ); | ||
134 | memset( p, 0, 2 * len * sizeof( leg_t )); | ||
135 | if( len == 2 ) | ||
136 | { | ||
137 | mul( p, a, b, len ); | ||
138 | printf( "L%d p:", len ); dump_int( p, len*2, 1 ); | ||
139 | } | ||
140 | else | ||
141 | { | ||
142 | leg_t scratch[len]; | ||
143 | int sign = 0, n = len / 2; | ||
144 | |||
145 | if( mp_cmp( a, a + n, n ) > 0 ) | ||
146 | sub_mod( scratch, a, a + n, n ); | ||
147 | else | ||
148 | { | ||
149 | sub_mod( scratch, a + n, a, n ); | ||
150 | sign = 1; | ||
151 | } | ||
152 | |||
153 | if( mp_cmp( b, b + n, n ) > 0 ) | ||
154 | sub_mod( scratch + n, b, b + n, n ); | ||
155 | else | ||
156 | { | ||
157 | sub_mod( scratch + n, b + n, b, n ); | ||
158 | sign ^= 1; | ||
159 | } | ||
160 | |||
161 | printf( "L%d (x0-x1)*(y1-y0)=", len ); dump_int( scratch, n, 0 ); putchar('*'); dump_int( scratch + n, n, 0 ); printf( "(sign = %d)\n", sign ); | ||
162 | |||
163 | kara2( p + n, scratch, scratch + n, n ); | ||
164 | if( !sign ) | ||
165 | negate( p, len + n ); | ||
166 | |||
167 | printf( "L%d psub:", len ); dump_int( p, len*2, 1 ); | ||
168 | |||
169 | kara2( scratch, a + n, b + n, n ); | ||
170 | adm( p + len, scratch, len, len ); | ||
171 | adm( p + n, scratch, len, n ); | ||
172 | |||
173 | kara2( scratch, a, b, n ); | ||
174 | adm( p + n, scratch, len, n ); | ||
175 | adm( p, scratch, len, 0 ); | ||
176 | printf( "L%d p:", len ); dump_int( p, len*2, 1 ); | ||
177 | } | ||
178 | } | ||
179 | |||
180 | void kara( uint8_t* p, uint8_t *a, uint8_t *b, int len ) | ||
181 | { | ||
182 | int i; | ||
183 | memset( p, 0, 2 * len * sizeof(leg_t ) ); | ||
184 | if( len <= 3 ) | ||
185 | { | ||
186 | printf( "L%d karaend: a*b ", len ); dump_int( a, len, 0 ); putchar('*'); dump_int( b, len, 1 ); | ||
187 | mul( p, a, b, len ); | ||
188 | printf( "L%d result:\n", len ); dump_int( p, len*2, 1 ); | ||
189 | } | ||
190 | else { | ||
191 | int n = len / 2; | ||
192 | int l = ( len + 1 ) / 2; | ||
193 | int o = l - n; | ||
194 | int sh = l + 1; | ||
195 | uint8_t scratch[2*sh]; | ||
196 | uint8_t *midpoint = p + len * 2 - n; | ||
197 | |||
198 | printf( "L%d karaloop: a*b ", len ); dump_int( a, len, 0 ); putchar('*'); dump_int( b, len, 1 ); | ||
199 | printf( "Values n: %d l: %d o: %d sh: %d\n", n, l, o, sh ); | ||
200 | |||
201 | memset( scratch, 0, sizeof(scratch)); | ||
202 | memcpy( scratch + 1 , a, l ); | ||
203 | memcpy( scratch + 1 + sh , b, l ); | ||
204 | |||
205 | /* P3 = (a1+a2)*(b1+b2) */ | ||
206 | add( scratch + 1 + o, a + l, n ); | ||
207 | add( scratch + sh + 1 + o, b + l, n ); | ||
208 | |||
209 | printf( "L%d a1+a2: \n", len ); dump_int( scratch, sh, 1 ); | ||
210 | printf( "L%d b1+b2: \n", len ); dump_int( scratch+sh, sh, 1 ); | ||
211 | |||
212 | kara( midpoint - 2 * sh, scratch, scratch + sh, sh ); | ||
213 | |||
214 | printf( "L%d (a1+a2)*(b1+b2)<<n: \n", len ); dump_int( p, len*2, 1 ); | ||
215 | |||
216 | /* P2 = a2*b2 */ | ||
217 | kara( scratch, a + l, b + l, n ); | ||
218 | |||
219 | printf( "L%d a2*b2: \n", len ); dump_int( scratch, n*2, 1 ); | ||
220 | |||
221 | /* P = ( P3 - P2 ) * 2^n + P2 ... and later: + P1 * 2^(2*n) - P1 * 2^n*/ | ||
222 | add( midpoint - n, scratch, n * 2 ); | ||
223 | |||
224 | printf( "L%d p + p2: \n", len ); dump_int( p, len*2, 1 ); | ||
225 | |||
226 | sub( midpoint - n * 2, scratch, n * 2 ); | ||
227 | |||
228 | printf( "L%d p - p2<<n: \n", len ); dump_int( p, len*2, 1 ); | ||
229 | |||
230 | /* now + P1 * 2^(2*n) - P1 * 2^n*/ | ||
231 | kara( scratch, a, b, l ); | ||
232 | |||
233 | printf( "L%d a1*b1: \n", len ); dump_int( scratch, l*2, 1 ); | ||
234 | |||
235 | sub( midpoint - l * 2, scratch, l * 2 ); | ||
236 | |||
237 | printf( "L%d p - p1<<n: \n", len ); dump_int( p, len*2, 1 ); | ||
238 | |||
239 | add( midpoint - n - l * 2, scratch, l * 2 ); | ||
240 | |||
241 | printf( "L%d p + p1<<n^2: \n", len ); dump_int( p, len*2, 1 ); | ||
242 | } | ||
243 | } | ||
244 | |||
245 | #undef putchar | ||
246 | int main() { | ||
247 | uint8_t a[] = { 0x99,0xf3,0x0a,0x3d,0x8a,0xdd,0x2d,0x1f }; | ||
248 | uint8_t b[] = { 0xf1,0x48,0x71,0x4b,0xf7,0xb1,0x9a,0x30 }; | ||
249 | |||
250 | // uint8_t a[g_legs], b[g_legs]; | ||
251 | // uint8_t p[2*g_legs]; | ||
252 | |||
253 | /* | ||
254 | static const leg_t b[] = { | ||
255 | 0x3D, 0xA9, 0x76, 0x59, 0xE2, 0x80, 0xDB, 0x0B, 0xE6, 0x5B, 0xCC, 0x3A, 0xB7, 0x8F, 0xDA, 0xA9, 0xB7, 0xB7, 0x68, 0xC8, 0x99, 0x31, 0xD7, 0x8D, 0xF8, 0xB1, 0x17, 0x25, 0x33, 0x9E, 0xBC, 0x93, | ||
256 | 0xAA, 0x7F, 0xBD, 0x95, 0x62, 0x05, 0x9F, 0x1F, 0x4C, 0x2D, 0xE6, 0x7D, 0xAD, 0x47, 0x52, 0x7E, 0x52, 0x6A, 0x65, 0x3A, 0x7A, 0x67, 0x4B, 0xD5, 0x54, 0x01, 0xEA, 0x4F, 0x3E, 0xD7, 0x3A, 0x2F, | ||
257 | 0x70, 0xB5, 0x6F, 0x52, 0x7F, 0x6F, 0x60, 0x4F, 0xF3, 0xE5, 0x6C, 0xC2, 0xBD, 0x9F, 0x04, 0x8C, 0xFE, 0xA8, 0x0D, 0x9A, 0x6E, 0xC9, 0xFC, 0xD3, 0x3A, 0xD3, 0x6F, 0xD8, 0x22, 0xC3, 0x9F, 0x34, | ||
258 | 0x91, 0xF3, 0x0C, 0x52, 0xF7, 0x98, 0xDA, 0x6A, 0x18, 0xC3, 0xDC, 0xE2, 0x01, 0x88, 0xD8, 0x4C, 0xCB, 0x22, 0x51, 0x76, 0x25, 0x9E, 0x08, 0x0F, 0xA8, 0x9D, 0x1D, 0xCD, 0x9D, 0x38, 0x1C, 0xC5, | ||
259 | 0xBE, 0xAC, 0xD4, 0x6F, 0x3C, 0xDD, 0x11, 0x96, 0x91, 0xA4, 0xF5, 0x57, 0x29, 0x29, 0xB9, 0x0C, 0x67, 0xDE, 0x8F, 0xA0, 0x23, 0x86, 0x47, 0x14, 0xC2, 0x8A, 0x61, 0xD4, 0x74, 0x11, 0x40, 0x2D, | ||
260 | 0x6C, 0x09, 0x06, 0x0D, 0x41, 0x05, 0x86, 0x39, 0xE4, 0xFC, 0xCF, 0x1D, 0x63, 0x8F, 0x45, 0x66, 0x9E, 0x10, 0xFD, 0xE2, 0x8E, 0x54, 0x80, 0x6B, 0xF1, 0xD2, 0x7D, 0x0B, 0x5C, 0x7D, 0xC9, 0xC2, | ||
261 | 0xB6, 0x16, 0xD6, 0xFA, 0x8B, 0xE2, 0xC9, 0x1D, 0xE8, 0x10, 0x54, 0x64, 0xE9, 0xF8, 0x0A, 0x5F, 0x34, 0x72, 0x08, 0x69, 0x90, 0xDA, 0xCF, 0x1A, 0x7E, 0x2C, 0x75, 0xA5, 0x8E, 0x25, 0xF1, 0x42, | ||
262 | 0x8F, 0xB4, 0x83, 0x2E, 0xF8, 0x27, 0xDE, 0x84, 0xCA, 0x06, 0xDA, 0x91, 0xC2, 0xB3, 0xE7, 0xE2, 0x4F, 0x02, 0x41, 0x93, 0x78, 0x7A, 0x82, 0x78, 0x8B, 0xD7, 0x05, 0x62, 0xDA, 0x60, 0xE3, 0x92, | ||
263 | 0x46, 0xBD, 0xB7, 0x33, 0x6E, 0x84, 0x52, 0xD9, 0x1D, 0x7D, 0x37, 0xA2, 0x3F, 0xB8, 0xCF, 0x61, 0xB1, 0x8A, 0x9E, 0xF1, 0x50, 0xC8, 0x95, 0x3A, 0xC9, 0x39, 0x19, 0xD1, 0x2A, 0x4B, 0x1A, 0x67, | ||
264 | 0xFD, 0xC6, 0x5A, 0x26, 0x9B, 0x51, 0xC1, 0xEF, 0x09, 0x95, 0x48, 0x43, 0x20, 0xE7, 0x39, 0xF4, 0x6C, 0x79, 0x51, 0xA5, 0x23, 0xCE, 0xF7, 0x85, 0xE6, 0x0C, 0x6E, 0xFD, 0xFF, 0xB7, 0xA9, 0xA9, | ||
265 | 0x5F, 0x66, 0x61, 0x46, 0xCB, 0x44, 0x1F, 0x59, 0xAE, 0x01, 0xE0, 0xF3, 0x63, 0xA9, 0x31, 0x5D, 0x04, 0xBA, 0x04, 0x4A, 0xEB, 0x4E, 0xEF, 0xD4, 0x85, 0x63, 0x21, 0x5F, 0x72, 0xC8, 0xD9, 0x89, | ||
266 | 0x62, 0xD2, 0x18, 0x77, 0x12, 0x96, 0xEF, 0x6A, 0x20, 0xBD, 0x72, 0xB9, 0xB2, 0x1E, 0x6B, 0x3D, 0x53, 0xC4, 0x4F, 0xAB, 0x73, 0x48, 0x10, 0xF7, 0xD2, 0x03, 0xA9, 0xE0, 0xD7, 0xCE, 0x25, 0xD0, | ||
267 | 0x2E, 0x52, 0x98, 0x9E, 0xCC, 0xF8, 0x5F, 0x34, 0x91, 0x2A, 0x04, 0x91, 0x3E, 0x9E, 0xBD, 0x87, 0x82, 0x67, 0x53, 0x7D, 0x4A, 0x61, 0x2A, 0x18, 0x51, 0xE7, 0x5D, 0x99, 0x98, 0xF0, 0x01, 0xDB, | ||
268 | 0xC9, 0xC7, 0x7F, 0x0C, 0x35, 0x2D, 0x40, 0x8C, 0xA7, 0x96, 0xD1, 0x82, 0x04, 0xA6, 0x36, 0xF7, 0xE4, 0x40, 0x40, 0x92, 0x1C, 0x1E, 0x46, 0x7C, 0x52, 0x4E, 0x7C, 0x7A, 0x7E, 0xD3, 0x6C, 0x41, | ||
269 | 0x2A, 0x43, 0x4C, 0xEB, 0x23, 0x0B, 0x2D, 0xFE, 0x35, 0x49, 0xC5, 0x77, 0x7A, 0x17, 0xFB, 0x04, 0xB8, 0x50, 0xDE, 0x95, 0xD9, 0x7A, 0xC4, 0x0A, 0x55, 0xEA, 0x6F, 0x75, 0x41, 0xC4, 0xF8, 0x2B, | ||
270 | 0x37, 0xBF, 0x90, 0xFE, 0x52, 0x07, 0x4F, 0x19, 0xFA, 0x8F, 0x75, 0xF0, 0x06, 0x7E, 0x82, 0xB1, 0x8A, 0x1A, 0xC0, 0x24, 0xB3, 0x0E, 0x9B, 0x12, 0xC1, 0x4A, 0xB0, 0xDD, 0xCC, 0x03, 0xAA, 0x20 | ||
271 | }; | ||
272 | uint8_t a[sizeof(b)]; | ||
273 | */ | ||
274 | |||
275 | uint8_t p[2*sizeof(b)]; | ||
276 | int i; | ||
277 | |||
278 | srandomdev(); | ||
279 | |||
280 | // while( 1 ) | ||
281 | { | ||
282 | // memset( a, 0, sizeof(a) ); | ||
283 | // a[sizeof(b)-1]=2; | ||
284 | /* | ||
285 | memset( b, 0, sizeof(b) ); | ||
286 | a[g_legs-1] = random() & 0xff; | ||
287 | for( i=g_legs/2; i<sizeof(a); ++i ) | ||
288 | { | ||
289 | b[i] = random() & 0xff; | ||
290 | } | ||
291 | */ | ||
292 | kara2( p, a, b, sizeof(b)); | ||
293 | |||
294 | dump_int2( a, sizeof(a), 0 ); | ||
295 | putchar( '*' ); | ||
296 | dump_int2( b, sizeof(b), 0 ); | ||
297 | putchar( '-' ); | ||
298 | dump_int2( p, sizeof(p), 1 ); | ||
299 | } | ||
300 | return 0; | ||
301 | } | ||