bigint/bigint.c

   1 /* bigint.c */
   2 /*
   3     This file is part of the ARM-Crypto-Lib.
   4     Copyright (C) 2006-2015 Daniel Otte (bg@nerilex.org)
   5
   6     This program is free software: you can redistribute it and/or modify
   7     it under the terms of the GNU General Public License as published by
   8     the Free Software Foundation, either version 3 of the License, or
   9     (at your option) any later version.
  10
  11     This program is distributed in the hope that it will be useful,
  12     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14     GNU General Public License for more details.
  15
  16     You should have received a copy of the GNU General Public License
  17     along with this program.  If not, see <http://www.gnu.org/licenses/>.
  18 */
  19 /**
  20  * \file                bigint.c
  21  * \author              Daniel Otte
  22  * \date                2010-02-22
  23  *
  24  * \license         GPLv3 or later
  25  *
  26  */
  27
  28
  29 #define STRING2(x) #x
  30 #define STRING(x) STRING2(x)
  31 #define STR_LINE STRING(__LINE__)
  32
  33 #include "bigint.h"
  34 #include <string.h>
  35 #include <stdio.h>
  36
  37 #define PREFERE_HEAP_SPACE 1
  38
  39 #if PREFERE_HEAP_SPACE
  40 #include <stdlib.h>
  41
  42 #define ALLOC_BIGINT_WORDS(var,words)                                      \
  43     bigint_word_t *(var) = malloc((words) * sizeof(bigint_word_t));        \
  44     if (!(var)) {                                                          \
  45         puts_P(PSTR("\n\nDBG: OOM ERROR (in arithmeics)!\n"));             \
  46         uart0_flush();                                                     \
  47         for(;;)                                                            \
  48             ;                                                              \
  49     }
  50
  51 #define FREE(x) free(x)
  52
  53 #else
  54
  55 #define ALLOC_BIGINT_WORDS(var,words) bigint_word_t var[words]
  56 #define FREE(x)
  57
  58 #endif
  59
  60 #define DEBUG 1
  61
  62 #if DEBUG
  63 #include "cli.h"
  64 #include "uart_i.h"
  65 #include "bigint_io.h"
  66 #endif
  67
  68 #ifndef MAX
  69  #define MAX(a,b) (((a)>(b))?(a):(b))
  70 #endif
  71
  72 #ifndef MIN
  73  #define MIN(a,b) (((a) < (b)) ? (a) : (b))
  74 #endif
  75
  76 #define SET_FBS(a, v) do {(a)->info &= ~BIGINT_FBS_MASK; (a)->info |= (v);} while(0)
  77 #define GET_FBS(a)   ((a)->info & BIGINT_FBS_MASK)
  78 #define SET_NEG(a)   (a)->info |= BIGINT_NEG_MASK
  79 #define SET_POS(a)   (a)->info &= ~BIGINT_NEG_MASK
  80 #define XCHG(a,b)    do{(a) ^= (b); (b) ^= (a); (a) ^= (b);} while(0)
  81 #define XCHG_PTR(a,b)    do{ a = (void*)(((intptr_t)(a)) ^ ((intptr_t)(b))); \
  82                                  b = (void*)(((intptr_t)(a)) ^ ((intptr_t)(b))); \
  83                                  a = (void*)(((intptr_t)(a)) ^ ((intptr_t)(b)));} while(0)
  84
  85 #define GET_SIGN(a) ((a)->info & BIGINT_NEG_MASK)
  86
  87 /******************************************************************************/
  88 void bigint_adjust(bigint_t *a){
  89         while (a->length_W != 0 && a->wordv[a->length_W - 1] == 0) {
  90                 a->length_W--;
  91         }
  92         if (a->length_W == 0) {
  93                 a->info=0;
  94                 return;
  95         }
  96         bigint_word_t t;
  97         uint8_t i = BIGINT_WORD_SIZE - 1;
  98         t = a->wordv[a->length_W - 1];
  99         while ((t & (((bigint_word_t)1) << (BIGINT_WORD_SIZE - 1))) == 0 && i) {
 100                 t <<= 1;
 101                 i--;
 102         }
 103         SET_FBS(a, i);
 104 }
 105
 106 /******************************************************************************/
 107
 108 bigint_length_t bigint_length_b(const bigint_t *a){
 109         if(!a->length_W || a->length_W==0){
 110                 return 0;
 111         }
 112         return (a->length_W-1) * BIGINT_WORD_SIZE + GET_FBS(a);
 113 }
 114
 115 /******************************************************************************/
 116
 117 bigint_length_t bigint_length_B(const bigint_t *a){
 118         return a->length_W * sizeof(bigint_word_t);
 119 }
 120
 121 /******************************************************************************/
 122
 123 int32_t bigint_get_first_set_bit(const bigint_t *a){
 124         if(a->length_W == 0) {
 125                 return -1;
 126         }
 127         return (a->length_W-1) * sizeof(bigint_word_t) * CHAR_BIT + GET_FBS(a);
 128 }
 129
 130
 131 /******************************************************************************/
 132
 133 int32_t bigint_get_last_set_bit(const bigint_t *a){
 134         uint32_t r = 0;
 135         uint8_t b = 0;
 136         bigint_word_t x = 1;
 137         if (a->length_W == 0) {
 138                 return -1;
 139         }
 140         while (a->wordv[r] == 0 && r < a->length_W) {
 141                 ++r;
 142         }
 143         if (a->wordv[r] == 0) {
 144                 return (uint32_t)(-1);
 145         }
 146         while ((x&a->wordv[r])==0) {
 147                 ++b;
 148                 x <<= 1;
 149         }
 150         return r * BIGINT_WORD_SIZE + b;
 151 }
 152
 153 /******************************************************************************/
 154
 155 void bigint_copy(bigint_t *dest, const bigint_t *src){
 156     if(dest->wordv != src->wordv){
 157             memcpy(dest->wordv, src->wordv, src->length_W * sizeof(bigint_word_t));
 158     }
 159     dest->length_W = src->length_W;
 160         dest->info = src->info;
 161 }
 162
 163 /******************************************************************************/
 164
 165 /* this should be implemented in assembly */
 166 void bigint_add_u(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 167         bigint_length_t i;
 168         bigint_wordplus_t t = 0LL;
 169         if (a->length_W < b->length_W) {
 170                 XCHG_PTR(a, b);
 171         }
 172         for(i = 0; i < b->length_W; ++i) {
 173                 t += a->wordv[i];
 174                 t += b->wordv[i];
 175                 dest->wordv[i] = (bigint_word_t)t;
 176                 t >>= BIGINT_WORD_SIZE;
 177         }
 178         for(; i < a->length_W; ++i){
 179                 t += a->wordv[i];
 180                 dest->wordv[i] = (bigint_word_t)t;
 181                 t >>= BIGINT_WORD_SIZE;
 182         }
 183         if(t){
 184                 dest->wordv[i++] = (bigint_word_t)t;
 185         }
 186         dest->length_W = i;
 187         SET_POS(dest);
 188         bigint_adjust(dest);
 189 }
 190
 191 /******************************************************************************/
 192
 193 /* this should be implemented in assembly */
 194 void bigint_add_scale_u(bigint_t *dest, const bigint_t *a, bigint_length_t scale){
 195         if(a->length_W == 0){
 196                 return;
 197         }
 198         if(scale == 0){
 199                 bigint_add_u(dest, dest, a);
 200                 return;
 201         }
 202         bigint_t x;
 203 #if BIGINT_WORD_SIZE == 8
 204         memset(dest->wordv + dest->length_W, 0, MAX(dest->length_W, a->length_W + scale) - dest->length_W);
 205         x.wordv = dest->wordv + scale;
 206         x.length_W = dest->length_W - scale;
 207         if((int16_t)x.length_W < 0){
 208                 x.length_W = 0;
 209                 x.info = 0;
 210         } else {
 211                 x.info = dest->info;
 212         }
 213         bigint_add_u(&x, &x, a);
 214         dest->length_W = x.length_W + scale;
 215         dest->info = 0;
 216         bigint_adjust(dest);
 217 #else
 218         bigint_t s;
 219         bigint_length_t word_shift = scale / sizeof(bigint_word_t), byte_shift = scale % sizeof(bigint_word_t);
 220         bigint_word_t bv[a->length_W + 1];
 221         s.wordv = bv;
 222         bv[0] = bv[a->length_W] = 0;
 223         memcpy((uint8_t*)bv + byte_shift, a->wordv, a->length_W * sizeof(bigint_word_t));
 224         s.length_W = a->length_W + 1;
 225         bigint_adjust(&s);
 226         memset(dest->wordv + dest->length_W, 0, (MAX(dest->length_W, s.length_W + word_shift) - dest->length_W) * sizeof(bigint_word_t));
 227         x.wordv = dest->wordv + word_shift;
 228         x.length_W = dest->length_W - word_shift;
 229         if((int16_t)x.length_W < 0){
 230                 x.length_W = 0;
 231                 x.info = 0;
 232         }else{
 233                 x.info = dest->info;
 234         }
 235         bigint_add_u(&x, &x, &s);
 236         dest->length_W = x.length_W + word_shift;
 237         dest->info = 0;
 238         bigint_adjust(dest);
 239 #endif
 240 }
 241
 242 /******************************************************************************/
 243
 244 /* this should be implemented in assembly */
 245 void bigint_sub_u(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 246         int8_t borrow = 0;
 247         int8_t  r;
 248         bigint_wordplus_signed_t t = 0;
 249         bigint_length_t i;
 250         if(b->length_W == 0){
 251                 bigint_copy(dest, a);
 252                 SET_POS(dest);
 253                 return;
 254         }
 255         if(a->length_W == 0){
 256                 bigint_copy(dest, b);
 257                 SET_NEG(dest);
 258                 return;
 259         }
 260     r = bigint_cmp_u(a,b);
 261     if(r == 0){
 262         bigint_set_zero(dest);
 263         return;
 264     }
 265         if(r < 0){
 266                 bigint_sub_u(dest, b, a);
 267                 SET_NEG(dest);
 268                 return;
 269         }
 270         for(i = 0; i < a->length_W; ++i){
 271                 t = a->wordv[i];
 272                 if(i < b->length_W){
 273                         t -= b->wordv[i];
 274                 }
 275                 t -= borrow;
 276                 dest->wordv[i] = (bigint_word_t)t;
 277                 borrow = t < 0 ? 1 : 0;
 278         }
 279         SET_POS(dest);
 280         dest->length_W = i;
 281         bigint_adjust(dest);
 282 }
 283
 284 /******************************************************************************/
 285
 286 int8_t bigint_cmp_u(const bigint_t *a, const bigint_t *b){
 287         if(a->length_W > b->length_W){
 288                 return 1;
 289         }
 290         if(a->length_W < b->length_W){
 291                 return -1;
 292         }
 293         if(a->length_W == 0){
 294                 return 0;
 295         }
 296         bigint_length_t i;
 297         i = a->length_W - 1;
 298         do{
 299                 if(a->wordv[i] != b->wordv[i]){
 300                         if(a->wordv[i] > b->wordv[i]){
 301                                 return 1;
 302                         }else{
 303                                 return -1;
 304                         }
 305                 }
 306         }while(i--);
 307         return 0;
 308 }
 309
 310 /******************************************************************************/
 311
 312 void bigint_add_s(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 313         uint8_t s;
 314         s  = GET_SIGN(a)?2:0;
 315         s |= GET_SIGN(b)?1:0;
 316         switch(s){
 317                 case 0: /* both positive */
 318                         bigint_add_u(dest, a,b);
 319                         SET_POS(dest);
 320                         break;
 321                 case 1: /* a positive, b negative */
 322                         bigint_sub_u(dest, a, b);
 323                         break;
 324                 case 2: /* a negative, b positive */
 325                         bigint_sub_u(dest, b, a);
 326                         break;
 327                 case 3: /* both negative */
 328                         bigint_add_u(dest, a, b);
 329                         SET_NEG(dest);
 330                         break;
 331                 default: /* how can this happen?*/
 332                         break;
 333         }
 334 }
 335
 336 /******************************************************************************/
 337
 338 void bigint_sub_s(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 339         uint8_t s;
 340         s  = GET_SIGN(a)?2:0;
 341         s |= GET_SIGN(b)?1:0;
 342         switch(s){
 343                 case 0: /* both positive */
 344                         bigint_sub_u(dest, a,b);
 345                         break;
 346                 case 1: /* a positive, b negative */
 347                         bigint_add_u(dest, a, b);
 348                         SET_POS(dest);
 349                         break;
 350                 case 2: /* a negative, b positive */
 351                         bigint_add_u(dest, a, b);
 352                         SET_NEG(dest);
 353                         break;
 354                 case 3: /* both negative */
 355                         bigint_sub_u(dest, b, a);
 356                         break;
 357                 default: /* how can this happen?*/
 358                                         break;
 359         }
 360
 361 }
 362
 363 /******************************************************************************/
 364
 365 int8_t bigint_cmp_s(const bigint_t *a, const bigint_t *b){
 366         uint8_t s;
 367         if(a->length_W==0 && b->length_W==0){
 368                 return 0;
 369         }
 370         s  = GET_SIGN(a)?2:0;
 371         s |= GET_SIGN(b)?1:0;
 372         switch(s){
 373                 case 0: /* both positive */
 374                         return bigint_cmp_u(a, b);
 375                         break;
 376                 case 1: /* a positive, b negative */
 377                         return 1;
 378                         break;
 379                 case 2: /* a negative, b positive */
 380                         return -1;
 381                         break;
 382                 case 3: /* both negative */
 383                         return bigint_cmp_u(b, a);
 384                         break;
 385                 default: /* how can this happen?*/
 386                                         break;
 387         }
 388         return 0; /* just to satisfy the compiler */
 389 }
 390
 391 /******************************************************************************/
 392
 393 void bigint_shiftleft_bits(bigint_t *a, uint8_t shift) {
 394     bigint_length_t i;
 395     bigint_wordplus_t t = 0;
 396
 397     for (i = 0; i < a->length_W; ++i) {
 398         t |= ((bigint_wordplus_t)a->wordv[i]) << shift;
 399         a->wordv[i] = (bigint_word_t)t;
 400         t >>= BIGINT_WORD_SIZE;
 401     }
 402     if (t) {
 403         a->wordv[a->length_W++] = (bigint_word_t)t;
 404     }
 405     bigint_adjust(a);
 406 }
 407
 408 /******************************************************************************/
 409
 410 void bigint_shiftleft(bigint_t *a, bigint_length_t shift){
 411         bigint_length_t byteshift;
 412         uint8_t bitshift;
 413
 414         if (a->length_W == 0 || shift == 0) {
 415                 return;
 416         }
 417         byteshift = shift / 8;
 418         bitshift = shift & 7;
 419     if (byteshift % sizeof(bigint_word_t)) {
 420         a->wordv[a->length_W + byteshift / sizeof(bigint_t)] = 0;
 421     }
 422         if (byteshift) {
 423                 memmove(((uint8_t*)a->wordv) + byteshift, a->wordv, a->length_W * sizeof(bigint_word_t));
 424                 memset(a->wordv, 0, byteshift);
 425         a->length_W += (byteshift + sizeof(bigint_word_t) - 1) / sizeof(bigint_word_t);
 426         }
 427
 428         if (bitshift == 0) {
 429             bigint_adjust(a);
 430         } else {
 431             bigint_shiftleft_bits(a, bitshift);
 432         }
 433 }
 434
 435 /******************************************************************************/
 436
 437 void bigint_shiftright_bits(bigint_t *a, uint8_t shift){
 438     bigint_length_t i;
 439     bigint_wordplus_t t = 0;
 440
 441     i = a->length_W;
 442     while (i--) {
 443         t |= ((bigint_wordplus_t)(a->wordv[i])) << (BIGINT_WORD_SIZE - shift);
 444         a->wordv[i] = (bigint_word_t)(t >> BIGINT_WORD_SIZE);
 445         t <<= BIGINT_WORD_SIZE;
 446     }
 447
 448     bigint_adjust(a);
 449 }
 450
 451 /******************************************************************************/
 452
 453 void bigint_shiftright_1bit(bigint_t *a){
 454     bigint_length_t i;
 455     bigint_word_t t1 = 0, t2;
 456
 457     i = a->length_W;
 458     while (i--) {
 459         t2 = a->wordv[i] & 1;
 460         a->wordv[i] >>= 1;
 461         a->wordv[i] |= t1;
 462         t1 = t2 << (BIGINT_WORD_SIZE - 1);
 463     }
 464     bigint_adjust(a);
 465 }
 466
 467 /******************************************************************************/
 468
 469 void bigint_shiftright_1word(bigint_t *a){
 470     if (a->length_W == 0) {
 471         return;
 472     }
 473     memmove(a->wordv, &a->wordv[1], (a->length_W - 1) * sizeof(bigint_word_t));
 474     a->length_W--;
 475 }
 476
 477 /******************************************************************************/
 478
 479 void bigint_shiftright(bigint_t *a, bigint_length_t shift){
 480         bigint_length_t byteshift = shift / 8;
 481         uint8_t bitshift = shift & 7;
 482
 483         if (a->length_W == 0) {
 484             return;
 485         }
 486
 487         if(bigint_get_first_set_bit(a) < shift){ /* we would shift out more than we have */
 488                 bigint_set_zero(a);
 489                 return;
 490         }
 491
 492         if(byteshift){
 493                 memmove(a->wordv, (uint8_t*)a->wordv + byteshift, a->length_W * sizeof(bigint_word_t) - byteshift);
 494                 memset((uint8_t*)&a->wordv[a->length_W] - byteshift, 0, byteshift);
 495             a->length_W -= byteshift / sizeof(bigint_word_t);
 496             bigint_adjust(a);
 497         }
 498
 499     if(bitshift != 0){
 500         bigint_shiftright_bits(a, bitshift);
 501         }
 502 }
 503
 504 /******************************************************************************/
 505
 506 void bigint_xor(bigint_t *dest, const bigint_t *a){
 507         bigint_length_t i;
 508         for(i=0; i<a->length_W; ++i){
 509                 dest->wordv[i] ^= a->wordv[i];
 510         }
 511         bigint_adjust(dest);
 512 }
 513
 514 /******************************************************************************/
 515
 516 void bigint_set_zero(bigint_t *a){
 517         a->length_W = 0;
 518 }
 519
 520 /******************************************************************************/
 521
 522 /* using the Karatsuba-Algorithm */
 523 /* x*y = (xh*yh)*b**2n + ((xh+xl)*(yh+yl) - xh*yh - xl*yl)*b**n + yh*yl */
 524 void bigint_mul_u(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 525         if (a->length_W == 0 || b->length_W == 0) {
 526                 bigint_set_zero(dest);
 527                 return;
 528         }
 529         if (dest == a || dest == b) {
 530                 bigint_t d;
 531                 bigint_word_t d_b[a->length_W + b->length_W];
 532                 d.wordv = d_b;
 533                 bigint_mul_u(&d, a, b);
 534                 bigint_copy(dest, &d);
 535                 return;
 536         }
 537         if (a->length_W == 1 || b->length_W == 1) {
 538                 if (a->length_W != 1) {
 539                         XCHG_PTR(a,b);
 540                 }
 541                 bigint_wordplus_t t = 0;
 542                 bigint_length_t i;
 543                 bigint_word_t x = a->wordv[0];
 544                 for (i = 0; i < b->length_W; ++i) {
 545                         t += ((bigint_wordplus_t)b->wordv[i]) * ((bigint_wordplus_t)x);
 546                         dest->wordv[i] = (bigint_word_t)t;
 547                         t >>= BIGINT_WORD_SIZE;
 548                 }
 549                 dest->length_W = i;
 550                 if (t) {
 551                     dest->wordv[i] = (bigint_word_t)t;
 552                     dest->length_W += 1;
 553                 }
 554                 dest->info = 0;
 555                 bigint_adjust(dest);
 556                 return;
 557         }
 558         if (a->length_W * sizeof(bigint_word_t) <= 4 && b->length_W * sizeof(bigint_word_t) <= 4) {
 559                 uint32_t p = 0, q = 0;
 560                 uint64_t r;
 561                 memcpy(&p, a->wordv, a->length_W * sizeof(bigint_word_t));
 562                 memcpy(&q, b->wordv, b->length_W * sizeof(bigint_word_t));
 563                 r = (uint64_t)p * (uint64_t)q;
 564                 memcpy(dest->wordv, &r, (dest->length_W = a->length_W + b->length_W) * sizeof(bigint_word_t));
 565                 bigint_adjust(dest);
 566                 return;
 567         }
 568         /* split a in xh & xl; split b in yh & yl */
 569         const bigint_length_t n = (MAX(a->length_W, b->length_W) + 1) / 2;
 570         bigint_t xl, xh, yl, yh;
 571         xl.wordv = a->wordv;
 572         yl.wordv = b->wordv;
 573         if (a->length_W <= n) {
 574                 bigint_set_zero(&xh);
 575                 xl.length_W = a->length_W;
 576                 xl.info = a->info;
 577         } else {
 578                 xl.length_W = n;
 579                 xl.info = 0;
 580                 bigint_adjust(&xl);
 581                 xh.wordv = &(a->wordv[n]);
 582                 xh.length_W = a->length_W-n;
 583                 xh.info = a->info;
 584         }
 585         if (b->length_W <= n) {
 586                 bigint_set_zero(&yh);
 587                 yl.length_W = b->length_W;
 588                 yl.info = b->info;
 589         } else {
 590                 yl.length_W = n;
 591                 yl.info = 0;
 592                 bigint_adjust(&yl);
 593                 yh.wordv = &(b->wordv[n]);
 594                 yh.length_W = b->length_W-n;
 595                 yh.info = b->info;
 596         }
 597         /* now we have split up a and b */
 598         /* remember we want to do:
 599          * x*y = (xh * b**n + xl) * (yh * b**n + yl)
 600          *     = (xh * yh) * b**2n + xh * b**n * yl + yh * b**n * xl + xl * yl
 601          *     = (xh * yh) * b**2n + (xh * yl + yh * xl) * b**n + xl *yl
 602          *     // xh * yl + yh * xl = (xh + yh) * (xl + yl) - xh * yh - xl * yl
 603          * x*y = (xh * yh) * b**2n + ((xh+xl)*(yh+yl) - xh*yh - xl*yl)*b**n + xl*yl
 604          *          5              9     2   4   3    7   5   6   1         8   1
 605          */
 606         ALLOC_BIGINT_WORDS(tmp_w, 2 * n + 2);
 607         ALLOC_BIGINT_WORDS(m_w, 2 * n + 2);
 608         bigint_t tmp, tmp2, m;
 609         tmp.wordv = tmp_w;
 610         tmp2.wordv = &(tmp_w[n + 1]);
 611         m.wordv = m_w;
 612
 613         bigint_mul_u(dest, &xl, &yl);  /* 1: dest <= xl*yl     */
 614         bigint_add_u(&tmp2, &xh, &xl); /* 2: tmp2 <= xh+xl     */
 615         bigint_add_u(&tmp, &yh, &yl);  /* 3: tmp  <= yh+yl     */
 616         bigint_mul_u(&m, &tmp2, &tmp); /* 4: m    <= tmp2*tmp  */
 617         bigint_mul_u(&tmp, &xh, &yh);  /* 5: tmp  <= xh*yh     */
 618         bigint_sub_u(&m, &m, dest);    /* 6: m    <= m-dest    */
 619     bigint_sub_u(&m, &m, &tmp);    /* 7: m    <= m-tmp     */
 620         bigint_add_scale_u(dest, &m, n * sizeof(bigint_word_t));       /* 8: dest <= dest+m**n*/
 621         bigint_add_scale_u(dest, &tmp, 2 * n * sizeof(bigint_word_t)); /* 9: dest <= dest+tmp**(2*n) */
 622         FREE(m_w);
 623         FREE(tmp_w);
 624 }
 625
 626 /******************************************************************************/
 627
 628 void bigint_mul_s(bigint_t *dest, const bigint_t *a, const bigint_t *b){
 629         uint8_t s;
 630         s  = GET_SIGN(a)?2:0;
 631         s |= GET_SIGN(b)?1:0;
 632         switch(s){
 633                 case 0: /* both positive */
 634                         bigint_mul_u(dest, a,b);
 635                         SET_POS(dest);
 636                         break;
 637                 case 1: /* a positive, b negative */
 638                         bigint_mul_u(dest, a,b);
 639                         SET_NEG(dest);
 640                         break;
 641                 case 2: /* a negative, b positive */
 642                         bigint_mul_u(dest, a,b);
 643                         SET_NEG(dest);
 644                         break;
 645                 case 3: /* both negative */
 646                         bigint_mul_u(dest, a,b);
 647                         SET_POS(dest);
 648                         break;
 649                 default: /* how can this happen?*/
 650                         break;
 651         }
 652 }
 653
 654 /******************************************************************************/
 655
 656 void bigint_square(bigint_t *dest, const bigint_t *a) {
 657     union __attribute__((packed)) {
 658         bigint_word_t u[2];
 659         bigint_wordplus_t uv;
 660     } acc;
 661     bigint_word_t q, c1, c2;
 662     bigint_length_t i, j;
 663
 664     if (a->length_W * sizeof(bigint_word_t) <= 4) {
 665         uint64_t r = 0;
 666         memcpy(&r, a->wordv, a->length_W * sizeof(bigint_word_t));
 667         r = r * r;
 668         memcpy(dest->wordv, &r, 2 * a->length_W * sizeof(bigint_word_t));
 669         SET_POS(dest);
 670         dest->length_W = 2 * a->length_W;
 671         bigint_adjust(dest);
 672         return;
 673     }
 674
 675     if (dest->wordv == a->wordv) {
 676         bigint_t d;
 677         ALLOC_BIGINT_WORDS(d_w, a->length_W * 2);
 678         d.wordv = d_w;
 679         bigint_square(&d, a);
 680         bigint_copy(dest, &d);
 681         FREE(d_w);
 682         return;
 683     }
 684
 685     memset(dest->wordv, 0, a->length_W * 2 * sizeof(bigint_word_t));
 686
 687     for (i = 0; i < a->length_W; ++i) {
 688         acc.uv = (bigint_wordplus_t)a->wordv[i] * (bigint_wordplus_t)a->wordv[i] + (bigint_wordplus_t)dest->wordv[2 * i];
 689         dest->wordv[2 * i] = acc.u[0];
 690         c1 = acc.u[1];
 691         c2 = 0;
 692         for (j = i + 1; j < a->length_W; ++j) {
 693             acc.uv = (bigint_wordplus_t)a->wordv[i] * (bigint_wordplus_t)a->wordv[j];
 694             q = acc.u[1] >> (BIGINT_WORD_SIZE - 1);
 695             acc.uv <<= 1;
 696             acc.uv += dest->wordv[i + j];
 697             q += (acc.uv < dest->wordv[i + j]);
 698             acc.uv += c1;
 699             q += (acc.uv < c1);
 700             dest->wordv[i + j] = acc.u[0];
 701             c1 = (bigint_wordplus_t)acc.u[1] + c2;
 702             c2 = q + (c1 < c2);
 703         }
 704         dest->wordv[i + a->length_W] += c1;
 705         if (i < a->length_W - 1) {
 706             dest->wordv[i + a->length_W + 1] = c2 + (dest->wordv[i + a->length_W] < c1);
 707         }
 708     }
 709     dest->info = 0;
 710     dest->length_W = 2 * a->length_W;
 711     bigint_adjust(dest);
 712 }
 713
 714 /******************************************************************************/
 715
 716 void bigint_sub_u_bitscale(bigint_t *a, const bigint_t *b, bigint_length_t bitscale){
 717         bigint_t tmp, x;
 718         const bigint_length_t word_shift = bitscale / BIGINT_WORD_SIZE;
 719
 720         if (a->length_W < b->length_W + word_shift) {
 721 #if DEBUG
 722                 cli_putstr("\r\nDBG: *bang*\r\n");
 723 #endif
 724                 bigint_set_zero(a);
 725                 return;
 726         }
 727     ALLOC_BIGINT_WORDS(tmp_w, b->length_W + 1);
 728         tmp.wordv = tmp_w;
 729         bigint_copy(&tmp, b);
 730         bigint_shiftleft_bits(&tmp, bitscale % BIGINT_WORD_SIZE);
 731
 732         x.info = a->info;
 733         x.wordv = &(a->wordv[word_shift]);
 734         x.length_W = a->length_W - word_shift;
 735
 736         bigint_sub_u(&x, &x, &tmp);
 737         FREE(tmp_w);
 738         bigint_adjust(a);
 739         return;
 740 }
 741
 742 /******************************************************************************/
 743
 744 void bigint_reduce(bigint_t *a, const bigint_t *r){
 745         uint8_t rfbs = GET_FBS(r);
 746         if (r->length_W == 0 || a->length_W == 0) {
 747                 return;
 748         }
 749
 750         if (bigint_length_b(a) + 3 > bigint_length_b(r)) {
 751         if ((r->length_W * sizeof(bigint_word_t) <= 4) && (a->length_W * sizeof(bigint_word_t) <= 4)) {
 752             uint32_t p = 0, q = 0;
 753             memcpy(&p, a->wordv, a->length_W * sizeof(bigint_word_t));
 754             memcpy(&q, r->wordv, r->length_W * sizeof(bigint_word_t));
 755             p %= q;
 756             memcpy(a->wordv, &p, a->length_W * sizeof(bigint_word_t));
 757             a->length_W = r->length_W;
 758             bigint_adjust(a);
 759             return;
 760         }
 761         unsigned shift;
 762         while (a->length_W > r->length_W) {
 763             shift = (a->length_W - r->length_W) * CHAR_BIT * sizeof(bigint_word_t) + GET_FBS(a) - rfbs - 1;
 764             bigint_sub_u_bitscale(a, r, shift);
 765         }
 766         while ((GET_FBS(a) > rfbs) && (a->length_W == r->length_W)) {
 767             shift = GET_FBS(a) - rfbs - 1;
 768             bigint_sub_u_bitscale(a, r, shift);
 769         }
 770         }
 771         while (bigint_cmp_u(a, r) >= 0) {
 772                 bigint_sub_u(a, a, r);
 773         }
 774         bigint_adjust(a);
 775 }
 776
 777 /******************************************************************************/
 778
 779 /* calculate dest = a**exp % r */
 780 /* using square&multiply */
 781 void bigint_expmod_u_sam(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r){
 782         if (a->length_W == 0) {
 783             bigint_set_zero(dest);
 784                 return;
 785         }
 786
 787     if(exp->length_W == 0){
 788         dest->info = 0;
 789         dest->length_W = 1;
 790         dest->wordv[0] = 1;
 791         return;
 792     }
 793
 794     bigint_t res, base;
 795     bigint_word_t t;
 796         bigint_length_t i;
 797         uint8_t j;
 798
 799         ALLOC_BIGINT_WORDS(base_w, MAX(a->length_W, r->length_W));
 800         ALLOC_BIGINT_WORDS(res_w, r->length_W * 2);
 801
 802         res.wordv = res_w;
 803         base.wordv = base_w;
 804         bigint_copy(&base, a);
 805         bigint_reduce(&base, r);
 806         res.wordv[0] = 1;
 807         res.length_W = 1;
 808         res.info = 0;
 809         bigint_adjust(&res);
 810         bigint_copy(&res, &base);
 811         uint8_t flag = 0;
 812         t = exp->wordv[exp->length_W - 1];
 813         for (i = exp->length_W; i > 0; --i) {
 814                 t = exp->wordv[i - 1];
 815                 for (j = BIGINT_WORD_SIZE; j > 0; --j) {
 816                         if (!flag) {
 817                                 if (t & (((bigint_word_t)1) << (BIGINT_WORD_SIZE - 1))) {
 818                                         flag = 1;
 819                                 }
 820                         } else {
 821                                 bigint_square(&res, &res);
 822                                 bigint_reduce(&res, r);
 823                                 if(t & (((bigint_word_t)1) << (BIGINT_WORD_SIZE - 1))){
 824                                         bigint_mul_u(&res, &res, &base);
 825                                         bigint_reduce(&res, r);
 826                                 }
 827                         }
 828                         t <<= 1;
 829                 }
 830         }
 831
 832         SET_POS(&res);
 833         bigint_copy(dest, &res);
 834         FREE(res_w);
 835         FREE(base_w);
 836 }
 837
 838 /******************************************************************************/
 839 /* gcd <-- gcd(x,y) a*x+b*y=gcd */
 840 void bigint_gcdext(bigint_t *gcd, bigint_t *a, bigint_t *b, const bigint_t *x, const bigint_t *y){
 841          bigint_length_t i = 0;
 842          if(x->length_W == 0 || y->length_W == 0){
 843              if(gcd){
 844                  bigint_set_zero(gcd);
 845              }
 846              if(a){
 847                  bigint_set_zero(a);
 848              }
 849          if(b){
 850              bigint_set_zero(b);
 851          }
 852                  return;
 853          }
 854          if(x->length_W == 1 && x->wordv[0] == 1){
 855              if(gcd){
 856              gcd->length_W = 1;
 857              gcd->wordv[0] = 1;
 858              gcd->info = 0;
 859              }
 860                  if(a){
 861                          a->length_W = 1;
 862                          a->wordv[0] = 1;
 863                          SET_POS(a);
 864                          bigint_adjust(a);
 865                  }
 866                  if(b){
 867                          bigint_set_zero(b);
 868                  }
 869                  return;
 870          }
 871          if(y->length_W == 1 && y->wordv[0] == 1){
 872                  if(gcd){
 873              gcd->length_W = 1;
 874              gcd->wordv[0] = 1;
 875              gcd->info = 0;
 876                  }
 877                  if(b){
 878                          b->length_W = 1;
 879                          b->wordv[0] = 1;
 880                          SET_POS(b);
 881                          bigint_adjust(b);
 882                  }
 883                  if(a){
 884                          bigint_set_zero(a);
 885                  }
 886                  return;
 887          }
 888
 889          while(x->wordv[i] == 0 && y->wordv[i] == 0){
 890                  ++i;
 891          }
 892
 893          ALLOC_BIGINT_WORDS(g_w, i + 2);
 894          ALLOC_BIGINT_WORDS(x_w, x->length_W - i);
 895          ALLOC_BIGINT_WORDS(y_w, y->length_W - i);
 896          ALLOC_BIGINT_WORDS(u_w, x->length_W - i);
 897      ALLOC_BIGINT_WORDS(v_w, y->length_W - i);
 898      ALLOC_BIGINT_WORDS(a_w, y->length_W + 2);
 899      ALLOC_BIGINT_WORDS(c_w, y->length_W + 2);
 900      ALLOC_BIGINT_WORDS(b_w, x->length_W + 2);
 901      ALLOC_BIGINT_WORDS(d_w, x->length_W + 2);
 902
 903          bigint_t g, x_, y_, u, v, a_, b_, c_, d_;
 904
 905          g.wordv = g_w;
 906          x_.wordv = x_w;
 907          y_.wordv = y_w;
 908          memset(g_w, 0, i * sizeof(bigint_word_t));
 909          g_w[i] = 1;
 910          g.length_W = i + 1;
 911          g.info = 0;
 912          x_.info = y_.info = 0;
 913          x_.length_W = x->length_W - i;
 914          y_.length_W = y->length_W - i;
 915          memcpy(x_.wordv, x->wordv + i, x_.length_W * sizeof(bigint_word_t));
 916          memcpy(y_.wordv, y->wordv + i, y_.length_W * sizeof(bigint_word_t));
 917
 918          for(i = 0; (x_.wordv[0] & ((bigint_word_t)1 << i)) == 0 && (y_.wordv[0] & ((bigint_word_t)1 << i)) == 0; ++i)
 919              ;
 920
 921          bigint_adjust(&x_);
 922          bigint_adjust(&y_);
 923
 924          if(i){
 925                  bigint_shiftleft_bits(&g, i);
 926                  bigint_shiftright_bits(&x_, i);
 927                  bigint_shiftright_bits(&y_, i);
 928          }
 929
 930          u.wordv = u_w;
 931          v.wordv = v_w;
 932          a_.wordv = a_w;
 933          b_.wordv = b_w;
 934          c_.wordv = c_w;
 935          d_.wordv = d_w;
 936
 937          bigint_copy(&u, &x_);
 938          bigint_copy(&v, &y_);
 939          a_.wordv[0] = 1;
 940          a_.length_W = 1;
 941          a_.info = 0;
 942          d_.wordv[0] = 1;
 943          d_.length_W = 1;
 944          d_.info = 0;
 945          bigint_set_zero(&b_);
 946          bigint_set_zero(&c_);
 947          do {
 948                  while ((u.wordv[0] & 1) == 0) {
 949                          bigint_shiftright_1bit(&u);
 950                          if((a_.wordv[0] & 1) || (b_.wordv[0] & 1)){
 951                                  bigint_add_s(&a_, &a_, &y_);
 952                                  bigint_sub_s(&b_, &b_, &x_);
 953                          }
 954                          bigint_shiftright_1bit(&a_);
 955                          bigint_shiftright_1bit(&b_);
 956                  }
 957                  while ((v.wordv[0] & 1) == 0) {
 958                      bigint_shiftright_1bit(&v);
 959                          if((c_.wordv[0] & 1) || (d_.wordv[0] & 1)){
 960                                  bigint_add_s(&c_, &c_, &y_);
 961                                  bigint_sub_s(&d_, &d_, &x_);
 962                          }
 963                          bigint_shiftright_1bit(&c_);
 964                          bigint_shiftright_1bit(&d_);
 965                  }
 966                  if(bigint_cmp_u(&u, &v) >= 0){
 967                         bigint_sub_u(&u, &u, &v);
 968                         bigint_sub_s(&a_, &a_, &c_);
 969                         bigint_sub_s(&b_, &b_, &d_);
 970                  }else{
 971                         bigint_sub_u(&v, &v, &u);
 972                         bigint_sub_s(&c_, &c_, &a_);
 973                         bigint_sub_s(&d_, &d_, &b_);
 974                  }
 975          } while(u.length_W);
 976          if(gcd){
 977                  bigint_mul_s(gcd, &v, &g);
 978          }
 979          if(a){
 980                 bigint_copy(a, &c_);
 981          }
 982          if(b){
 983                  bigint_copy(b, &d_);
 984          }
 985
 986          FREE(d_w);
 987          FREE(b_w);
 988          FREE(c_w);
 989          FREE(a_w);
 990      FREE(v_w);
 991      FREE(u_w);
 992      FREE(y_w);
 993      FREE(x_w);
 994      FREE(g_w);
 995 }
 996
 997 /******************************************************************************/
 998
 999 void bigint_inverse(bigint_t *dest, const bigint_t *a, const bigint_t *m){
1000         bigint_gcdext(NULL, dest, NULL, a, m);
1001         while(dest->info&BIGINT_NEG_MASK){
1002                 bigint_add_s(dest, dest, m);
1003         }
1004 }
1005
1006 /******************************************************************************/
1007
1008 void bigint_changeendianess(bigint_t *a){
1009         uint8_t t, *p, *q;
1010         p = (uint8_t*)a->wordv;
1011         q = p + a->length_W * sizeof(bigint_word_t) - 1;
1012         while(p < q){
1013                 t = *p;
1014                 *p = *q;
1015                 *q = t;
1016                 ++p; --q;
1017         }
1018 }
1019
1020 /******************************************************************************/
1021
1022 void bigint_mul_word_u(bigint_t *a, bigint_word_t b){
1023     bigint_wordplus_t c0 = 0, c1 = 0;
1024     bigint_length_t i;
1025
1026     if(b == 0){
1027         bigint_set_zero(a);
1028         return;
1029     }
1030
1031     for(i = 0; i < a->length_W; ++i){
1032         c1 = ((bigint_wordplus_t)(a->wordv[i])) * ((bigint_wordplus_t)b);
1033         c1 += c0;
1034         a->wordv[i] = (bigint_word_t)c1;
1035         c0 = c1 >> BIGINT_WORD_SIZE;
1036     }
1037     if(c0){
1038         a->wordv[a->length_W] = (bigint_word_t)c0;
1039         a->length_W += 1;
1040     }
1041     bigint_adjust(a);
1042 }
1043
1044 /******************************************************************************/
1045
1046 void bigint_clip(bigint_t *dest, bigint_length_t length_W){
1047     if(dest->length_W > length_W){
1048         dest->length_W = length_W;
1049     }
1050     bigint_adjust(dest);
1051 }
1052
1053 /******************************************************************************/
1054 /*
1055  * m_ = m * m'[0]
1056  * dest = (a * b) % m (?)
1057  */
1058
1059 void bigint_mont_mul(bigint_t *dest, const bigint_t *a, const bigint_t *b, const bigint_t *m, const bigint_t *m_){
1060     const bigint_length_t s = MAX(MAX(a->length_W, b->length_W), m->length_W);
1061     bigint_t u, t;
1062
1063     bigint_length_t i;
1064
1065     if (a->length_W == 0 || b->length_W == 0) {
1066         bigint_set_zero(dest);
1067         return;
1068     }
1069     ALLOC_BIGINT_WORDS(u_w, s + 2);
1070     ALLOC_BIGINT_WORDS(t_w, s + 2);
1071     u.wordv = u_w;
1072     u.info = 0;
1073     u.length_W = 0;
1074     t.wordv = t_w;
1075     for (i = 0; i < a->length_W; ++i) {
1076         bigint_copy(&t, b);
1077         bigint_mul_word_u(&t, a->wordv[i]);
1078         bigint_add_u(&u, &u, &t);
1079         bigint_copy(&t, m_);
1080         if (u.length_W != 0) {
1081             bigint_mul_word_u(&t, u.wordv[0]);
1082             bigint_add_u(&u, &u, &t);
1083         }
1084         bigint_shiftright_1word(&u);
1085     }
1086     for (; i < s; ++i) {
1087         bigint_copy(&t, m_);
1088         if (u.length_W != 0) {
1089             bigint_mul_word_u(&t, u.wordv[0]);
1090             bigint_add_u(&u, &u, &t);
1091         }
1092         bigint_shiftright_1word(&u);
1093     }
1094     bigint_reduce(&u, m);
1095     bigint_copy(dest, &u);
1096     FREE(t_w);
1097     FREE(u_w);
1098 }
1099
1100 /******************************************************************************/
1101
1102 void bigint_mont_red(bigint_t *dest, const bigint_t *a, const bigint_t *m, const bigint_t *m_){
1103     bigint_t u, t;
1104     bigint_length_t i, s = MAX(a->length_W, MAX(m->length_W, m_->length_W));
1105
1106     if (a->length_W == 0) {
1107         bigint_set_zero(dest);
1108         return;
1109     }
1110
1111     ALLOC_BIGINT_WORDS(u_w, s + 2);
1112     ALLOC_BIGINT_WORDS(t_w, s + 2);
1113     t.wordv = t_w;
1114     u.wordv = u_w;
1115     bigint_copy(&u, a);
1116     for (i = 0; i < m->length_W; ++i) {
1117         bigint_copy(&t, m_);
1118         if (u.length_W != 0) {
1119             bigint_mul_word_u(&t, u.wordv[0]);
1120             bigint_add_u(&u, &u, &t);
1121         }
1122         bigint_shiftright_1word(&u);
1123     }
1124     bigint_reduce(&u, m);
1125     bigint_copy(dest, &u);
1126     FREE(t_w);
1127     FREE(u_w);
1128 }
1129
1130 /******************************************************************************/
1131 /*
1132  * m_ = m * (- m0^-1 (mod 2^W))
1133  */
1134 void bigint_mont_gen_m_(bigint_t* dest, const bigint_t* m){
1135     bigint_word_t x_w[2], m_w_0[1];
1136     bigint_t x, m_0;
1137     if (m->length_W == 0) {
1138         bigint_set_zero(dest);
1139         return;
1140     }
1141     if ((m->wordv[0] & 1) == 0) {
1142 #if DEBUG
1143         printf_P(PSTR("ERROR: m must not be even, m = "));
1144         bigint_print_hex(m);
1145         putchar('\n');
1146         uart0_flush();
1147 #endif
1148         return;
1149     }
1150     x.wordv = x_w;
1151     x.info = 0;
1152     x.length_W = 2;
1153     x_w[0] = 0;
1154     x_w[1] = 1;
1155     m_0.wordv = m_w_0;
1156     m_0.info = 0;
1157     m_0.length_W = 1;
1158     m_0.wordv[0] = m->wordv[0];
1159     bigint_adjust(&x);
1160     bigint_adjust(&m_0);
1161     bigint_inverse(dest, &m_0, &x);
1162     bigint_sub_s(&x, &x, dest);
1163     bigint_copy(dest, m);
1164     bigint_mul_word_u(dest, x.wordv[0]);
1165
1166 }
1167
1168 /******************************************************************************/
1169
1170 /*
1171  * dest = a * R mod m
1172  */
1173 void bigint_mont_trans(bigint_t *dest, const bigint_t *a, const bigint_t *m){
1174     bigint_t t;
1175
1176     ALLOC_BIGINT_WORDS(t_w, a->length_W + m->length_W);
1177     t.wordv = t_w;
1178     memset(t_w, 0, m->length_W * sizeof(bigint_word_t));
1179     memcpy(&t_w[m->length_W], a->wordv, a->length_W * sizeof(bigint_word_t));
1180     t.info = a->info;
1181     t.length_W = a->length_W + m->length_W;
1182     bigint_reduce(&t, m);
1183     bigint_copy(dest, &t);
1184     FREE(t_w);
1185 }
1186
1187 /******************************************************************************/
1188
1189 /* calculate dest = a**exp % r */
1190 /* using square&multiply */
1191 void bigint_expmod_u_mont_accel(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r, const bigint_t *m_){
1192     if(r->length_W == 0) {
1193         return;
1194     }
1195
1196     bigint_t res, ax;
1197     bigint_word_t t;
1198     bigint_length_t i;
1199     uint8_t j;
1200
1201     if (exp->length_W == 0) {
1202         dest->length_W = 1;
1203         dest->info = 0;
1204         dest->wordv[0] = 1;
1205         return;
1206     }
1207
1208     ALLOC_BIGINT_WORDS(res_w, r->length_W * 2);
1209     ALLOC_BIGINT_WORDS(ax_w, MAX(r->length_W, a->length_W));
1210
1211     res.wordv = res_w;
1212     ax.wordv = ax_w;
1213
1214     res.wordv[0] = 1;
1215     res.length_W = 1;
1216     res.info = 0;
1217
1218     bigint_copy(&ax, a);
1219     bigint_reduce(&ax, r);
1220
1221     bigint_mont_trans(&ax, &ax, r);
1222     bigint_mont_trans(&res, &res, r);
1223
1224     uint8_t flag = 0;
1225     t = exp->wordv[exp->length_W - 1];
1226     for (i = exp->length_W; i > 0; --i) {
1227         t = exp->wordv[i - 1];
1228         for(j = BIGINT_WORD_SIZE; j > 0; --j){
1229             if (!flag) {
1230                 if(t & (((bigint_word_t)1) << (BIGINT_WORD_SIZE - 1))){
1231                     flag = 1;
1232                 }
1233             }
1234             if (flag) {
1235                 bigint_square(&res, &res);
1236                 bigint_mont_red(&res, &res, r, m_);
1237                 if (t & (((bigint_word_t)1) << (BIGINT_WORD_SIZE - 1))) {
1238                     bigint_mont_mul(&res, &res, &ax, r, m_);
1239                 }
1240             }
1241             t <<= 1;
1242         }
1243     }
1244     SET_POS(&res);
1245     bigint_mont_red(dest, &res, r, m_);
1246     FREE(ax_w);
1247     FREE(res_w);
1248 }
1249
1250 /******************************************************************************/
1251
1252 void bigint_expmod_u_mont_sam(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r){
1253     if(r->length_W == 0) {
1254         return;
1255     }
1256     if(a->length_W == 0) {
1257         bigint_set_zero(dest);
1258         return;
1259     }
1260     bigint_t m_;
1261     bigint_word_t m_w_[r->length_W + 1];
1262     m_.wordv = m_w_;
1263     bigint_mont_gen_m_(&m_, r);
1264     bigint_expmod_u_mont_accel(dest, a, exp, r,&m_);
1265 }
1266
1267 /******************************************************************************/
1268
1269 void bigint_expmod_u(bigint_t *dest, const bigint_t *a, const bigint_t *exp, const bigint_t *r){
1270 #if 0
1271     printf("\nDBG: expmod_u (a ** e %% m) <%s %s %d>\n\ta: ", __FILE__, __func__, __LINE__);
1272     bigint_print_hex(a);
1273     printf("\n\te: ");
1274     bigint_print_hex(exp);
1275     printf("\n\tm: ");
1276     bigint_print_hex(r);
1277 #endif
1278     if (r->wordv[0] & 1) {
1279         bigint_expmod_u_mont_sam(dest, a, exp, r);
1280     } else {
1281         bigint_expmod_u_sam(dest, a, exp, r);
1282     }
1283 }
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299