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@*/

Generated on Wed Sep 4 12:49:48 2002 for rpm by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002