utils/zenutils/libraries/beecrypt-4.1.2/beecrypt/mp.c

   1 /*
   2  * Copyright (c) 2002, 2003 Bob Deblier
   3  *
   4  * This library is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU Lesser General Public
   6  * License as published by the Free Software Foundation; either
   7  * version 2.1 of the License, or (at your option) any later version.
   8  *
   9  * This library is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12  * Lesser General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU Lesser General Public
  15  * License along with this library; if not, write to the Free Software
  16  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  17  *
  18  */
  19
  20 /*!\file mp.c
  21  * \brief Multi-precision integer routines.
  22  * \author Bob Deblier <bob.deblier@pandora.be>
  23  * \ingroup MP_m
  24  */
  25
  26 #define BEECRYPT_DLL_EXPORT
  27
  28 #if HAVE_CONFIG_H
  29 # include "config.h"
  30 #endif
  31
  32 #include "beecrypt/mp.h"
  33 #include "beecrypt/mpopt.h"
  34
  35 #ifndef ASM_MPZERO
  36 void mpzero(size_t size, mpw* data)
  37 {
  38         while (size--)
  39                 *(data++) = 0;
  40 }
  41 #endif
  42
  43 #ifndef ASM_MPFILL
  44 void mpfill(size_t size, mpw* data, mpw fill)
  45 {
  46         while (size--)
  47                 *(data++) = fill;
  48 }
  49 #endif
  50
  51 #ifndef ASM_MPODD
  52 int mpodd(size_t size, const mpw* data)
  53 {
  54         return (int)(data[size-1] & 0x1);
  55 }
  56 #endif
  57
  58 #ifndef ASM_MPEVEN
  59 int mpeven(size_t size, const mpw* data)
  60 {
  61         return !(int)(data[size-1] & 0x1);
  62 }
  63 #endif
  64
  65 #ifndef ASM_MPZ
  66 int mpz(size_t size, const mpw* data)
  67 {
  68         while (size--)
  69                 if (*(data++))
  70                         return 0;
  71         return 1;
  72 }
  73 #endif
  74
  75 #ifndef ASM_MPNZ
  76 int mpnz(size_t size, const mpw* data)
  77 {
  78         while (size--)
  79                 if (*(data++))
  80                         return 1;
  81         return 0;
  82 }
  83 #endif
  84
  85 #ifndef ASM_MPEQ
  86 int mpeq(size_t size, const mpw* xdata, const mpw* ydata)
  87 {
  88         while (size--)
  89         {
  90                 if (*xdata == *ydata)
  91                 {
  92                         xdata++;
  93                         ydata++;
  94                 }
  95                 else
  96                         return 0;
  97         }
  98         return 1;
  99 }
 100 #endif
 101
 102 #ifndef ASM_MPEQX
 103 int mpeqx(size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata)
 104 {
 105         if (xsize > ysize)
 106         {
 107                 register size_t diff = xsize - ysize;
 108                 return mpeq(ysize, xdata+diff, ydata) && mpz(diff, xdata);
 109         }
 110         else if (xsize < ysize)
 111         {
 112                 register size_t diff = ysize - xsize;
 113                 return mpeq(xsize, ydata+diff, xdata) && mpz(diff, ydata);
 114         }
 115         else
 116                 return mpeq(xsize, xdata, ydata);
 117 }
 118 #endif
 119
 120 #ifndef ASM_MPNE
 121 int mpne(size_t size, const mpw* xdata, const mpw* ydata)
 122 {
 123         while (size--)
 124         {
 125                 if (*xdata == *ydata)
 126                 {
 127                         xdata++;
 128                         ydata++;
 129                 }
 130                 else
 131                         return 1;
 132         }
 133         return 0;
 134 }
 135 #endif
 136
 137 #ifndef ASM_MPNEX
 138 int mpnex(size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata)
 139 {
 140         if (xsize > ysize)
 141         {
 142                 register size_t diff = xsize - ysize;
 143                 return mpnz(diff, xdata) || mpne(ysize, xdata+diff, ydata);
 144         }
 145         else if (xsize < ysize)
 146         {
 147                 register size_t diff = ysize - xsize;
 148                 return mpnz(diff, ydata) || mpne(xsize, ydata+diff, xdata);
 149         }
 150         else
 151                 return mpne(xsize, xdata, ydata);
 152 }
 153 #endif
 154
 155 #ifndef ASM_MPGT
 156 int mpgt(size_t size, const mpw* xdata, const mpw* ydata)
 157 {
 158         while (size--)
 159         {
 160                 if (*xdata < *ydata)
 161                         return 0;
 162                 if (*xdata > *ydata)
 163                         return 1;
 164                 xdata++; ydata++;
 165         }
 166         return 0;
 167 }
 168 #endif
 169
 170 #ifndef ASM_MPGTX
 171 int mpgtx(size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata)
 172 {
 173         if (xsize > ysize)
 174         {
 175                 register size_t diff = xsize - ysize;
 176                 return mpnz(diff, xdata) || mpgt(ysize, xdata + diff, ydata);
 177         }
 178         else if (xsize < ysize)
 179         {
 180                 register size_t diff = ysize - xsize;
 181                 return mpz(diff, ydata) && mpgt(xsize, xdata, ydata + diff);
 182         }
 183         else
 184                 return mpgt(xsize, xdata, ydata);
 185 }
 186 #endif
 187
 188 #ifndef ASM_MPLT
 189 int mplt(size_t size, const mpw* xdata, const mpw* ydata)
 190 {
 191         while (size--)
 192         {
 193                 if (*xdata > *ydata)
 194                         return 0;
 195                 if (*xdata < *ydata)
 196                         return 1;
 197                 xdata++; ydata++;
 198         }
 199         return 0;
 200 }
 201 #endif
 202
 203 #ifndef ASM_MPLTX
 204 int mpltx(size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata)
 205 {
 206         if (xsize > ysize)
 207         {
 208                 register size_t diff = xsize - ysize;
 209                 return mpz(diff, xdata) && mplt(ysize, xdata+diff, ydata);
 210         }
 211         else if (xsize < ysize)
 212         {
 213                 register size_t diff = ysize - xsize;
 214                 return mpnz(diff, ydata) || mplt(xsize, xdata, ydata+diff);
 215         }
 216         else
 217                 return mplt(xsize, xdata, ydata);
 218 }
 219 #endif
 220
 221 #ifndef ASM_MPGE
 222 int mpge(size_t size, const mpw* xdata, const mpw* ydata)
 223 {
 224         while (size--)
 225         {
 226                 if (*xdata < *ydata)
 227                         return 0;
 228                 if (*xdata > *ydata)
 229                         return 1;
 230                 xdata++; ydata++;
 231         }
 232         return 1;
 233 }
 234 #endif
 235
 236 #ifndef ASM_MPGEX
 237 int mpgex(size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata)
 238 {
 239         if (xsize > ysize)
 240         {
 241                 register size_t diff = xsize - ysize;
 242                 return mpnz(diff, xdata) || mpge(ysize, xdata+diff, ydata);
 243         }
 244         else if (xsize < ysize)
 245         {
 246                 register size_t diff = ysize - xsize;
 247                 return mpz(diff, ydata) && mpge(xsize, xdata, ydata+diff);
 248         }
 249         else
 250                 return mpge(xsize, xdata, ydata);
 251 }
 252 #endif
 253
 254 #ifndef ASM_MPLE
 255 int mple(size_t size, const mpw* xdata, const mpw* ydata)
 256 {
 257         while (size--)
 258         {
 259                 if (*xdata < *ydata)
 260                         return 1;
 261                 if (*xdata > *ydata)
 262                         return 0;
 263                 xdata++; ydata++;
 264         }
 265         return 1;
 266 }
 267 #endif
 268
 269 #ifndef ASM_MPLEX
 270 int mplex(size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata)
 271 {
 272         if (xsize > ysize)
 273         {
 274                 register size_t diff = xsize - ysize;
 275                 return mpz(diff, xdata) && mple(ysize, xdata+ diff, ydata);
 276         }
 277         else if (xsize < ysize)
 278         {
 279                 register size_t diff = ysize - xsize;
 280                 return mpnz(diff, ydata) || mple(xsize, xdata, ydata+diff);
 281         }
 282         else
 283                 return mple(xsize, xdata, ydata);
 284 }
 285 #endif
 286
 287 #ifndef ASM_MPISONE
 288 int mpisone(size_t size, const mpw* data)
 289 {
 290         data += size;
 291         if (*(--data) == 1)
 292         {
 293                 while (--size)
 294                         if (*(--data))
 295                                 return 0;
 296                 return 1;
 297         }
 298         return 0;
 299 }
 300 #endif
 301
 302 #ifndef ASM_MPISTWO
 303 int mpistwo(size_t size, const mpw* data)
 304 {
 305         data += size;
 306         if (*(--data) == 2)
 307         {
 308                 while (--size)
 309                         if (*(--data))
 310                                 return 0;
 311                 return 1;
 312         }
 313         return 0;
 314 }
 315 #endif
 316
 317 #ifndef ASM_MPEQMONE
 318 int mpeqmone(size_t size, const mpw* xdata, const mpw* ydata)
 319 {
 320     xdata += size;
 321     ydata += size;
 322
 323     if (*(--xdata)+1 == *(--ydata))
 324     {
 325         while (--size)
 326             if (*(--xdata) != *(--ydata))
 327                 return 0;
 328         return 1;
 329     }
 330     return 0;
 331 }
 332 #endif
 333
 334 #ifndef ASM_MPLEONE
 335 int mpleone(size_t size, const mpw* data)
 336 {
 337         data += size;
 338         if (*(--data) > 1)
 339                 return 0;
 340         else
 341         {
 342                 while (--size)
 343                         if (*(--data))
 344                                 return 0;
 345                 return 1;
 346         }
 347 }
 348 #endif
 349
 350 #ifndef ASM_MPMSBSET
 351 int mpmsbset(size_t size, const mpw* data)
 352 {
 353         return (int)((*data) >> (MP_WBITS-1));
 354 }
 355 #endif
 356
 357 #ifndef ASM_MPLSBSET
 358 int mplsbset(size_t size, const mpw* data)
 359 {
 360     return (int)(data[size-1] & 0x1);
 361 }
 362 #endif
 363
 364 #ifndef ASM_MPSETMSB
 365 void mpsetmsb(size_t size, mpw* data)
 366 {
 367         *data |= MP_MSBMASK;
 368 }
 369 #endif
 370
 371 #ifndef ASM_MPSETLSB
 372 void mpsetlsb(size_t size, mpw* data)
 373 {
 374         data[size-1] |= MP_LSBMASK;
 375 }
 376 #endif
 377
 378 #ifndef ASM_MPCLRMSB
 379 void mpclrmsb(size_t size, mpw* data)
 380 {
 381         *data &= ~ MP_MSBMASK;
 382 }
 383 #endif
 384
 385 #ifndef ASM_MPCLRLSB
 386 void mpclrlsb(size_t size, mpw* data)
 387 {
 388     data[size-1] &= ~ MP_LSBMASK;
 389 }
 390 #endif
 391
 392 #ifndef ASM_MPAND
 393 void mpand(size_t size, mpw* xdata, const mpw* ydata)
 394 {
 395         while (size--)
 396                 xdata[size] &= ydata[size];
 397 }
 398 #endif
 399
 400 #ifndef ASM_MPOR
 401 void mpor(size_t size, mpw* xdata, const mpw* ydata)
 402 {
 403         while (size--)
 404                 xdata[size] |= ydata[size];
 405 }
 406 #endif
 407
 408 #ifndef ASM_MPXOR
 409 void mpxor(size_t size, mpw* xdata, const mpw* ydata)
 410 {
 411         while (size--)
 412                 xdata[size] ^= ydata[size];
 413 }
 414 #endif
 415
 416 #ifndef ASM_MPNOT
 417 void mpnot(size_t size, mpw* data)
 418 {
 419         while (size--)
 420                 data[size] = ~data[size];
 421 }
 422 #endif
 423
 424 #ifndef ASM_MPSETW
 425 void mpsetw(size_t size, mpw* xdata, mpw y)
 426 {
 427         while (--size)
 428                 *(xdata++) = 0;
 429         *(xdata++) = y;
 430 }
 431 #endif
 432
 433 #ifndef ASM_MPSETX
 434 void mpsetx(size_t xsize, mpw* xdata, size_t ysize, const mpw* ydata)
 435 {
 436         while (xsize > ysize)
 437         {
 438                 xsize--;
 439                 *(xdata++) = 0;
 440         }
 441         while (ysize > xsize)
 442         {
 443                 ysize--;
 444                 ydata++;
 445         }
 446         while (xsize--)
 447                 *(xdata++) = *(ydata++);
 448 }
 449 #endif
 450
 451 #ifndef ASM_MPADDW
 452 int mpaddw(size_t size, mpw* xdata, mpw y)
 453 {
 454         register mpw load, temp;
 455         register int carry = 0;
 456
 457         xdata += size-1;
 458
 459         load = *xdata;
 460         temp = load + y;
 461         *(xdata--) = temp;
 462         carry = (load > temp);
 463
 464         while (--size && carry)
 465         {
 466                 load = *xdata;
 467                 temp = load + 1;
 468                 *(xdata--) = temp;
 469                 carry = (load > temp);
 470         }
 471         return carry;
 472 }
 473 #endif
 474
 475 #ifndef ASM_MPADD
 476 int mpadd(size_t size, mpw* xdata, const mpw* ydata)
 477 {
 478         register mpw load, temp;
 479         register int carry = 0;
 480
 481         xdata += size-1;
 482         ydata += size-1;
 483
 484         while (size--)
 485         {
 486                 temp = *(ydata--);
 487                 load = *xdata;
 488                 temp = carry ? (load + temp + 1) : (load + temp);
 489                 *(xdata--) = temp;
 490                 carry = carry ? (load >= temp) : (load > temp);
 491         }
 492         return carry;
 493 }
 494 #endif
 495
 496 #ifndef ASM_MPADDX
 497 int mpaddx(size_t xsize, mpw* xdata, size_t ysize, const mpw* ydata)
 498 {
 499         if (xsize > ysize)
 500         {
 501                 register size_t diff = xsize - ysize;
 502                 return mpaddw(diff, xdata, (mpw) mpadd(ysize, xdata+diff, ydata));
 503         }
 504         else
 505         {
 506                 register size_t diff = ysize - xsize;
 507                 return mpadd(xsize, xdata, ydata+diff);
 508         }
 509 }
 510 #endif
 511
 512 #ifndef ASM_MPSUBW
 513 int mpsubw(size_t size, mpw* xdata, mpw y)
 514 {
 515         register mpw load, temp;
 516         register int carry = 0;
 517
 518         xdata += size-1;
 519
 520         load = *xdata;
 521         temp = load - y;
 522         *(xdata--) = temp;
 523         carry = (load < temp);
 524
 525         while (--size && carry)
 526         {
 527                 load = *xdata;
 528                 temp = load - 1;
 529                 *(xdata--) = temp;
 530                 carry = (load < temp);
 531         }
 532         return carry;
 533 }
 534 #endif
 535
 536 #ifndef ASM_MPSUB
 537 int mpsub(size_t size, mpw* xdata, const mpw* ydata)
 538 {
 539         register mpw load, temp;
 540         register int carry = 0;
 541
 542         xdata += size-1;
 543         ydata += size-1;
 544
 545         while (size--)
 546         {
 547                 temp = *(ydata--);
 548                 load = *xdata;
 549                 temp = carry ? (load - temp - 1) : (load - temp);
 550                 *(xdata--) = temp;
 551                 carry = carry ? (load <= temp) : (load < temp);
 552         }
 553         return carry;
 554 }
 555 #endif
 556
 557 #ifndef ASM_MPSUBX
 558 int mpsubx(size_t xsize, mpw* xdata, size_t ysize, const mpw* ydata)
 559 {
 560         if (xsize > ysize)
 561         {
 562                 register size_t diff = xsize - ysize;
 563                 return mpsubw(diff, xdata, (mpw) mpsub(ysize, xdata+diff, ydata));
 564         }
 565         else
 566         {
 567                 register size_t diff = ysize - xsize;
 568                 return mpsub(xsize, xdata, ydata+diff);
 569         }
 570 }
 571 #endif
 572
 573 #ifndef ASM_MPNEG
 574 void mpneg(size_t size, mpw* data)
 575 {
 576         mpnot(size, data);
 577         mpaddw(size, data, 1);
 578 }
 579 #endif
 580
 581 #ifndef ASM_MPSETMUL
 582 mpw mpsetmul(size_t size, mpw* result, const mpw* data, mpw y)
 583 {
 584         #if HAVE_MPDW
 585         register mpdw temp;
 586         register mpw  carry = 0;
 587
 588         data += size;
 589         result += size;
 590
 591         while (size--)
 592         {
 593                 temp = *(--data);
 594                 temp *= y;
 595                 temp += carry;
 596                 *(--result) = (mpw) temp;
 597                 carry = (mpw)(temp >> MP_WBITS);
 598         }
 599         #else
 600         register mpw temp, load, carry = 0;
 601         register mphw ylo, yhi;
 602
 603         ylo = (mphw)  y;
 604         yhi = (mphw) (y >> MP_HWBITS);
 605
 606         data += size;
 607         result += size;
 608
 609         while (size--)
 610         {
 611                 register mphw xlo, xhi;
 612                 register mpw rlo, rhi;
 613
 614                 xlo = (mphw) (temp = *(--data));
 615                 xhi = (mphw) (temp >> MP_HWBITS);
 616
 617                 rlo = (mpw) xlo * ylo;
 618                 rhi = (mpw) xhi * yhi;
 619                 load = rlo;
 620                 temp = (mpw) xhi * ylo;
 621                 rlo += (temp << MP_HWBITS);
 622                 rhi += (temp >> MP_HWBITS) + (load > rlo);
 623                 load = rlo;
 624                 temp = (mpw) xlo * yhi;
 625                 rlo += (temp << MP_HWBITS);
 626                 rhi += (temp >> MP_HWBITS) + (load > rlo);
 627                 load = rlo;
 628                 temp = rlo + carry;
 629                 carry = rhi + (load > temp);
 630                 *(--result) = temp;
 631         }
 632         #endif
 633         return carry;
 634 }
 635 #endif
 636
 637 #ifndef ASM_MPADDMUL
 638 mpw mpaddmul(size_t size, mpw* result, const mpw* data, mpw y)
 639 {
 640         #if HAVE_MPDW
 641         register mpdw temp;
 642         register mpw  carry = 0;
 643
 644         data += size;
 645         result += size;
 646
 647         while (size--)
 648         {
 649                 temp = *(--data);
 650                 temp *= y;
 651                 temp += carry;
 652                 temp += *(--result);
 653                 *result = (mpw) temp;
 654                 carry = (mpw)(temp >> MP_WBITS);
 655         }
 656         #else
 657         register mpw temp, load, carry = 0;
 658         register mphw ylo, yhi;
 659
 660         ylo = (mphw)  y;
 661         yhi = (mphw) (y >> MP_HWBITS);
 662
 663         data += size;
 664         result += size;
 665
 666         while (size--)
 667         {
 668                 register mphw xlo, xhi;
 669                 register mpw rlo, rhi;
 670
 671                 xlo = (mphw) (temp = *(--data));
 672                 xhi = (mphw) (temp >> MP_HWBITS);
 673
 674                 rlo = (mpw) xlo * ylo;
 675                 rhi = (mpw) xhi * yhi;
 676                 load = rlo;
 677                 temp = (mpw) xhi * ylo;
 678                 rlo += (temp << MP_HWBITS);
 679                 rhi += (temp >> MP_HWBITS) + (load > rlo);
 680                 load = rlo;
 681                 temp = (mpw) xlo * yhi;
 682                 rlo += (temp << MP_HWBITS);
 683                 rhi += (temp >> MP_HWBITS) + (load > rlo);
 684                 load = rlo;
 685                 rlo += carry;
 686                 temp = (load > rlo);
 687                 load = rhi;
 688                 rhi += temp;
 689                 carry = (load > rhi);
 690                 load = rlo;
 691                 rlo += *(--result);
 692                 *result = rlo;
 693                 carry += rhi + (load > rlo);
 694         }
 695         #endif
 696         return carry;
 697 }
 698 #endif
 699
 700 #ifndef ASM_MPMUL
 701 void mpmul(mpw* result, size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata)
 702 {
 703         /* preferred passing of parameters is x the larger of the two numbers */
 704         if (xsize >= ysize)
 705         {
 706                 register mpw rc;
 707
 708                 result += ysize;
 709                 ydata += ysize;
 710
 711                 rc = mpsetmul(xsize, result, xdata, *(--ydata));
 712                 *(--result) = rc;
 713
 714                 while (--ysize)
 715                 {
 716                         rc = mpaddmul(xsize, result, xdata, *(--ydata));
 717                         *(--result) = rc;
 718                 }
 719         }
 720         else
 721         {
 722                 register mpw rc;
 723
 724                 result += xsize;
 725                 xdata += xsize;
 726
 727                 rc = mpsetmul(ysize, result, ydata, *(--xdata));
 728                 *(--result) = rc;
 729
 730                 while (--xsize)
 731                 {
 732                         rc = mpaddmul(ysize, result, ydata, *(--xdata));
 733                         *(--result) = rc;
 734                 }
 735         }
 736 }
 737 #endif
 738
 739 #ifndef ASM_MPADDSQRTRC
 740 void mpaddsqrtrc(size_t size, mpw* result, const mpw* data)
 741 {
 742         #if HAVE_MPDW
 743         register mpdw temp;
 744         register mpw  load, carry = 0;
 745
 746         result += (size << 1);
 747
 748         while (size--)
 749         {
 750                 temp = load = data[size];
 751                 temp *= load;
 752                 temp += carry;
 753                 temp += *(--result);
 754                 *result = (mpw) temp;
 755                 temp >>= MP_WBITS;
 756                 temp += *(--result);
 757                 *result = (mpw) temp;
 758                 carry = (mpw)(temp >> MP_WBITS);
 759         }
 760         #else
 761         register mpw temp, load, carry = 0;
 762
 763         result += (size << 1);
 764
 765         while (size--)
 766         {
 767                 register mphw xlo, xhi;
 768                 register mpw rlo, rhi;
 769
 770                 xlo = (mphw) (temp = data[size]);
 771                 xhi = (mphw) (temp >> MP_HWBITS);
 772
 773                 rlo = (mpw) xlo * xlo;
 774                 rhi = (mpw) xhi * xhi;
 775                 temp = (mpw) xhi * xlo;
 776                 load = rlo;
 777                 rlo += (temp << MP_HWBITS);
 778                 rhi += (temp >> MP_HWBITS) + (load > rlo);
 779                 load = rlo;
 780                 rlo += (temp << MP_HWBITS);
 781                 rhi += (temp >> MP_HWBITS) + (load > rlo);
 782                 load = rlo;
 783                 rlo += carry;
 784                 rhi += (load > rlo);
 785                 load = rlo;
 786                 rlo += *(--result);
 787                 *result = rlo;
 788                 temp = (load > rlo);
 789                 load = rhi;
 790                 rhi += temp;
 791                 carry = (load > rhi);
 792                 load = rhi;
 793                 rhi += *(--result);
 794                 *result = rhi;
 795                 carry += (load > rhi);
 796         }
 797         #endif
 798 }
 799 #endif
 800
 801 #ifndef ASM_MPSQR
 802 void mpsqr(mpw* result, size_t size, const mpw* data)
 803 {
 804         register mpw rc;
 805         register size_t n = size-1;
 806
 807         result += size;
 808         result[n] = 0;
 809
 810         if (n)
 811         {
 812                 rc = mpsetmul(n, result, data, data[n]);
 813                 *(--result) = rc;
 814                 while (--n)
 815                 {
 816                         rc = mpaddmul(n, result, data, data[n]);
 817                         *(--result) = rc;
 818                 }
 819         }
 820
 821         *(--result) = 0;
 822
 823         mpmultwo(size << 1, result);
 824
 825         mpaddsqrtrc(size, result, data);
 826 }
 827 #endif
 828
 829 #ifndef ASM_MPSIZE
 830 size_t mpsize(size_t size, const mpw* data)
 831 {
 832         while (size)
 833         {
 834                 if (*data)
 835                         return size;
 836                 data++;
 837                 size--;
 838         }
 839         return 0;
 840 }
 841 #endif
 842
 843 #ifndef ASM_MPBITS
 844 size_t mpbits(size_t size, const mpw* data)
 845 {
 846         return MP_WORDS_TO_BITS(size) - mpmszcnt(size, data);
 847 }
 848 #endif
 849
 850 #ifndef ASM_MPNORM
 851 size_t mpnorm(size_t size, mpw* data)
 852 {
 853         register size_t shift = mpmszcnt(size, data);
 854         mplshift(size, data, shift);
 855         return shift;
 856 }
 857 #endif
 858
 859 #ifndef ASM_MPDIVTWO
 860 void mpdivtwo(size_t size, mpw* data)
 861 {
 862         register mpw temp, carry = 0;
 863
 864         while (size--)
 865         {
 866                 temp = *data;
 867                 *(data++) = (temp >> 1) | carry;
 868                 carry = (temp << (MP_WBITS-1));
 869         }
 870 }
 871 #endif
 872
 873 #ifndef ASM_MPSDIVTWO
 874 void mpsdivtwo(size_t size, mpw* data)
 875 {
 876         int carry = mpmsbset(size, data);
 877         mpdivtwo(size, data);
 878         if (carry)
 879                 mpsetmsb(size, data);
 880 }
 881 #endif
 882
 883 #ifndef ASM_MPMULTWO
 884 int mpmultwo(size_t size, mpw* data)
 885 {
 886         register mpw temp, carry = 0;
 887
 888         data += size;
 889         while (size--)
 890         {
 891                 temp = *(--data);
 892                 *data = (temp << 1) | carry;
 893                 carry = (temp >> (MP_WBITS-1));
 894         }
 895         return (int) carry;
 896 }
 897 #endif
 898
 899 #ifndef ASM_MPMSZCNT
 900 size_t mpmszcnt(size_t size, const mpw* data)
 901 {
 902         register size_t zbits = 0;
 903         register size_t i = 0;
 904
 905         while (i < size)
 906         {
 907                 register mpw temp = data[i++];
 908                 if (temp)
 909                 {
 910                         while (!(temp & MP_MSBMASK))
 911                         {
 912                                 zbits++;
 913                                 temp <<= 1;
 914                         }
 915                         break;
 916                 }
 917                 else
 918                         zbits += MP_WBITS;
 919         }
 920         return zbits;
 921 }
 922 #endif
 923
 924 #ifndef ASM_MPLSZCNT
 925 size_t mplszcnt(size_t size, const mpw* data)
 926 {
 927         register size_t zbits = 0;
 928
 929         while (size--)
 930         {
 931                 register mpw temp = data[size];
 932                 if (temp)
 933                 {
 934                         while (!(temp & MP_LSBMASK))
 935                         {
 936                                 zbits++;
 937                                 temp >>= 1;
 938                         }
 939                         break;
 940                 }
 941                 else
 942                         zbits += MP_WBITS;
 943         }
 944         return zbits;
 945 }
 946 #endif
 947
 948 #ifndef ASM_MPLSHIFT
 949 void mplshift(size_t size, mpw* data, size_t count)
 950 {
 951         register size_t words = MP_BITS_TO_WORDS(count);
 952
 953         if (words < size)
 954         {
 955                 register short lbits = (short) (count & (MP_WBITS-1));
 956
 957                 /* first do the shifting, then do the moving */
 958                 if (lbits)
 959                 {
 960                         register mpw temp, carry = 0;
 961                         register short rbits = MP_WBITS - lbits;
 962                         register size_t i = size;
 963
 964                         while (i > words)
 965                         {
 966                                 temp = data[--i];
 967                                 data[i] = (temp << lbits) | carry;
 968                                 carry = (temp >> rbits);
 969                         }
 970                 }
 971                 if (words)
 972                 {
 973                         mpmove(size-words, data, data+words);
 974                         mpzero(words, data+size-words);
 975                 }
 976         }
 977         else
 978                 mpzero(size, data);
 979 }
 980 #endif
 981
 982 #ifndef ASM_MPRSHIFT
 983 void mprshift(size_t size, mpw* data, size_t count)
 984 {
 985         register size_t words = MP_BITS_TO_WORDS(count);
 986
 987         if (words < size)
 988         {
 989                 register short rbits = (short) (count & (MP_WBITS-1));
 990
 991                 /* first do the shifting, then do the moving */
 992                 if (rbits)
 993                 {
 994                         register mpw temp, carry = 0;
 995                         register short lbits = MP_WBITS - rbits;
 996                         register size_t i = 0;
 997
 998                         while (i < size-words)
 999                         {
1000                                 temp = data[i];
1001                                 data[i++] = (temp >> rbits) | carry;
1002                                 carry = (temp << lbits);
1003                         }
1004                 }
1005                 if (words)
1006                 {
1007                         mpmove(size-words, data+words, data);
1008                         mpzero(words, data);
1009                 }
1010         }
1011         else
1012                 mpzero(size, data);
1013 }
1014 #endif
1015
1016 #ifndef ASM_MPRSHIFTLSZ
1017 size_t mprshiftlsz(size_t size, mpw* data)
1018 {
1019         register mpw* slide = data+size-1;
1020         register size_t  zwords = 0; /* counter for 'all zero bit' words */
1021         register short   lbits, rbits = 0; /* counter for 'least significant zero' bits */
1022         register mpw  temp, carry = 0;
1023
1024         data = slide;
1025
1026         /* count 'all zero' words and move src pointer */
1027         while (size--)
1028         {
1029                 /* test if we have a non-zero word */
1030                 if ((carry = *(slide--)))
1031                 {
1032                         /* count 'least signification zero bits and set zbits counter */
1033                         while (!(carry & MP_LSBMASK))
1034                         {
1035                                 carry >>= 1;
1036                                 rbits++;
1037                         }
1038                         break;
1039                 }
1040                 zwords++;
1041         }
1042
1043         if ((rbits == 0) && (zwords == 0))
1044                 return 0;
1045
1046         /* prepare right-shifting of data */
1047         lbits = MP_WBITS - rbits;
1048
1049         /* shift data */
1050         while (size--)
1051         {
1052                 temp = *(slide--);
1053                 *(data--) = (temp << lbits) | carry;
1054                 carry = (temp >> rbits);
1055         }
1056
1057         /* store the final carry */
1058         *(data--) = carry;
1059
1060         /* store the return value in size */
1061         size = MP_WORDS_TO_BITS(zwords) + rbits;
1062
1063         /* zero the (zwords) most significant words */
1064         while (zwords--)
1065                 *(data--) = 0;
1066
1067         return size;
1068 }
1069 #endif
1070
1071 /* try an alternate version here, with descending sizes */
1072 /* also integrate lszcnt and rshift properly into one function */
1073 #ifndef ASM_MPGCD_W
1074 /*
1075  * mpgcd_w
1076  *  need workspace of (size) words
1077  */
1078 void mpgcd_w(size_t size, const mpw* xdata, const mpw* ydata, mpw* result, mpw* wksp)
1079 {
1080         register size_t shift, temp;
1081
1082         if (mpge(size, xdata, ydata))
1083         {
1084                 mpcopy(size, wksp, xdata);
1085                 mpcopy(size, result, ydata);
1086         }
1087         else
1088         {
1089                 mpcopy(size, wksp, ydata);
1090                 mpcopy(size, result, xdata);
1091         }
1092
1093         /* get the smallest returned values, and set shift to that */
1094
1095         shift = mprshiftlsz(size, wksp);
1096         temp = mprshiftlsz(size, result);
1097
1098         if (shift > temp)
1099                 shift = temp;
1100
1101         while (mpnz(size, wksp))
1102         {
1103                 mprshiftlsz(size, wksp);
1104                 mprshiftlsz(size, result);
1105
1106                 if (mpge(size, wksp, result))
1107                         mpsub(size, wksp, result);
1108                 else
1109                         mpsub(size, result, wksp);
1110
1111                 /* slide past zero words in both operands by increasing pointers and decreasing size */
1112                 if ((*wksp == 0) && (*result == 0))
1113                 {
1114                         size--;
1115                         wksp++;
1116                         result++;
1117                 }
1118         }
1119
1120         /* figure out if we need to slide the result pointer back */
1121         if ((temp = MP_BITS_TO_WORDS(shift)))
1122         {
1123                 size += temp;
1124                 result -= temp;
1125         }
1126
1127         mplshift(size, result, shift);
1128 }
1129 #endif
1130
1131 #ifndef ASM_MPEXTGCD_W
1132 /* needs workspace of (6*size+6) words */
1133 /* used to compute the modular inverse */
1134 int mpextgcd_w(size_t size, const mpw* xdata, const mpw* ydata, mpw* result, mpw* wksp)
1135 {
1136         /*
1137          * For computing a modular inverse, pass the modulus as xdata and the number
1138          * to be inverted as ydata.
1139          *
1140          * Fact: if a element of Zn, then a is invertible if and only if gcd(a,n) = 1
1141          * Hence: if n is even, then a must be odd, otherwise the gcd(a,n) >= 2
1142          *
1143          * The calling routine must guarantee this condition.
1144          */
1145
1146         register size_t sizep = size+1;
1147         register int full;
1148
1149         mpw* udata = wksp;
1150         mpw* vdata = udata+sizep;
1151         mpw* adata = vdata+sizep;
1152         mpw* bdata = adata+sizep;
1153         mpw* cdata = bdata+sizep;
1154         mpw* ddata = cdata+sizep;
1155
1156         mpsetx(sizep, udata, size, xdata);
1157         mpsetx(sizep, vdata, size, ydata);
1158         mpzero(sizep, bdata);
1159         mpsetw(sizep, ddata, 1);
1160
1161         if ((full = mpeven(sizep, udata)))
1162         {
1163                 mpsetw(sizep, adata, 1);
1164                 mpzero(sizep, cdata);
1165         }
1166
1167         while (1)
1168         {
1169                 while (mpeven(sizep, udata))
1170                 {
1171                         mpdivtwo(sizep, udata);
1172
1173                         if (mpodd(sizep, bdata) || (full && mpodd(sizep, adata)))
1174                         {
1175                                 if (full) mpaddx(sizep, adata, size, ydata);
1176                                 mpsubx(sizep, bdata, size, xdata);
1177                         }
1178
1179                         if (full) mpsdivtwo(sizep, adata);
1180                         mpsdivtwo(sizep, bdata);
1181                 }
1182                 while (mpeven(sizep, vdata))
1183                 {
1184                         mpdivtwo(sizep, vdata);
1185
1186                         if (mpodd(sizep, ddata) || (full && mpodd(sizep, cdata)))
1187                         {
1188                                 if (full) mpaddx(sizep, cdata, size, ydata);
1189                                 mpsubx(sizep, ddata, size, xdata);
1190                         }
1191
1192                         if (full) mpsdivtwo(sizep, cdata);
1193                         mpsdivtwo(sizep, ddata);
1194                 }
1195                 if (mpge(sizep, udata, vdata))
1196                 {
1197                         mpsub(sizep, udata, vdata);
1198                         if (full) mpsub(sizep, adata, cdata);
1199                         mpsub(sizep, bdata, ddata);
1200                 }
1201                 else
1202                 {
1203                         mpsub(sizep, vdata, udata);
1204                         if (full) mpsub(sizep, cdata, adata);
1205                         mpsub(sizep, ddata, bdata);
1206                 }
1207                 if (mpz(sizep, udata))
1208                 {
1209                         if (mpisone(sizep, vdata))
1210                         {
1211                                 if (result)
1212                                 {
1213                                         if (*ddata & MP_MSBMASK)
1214                                         {
1215                                                 /* keep adding the modulus until we get a carry */
1216                                                 while (!mpaddx(sizep, ddata, size, xdata));
1217                                         }
1218                                         else
1219                                         {
1220                                                 /* in some computations, d ends up > x, hence:
1221                                                  * keep subtracting n from d until d < x
1222                                                  */
1223                                                 while (mpgtx(sizep, ddata, size, xdata))
1224                                                         mpsubx(sizep, ddata, size, xdata);
1225                                         }
1226                                         mpsetx(size, result, sizep, ddata);
1227                                 }
1228                                 return 1;
1229                         }
1230                         return 0;
1231                 }
1232         }
1233 }
1234 #endif
1235
1236 #ifndef ASM_MPPNDIV
1237 mpw mppndiv(mpw xhi, mpw xlo, mpw y)
1238 {
1239         register mpw result = 0;
1240         register short count = MP_WBITS;
1241         register int carry = 0;
1242
1243         while (count--)
1244         {
1245                 if (carry | (xhi >= y))
1246                 {
1247                         xhi -= y;
1248                         result++;
1249                 }
1250                 carry = (xhi >> (MP_WBITS-1));
1251                 xhi <<= 1;
1252                 xhi |= (xlo >> (MP_WBITS-1));
1253                 xlo <<= 1;
1254                 result <<= 1;
1255         }
1256         if (carry | (xhi >= y))
1257         {
1258                 xhi -= y;
1259                 result++;
1260         }
1261         return result;
1262 }
1263 #endif
1264
1265 #ifndef ASM_MPMOD
1266 void mpmod(mpw* result, size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata, mpw* workspace)
1267 {
1268         /* result size xsize, workspace size 2*ysize+1 */
1269         mpw q, msw;
1270         mpw* rdata = result;
1271         mpw* ynorm = workspace+ysize+1;
1272         size_t shift, qsize = xsize-ysize;
1273
1274         mpcopy(ysize, ynorm, ydata);
1275         shift = mpnorm(ysize, ynorm);
1276         msw = *ynorm;
1277         mpcopy(xsize, rdata, xdata);
1278         if (mpge(ysize, rdata, ynorm))
1279                 mpsub(ysize, rdata, ynorm);
1280
1281         while (qsize--)
1282         {
1283                 q = mppndiv(rdata[0], rdata[1], msw);
1284
1285                 *workspace = mpsetmul(ysize, workspace+1, ynorm, q);
1286
1287                 while (mplt(ysize+1, rdata, workspace))
1288                 {
1289                         mpsubx(ysize+1, workspace, ysize, ynorm);
1290                         q--;
1291                 }
1292                 mpsub(ysize+1, rdata, workspace);
1293                 rdata++;
1294         }
1295         /* de-normalization steps */
1296         while (shift--)
1297         {
1298                 mpdivtwo(ysize, ynorm);
1299                 if (mpge(ysize, rdata, ynorm))
1300                         mpsub(ysize, rdata, ynorm);
1301         }
1302 }
1303 #endif
1304
1305 #ifndef ASM_MPNDIVMOD
1306 void mpndivmod(mpw* result, size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata, register mpw* workspace)
1307 {
1308         /* result must be xsize+1 in length */
1309         /* workspace must be ysize+1 in length */
1310         /* expect ydata to be normalized */
1311         mpw q;
1312         mpw msw = *ydata;
1313         size_t qsize = xsize-ysize;
1314
1315         *result = (mpge(ysize, xdata, ydata) ? 1 : 0);
1316         mpcopy(xsize, result+1, xdata);
1317
1318         if (*result)
1319                 (void) mpsub(ysize, result+1, ydata);
1320
1321         result++;
1322
1323         while (qsize--)
1324         {
1325                 q = mppndiv(result[0], result[1], msw);
1326
1327                 *workspace = mpsetmul(ysize, workspace+1, ydata, q);
1328
1329                 while (mplt(ysize+1, result, workspace))
1330                 {
1331                         mpsubx(ysize+1, workspace, ysize, ydata);
1332                         q--;
1333                 }
1334                 mpsub(ysize+1, result, workspace);
1335                 *(result++) = q;
1336         }
1337 }
1338 #endif
1339
1340 void mpprint(size_t size, const mpw* data)
1341 {
1342         mpfprint(stdout, size, data);
1343 }
1344
1345 void mpprintln(size_t size, const mpw* data)
1346 {
1347         mpfprintln(stdout, size, data);
1348 }
1349
1350 void mpfprint(FILE* f, size_t size, const mpw* data)
1351 {
1352         if (data == (mpw*) 0)
1353                 return;
1354
1355         if (f == (FILE*) 0)
1356                 return;
1357
1358         while (size--)
1359         {
1360                 #if (MP_WBITS == 32)
1361                 fprintf(f, "%08x", (unsigned) *(data++));
1362                 #elif (MP_WBITS == 64)
1363                 # if WIN32
1364                 fprintf(f, "%016I64x", *(data++));
1365                 # elif SIZEOF_UNSIGNED_LONG == 8
1366                 fprintf(f, "%016lx", *(data++));
1367                 # else
1368                 fprintf(f, "%016llx", *(data++));
1369                 # endif
1370                 #else
1371                 # error
1372                 #endif
1373         }
1374         fflush(f);
1375 }
1376
1377 void mpfprintln(FILE* f, size_t size, const mpw* data)
1378 {
1379         if (data == (mpw*) 0)
1380                 return;
1381
1382         if (f == (FILE*) 0)
1383                 return;
1384
1385         while (size--)
1386         {
1387                 #if (MP_WBITS == 32)
1388                 fprintf(f, "%08x", *(data++));
1389                 #elif (MP_WBITS == 64)
1390                 # if WIN32
1391                 fprintf(f, "%016I64x", *(data++));
1392                 # elif SIZEOF_UNSIGNED_LONG == 8
1393                 fprintf(f, "%016lx", *(data++));
1394                 # else
1395                 fprintf(f, "%016llx", *(data++));
1396                 # endif
1397                 #else
1398                 # error
1399                 #endif
1400         }
1401         fprintf(f, "\n");
1402         fflush(f);
1403 }
1404
1405 int i2osp(byte *osdata, size_t ossize, const mpw* idata, size_t isize)
1406 {
1407         #if WORDS_BIGENDIAN
1408         size_t max_bytes = MP_WORDS_TO_BYTES(isize);
1409         #endif
1410         size_t significant_bytes = (mpbits(isize, idata) + 7) >> 3;
1411
1412         /* verify that ossize is large enough to contain the significant bytes */
1413         if (ossize >= significant_bytes)
1414         {
1415                 /* looking good; check if we have more space than significant bytes */
1416                 if (ossize > significant_bytes)
1417                 {       /* fill most significant bytes with zero */
1418                         memset(osdata, 0, ossize - significant_bytes);
1419                         osdata += ossize - significant_bytes;
1420                 }
1421                 if (significant_bytes)
1422                 {       /* fill remaining bytes with endian-adjusted data */
1423                         #if !WORDS_BIGENDIAN
1424                         mpw w = idata[--isize];
1425                         byte shift = 0;
1426
1427                         /* fill right-to-left; much easier than left-to-right */
1428                         do
1429                         {
1430                                 osdata[--significant_bytes] = (byte)(w >> shift);
1431                                 shift += 8;
1432                                 if (shift == MP_WBITS)
1433                                 {
1434                                         shift = 0;
1435                                         w = idata[--isize];
1436                                 }
1437                         } while (significant_bytes);
1438                         #else
1439                         /* just copy data past zero bytes */
1440                         memcpy(osdata, ((byte*) idata) + (max_bytes - significant_bytes), significant_bytes);
1441                         #endif
1442                 }
1443                 return 0;
1444         }
1445         return -1;
1446 }
1447
1448 int os2ip(mpw* idata, size_t isize, const byte* osdata, size_t ossize)
1449 {
1450         size_t required;
1451
1452         /* skip non-significant leading zero bytes */
1453         while (!(*osdata) && ossize)
1454         {
1455                 osdata++;
1456                 ossize--;
1457         }
1458
1459         required = MP_BYTES_TO_WORDS(ossize + MP_WBYTES - 1);
1460
1461         if (isize >= required)
1462         {
1463                 /* yes, we have enough space and can proceed */
1464                 mpw w = 0;
1465                 /* adjust counter so that the loop will start by skipping the proper
1466                  * amount of leading bytes in the first significant word
1467                  */
1468                 byte b = (ossize % MP_WBYTES);
1469
1470                 if (isize > required)
1471                 {       /* fill initials words with zero */
1472                         mpzero(isize-required, idata);
1473                         idata += isize-required;
1474                 }
1475
1476                 if (b == 0)
1477                         b = MP_WBYTES;
1478
1479                 while (ossize--)
1480                 {
1481                         w <<= 8;
1482                         w |= *(osdata++);
1483                         b--;
1484
1485                         if (b == 0)
1486                         {
1487                                 *(idata++) = w;
1488                                 w = 0;
1489                                 b = MP_WBYTES;
1490                         }
1491                 }
1492
1493                 return 0;
1494         }
1495         return -1;
1496 }
1497
1498 int hs2ip(mpw* idata, size_t isize, const char* hsdata, size_t hssize)
1499 {
1500         size_t required = MP_NIBBLES_TO_WORDS(hssize + MP_WNIBBLES - 1);
1501
1502         if (isize >= required)
1503         {
1504                 register size_t i;
1505
1506
1507                 if (isize > required)
1508                 {       /* fill initial words with zero */
1509                         for (i = required; i < isize; i++)
1510                                 *(idata++) = 0;
1511                 }
1512                 while (hssize)
1513                 {
1514                         register mpw w = 0;
1515                         register size_t chunk = hssize & (MP_WNIBBLES - 1);
1516                         register char ch;
1517
1518                         if (chunk == 0) chunk = MP_WNIBBLES;
1519
1520                         for (i = 0; i < chunk; i++)
1521                         {
1522                                 ch = *(hsdata++);
1523                                 w <<= 4;
1524                                 if (ch >= '0' && ch <= '9')
1525                                         w += (ch - '0');
1526                                 else if (ch >= 'A' && ch <= 'F')
1527                                         w += (ch - 'A') + 10;
1528                                 else if (ch >= 'a' && ch <= 'f')
1529                                         w += (ch - 'a') + 10;
1530                         }
1531                         *(idata++) = w;
1532                         hssize -= chunk;
1533                 }
1534                 return 0;
1535         }
1536         return -1;
1537 }