Main Page Modules Data Structures File List Data Fields Globals Related Pages

# beecrypt/mp32barrett.c

Go to the documentation of this file.

    00001 /*@-sizeoftype -type@*/
    00013 /*
    00014  * Copyright (c) 1997, 1998, 1999, 2000, 2001 Virtual Unlimited B.V.
    00015  *
    00016  * Author: Bob Deblier <bob@virtualunlimited.com>
    00017  *
    00018  * This library is free software; you can redistribute it and/or
    00019  * modify it under the terms of the GNU Lesser General Public
    00020  * License as published by the Free Software Foundation; either
    00021  * version 2.1 of the License, or (at your option) any later version.
    00022  *
    00023  * This library is distributed in the hope that it will be useful,
    00024  * but WITHOUT ANY WARRANTY; without even the implied warranty of
    00025  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    00026  * Lesser General Public License for more details.
    00027  *
    00028  * You should have received a copy of the GNU Lesser General Public
    00029  * License along with this library; if not, write to the Free Software
    00030  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
    00031  *
    00032  */
    00033
    00034 #include "system.h"
    00035 #include "mp32.h"
    00036 #include "mp32prime.h"
    00037 #include "mp32barrett.h"
    00038 #include "debug.h"
    00039
    00043 void mp32bzero(mp32barrett* b)
    00044 {
    00045         b->size = 0;
    00046         b->modl = (uint32*) 0;
    00047         b->mu = (uint32*) 0;
    00048 }
    00049
    00050 /*@-nullstate@*/        /* b->modl may be null @*/
    00055 void mp32binit(mp32barrett* b, uint32 size)
    00056 {
    00057         b->size = size;
    00058         if (b->modl)
    00059                 free(b->modl);
    00060         b->modl = (uint32*) calloc(2*size+1, sizeof(uint32));
    00061
    00062         if (b->modl != (uint32*) 0)
    00063                 b->mu = b->modl+size;
    00064         else
    00065                 b->mu = (uint32*) 0;
    00066 }
    00067 /*@=nullstate@*/
    00068
    00072 void mp32bfree(mp32barrett* b)
    00073 {
    00074         if (b->modl != (uint32*) 0)
    00075         {
    00076                 free(b->modl);
    00077                 b->modl = (uint32*) 0;
    00078                 b->mu = (uint32*) 0;
    00079         }
    00080         b->size = 0;
    00081 }
    00082
    00083 /*@-boundswrite@*/
    00084 /*@-nullstate -compdef @*/      /* b->modl may be null @*/
    00085 void mp32bcopy(mp32barrett* b, const mp32barrett* copy)
    00086 {
    00087         register uint32 size = copy->size;
    00088
    00089         if (size)
    00090         {
    00091                 if (b->modl)
    00092                 {
    00093                         if (b->size != size)
    00094                                 b->modl = (uint32*) realloc(b->modl, (2*size+1) * sizeof(uint32));
    00095                 }
    00096                 else
    00097                         b->modl = (uint32*) malloc((2*size+1) * sizeof(uint32));
    00098
    00099                 if (b->modl)
    00100                 {
    00101                         b->size = size;
    00102                         b->mu = b->modl+copy->size;
    00103                         mp32copy(2*size+1, b->modl, copy->modl);
    00104                 }
    00105                 else
    00106                 {
    00107                         b->size = 0;
    00108                         b->mu = (uint32*) 0;
    00109                 }
    00110         }
    00111         else if (b->modl)
    00112         {
    00113                 free(b->modl);
    00114                 b->size = 0;
    00115                 b->modl = (uint32*) 0;
    00116                 b->mu = (uint32*) 0;
    00117         }
    00118         else
    00119                 {};
    00120 }
    00121 /*@=nullstate =compdef @*/
    00122 /*@=boundswrite@*/
    00123
    00124 /*@-boundswrite@*/
    00125 /*@-nullstate -compdef @*/      /* b->modl may be null @*/
    00129 void mp32bset(mp32barrett* b, uint32 size, const uint32 *data)
    00130 {
    00131         if (size > 0)
    00132         {
    00133                 if (b->modl)
    00134                 {
    00135                         if (b->size != size)
    00136                                 b->modl = (uint32*) realloc(b->modl, (2*size+1) * sizeof(uint32));
    00137                 }
    00138                 else
    00139                         b->modl = (uint32*) malloc((2*size+1) * sizeof(uint32));
    00140
    00141                 if (b->modl)
    00142                 {
    00143                         uint32* temp = (uint32*) malloc((6*size+4) * sizeof(uint32));
    00144
    00145                         b->size = size;
    00146                         b->mu = b->modl+size;
    00147                         mp32copy(size, b->modl, data);
    00148                         /*@-nullpass@*/         /* temp may be NULL */
    00149                         mp32bmu_w(b, temp);
    00150
    00151                         free(temp);
    00152                         /*@=nullpass@*/
    00153                 }
    00154                 else
    00155                 {
    00156                         b->size = 0;
    00157                         b->mu = (uint32*) 0;
    00158                 }
    00159         }
    00160 }
    00161 /*@=nullstate =compdef @*/
    00162 /*@=boundswrite@*/
    00163
    00164 /*@-boundswrite@*/
    00165 /*@-nullstate -compdef @*/      /* b->modl may be null @*/
    00166 void mp32bsethex(mp32barrett* b, const char* hex)
    00167 {
    00168         uint32 length = strlen(hex);
    00169         uint32 size = (length+7) >> 3;
    00170         uint8 rem = (uint8)(length & 0x7);
    00171
    00172         if (b->modl)
    00173         {
    00174                 if (b->size != size)
    00175                         b->modl = (uint32*) realloc(b->modl, (2*size+1) * sizeof(uint32));
    00176         }
    00177         else
    00178                 b->modl = (uint32*) malloc((2*size+1) * sizeof(uint32));
    00179
    00180         if (b->modl != (uint32*) 0)
    00181         {
    00182                 register uint32  val = 0;
    00183                 register uint32* dst = b->modl;
    00184                 register uint32* temp = (uint32*) malloc((6*size+4) * sizeof(uint32));
    00185                 register char ch;
    00186
    00187                 b->size = size;
    00188                 b->mu = b->modl+size;
    00189
    00190                 while (length-- > 0)
    00191                 {
    00192                         ch = *(hex++);
    00193                         val <<= 4;
    00194                         if (ch >= '0' && ch <= '9')
    00195                                 val += (ch - '0');
    00196                         else if (ch >= 'A' && ch <= 'F')
    00197                                 val += (ch - 'A') + 10;
    00198                         else if (ch >= 'a' && ch <= 'f')
    00199                                 val += (ch - 'a') + 10;
    00200                         else
    00201                                 {};
    00202
    00203                         if ((length & 0x7) == 0)
    00204                         {
    00205                                 *(dst++) = val;
    00206                                 val = 0;
    00207                         }
    00208                 }
    00209                 if (rem != 0)
    00210                         *dst = val;
    00211
    00212                 /*@-nullpass@*/         /* temp may be NULL */
    00213                 mp32bmu_w(b, temp);
    00214
    00215                 free(temp);
    00216                 /*@=nullpass@*/
    00217         }
    00218         else
    00219         {
    00220                 b->size = 0;
    00221                 b->mu = 0;
    00222         }
    00223 }
    00224 /*@=nullstate =compdef @*/
    00225 /*@=boundswrite@*/
    00226
    00231 /*@-boundswrite@*/
    00232 void mp32bmu_w(mp32barrett* b, uint32* wksp)
    00233 {
    00234         register uint32  size = b->size;
    00235         register uint32* divmod = wksp;
    00236         register uint32* dividend = divmod+(size*2+2);
    00237         register uint32* workspace = dividend+(size*2+1);
    00238         register uint32  shift;
    00239
    00240         /* normalize modulus before division */
    00241         shift = mp32norm(size, b->modl);
    00242         /* make the dividend, initialize first word to 1 (shifted); the rest is zero */
    00243         *dividend = (uint32) (1 << shift);
    00244         mp32zero(size*2, dividend+1);
    00245         mp32ndivmod(divmod, size*2+1, dividend, size, b->modl, workspace);
    00246         /*@-nullpass@*/ /* b->mu may be NULL */
    00247         mp32copy(size+1, b->mu, divmod+1);
    00248         /*@=nullpass@*/
    00249         /* de-normalize */
    00250         mp32rshift(size, b->modl, shift);
    00251 }
    00252 /*@=boundswrite@*/
    00253
    00258 /*@-boundswrite@*/
    00259 void mp32brnd_w(const mp32barrett* b, randomGeneratorContext* rc, uint32* result, uint32* wksp)
    00260 {
    00261         uint32 msz = mp32mszcnt(b->size, b->modl);
    00262
    00263         mp32copy(b->size, wksp, b->modl);
    00264         (void) mp32subw(b->size, wksp, 1);
    00265
    00266         do
    00267         {
    00268                 /*@-noeffectuncon@*/ /* LCL: ??? */
    00269                 (void) rc->rng->next(rc->param, result, b->size);
    00270                 /*@=noeffectuncon@*/
    00271
    00272                 /*@-shiftimplementation -usedef@*/
    00273                 result[0] &= (0xffffffff >> msz);
    00274                 /*@=shiftimplementation =usedef@*/
    00275
    00276                 while (mp32ge(b->size, result, wksp))
    00277                         (void) mp32sub(b->size, result, wksp);
    00278         } while (mp32leone(b->size, result));
    00279 }
    00280 /*@=boundswrite@*/
    00281
    00286 /*@-boundswrite@*/
    00287 void mp32brndodd_w(const mp32barrett* b, randomGeneratorContext* rc, uint32* result, uint32* wksp)
    00288 {
    00289         uint32 msz = mp32mszcnt(b->size, b->modl);
    00290
    00291         mp32copy(b->size, wksp, b->modl);
    00292         (void) mp32subw(b->size, wksp, 1);
    00293
    00294         do
    00295         {
    00296                 /*@-noeffectuncon@*/ /* LCL: ??? */
    00297                 (void) rc->rng->next(rc->param, result, b->size);
    00298                 /*@=noeffectuncon@*/
    00299
    00300                 /*@-shiftimplementation -usedef@*/
    00301                 result[0] &= (0xffffffff >> msz);
    00302                 /*@=shiftimplementation =usedef@*/
    00303                 mp32setlsb(b->size, result);
    00304
    00305                 while (mp32ge(b->size, result, wksp))
    00306                 {
    00307                         (void) mp32sub(b->size, result, wksp);
    00308                         mp32setlsb(b->size, result);
    00309                 }
    00310         } while (mp32leone(b->size, result));
    00311 }
    00312 /*@=boundswrite@*/
    00313
    00318 void mp32brndinv_w(const mp32barrett* b, randomGeneratorContext* rc, uint32* result, uint32* inverse, uint32* wksp)
    00319 {
    00320         register uint32  size = b->size;
    00321
    00322         do
    00323         {
    00324                 if (mp32even(size, b->modl))
    00325                         mp32brndodd_w(b, rc, result, wksp);
    00326                 else
    00327                         mp32brnd_w(b, rc, result, wksp);
    00328
    00329         } while (mp32binv_w(b, size, result, inverse, wksp) == 0);
    00330 }
    00331
    00336 /*@-boundswrite@*/
    00337 void mp32bmod_w(const mp32barrett* b, const uint32* xdata, uint32* result, uint32* wksp)
    00338 {
    00339         register uint32 rc;
    00340         register uint32 sp = 2;
    00341         register const uint32* src = xdata+b->size+1;
    00342         register       uint32* dst = wksp +b->size+1;
    00343
    00344         /*@-nullpass@*/ /* b->mu may be NULL */
    00345         rc = mp32setmul(sp, dst, b->mu, *(--src));
    00346         *(--dst) = rc;
    00347
    00348         while (sp <= b->size)
    00349         {
    00350                 sp++;
    00351                 if ((rc = *(--src)))
    00352                 {
    00353                         rc = mp32addmul(sp, dst, b->mu, rc);
    00354                         *(--dst) = rc;
    00355                 }
    00356                 else
    00357                         *(--dst) = 0;
    00358         }
    00359         if ((rc = *(--src)))
    00360         {
    00361                 rc = mp32addmul(sp, dst, b->mu, rc);
    00362                 *(--dst) = rc;
    00363         }
    00364         else
    00365                 *(--dst) = 0;
    00366         /*@=nullpass@*/
    00367
    00368         /* q3 is one word larger than b->modl */
    00369         /* r2 is (2*size+1) words, of which we only needs the (size+1) lsw's */
    00370
    00371         sp = b->size;
    00372         rc = 0;
    00373
    00374         dst = wksp+b->size+1;
    00375         src = dst;
    00376
    00377         /*@-evalorder@*/ /* --src side effect, dst/src aliases */
    00378         *dst = mp32setmul(sp, dst+1, b->modl, *(--src));
    00379         /*@=evalorder@*/
    00380
    00381         while (sp > 0)
    00382         {
    00383                 (void) mp32addmul(sp--, dst, b->modl+(rc++), *(--src));
    00384         }
    00385
    00386         mp32setx(b->size+1, wksp, b->size*2, xdata);
    00387         (void) mp32sub(b->size+1, wksp, wksp+b->size+1);
    00388
    00389         while (mp32gex(b->size+1, wksp, b->size, b->modl))
    00390         {
    00391                 (void) mp32subx(b->size+1, wksp, b->size, b->modl);
    00392         }
    00393         mp32copy(b->size, result, wksp+1);
    00394 }
    00395 /*@=boundswrite@*/
    00396
    00400 /*@-boundswrite@*/
    00401 void mp32bsubone(const mp32barrett* b, uint32* result)
    00402 {
    00403         register uint32 size = b->size;
    00404
    00405         mp32copy(size, result, b->modl);
    00406         (void) mp32subw(size, result, 1);
    00407 }
    00408 /*@=boundswrite@*/
    00409
    00413 /*@-boundswrite@*/
    00414 void mp32bneg(const mp32barrett* b, const uint32* xdata, uint32* result)
    00415 {
    00416         register uint32  size = b->size;
    00417
    00418         mp32copy(size, result, xdata);
    00419         mp32neg(size, result);
    00420         (void) mp32add(size, result, b->modl);
    00421 }
    00422 /*@=boundswrite@*/
    00423
    00428 void mp32baddmod_w(const mp32barrett* b, uint32 xsize, const uint32* xdata, uint32 ysize, const uint32* ydata, uint32* result, uint32* wksp)
    00429 {
    00430         /* xsize and ysize must be less than or equal to b->size */
    00431         register uint32  size = b->size;
    00432         register uint32* temp = wksp + size*2+2;
    00433
    00434         mp32setx(2*size, temp, xsize, xdata);
    00435         (void) mp32addx(2*size, temp, ysize, ydata);
    00436
    00437         mp32bmod_w(b, temp, result, wksp);
    00438 }
    00439
    00444 void mp32bsubmod_w(const mp32barrett* b, uint32 xsize, const uint32* xdata, uint32 ysize, const uint32* ydata, uint32* result, uint32* wksp)
    00445 {
    00446         /* xsize and ysize must be less than or equal to b->size */
    00447         register uint32  size = b->size;
    00448         register uint32* temp = wksp + size*2+2;
    00449
    00450         mp32setx(2*size, temp, xsize, xdata);
    00451         if (mp32subx(2*size, temp, ysize, ydata)) /* if there's carry, i.e. the result would be negative, add the modulus */
    00452                 (void) mp32addx(2*size, temp, size, b->modl);
    00453
    00454         mp32bmod_w(b, temp, result, wksp);
    00455 }
    00456
    00461 void mp32bmulmod_w(const mp32barrett* b, uint32 xsize, const uint32* xdata, uint32 ysize, const uint32* ydata, uint32* result, uint32* wksp)
    00462 {
    00463         /* xsize and ysize must be <= b->size */
    00464         register uint32  size = b->size;
    00465         register uint32* temp = wksp + size*2+2;
    00466         register uint32  fill = size*2-xsize-ysize;
    00467
    00468         if (fill)
    00469                 mp32zero(fill, temp);
    00470
    00471         mp32mul(temp+fill, xsize, xdata, ysize, ydata);
    00472         /*@-compdef@*/  /* *temp undefined */
    00473         mp32bmod_w(b, temp, result, wksp);
    00474         /*@=compdef@*/
    00475 }
    00476
    00481 void mp32bsqrmod_w(const mp32barrett* b, uint32 xsize, const uint32* xdata, uint32* result, uint32* wksp)
    00482 {
    00483         /* xsize must be <= b->size */
    00484         register uint32  size = b->size;
    00485         register uint32* temp = wksp + size*2+2;
    00486         register uint32  fill = 2*(size-xsize);
    00487
    00488         if (fill)
    00489                 mp32zero(fill, temp);
    00490
    00491         mp32sqr(temp+fill, xsize, xdata);
    00492         /*@-compdef@*/  /* *temp undefined */
    00493         mp32bmod_w(b, temp, result, wksp);
    00494         /*@=compdef@*/
    00495 }
    00496
    00534 static void mp32bslide_w(const mp32barrett* b, const uint32 xsize, const uint32* xdata, /*@out@*/ uint32* slide, /*@out@*/ uint32* wksp)
    00535         /*@modifies slide, wksp @*/
    00536 {
    00537         register uint32 size = b->size;
    00538         mp32bsqrmod_w(b, xsize, xdata,                     slide       , wksp); /* x^2 mod b, temp */
    00539         mp32bmulmod_w(b, xsize, xdata, size, slide       , slide+size  , wksp); /* x^3 mod b */
    00540         mp32bmulmod_w(b,  size, slide, size, slide+size  , slide+2*size, wksp); /* x^5 mod b */
    00541         mp32bmulmod_w(b,  size, slide, size, slide+2*size, slide+3*size, wksp); /* x^7 mod b */
    00542         mp32bmulmod_w(b,  size, slide, size, slide+3*size, slide+4*size, wksp); /* x^9 mod b */
    00543         mp32bmulmod_w(b,  size, slide, size, slide+4*size, slide+5*size, wksp); /* x^11 mod b */
    00544         mp32bmulmod_w(b,  size, slide, size, slide+5*size, slide+6*size, wksp); /* x^13 mod b */
    00545         mp32bmulmod_w(b,  size, slide, size, slide+6*size, slide+7*size, wksp); /* x^15 mod b */
    00546         mp32setx(size, slide, xsize, xdata);                                    /* x^1 mod b */
    00547 }
    00548
    00549 /*@observer@*/ /*@unchecked@*/
    00550 static byte mp32bslide_presq[16] =
    00551 { 0, 1, 1, 2, 1, 3, 2, 3, 1, 4, 3, 4, 2, 4, 3, 4 };
    00552
    00553 /*@observer@*/ /*@unchecked@*/
    00554 static byte mp32bslide_mulg[16] =
    00555 { 0, 0, 0, 1, 0, 2, 1, 3, 0, 4, 2, 5, 1, 6, 3, 7 };
    00556
    00557 /*@observer@*/ /*@unchecked@*/
    00558 static byte mp32bslide_postsq[16] =
    00559 { 0, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0 };
    00560
    00565 /*@-boundsread@*/
    00566 void mp32bpowmod_w(const mp32barrett* b, uint32 xsize, const uint32* xdata, uint32 psize, const uint32* pdata, uint32* result, uint32* wksp)
    00567 {
    00568         /*
    00569          * Modular exponention
    00570          *
    00571          * Uses sliding window exponentiation; needs extra storage: if K=3, needs 8*size, if K=4, needs 16*size
    00572          *
    00573          */
    00574
    00575         /* K == 4 for the first try */
    00576
    00577         uint32  size = b->size;
    00578         uint32  temp = 0;
    00579
    00580         while (psize)
    00581         {
    00582                 if ((temp = *(pdata++))) /* break when first non-zero word found */
    00583                         break;
    00584                 psize--;
    00585         }
    00586
    00587         /* if temp is still zero, then we're trying to raise x to power zero, and result stays one */
    00588         if (temp)
    00589         {
    00590                 uint32* slide = (uint32*) malloc((8*size)*sizeof(uint32));
    00591
    00592                 /*@-nullpass@*/         /* slide may be NULL */
    00593                 mp32bslide_w(b, xsize, xdata, slide, wksp);
    00594
    00595                 /*@-internalglobs -mods@*/ /* noisy */
    00596                 mp32bpowmodsld_w(b, slide, psize, pdata-1, result, wksp);
    00597                 /*@=internalglobs =mods@*/
    00598
    00599                 free(slide);
    00600                 /*@=nullpass@*/
    00601         }
    00602 }
    00603 /*@=boundsread@*/
    00604
    00605 /*@-boundsread@*/
    00606 void mp32bpowmodsld_w(const mp32barrett* b, const uint32* slide, uint32 psize, const uint32* pdata, uint32* result, uint32* wksp)
    00607 {
    00608         /*
    00609          * Modular exponentiation with precomputed sliding window table, so no x is required
    00610          *
    00611          */
    00612
    00613         uint32 size = b->size;
    00614         uint32 temp = 0;
    00615
    00616         mp32setw(size, result, 1);
    00617
    00618         while (psize)
    00619         {
    00620                 if ((temp = *(pdata++))) /* break when first non-zero word found in power */
    00621                         break;
    00622                 psize--;
    00623         }
    00624
    00625         /*@+charindex@*/
    00626         /* if temp is still zero, then we're trying to raise x to power zero, and result stays one */
    00627         if (temp)
    00628         {
    00629                 uint8 l = 0, n = 0, count = 32;
    00630
    00631                 /* first skip bits until we reach a one */
    00632                 while (count != 0)
    00633                 {
    00634                         if (temp & 0x80000000)
    00635                                 break;
    00636                         temp <<= 1;
    00637                         count--;
    00638                 }
    00639
    00640                 while (psize)
    00641                 {
    00642                         while (count != 0)
    00643                         {
    00644                                 uint8 bit = (temp & 0x80000000) != 0;
    00645
    00646                                 n <<= 1;
    00647                                 n += bit;
    00648
    00649                                 if (n != 0)
    00650                                 {
    00651                                         if (l != 0)
    00652                                                 l++;
    00653                                         else if (bit != 0)
    00654                                                 l = 1;
    00655                                         else
    00656                                                 {};
    00657
    00658                                         if (l == 4)
    00659                                         {
    00660                                                 uint8 s = mp32bslide_presq[n];
    00661
    00662                                                 while (s--)
    00663                                                         mp32bsqrmod_w(b, size, result, result, wksp);
    00664
    00665                                                 mp32bmulmod_w(b, size, result, size, slide+mp32bslide_mulg[n]*size, result, wksp);
    00666
    00667                                                 s = mp32bslide_postsq[n];
    00668
    00669                                                 while (s--)
    00670                                                         mp32bsqrmod_w(b, size, result, result, wksp);
    00671
    00672                                                 l = n = 0;
    00673                                         }
    00674                                 }
    00675                                 else
    00676                                         mp32bsqrmod_w(b, size, result, result, wksp);
    00677
    00678                                 temp <<= 1;
    00679                                 count--;
    00680                         }
    00681                         if (--psize)
    00682                         {
    00683                                 count = 32;
    00684                                 temp = *(pdata++);
    00685                         }
    00686                 }
    00687
    00688                 if (n != 0)
    00689                 {
    00690                         uint8 s = mp32bslide_presq[n];
    00691                         while (s--)
    00692                                 mp32bsqrmod_w(b, size, result, result, wksp);
    00693
    00694                         mp32bmulmod_w(b, size, result, size, slide+mp32bslide_mulg[n]*size, result, wksp);
    00695
    00696                         s = mp32bslide_postsq[n];
    00697
    00698                         while (s--)
    00699                                 mp32bsqrmod_w(b, size, result, result, wksp);
    00700                 }
    00701         }
    00702         /*@=charindex@*/
    00703 }
    00704 /*@=boundsread@*/
    00705
    00710 /*@-boundsread@*/
    00711 void mp32btwopowmod_w(const mp32barrett* b, uint32 psize, const uint32* pdata, uint32* result, uint32* wksp)
    00712 {
    00713         /*
    00714          * Modular exponention, 2^p mod modulus, special optimization
    00715          *
    00716          * Uses left-to-right exponentiation; needs no extra storage
    00717          *
    00718          */
    00719
    00720         /* this routine calls mp32bmod, which needs (size*2+2), this routine needs (size*2) for sdata */
    00721
    00722         register uint32 size = b->size;
    00723         register uint32 temp = 0;
    00724
    00725         mp32setw(size, result, 1);
    00726
    00727         while (psize)
    00728         {
    00729                 if ((temp = *(pdata++))) /* break when first non-zero word found */
    00730                         break;
    00731                 psize--;
    00732         }
    00733
    00734         /* if temp is still zero, then we're trying to raise x to power zero, and result stays one */
    00735         if (temp)
    00736         {
    00737                 register int count = 32;
    00738
    00739                 /* first skip bits until we reach a one */
    00740                 while (count)
    00741                 {
    00742                         if (temp & 0x80000000)
    00743                                 break;
    00744                         temp <<= 1;
    00745                         count--;
    00746                 }
    00747
    00748                 while (psize--)
    00749                 {
    00750                         while (count)
    00751                         {
    00752                                 /* always square */
    00753                                 mp32bsqrmod_w(b, size, result, result, wksp);
    00754
    00755                                 /* multiply by two if bit is 1 */
    00756                                 if (temp & 0x80000000)
    00757                                 {
    00758                                         if (mp32add(size, result, result) || mp32ge(size, result, b->modl))
    00759                                         {
    00760                                                 /* there was carry, or the result is greater than the modulus, so we need to adjust */
    00761                                                 (void) mp32sub(size, result, b->modl);
    00762                                         }
    00763                                 }
    00764
    00765                                 temp <<= 1;
    00766                                 count--;
    00767                         }
    00768                         count = 32;
    00769                         temp = *(pdata++);
    00770                 }
    00771         }
    00772 }
    00773 /*@=boundsread@*/
    00774
    00775 #ifdef  DYING
    00776
    00781 int mp32binv_w(const mp32barrett* b, uint32 xsize, const uint32* xdata, uint32* result, uint32* wksp)
    00782 {
    00783         /*
    00784          * Fact: if a element of Zn, then a is invertible if and only if gcd(a,n) = 1
    00785          * Hence: if b->modl is even, then x must be odd, otherwise the gcd(x,n) >= 2
    00786          *
    00787          * The calling routine must guarantee this condition.
    00788          */
    00789
    00790         register uint32  size = b->size;
    00791
    00792         uint32* udata = wksp;
    00793         uint32* vdata = udata+size+1;
    00794         uint32* adata = vdata+size+1;
    00795         uint32* bdata = adata+size+1;
    00796         uint32* cdata = bdata+size+1;
    00797         uint32* ddata = cdata+size+1;
    00798
    00799         mp32setx(size+1, udata, size, b->modl);
    00800         mp32setx(size+1, vdata, xsize, xdata);
    00801         mp32zero(size+1, bdata);
    00802         mp32setw(size+1, ddata, 1);
    00803
    00804         if (mp32odd(size, b->modl))
    00805         {
    00806                 /* use simplified binary extended gcd algorithm */
    00807
    00808                 while (1)
    00809                 {
    00810                         while (mp32even(size+1, udata))
    00811                         {
    00812                                 mp32divtwo(size+1, udata);
    00813
    00814                                 if (mp32odd(size+1, bdata))
    00815                                         (void) mp32subx(size+1, bdata, size, b->modl);
    00816
    00817                                 mp32sdivtwo(size+1, bdata);
    00818                         }
    00819                         while (mp32even(size+1, vdata))
    00820                         {
    00821                                 mp32divtwo(size+1, vdata);
    00822
    00823                                 if (mp32odd(size+1, ddata))
    00824                                         (void) mp32subx(size+1, ddata, size, b->modl);
    00825
    00826                                 mp32sdivtwo(size+1, ddata);
    00827                         }
    00828                         if (mp32ge(size+1, udata, vdata))
    00829                         {
    00830                                 (void) mp32sub(size+1, udata, vdata);
    00831                                 (void) mp32sub(size+1, bdata, ddata);
    00832                         }
    00833                         else
    00834                         {
    00835                                 (void) mp32sub(size+1, vdata, udata);
    00836                                 (void) mp32sub(size+1, ddata, bdata);
    00837                         }
    00838
    00839                         if (mp32z(size+1, udata))
    00840                         {
    00841                                 if (mp32isone(size+1, vdata))
    00842                                 {
    00843                                         if (result)
    00844                                         {
    00845                                                 mp32setx(size, result, size+1, ddata);
    00846                                                 /*@-usedef@*/
    00847                                                 if (*ddata & 0x80000000)
    00848                                                 {
    00849                                                         /* keep adding the modulus until we get a carry */
    00850                                                         while (!mp32add(size, result, b->modl));
    00851                                                 }
    00852                                                 /*@=usedef@*/
    00853                                         }
    00854                                         return 1;
    00855                                 }
    00856                                 return 0;
    00857                         }
    00858                 }
    00859         }
    00860         else
    00861         {
    00862                 /* use full binary extended gcd algorithm */
    00863                 mp32setw(size+1, adata, 1);
    00864                 mp32zero(size+1, cdata);
    00865
    00866                 while (1)
    00867                 {
    00868                         while (mp32even(size+1, udata))
    00869                         {
    00870                                 mp32divtwo(size+1, udata);
    00871
    00872                                 if (mp32odd(size+1, adata) || mp32odd(size+1, bdata))
    00873                                 {
    00874                                         (void) mp32addx(size+1, adata, xsize, xdata);
    00875                                         (void) mp32subx(size+1, bdata, size, b->modl);
    00876                                 }
    00877
    00878                                 mp32sdivtwo(size+1, adata);
    00879                                 mp32sdivtwo(size+1, bdata);
    00880                         }
    00881                         while (mp32even(size+1, vdata))
    00882                         {
    00883                                 mp32divtwo(size+1, vdata);
    00884
    00885                                 if (mp32odd(size+1, cdata) || mp32odd(size+1, ddata))
    00886                                 {
    00887                                         (void) mp32addx(size+1, cdata, xsize, xdata);
    00888                                         (void) mp32subx(size+1, ddata, size, b->modl);
    00889                                 }
    00890
    00891                                 mp32sdivtwo(size+1, cdata);
    00892                                 mp32sdivtwo(size+1, ddata);
    00893                         }
    00894                         if (mp32ge(size+1, udata, vdata))
    00895                         {
    00896                                 (void) mp32sub(size+1, udata, vdata);
    00897                                 (void) mp32sub(size+1, adata, cdata);
    00898                                 (void) mp32sub(size+1, bdata, ddata);
    00899                         }
    00900                         else
    00901                         {
    00902                                 (void) mp32sub(size+1, vdata, udata);
    00903                                 (void) mp32sub(size+1, cdata, adata);
    00904                                 (void) mp32sub(size+1, ddata, bdata);
    00905                         }
    00906
    00907                         if (mp32z(size+1, udata))
    00908                         {
    00909                                 if (mp32isone(size+1, vdata))
    00910                                 {
    00911                                         if (result)
    00912                                         {
    00913                                                 mp32setx(size, result, size+1, ddata);
    00914                                                 /*@-usedef@*/
    00915                                                 if (*ddata & 0x80000000)
    00916                                                 {
    00917                                                         /* keep adding the modulus until we get a carry */
    00918                                                         while (!mp32add(size, result, b->modl));
    00919                                                 }
    00920                                                 /*@=usedef@*/
    00921                                         }
    00922                                         return 1;
    00923                                 }
    00924                                 return 0;
    00925                         }
    00926                 }
    00927         }
    00928 }
    00929 #else
    00930
    00931 /*@unchecked@*/
    00932 static int _debug = 0;
    00933
    00934 #undef  FULL_BINARY_EXTENDED_GCD
    00935
    00939 /*@-boundsread@*/
    00940 int mp32binv_w(const mp32barrett* b, uint32 xsize, const uint32* xdata, uint32* result, uint32* wksp)
    00941 {
    00942         uint32  ysize = b->size+1;
    00943         int k;
    00944         uint32* u  = wksp;
    00945         uint32* v  =  u+ysize;
    00946         uint32* u1 =  v+ysize;
    00947         uint32* v1 = u1+ysize;
    00948         uint32* t1 = v1+ysize;
    00949         uint32* u3 = t1+ysize;
    00950         uint32* v3 = u3+ysize;
    00951         uint32* t3 = v3+ysize;
    00952
    00953 #ifdef  FULL_BINARY_EXTENDED_GCD
    00954         uint32* u2 = t3+ysize;
    00955         uint32* v2 = u2+ysize;
    00956         uint32* t2 = v2+ysize;
    00957 #endif
    00958
    00959         mp32setx(ysize, u, xsize, xdata);
    00960         mp32setx(ysize, v, b->size, b->modl);
    00961
    00962         /* Y1. Find power of 2. */
    00963         for (k = 0; mp32even(ysize, u) && mp32even(ysize, v); k++) {
    00964                 mp32divtwo(ysize, u);
    00965                 mp32divtwo(ysize, v);
    00966         }
    00967
    00968         /* Y2. Initialize. */
    00969         mp32setw(ysize, u1, 1);
    00970         mp32setx(ysize, v1, ysize, v);
    00971         mp32setx(ysize, u3, ysize, u);
    00972         mp32setx(ysize, v3, ysize, v);
    00973
    00974 #ifdef  FULL_BINARY_EXTENDED_GCD
    00975         mp32zero(ysize, u2);
    00976         mp32setw(ysize, v2, 1);
    00977         (void) mp32sub(ysize, v2, u);
    00978 #endif
    00979
    00980 if (_debug < 0) {
    00981 /*@-modfilesys@*/
    00982 fprintf(stderr, "       u: "), mp32println(stderr, ysize, u);
    00983 fprintf(stderr, "       v: "), mp32println(stderr, ysize, v);
    00984 fprintf(stderr, "      u1: "), mp32println(stderr, ysize, u1);
    00985 #ifdef  FULL_BINARY_EXTENDED_GCD
    00986 fprintf(stderr, "      u2: "), mp32println(stderr, ysize, u2);
    00987 #endif
    00988 fprintf(stderr, "      u3: "), mp32println(stderr, ysize, u3);
    00989 fprintf(stderr, "      v1: "), mp32println(stderr, ysize, v1);
    00990 #ifdef  FULL_BINARY_EXTENDED_GCD
    00991 fprintf(stderr, "      v2: "), mp32println(stderr, ysize, v2);
    00992 #endif
    00993 fprintf(stderr, "      v3: "), mp32println(stderr, ysize, v3);
    00994 /*@=modfilesys@*/
    00995 }
    00996
    00997         if (mp32odd(ysize, u)) {
    00998                 mp32zero(ysize, t1);
    00999 #ifdef  FULL_BINARY_EXTENDED_GCD
    01000                 mp32zero(ysize, t2);
    01001                 mp32subw(ysize, t2, 1);
    01002 #endif
    01003                 mp32zero(ysize, t3);
    01004                 (void) mp32sub(ysize, t3, v);
    01005                 goto Y4;
    01006         } else {
    01007                 mp32setw(ysize, t1, 1);
    01008 #ifdef  FULL_BINARY_EXTENDED_GCD
    01009                 mp32zero(ysize, t2);
    01010 #endif
    01011                 mp32setx(ysize, t3, ysize, u);
    01012         }
    01013
    01014         do {
    01015             do {
    01016 #ifdef  FULL_BINARY_EXTENDED_GCD
    01017                 if (mp32odd(ysize, t1) ||  mp32odd(ysize, t2)) {
    01018                         (void) mp32add(ysize, t1, v);
    01019                         (void) mp32sub(ysize, t2, u);
    01020                 }
    01021 #else
    01022                 /* XXX this assumes v is odd, true for DSA inversion. */
    01023                 if (mp32odd(ysize, t1))
    01024                         (void) mp32add(ysize, t1, v);
    01025 #endif
    01026
    01027                 mp32sdivtwo(ysize, t1);
    01028 #ifdef  FULL_BINARY_EXTENDED_GCD
    01029                 mp32sdivtwo(ysize, t2);
    01030 #endif
    01031                 mp32sdivtwo(ysize, t3);
    01032 Y4:
    01033 if (_debug < 0) {
    01034 /*@-modfilesys@*/
    01035 fprintf(stderr, "-->Y4 t3: "), mp32println(stderr, ysize, t3);
    01036 #ifdef  FULL_BINARY_EXTENDED_GCD
    01037 fprintf(stderr, "      t2: "), mp32println(stderr, ysize, t2);
    01038 #endif
    01039 fprintf(stderr, "      t1: "), mp32println(stderr, ysize, t1);
    01040 /*@=modfilesys@*/
    01041 }
    01042             } while (mp32even(ysize, t3));
    01043
    01044             /* Y5. Reset max(u3,v3). */
    01045             if (!(*t3 & 0x80000000)) {
    01046                 mp32setx(ysize, u1, ysize, t1);
    01047 #ifdef  FULL_BINARY_EXTENDED_GCD
    01048                 mp32setx(ysize, u2, ysize, t2);
    01049 #endif
    01050                 mp32setx(ysize, u3, ysize, t3);
    01051 if (_debug < 0) {
    01052 /*@-modfilesys@*/
    01053 fprintf(stderr, "-->Y5 u1: "), mp32println(stderr, ysize, u1);
    01054 #ifdef  FULL_BINARY_EXTENDED_GCD
    01055 fprintf(stderr, "      u2: "), mp32println(stderr, ysize, u2);
    01056 #endif
    01057 fprintf(stderr, "      u3: "), mp32println(stderr, ysize, u3);
    01058 /*@=modfilesys@*/
    01059 }
    01060             } else {
    01061                 mp32setx(ysize, v1, ysize, v);
    01062                 (void) mp32sub(ysize, v1, t1);
    01063 #ifdef  FULL_BINARY_EXTENDED_GCD
    01064                 mp32setx(ysize, v2, ysize, u);
    01065                 mp32neg(ysize, v2);
    01066                 (void) mp32sub(ysize, v2, t2);
    01067 #endif
    01068                 mp32zero(ysize, v3);
    01069                 (void) mp32sub(ysize, v3, t3);
    01070 if (_debug < 0) {
    01071 /*@-modfilesys@*/
    01072 fprintf(stderr, "-->Y5 v1: "), mp32println(stderr, ysize, v1);
    01073 #ifdef  FULL_BINARY_EXTENDED_GCD
    01074 fprintf(stderr, "      v2: "), mp32println(stderr, ysize, v2);
    01075 #endif
    01076 fprintf(stderr, "      v3: "), mp32println(stderr, ysize, v3);
    01077 /*@=modfilesys@*/
    01078 }
    01079             }
    01080
    01081             /* Y6. Subtract. */
    01082             mp32setx(ysize, t1, ysize, u1);
    01083             (void) mp32sub(ysize, t1, v1);
    01084 #ifdef  FULL_BINARY_EXTENDED_GCD
    01085             mp32setx(ysize, t2, ysize, u2);
    01086             (void) mp32sub(ysize, t2, v2);
    01087 #endif
    01088             mp32setx(ysize, t3, ysize, u3);
    01089             (void) mp32sub(ysize, t3, v3);
    01090
    01091             if (*t1 & 0x80000000) {
    01092                 (void) mp32add(ysize, t1, v);
    01093 #ifdef  FULL_BINARY_EXTENDED_GCD
    01094                 (void) mp32sub(ysize, t2, u);
    01095 #endif
    01096             }
    01097
    01098 if (_debug < 0) {
    01099 /*@-modfilesys@*/
    01100 fprintf(stderr, "-->Y6 t1: "), mp32println(stderr, ysize, t1);
    01101 #ifdef  FULL_BINARY_EXTENDED_GCD
    01102 fprintf(stderr, "      t2: "), mp32println(stderr, ysize, t2);
    01103 #endif
    01104 fprintf(stderr, "      t3: "), mp32println(stderr, ysize, t3);
    01105 /*@=modfilesys@*/
    01106 }
    01107
    01108         } while (mp32nz(ysize, t3));
    01109
    01110         if (!mp32isone(ysize, u3) || !mp32isone(ysize, v3))
    01111                 return 0;
    01112
    01113         if (result) {
    01114                 while (--k > 0)
    01115                         (void) mp32add(ysize, u1, u1);
    01116                 mp32setx(b->size, result, ysize, u1);
    01117         }
    01118
    01119 if (_debug) {
    01120 /*@-modfilesys@*/
    01121 if (result)
    01122 fprintf(stderr, "=== EXIT: "), mp32println(stderr, b->size, result);
    01123 fprintf(stderr, "      u1: "), mp32println(stderr, ysize, u1);
    01124 #ifdef  FULL_BINARY_EXTENDED_GCD
    01125 fprintf(stderr, "      u2: "), mp32println(stderr, ysize, u2);
    01126 #endif
    01127 fprintf(stderr, "      u3: "), mp32println(stderr, ysize, u3);
    01128 fprintf(stderr, "      v1: "), mp32println(stderr, ysize, v1);
    01129 #ifdef  FULL_BINARY_EXTENDED_GCD
    01130 fprintf(stderr, "      v2: "), mp32println(stderr, ysize, v2);
    01131 #endif
    01132 fprintf(stderr, "      v3: "), mp32println(stderr, ysize, v3);
    01133 fprintf(stderr, "      t1: "), mp32println(stderr, ysize, t1);
    01134 #ifdef  FULL_BINARY_EXTENDED_GCD
    01135 fprintf(stderr, "      t2: "), mp32println(stderr, ysize, t2);
    01136 #endif
    01137 fprintf(stderr, "      t3: "), mp32println(stderr, ysize, t3);
    01138 /*@=modfilesys@*/
    01139 }
    01140
    01141         return 1;
    01142 }
    01143 /*@=boundsread@*/
    01144
    01145 #endif
    01146
    01150 /*@-boundsread@*/
    01151 int mp32bpprime_w(const mp32barrett* b, randomGeneratorContext* rc, int t, uint32* wksp)
    01152 {
    01153         /*
    01154          * This test works for candidate probable primes >= 3, which are also not small primes.
    01155          *
    01156          * It assumes that b->modl contains the candidate prime
    01157          *
    01158          */
    01159
    01160         uint32 size = b->size;
    01161
    01162         /* first test if modl is odd */
    01163
    01164         if (mp32odd(b->size, b->modl))
    01165         {
    01166                 /*
    01167                  * Small prime factor test:
    01168                  *
    01169                  * Tables in mp32spprod contain multi-precision integers with products of small primes
    01170                  * If the greatest common divisor of this product and the candidate is not one, then
    01171                  * the candidate has small prime factors, or is a small prime. Neither is acceptable when
    01172                  * we are looking for large probable primes =)
    01173                  *
    01174                  */
    01175
    01176                 if (size > SMALL_PRIMES_PRODUCT_MAX)
    01177                 {
    01178                         /*@-globs@*/
    01179                         mp32setx(size, wksp+size, SMALL_PRIMES_PRODUCT_MAX, mp32spprod[SMALL_PRIMES_PRODUCT_MAX-1]);
    01180                         /*@=globs@*/
    01181                         /*@-compdef@*/ /* LCL: wksp+size */
    01182                         mp32gcd_w(size, b->modl, wksp+size, wksp, wksp+2*size);
    01183                         /*@=compdef@*/
    01184                 }
    01185                 else
    01186                 {
    01187                         /*@-globs@*/
    01188                         mp32gcd_w(size, b->modl, mp32spprod[size-1], wksp, wksp+2*size);
    01189                         /*@=globs@*/
    01190                 }
    01191
    01192                 if (mp32isone(size, wksp))
    01193                 {
    01194                         return mp32pmilrab_w(b, rc, t, wksp);
    01195                 }
    01196         }
    01197
    01198         return 0;
    01199 }
    01200 /*@=boundsread@*/
    01201
    01202 void mp32bnrnd(const mp32barrett* b, randomGeneratorContext* rc, mp32number* result)
    01203 {
    01204         register uint32  size = b->size;
    01205         register uint32* temp = (uint32*) malloc(size * sizeof(uint32));
    01206
    01207         mp32nfree(result);
    01208         mp32nsize(result, size);
    01209         /*@-nullpass@*/         /* temp may be NULL */
    01210         /*@-usedef@*/           /* result->data unallocated? */
    01211         mp32brnd_w(b, rc, result->data, temp);
    01212         /*@=usedef@*/
    01213
    01214         free(temp);
    01215         /*@=nullpass@*/
    01216 }
    01217
    01218 void mp32bnmulmod(const mp32barrett* b, const mp32number* x, const mp32number* y, mp32number* result)
    01219 {
    01220         register uint32  size = b->size;
    01221         register uint32* temp = (uint32*) malloc((4*size+2) * sizeof(uint32));
    01222
    01223         /* xsize and ysize must be <= b->size */
    01224         register uint32  fill = 2*size-x->size-y->size;
    01225         /*@-nullptrarith@*/     /* temp may be NULL */
    01226         register uint32* opnd = temp+size*2+2;
    01227         /*@=nullptrarith@*/
    01228
    01229         mp32nfree(result);
    01230         mp32nsize(result, size);
    01231
    01232         if (fill)
    01233                 mp32zero(fill, opnd);
    01234
    01235         mp32mul(opnd+fill, x->size, x->data, y->size, y->data);
    01236         /*@-nullpass@*/         /* temp may be NULL */
    01237         /*@-usedef -compdef @*/ /* result->data unallocated? */
    01238         mp32bmod_w(b, opnd, result->data, temp);
    01239         /*@=usedef =compdef @*/
    01240
    01241         free(temp);
    01242         /*@=nullpass@*/
    01243 }
    01244
    01245 void mp32bnsqrmod(const mp32barrett* b, const mp32number* x, mp32number* result)
    01246 {
    01247         register uint32  size = b->size;
    01248         register uint32* temp = (uint32*) malloc(size * sizeof(uint32));
    01249
    01250         /* xsize must be <= b->size */
    01251         register uint32  fill = 2*(size-x->size);
    01252         /*@-nullptrarith@*/     /* temp may be NULL */
    01253         register uint32* opnd = temp + size*2+2;
    01254         /*@=nullptrarith@*/
    01255
    01256         mp32nfree(result);
    01257         mp32nsize(result, size);
    01258
    01259         if (fill)
    01260                 mp32zero(fill, opnd);
    01261
    01262         mp32sqr(opnd+fill, x->size, x->data);
    01263         /*@-nullpass@*/         /* temp may be NULL */
    01264         /*@-usedef -compdef @*/ /* result->data unallocated? */
    01265         mp32bmod_w(b, opnd, result->data, temp);
    01266         /*@=usedef =compdef @*/
    01267
    01268         free(temp);
    01269         /*@=nullpass@*/
    01270 }
    01271
    01272 void mp32bnpowmod(const mp32barrett* b, const mp32number* x, const mp32number* pow, mp32number* y)
    01273 {
    01274         register uint32  size = b->size;
    01275         register uint32* temp = (uint32*) malloc((4*size+2) * sizeof(uint32));
    01276
    01277         mp32nfree(y);
    01278         mp32nsize(y, size);
    01279
    01280         /*@-nullpass@*/         /* temp may be NULL */
    01281         mp32bpowmod_w(b, x->size, x->data, pow->size, pow->data, y->data, temp);
    01282
    01283         free(temp);
    01284         /*@=nullpass@*/
    01285 }
    01286
    01287 void mp32bnpowmodsld(const mp32barrett* b, const uint32* slide, const mp32number* pow, mp32number* y)
    01288 {
    01289         register uint32  size = b->size;
    01290         register uint32* temp = (uint32*) malloc((4*size+2) * sizeof(uint32));
    01291
    01292         mp32nfree(y);
    01293         mp32nsize(y, size);
    01294
    01295         /*@-nullpass@*/         /* temp may be NULL */
    01296         /*@-internalglobs -mods@*/ /* noisy */
    01297         mp32bpowmodsld_w(b, slide, pow->size, pow->data, y->data, temp);
    01298         /*@=internalglobs =mods@*/
    01299
    01300         free(temp);
    01301         /*@=nullpass@*/
    01302 }
    01303 /*@=sizeoftype =type@*/