mpi.c

00001 /**
00002  * @file mpi.c
00003  * @brief MPI (Multiple Precision Integer Arithmetic)
00004  *
00005  * @section License
00006  *
00007  * Copyright (C) 2010-2017 Oryx Embedded SARL. All rights reserved.
00008  *
00009  * This file is part of CycloneCrypto Open.
00010  *
00011  * This program is free software; you can redistribute it and/or
00012  * modify it under the terms of the GNU General Public License
00013  * as published by the Free Software Foundation; either version 2
00014  * of the License, or (at your option) any later version.
00015  *
00016  * This program is distributed in the hope that it will be useful,
00017  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00018  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00019  * GNU General Public License for more details.
00020  *
00021  * You should have received a copy of the GNU General Public License
00022  * along with this program; if not, write to the Free Software Foundation,
00023  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
00024  *
00025  * @author Oryx Embedded SARL (www.oryx-embedded.com)
00026  * @version 1.7.6
00027  **/
00028 
00029 //Switch to the appropriate trace level
00030 #define TRACE_LEVEL CRYPTO_TRACE_LEVEL
00031 
00032 //Dependencies
00033 #include <stdlib.h>
00034 #include <string.h>
00035 #include "crypto.h"
00036 #include "mpi.h"
00037 #include "debug.h"
00038 
00039 //Check crypto library configuration
00040 #if (MPI_SUPPORT == ENABLED)
00041 
00042 
00043 /**
00044  * @brief Initialize a multiple precision integer
00045  * @param[in,out] r Pointer to the multiple precision integer to be initialized
00046  **/
00047 
00048 void mpiInit(Mpi *r)
00049 {
00050    //Initialize structure
00051    r->sign = 1;
00052    r->size = 0;
00053    r->data = NULL;
00054 }
00055 
00056 
00057 /**
00058  * @brief Release a multiple precision integer
00059  * @param[in,out] r Pointer to the multiple precision integer to be freed
00060  **/
00061 
00062 void mpiFree(Mpi *r)
00063 {
00064    //Any memory previously allocated?
00065    if(r->data != NULL)
00066    {
00067       //Erase contents before releasing memory
00068       memset(r->data, 0, r->size * MPI_INT_SIZE);
00069       cryptoFreeMem(r->data);
00070    }
00071 
00072    //Set size to zero
00073    r->size = 0;
00074    r->data = NULL;
00075 }
00076 
00077 
00078 /**
00079  * @brief Adjust the size of multiple precision integer
00080  * @param[in,out] r A multiple precision integer whose size is to be increased
00081  * @param[in] size Desired size in words
00082  * @return Error code
00083  **/
00084 
00085 error_t mpiGrow(Mpi *r, uint_t size)
00086 {
00087    uint_t *data;
00088 
00089    //Ensure the parameter is valid
00090    size = MAX(size, 1);
00091 
00092    //Check the current size
00093    if(r->size >= size)
00094       return NO_ERROR;
00095 
00096    //Allocate a memory buffer
00097    data = cryptoAllocMem(size * MPI_INT_SIZE);
00098    //Failed to allocate memory?
00099    if(data == NULL)
00100       return ERROR_OUT_OF_MEMORY;
00101 
00102    //Clear buffer contents
00103    memset(data, 0, size * MPI_INT_SIZE);
00104 
00105    //Any data to copy?
00106    if(r->size > 0)
00107    {
00108       //Copy original data
00109       memcpy(data, r->data, r->size * MPI_INT_SIZE);
00110       //Free previously allocated memory
00111       cryptoFreeMem(r->data);
00112    }
00113 
00114    //Update the size of the multiple precision integer
00115    r->size = size;
00116    r->data = data;
00117 
00118    //Successful operation
00119    return NO_ERROR;
00120 }
00121 
00122 
00123 /**
00124  * @brief Get the actual length in words
00125  * @param[in] a Pointer to a multiple precision integer
00126  * @return The actual length in words
00127  **/
00128 
00129 uint_t mpiGetLength(const Mpi *a)
00130 {
00131    int_t i;
00132 
00133    //Check whether the specified multiple precision integer is empty
00134    if(a->size == 0)
00135       return 0;
00136 
00137    //Start from the most significant word
00138    for(i = a->size - 1; i >= 0; i--)
00139    {
00140       //Loop as long as the current word is zero
00141       if(a->data[i] != 0)
00142          break;
00143    }
00144 
00145    //Return the actual length
00146    return i + 1;
00147 }
00148 
00149 
00150 /**
00151  * @brief Get the actual length in bytes
00152  * @param[in] a Pointer to a multiple precision integer
00153  * @return The actual byte count
00154  **/
00155 
00156 uint_t mpiGetByteLength(const Mpi *a)
00157 {
00158    uint_t n;
00159    uint32_t m;
00160 
00161    //Check whether the specified multiple precision integer is empty
00162    if(a->size == 0)
00163       return 0;
00164 
00165    //Start from the most significant word
00166    for(n = a->size - 1; n > 0; n--)
00167    {
00168       //Loop as long as the current word is zero
00169       if(a->data[n] != 0)
00170          break;
00171    }
00172 
00173    //Get the current word
00174    m = a->data[n];
00175    //Convert the length to a byte count
00176    n *= MPI_INT_SIZE;
00177 
00178    //Adjust the byte count
00179    for(; m != 0; m >>= 8) n++;
00180 
00181    //Return the actual length in bytes
00182    return n;
00183 }
00184 
00185 
00186 /**
00187  * @brief Get the actual length in bits
00188  * @param[in] a Pointer to a multiple precision integer
00189  * @return The actual bit count
00190  **/
00191 
00192 uint_t mpiGetBitLength(const Mpi *a)
00193 {
00194    uint_t n;
00195    uint32_t m;
00196 
00197    //Check whether the specified multiple precision integer is empty
00198    if(a->size == 0)
00199       return 0;
00200 
00201    //Start from the most significant word
00202    for(n = a->size - 1; n > 0; n--)
00203    {
00204       //Loop as long as the current word is zero
00205       if(a->data[n] != 0)
00206          break;
00207    }
00208 
00209    //Get the current word
00210    m = a->data[n];
00211    //Convert the length to a bit count
00212    n *= MPI_INT_SIZE * 8;
00213 
00214    //Adjust the bit count
00215    for(; m != 0; m >>= 1) n++;
00216 
00217    //Return the actual length in bits
00218    return n;
00219 }
00220 
00221 
00222 /**
00223  * @brief Set the bit value at the specified index
00224  * @param[in] r Pointer to a multiple precision integer
00225  * @param[in] index Position of the bit to be written
00226  * @param[in] value Bit value
00227  * @return Error code
00228  **/
00229 
00230 error_t mpiSetBitValue(Mpi *r, uint_t index, uint_t value)
00231 {
00232    error_t error;
00233    uint_t n1;
00234    uint_t n2;
00235 
00236    //Retrieve the position of the bit to be written
00237    n1 = index / (MPI_INT_SIZE * 8);
00238    n2 = index % (MPI_INT_SIZE * 8);
00239 
00240    //Ajust the size of the multiple precision integer if necessary
00241    error = mpiGrow(r, (index + (MPI_INT_SIZE * 8) - 1) / (MPI_INT_SIZE * 8));
00242    //Failed to adjust the size?
00243    if(error)
00244       return error;
00245 
00246    //Set bit value
00247    if(value)
00248       r->data[n1] |= (1 << n2);
00249    else
00250       r->data[n1] &= ~(1 << n2);
00251 
00252    //No error to report
00253    return NO_ERROR;
00254 }
00255 
00256 
00257 /**
00258  * @brief Get the bit value at the specified index
00259  * @param[in] a Pointer to a multiple precision integer
00260  * @param[in] index Position where to read the bit
00261  * @return The actual bit value
00262  **/
00263 
00264 uint_t mpiGetBitValue(const Mpi *a, uint_t index)
00265 {
00266    uint_t n1;
00267    uint_t n2;
00268 
00269    //Retrieve the position of the bit to be read
00270    n1 = index / (MPI_INT_SIZE * 8);
00271    n2 = index % (MPI_INT_SIZE * 8);
00272 
00273    //Index out of range?
00274    if(n1 >= a->size)
00275       return 0;
00276 
00277    //Return the actual bit value
00278    return (a->data[n1] >> n2) & 0x01;
00279 }
00280 
00281 
00282 /**
00283  * @brief Compare two multiple precision integers
00284  * @param[in] a The first multiple precision integer to be compared
00285  * @param[in] b The second multiple precision integer to be compared
00286  * @return Comparison result
00287  **/
00288 
00289 int_t mpiComp(const Mpi *a, const Mpi *b)
00290 {
00291    uint_t m;
00292    uint_t n;
00293 
00294    //Determine the actual length of A and B
00295    m = mpiGetLength(a);
00296    n = mpiGetLength(b);
00297 
00298    //Compare lengths
00299    if(!m && !n)
00300       return 0;
00301    else if(m > n)
00302       return a->sign;
00303    else if(m < n)
00304       return -b->sign;
00305 
00306    //Compare signs
00307    if(a->sign > 0 && b->sign < 0)
00308       return 1;
00309    else if(a->sign < 0 && b->sign > 0)
00310       return -1;
00311 
00312    //Then compare values
00313    while(n--)
00314    {
00315       if(a->data[n] > b->data[n])
00316          return a->sign;
00317       else if(a->data[n] < b->data[n])
00318          return -a->sign;
00319    }
00320 
00321    //Multiple precision integers are equals
00322    return 0;
00323 }
00324 
00325 
00326 /**
00327  * @brief Compare a multiple precision integer with an integer
00328  * @param[in] a Multiple precision integer to be compared
00329  * @param[in] b Integer to be compared
00330  * @return Comparison result
00331  **/
00332 
00333 int_t mpiCompInt(const Mpi *a, int_t b)
00334 {
00335    uint_t value;
00336    Mpi t;
00337 
00338    //Initialize a temporary multiple precision integer
00339    value = (b >= 0) ? b : -b;
00340    t.sign = (b >= 0) ? 1 : -1;
00341    t.size = 1;
00342    t.data = &value;
00343 
00344    //Return comparison result
00345    return mpiComp(a, &t);
00346 }
00347 
00348 
00349 /**
00350  * @brief Compare the absolute value of two multiple precision integers
00351  * @param[in] a The first multiple precision integer to be compared
00352  * @param[in] b The second multiple precision integer to be compared
00353  * @return Comparison result
00354  **/
00355 
00356 int_t mpiCompAbs(const Mpi *a, const Mpi *b)
00357 {
00358    uint_t m;
00359    uint_t n;
00360 
00361    //Determine the actual length of A and B
00362    m = mpiGetLength(a);
00363    n = mpiGetLength(b);
00364 
00365    //Compare lengths
00366    if(!m && !n)
00367       return 0;
00368    else if(m > n)
00369       return 1;
00370    else if(m < n)
00371       return -1;
00372 
00373    //Then compare values
00374    while(n--)
00375    {
00376       if(a->data[n] > b->data[n])
00377          return 1;
00378       else if(a->data[n] < b->data[n])
00379          return -1;
00380    }
00381 
00382    //Operands are equals
00383    return 0;
00384 }
00385 
00386 
00387 /**
00388  * @brief Copy a multiple precision integer
00389  * @param[out] r Pointer to a multiple precision integer (destination)
00390  * @param[in] a Pointer to a multiple precision integer (source)
00391  * @return Error code
00392  **/
00393 
00394 error_t mpiCopy(Mpi *r, const Mpi *a)
00395 {
00396    error_t error;
00397    uint_t n;
00398 
00399    //R and A are the same instance?
00400    if(r == a)
00401       return NO_ERROR;
00402 
00403    //Determine the actual length of A
00404    n = mpiGetLength(a);
00405 
00406    //Ajust the size of the destination operand
00407    error = mpiGrow(r, n);
00408    //Any error to report?
00409    if(error)
00410       return error;
00411 
00412    //Clear the contents of the multiple precision integer
00413    memset(r->data, 0, r->size * MPI_INT_SIZE);
00414    //Let R = A
00415    memcpy(r->data, a->data, n * MPI_INT_SIZE);
00416    //Set the sign of R
00417    r->sign = a->sign;
00418 
00419    //Successful operation
00420    return NO_ERROR;
00421 }
00422 
00423 
00424 /**
00425  * @brief Set the value of a multiple precision integer
00426  * @param[out] r Pointer to a multiple precision integer
00427  * @param[in] a Value to be assigned to the multiple precision integer
00428  * @return Error code
00429  **/
00430 
00431 error_t mpiSetValue(Mpi *r, int_t a)
00432 {
00433    error_t error;
00434 
00435    //Ajust the size of the destination operand
00436    error = mpiGrow(r, 1);
00437    //Failed to adjust the size?
00438    if(error)
00439       return error;
00440 
00441    //Clear the contents of the multiple precision integer
00442    memset(r->data, 0, r->size * MPI_INT_SIZE);
00443    //Set the value or R
00444    r->data[0] = (a >= 0) ? a : -a;
00445    //Set the sign of R
00446    r->sign = (a >= 0) ? 1 : -1;
00447 
00448    //Successful operation
00449    return NO_ERROR;
00450 }
00451 
00452 
00453 /**
00454  * @brief Generate a random value
00455  * @param[out] r Pointer to a multiple precision integer
00456  * @param[in] length Desired length in bits
00457  * @param[in] prngAlgo PRNG algorithm
00458  * @param[in] prngContext Pointer to the PRNG context
00459  * @return Error code
00460  **/
00461 
00462 error_t mpiRand(Mpi *r, uint_t length, const PrngAlgo *prngAlgo, void *prngContext)
00463 {
00464    error_t error;
00465    uint_t m;
00466    uint_t n;
00467 
00468    //Compute the required length, in words
00469    n = (length + (MPI_INT_SIZE * 8) - 1) / (MPI_INT_SIZE * 8);
00470    //Number of bits in the most significant word
00471    m = length % (MPI_INT_SIZE * 8);
00472 
00473    //Ajust the size of the multiple precision integer if necessary
00474    error = mpiGrow(r, n);
00475    //Failed to adjust the size?
00476    if(error)
00477       return error;
00478 
00479    //Clear the contents of the multiple precision integer
00480    memset(r->data, 0, r->size * MPI_INT_SIZE);
00481    //Set the sign of R
00482    r->sign = 1;
00483 
00484    //Generate a random pattern
00485    error = prngAlgo->read(prngContext, (uint8_t *) r->data, n * MPI_INT_SIZE);
00486    //Any error to report?
00487    if(error)
00488       return error;
00489 
00490    //Remove the meaningless bits in the most significant word
00491    if(n > 0 && m > 0)
00492       r->data[n - 1] &= (1 << m) - 1;
00493 
00494    //Successful operation
00495    return NO_ERROR;
00496 }
00497 
00498 
00499 /**
00500  * @brief Octet string to integer conversion
00501  *
00502  * Converts an octet string to a non-negative integer
00503  *
00504  * @param[out] r Non-negative integer resulting from the conversion
00505  * @param[in] data Octet string to be converted
00506  * @param[in] length Length of the octet string
00507  * @return Error code
00508  **/
00509 
00510 error_t mpiReadRaw(Mpi *r, const uint8_t *data, uint_t length)
00511 {
00512    error_t error;
00513    uint_t i;
00514 
00515    //Skip leading zeroes
00516    while(length > 1 && *data == 0)
00517    {
00518       //Advance read pointer
00519       data++;
00520       length--;
00521    }
00522 
00523    //Ajust the size of the multiple precision integer
00524    error = mpiGrow(r, (length + MPI_INT_SIZE - 1) / MPI_INT_SIZE);
00525    //Failed to adjust the size?
00526    if(error)
00527       return error;
00528 
00529    //Clear the contents of the multiple precision integer
00530    memset(r->data, 0, r->size * MPI_INT_SIZE);
00531    //Set sign
00532    r->sign = 1;
00533 
00534    //Start from the least significant byte
00535    data += length - 1;
00536 
00537    //Copy data
00538    for(i = 0; i < length; i++, data--)
00539       r->data[i / MPI_INT_SIZE] |= *data << ((i % MPI_INT_SIZE) * 8);
00540 
00541    //The conversion succeeded
00542    return NO_ERROR;
00543 }
00544 
00545 
00546 /**
00547  * @brief Integer to octet string conversion
00548  *
00549  * Converts an integer to an octet string of a specified length
00550  *
00551  * @param[in] a Non-negative integer to be converted
00552  * @param[out] data Octet string resulting from the conversion
00553  * @param[in] length Intended length of the resulting octet string
00554  * @return Error code
00555  **/
00556 
00557 error_t mpiWriteRaw(const Mpi *a, uint8_t *data, uint_t length)
00558 {
00559    uint_t i;
00560 
00561    //Get the actual length in bytes
00562    uint_t n = mpiGetByteLength(a);
00563 
00564    //Make sure the output buffer is large enough
00565    if(n > length)
00566       return ERROR_INVALID_LENGTH;
00567 
00568    //Clear output buffer
00569    memset(data, 0, length);
00570 
00571    //Start from the least significant word
00572    data += length - 1;
00573 
00574    //Copy data
00575    for(i = 0; i < n; i++, data--)
00576       *data = a->data[i / MPI_INT_SIZE] >> ((i % MPI_INT_SIZE) * 8);
00577 
00578    //The conversion succeeded
00579    return NO_ERROR;
00580 }
00581 
00582 
00583 /**
00584  * @brief Multiple precision addition
00585  * @param[out] r Resulting integer R = A + B
00586  * @param[in] a First operand A
00587  * @param[in] b Second operand B
00588  * @return Error code
00589  **/
00590 
00591 error_t mpiAdd(Mpi *r, const Mpi *a, const Mpi *b)
00592 {
00593    error_t error;
00594    int_t sign;
00595 
00596    //Retrieve the sign of A
00597    sign = a->sign;
00598 
00599    //Both operands have the same sign?
00600    if(a->sign == b->sign)
00601    {
00602       //Perform addition
00603       error = mpiAddAbs(r, a, b);
00604       //Set the sign of the resulting number
00605       r->sign = sign;
00606    }
00607    //Operands have different signs?
00608    else
00609    {
00610       //Compare the absolute value of A and B
00611       if(mpiCompAbs(a, b) >= 0)
00612       {
00613          //Perform subtraction
00614          error = mpiSubAbs(r, a, b);
00615          //Set the sign of the resulting number
00616          r->sign = sign;
00617       }
00618       else
00619       {
00620          //Perform subtraction
00621          error = mpiSubAbs(r, b, a);
00622          //Set the sign of the resulting number
00623          r->sign = -sign;
00624       }
00625    }
00626 
00627    //Return status code
00628    return error;
00629 }
00630 
00631 
00632 /**
00633  * @brief Add an integer to a multiple precision integer
00634  * @param[out] r Resulting integer R = A + B
00635  * @param[in] a First operand A
00636  * @param[in] b Second operand B
00637  * @return Error code
00638  **/
00639 
00640 error_t mpiAddInt(Mpi *r, const Mpi *a, int_t b)
00641 {
00642    uint_t value;
00643    Mpi t;
00644 
00645    //Convert the second operand to a multiple precision integer
00646    value = (b >= 0) ? b : -b;
00647    t.sign = (b >= 0) ? 1 : -1;
00648    t.size = 1;
00649    t.data = &value;
00650 
00651    //Perform addition
00652    return mpiAdd(r, a ,&t);
00653 }
00654 
00655 
00656 /**
00657  * @brief Multiple precision subtraction
00658  * @param[out] r Resulting integer R = A - B
00659  * @param[in] a First operand A
00660  * @param[in] b Second operand B
00661  * @return Error code
00662  **/
00663 
00664 error_t mpiSub(Mpi *r, const Mpi *a, const Mpi *b)
00665 {
00666    error_t error;
00667    int_t sign;
00668 
00669    //Retrieve the sign of A
00670    sign = a->sign;
00671 
00672    //Both operands have the same sign?
00673    if(a->sign == b->sign)
00674    {
00675       //Compare the absolute value of A and B
00676       if(mpiCompAbs(a, b) >= 0)
00677       {
00678          //Perform subtraction
00679          error = mpiSubAbs(r, a, b);
00680          //Set the sign of the resulting number
00681          r->sign = sign;
00682       }
00683       else
00684       {
00685          //Perform subtraction
00686          error = mpiSubAbs(r, b, a);
00687          //Set the sign of the resulting number
00688          r->sign = -sign;
00689       }
00690    }
00691    //Operands have different signs?
00692    else
00693    {
00694       //Perform addition
00695       error = mpiAddAbs(r, a, b);
00696       //Set the sign of the resulting number
00697       r->sign = sign;
00698    }
00699 
00700    //Return status code
00701    return error;
00702 }
00703 
00704 
00705 /**
00706  * @brief Subtract an integer from a multiple precision integer
00707  * @param[out] r Resulting integer R = A - B
00708  * @param[in] a First operand A
00709  * @param[in] b Second operand B
00710  * @return Error code
00711  **/
00712 
00713 error_t mpiSubInt(Mpi *r, const Mpi *a, int_t b)
00714 {
00715    uint_t value;
00716    Mpi t;
00717 
00718    //Convert the second operand to a multiple precision integer
00719    value = (b >= 0) ? b : -b;
00720    t.sign = (b >= 0) ? 1 : -1;
00721    t.size = 1;
00722    t.data = &value;
00723 
00724    //Perform subtraction
00725    return mpiSub(r, a ,&t);
00726 }
00727 
00728 
00729 /**
00730  * @brief Helper routine for multiple precision addition
00731  * @param[out] r Resulting integer R = |A + B|
00732  * @param[in] a First operand A
00733  * @param[in] b Second operand B
00734  * @return Error code
00735  **/
00736 
00737 error_t mpiAddAbs(Mpi *r, const Mpi *a, const Mpi *b)
00738 {
00739    error_t error;
00740    uint_t i;
00741    uint_t n;
00742    uint_t c;
00743    uint_t d;
00744 
00745    //R and B are the same instance?
00746    if(r == b)
00747    {
00748       //Swap A and B
00749       const Mpi *t = a;
00750       a = b;
00751       b = t;
00752    }
00753    //R is neither A nor B?
00754    else if(r != a)
00755    {
00756       //Copy the first operand to R
00757       MPI_CHECK(mpiCopy(r, a));
00758    }
00759 
00760    //Determine the actual length of B
00761    n = mpiGetLength(b);
00762    //Extend the size of the destination register as needed
00763    MPI_CHECK(mpiGrow(r, n));
00764 
00765    //The result is always positive
00766    r->sign = 1;
00767    //Clear carry bit
00768    c = 0;
00769 
00770    //Add operands
00771    for(i = 0; i < n; i++)
00772    {
00773       //Add carry bit
00774       d = r->data[i] + c;
00775       //Update carry bit
00776       if(d != 0) c = 0;
00777       //Perform addition
00778       d += b->data[i];
00779       //Update carry bit
00780       if(d < b->data[i]) c = 1;
00781       //Save result
00782       r->data[i] = d;
00783    }
00784 
00785    //Loop as long as the carry bit is set
00786    for(i = n; c && i < r->size; i++)
00787    {
00788       //Add carry bit
00789       r->data[i] += c;
00790       //Update carry bit
00791       if(r->data[i] != 0) c = 0;
00792    }
00793 
00794    //Check the final carry bit
00795    if(c && n >= r->size)
00796    {
00797       //Extend the size of the destination register
00798       MPI_CHECK(mpiGrow(r, n + 1));
00799       //Add carry bit
00800       r->data[n] = 1;
00801    }
00802 
00803 end:
00804    //Return status code
00805    return error;
00806 }
00807 
00808 
00809 /**
00810  * @brief Helper routine for multiple precision subtraction
00811  * @param[out] r Resulting integer R = |A - B|
00812  * @param[in] a First operand A
00813  * @param[in] b Second operand B
00814  * @return Error code
00815  **/
00816 
00817 error_t mpiSubAbs(Mpi *r, const Mpi *a, const Mpi *b)
00818 {
00819    error_t error;
00820    uint_t c;
00821    uint_t d;
00822    uint_t i;
00823    uint_t m;
00824    uint_t n;
00825 
00826    //Check input parameters
00827    if(mpiCompAbs(a, b) < 0)
00828    {
00829       //Swap A and B if necessary
00830       const Mpi *t = b;
00831       a = b;
00832       b = t;
00833    }
00834 
00835    //Determine the actual length of A
00836    m = mpiGetLength(a);
00837    //Determine the actual length of B
00838    n = mpiGetLength(b);
00839 
00840    //Extend the size of the destination register as needed
00841    MPI_CHECK(mpiGrow(r, m));
00842 
00843    //The result is always positive
00844    r->sign = 1;
00845    //Clear carry bit
00846    c = 0;
00847 
00848    //Subtract operands
00849    for(i = 0; i < n; i++)
00850    {
00851       //Read first operand
00852       d = a->data[i];
00853 
00854       //Check the carry bit
00855       if(c)
00856       {
00857          //Update carry bit
00858          if(d != 0) c = 0;
00859          //Propagate carry bit
00860          d -= 1;
00861       }
00862 
00863       //Update carry bit
00864       if(d < b->data[i]) c = 1;
00865       //Perform subtraction
00866       r->data[i] = d - b->data[i];
00867    }
00868 
00869    //Loop as long as the carry bit is set
00870    for(i = n; c && i < m; i++)
00871    {
00872       //Update carry bit
00873       if(a->data[i] != 0) c = 0;
00874       //Propagate carry bit
00875       r->data[i] = a->data[i] - 1;
00876    }
00877 
00878    //R and A are not the same instance?
00879    if(r != a)
00880    {
00881       //Copy the remaining words
00882       for(; i < m; i++)
00883          r->data[i] = a->data[i];
00884 
00885       //Zero the upper part of R
00886       for(; i < r->size; i++)
00887          r->data[i] = 0;
00888    }
00889 
00890 end:
00891    //Return status code
00892    return error;
00893 }
00894 
00895 
00896 /**
00897  * @brief Left shift operation
00898  * @param[in,out] r The multiple precision integer to be shifted to the left
00899  * @param[in] n The number of bits to shift
00900  * @return Error code
00901  **/
00902 
00903 error_t mpiShiftLeft(Mpi *r, uint_t n)
00904 {
00905    error_t error;
00906    uint_t i;
00907 
00908    //Number of 32-bit words to shift
00909    uint_t n1 = n / (MPI_INT_SIZE * 8);
00910    //Number of bits to shift
00911    uint_t n2 = n % (MPI_INT_SIZE * 8);
00912 
00913    //Check parameters
00914    if(!r->size || !n)
00915       return NO_ERROR;
00916 
00917    //Increase the size of the multiple-precision number
00918    error = mpiGrow(r, r->size + (n + 31) / 32);
00919    //Check return code
00920    if(error)
00921       return error;
00922 
00923    //First, shift words
00924    if(n1 > 0)
00925    {
00926       //Process the most significant words
00927       for(i = r->size - 1; i >= n1; i--)
00928          r->data[i] = r->data[i - n1];
00929       //Fill the rest with zeroes
00930       for(i = 0; i < n1; i++)
00931          r->data[i] = 0;
00932    }
00933    //Then shift bits
00934    if(n2 > 0)
00935    {
00936       //Process the most significant words
00937       for(i = r->size - 1; i >= 1; i--)
00938          r->data[i] = (r->data[i] << n2) | (r->data[i - 1] >> (32 - n2));
00939       //The least significant word requires a special handling
00940       r->data[0] <<= n2;
00941    }
00942 
00943    //Shift operation is complete
00944    return NO_ERROR;
00945 }
00946 
00947 
00948 /**
00949  * @brief Right shift operation
00950  * @param[in,out] r The multiple precision integer to be shifted to the right
00951  * @param[in] n The number of bits to shift
00952  * @return Error code
00953  **/
00954 
00955 error_t mpiShiftRight(Mpi *r, uint_t n)
00956 {
00957    uint_t i;
00958    uint_t m;
00959 
00960    //Number of 32-bit words to shift
00961    uint_t n1 = n / (MPI_INT_SIZE * 8);
00962    //Number of bits to shift
00963    uint_t n2 = n % (MPI_INT_SIZE * 8);
00964 
00965    //Check parameters
00966    if(n1 >= r->size)
00967    {
00968       memset(r->data, 0, r->size * MPI_INT_SIZE);
00969       return NO_ERROR;
00970    }
00971 
00972    //First, shift words
00973    if(n1 > 0)
00974    {
00975       //Process the least significant words
00976       for(m = r->size - n1, i = 0; i < m; i++)
00977          r->data[i] = r->data[i + n1];
00978       //Fill the rest with zeroes
00979       for(i = m; i < r->size; i++)
00980          r->data[i] = 0;
00981    }
00982    //Then shift bits
00983    if(n2 > 0)
00984    {
00985       //Process the least significant words
00986       for(m = r->size - n1 - 1, i = 0; i < m; i++)
00987          r->data[i] = (r->data[i] >> n2) | (r->data[i + 1] << (32 - n2));
00988       //The most significant word requires a special handling
00989       r->data[m] >>= n2;
00990    }
00991 
00992    //Shift operation is complete
00993    return NO_ERROR;
00994 }
00995 
00996 
00997 /**
00998  * @brief Multiple precision multiplication
00999  * @param[out] r Resulting integer R = A * B
01000  * @param[in] a First operand A
01001  * @param[in] b Second operand B
01002  * @return Error code
01003  **/
01004 
01005 error_t mpiMul(Mpi *r, const Mpi *a, const Mpi *b)
01006 {
01007    error_t error;
01008    int_t i;
01009    int_t m;
01010    int_t n;
01011    Mpi ta;
01012    Mpi tb;
01013 
01014    //Initialize multiple precision integers
01015    mpiInit(&ta);
01016    mpiInit(&tb);
01017 
01018    //R and A are the same instance?
01019    if(r == a)
01020    {
01021       //Copy A to TA
01022       MPI_CHECK(mpiCopy(&ta, a));
01023       //Use TA instead of A
01024       a = &ta;
01025    }
01026 
01027    //R and B are the same instance?
01028    if(r == b)
01029    {
01030       //Copy B to TB
01031       MPI_CHECK(mpiCopy(&tb, b));
01032       //Use TB instead of B
01033       b = &tb;
01034    }
01035 
01036    //Determine the actual length of A and B
01037    m = mpiGetLength(a);
01038    n = mpiGetLength(b);
01039 
01040    //Adjust the size of R
01041    MPI_CHECK(mpiGrow(r, m + n));
01042    //Set the sign of R
01043    r->sign = (a->sign == b->sign) ? 1 : -1;
01044 
01045    //Clear the contents of the destination integer
01046    memset(r->data, 0, r->size * MPI_INT_SIZE);
01047 
01048    //Perform multiplication
01049    if(m < n)
01050    {
01051       for(i = 0; i < m; i++)
01052          mpiMulAccCore(&r->data[i], b->data, n, a->data[i]);
01053    }
01054    else
01055    {
01056       for(i = 0; i < n; i++)
01057          mpiMulAccCore(&r->data[i], a->data, m, b->data[i]);
01058    }
01059 
01060 end:
01061    //Release multiple precision integers
01062    mpiFree(&ta);
01063    mpiFree(&tb);
01064 
01065    //Return status code
01066    return error;
01067 }
01068 
01069 
01070 /**
01071  * @brief Multiply a multiple precision integer by an integer
01072  * @param[out] r Resulting integer R = A * B
01073  * @param[in] a First operand A
01074  * @param[in] b Second operand B
01075  * @return Error code
01076  **/
01077 
01078 error_t mpiMulInt(Mpi *r, const Mpi *a, int_t b)
01079 {
01080    uint_t value;
01081    Mpi t;
01082 
01083    //Convert the second operand to a multiple precision integer
01084    value = (b >= 0) ? b : -b;
01085    t.sign = (b >= 0) ? 1 : -1;
01086    t.size = 1;
01087    t.data = &value;
01088 
01089    //Perform multiplication
01090    return mpiMul(r, a, &t);
01091 }
01092 
01093 
01094 /**
01095  * @brief Multiple precision division
01096  * @param[out] q The quotient Q = A / B
01097  * @param[out] r The remainder R = A mod B
01098  * @param[in] a The dividend A
01099  * @param[in] b The divisor B
01100  * @return Error code
01101  **/
01102 
01103 error_t mpiDiv(Mpi *q, Mpi *r, const Mpi *a, const Mpi *b)
01104 {
01105    error_t error;
01106    uint_t m;
01107    uint_t n;
01108    Mpi c;
01109    Mpi d;
01110    Mpi e;
01111 
01112    //Check whether the divisor is equal to zero
01113    if(!mpiCompInt(b, 0))
01114       return ERROR_INVALID_PARAMETER;
01115 
01116    //Initialize multiple precision integers
01117    mpiInit(&c);
01118    mpiInit(&d);
01119    mpiInit(&e);
01120 
01121    MPI_CHECK(mpiCopy(&c, a));
01122    MPI_CHECK(mpiCopy(&d, b));
01123    MPI_CHECK(mpiSetValue(&e, 0));
01124 
01125    m = mpiGetBitLength(&c);
01126    n = mpiGetBitLength(&d);
01127 
01128    if(m > n)
01129       MPI_CHECK(mpiShiftLeft(&d, m - n));
01130 
01131    while(n++ <= m)
01132    {
01133       MPI_CHECK(mpiShiftLeft(&e, 1));
01134 
01135       if(mpiComp(&c, &d) >= 0)
01136       {
01137          MPI_CHECK(mpiSetBitValue(&e, 0, 1));
01138          MPI_CHECK(mpiSub(&c, &c, &d));
01139       }
01140 
01141       MPI_CHECK(mpiShiftRight(&d, 1));
01142    }
01143 
01144    if(q != NULL)
01145       MPI_CHECK(mpiCopy(q, &e));
01146 
01147    if(r != NULL)
01148       MPI_CHECK(mpiCopy(r, &c));
01149 
01150 end:
01151    //Release previously allocated memory
01152    mpiFree(&c);
01153    mpiFree(&d);
01154    mpiFree(&e);
01155 
01156    //Return status code
01157    return error;
01158 }
01159 
01160 
01161 /**
01162  * @brief Divide a multiple precision integer by an integer
01163  * @param[out] q The quotient Q = A / B
01164  * @param[out] r The remainder R = A mod B
01165  * @param[in] a The dividend A
01166  * @param[in] b The divisor B
01167  * @return Error code
01168  **/
01169 
01170 error_t mpiDivInt(Mpi *q, Mpi *r, const Mpi *a, int_t b)
01171 {
01172    uint_t value;
01173    Mpi t;
01174 
01175    //Convert the divisor to a multiple precision integer
01176    value = (b >= 0) ? b : -b;
01177    t.sign = (b >= 0) ? 1 : -1;
01178    t.size = 1;
01179    t.data = &value;
01180 
01181    //Perform division
01182    return mpiDiv(q, r, a, &t);
01183 }
01184 
01185 
01186 /**
01187  * @brief Modulo operation
01188  * @param[out] r Resulting integer R = A mod P
01189  * @param[in] a The multiple precision integer to be reduced
01190  * @param[in] p The modulus P
01191  * @return Error code
01192  **/
01193 
01194 error_t mpiMod(Mpi *r, const Mpi *a, const Mpi *p)
01195 {
01196    error_t error;
01197    int_t sign;
01198    uint_t m;
01199    uint_t n;
01200    Mpi c;
01201 
01202    //Make sure the modulus is positive
01203    if(mpiCompInt(p, 0) <= 0)
01204       return ERROR_INVALID_PARAMETER;
01205 
01206    //Initialize multiple precision integer
01207    mpiInit(&c);
01208 
01209    //Save the sign of A
01210    sign = a->sign;
01211    //Determine the actual length of A
01212    m = mpiGetBitLength(a);
01213    //Determine the actual length of P
01214    n = mpiGetBitLength(p);
01215 
01216    //Let R = A
01217    MPI_CHECK(mpiCopy(r, a));
01218 
01219    if(m >= n)
01220    {
01221       MPI_CHECK(mpiCopy(&c, p));
01222       MPI_CHECK(mpiShiftLeft(&c, m - n));
01223 
01224       while(mpiCompAbs(r, p) >= 0)
01225       {
01226          if(mpiCompAbs(r, &c) >= 0)
01227          {
01228             MPI_CHECK(mpiSubAbs(r, r, &c));
01229          }
01230 
01231          MPI_CHECK(mpiShiftRight(&c, 1));
01232       }
01233    }
01234 
01235    if(sign < 0)
01236    {
01237       MPI_CHECK(mpiSubAbs(r, p, r));
01238    }
01239 
01240 end:
01241    //Release previously allocated memory
01242    mpiFree(&c);
01243 
01244    //Return status code
01245    return error;
01246 }
01247 
01248 
01249 
01250 /**
01251  * @brief Modular addition
01252  * @param[out] r Resulting integer R = A + B mod P
01253  * @param[in] a The first operand A
01254  * @param[in] b The second operand B
01255  * @param[in] p The modulus P
01256  * @return Error code
01257  **/
01258 
01259 error_t mpiAddMod(Mpi *r, const Mpi *a, const Mpi *b, const Mpi *p)
01260 {
01261    error_t error;
01262 
01263    //Perform modular addition
01264    MPI_CHECK(mpiAdd(r, a, b));
01265    MPI_CHECK(mpiMod(r, r, p));
01266 
01267 end:
01268    //Return status code
01269    return error;
01270 }
01271 
01272 
01273 /**
01274  * @brief Modular subtraction
01275  * @param[out] r Resulting integer R = A - B mod P
01276  * @param[in] a The first operand A
01277  * @param[in] b The second operand B
01278  * @param[in] p The modulus P
01279  * @return Error code
01280  **/
01281 
01282 error_t mpiSubMod(Mpi *r, const Mpi *a, const Mpi *b, const Mpi *p)
01283 {
01284    error_t error;
01285 
01286    //Perform modular subtraction
01287    MPI_CHECK(mpiSub(r, a, b));
01288    MPI_CHECK(mpiMod(r, r, p));
01289 
01290 end:
01291    //Return status code
01292    return error;
01293 }
01294 
01295 
01296 /**
01297  * @brief Modular multiplication
01298  * @param[out] r Resulting integer R = A * B mod P
01299  * @param[in] a The first operand A
01300  * @param[in] b The second operand B
01301  * @param[in] p The modulus P
01302  * @return Error code
01303  **/
01304 
01305 error_t mpiMulMod(Mpi *r, const Mpi *a, const Mpi *b, const Mpi *p)
01306 {
01307    error_t error;
01308 
01309    //Perform modular multiplication
01310    MPI_CHECK(mpiMul(r, a, b));
01311    MPI_CHECK(mpiMod(r, r, p));
01312 
01313 end:
01314    //Return status code
01315    return error;
01316 }
01317 
01318 
01319 /**
01320  * @brief Modular inverse
01321  * @param[out] r Resulting integer R = A^-1 mod P
01322  * @param[in] a The multiple precision integer A
01323  * @param[in] p The modulus P
01324  * @return Error code
01325  **/
01326 
01327 error_t mpiInvMod(Mpi *r, const Mpi *a, const Mpi *p)
01328 {
01329    error_t error;
01330    Mpi b;
01331    Mpi c;
01332    Mpi q0;
01333    Mpi r0;
01334    Mpi t;
01335    Mpi u;
01336    Mpi v;
01337 
01338    //Initialize multiple precision integers
01339    mpiInit(&b);
01340    mpiInit(&c);
01341    mpiInit(&q0);
01342    mpiInit(&r0);
01343    mpiInit(&t);
01344    mpiInit(&u);
01345    mpiInit(&v);
01346 
01347    MPI_CHECK(mpiCopy(&b, p));
01348    MPI_CHECK(mpiCopy(&c, a));
01349    MPI_CHECK(mpiSetValue(&u, 0));
01350    MPI_CHECK(mpiSetValue(&v, 1));
01351 
01352    while(mpiCompInt(&c, 0) > 0)
01353    {
01354       MPI_CHECK(mpiDiv(&q0, &r0, &b, &c));
01355 
01356       MPI_CHECK(mpiCopy(&b, &c));
01357       MPI_CHECK(mpiCopy(&c, &r0));
01358 
01359       MPI_CHECK(mpiCopy(&t, &v));
01360       MPI_CHECK(mpiMul(&q0, &q0, &v));
01361       MPI_CHECK(mpiSub(&v, &u, &q0));
01362       MPI_CHECK(mpiCopy(&u, &t));
01363    }
01364 
01365    if(mpiCompInt(&b, 1))
01366    {
01367       MPI_CHECK(ERROR_FAILURE);
01368    }
01369 
01370    if(mpiCompInt(&u, 0) > 0)
01371    {
01372       MPI_CHECK(mpiCopy(r, &u));
01373    }
01374    else
01375    {
01376       MPI_CHECK(mpiAdd(r, &u, p));
01377    }
01378 
01379 end:
01380    //Release previously allocated memory
01381    mpiFree(&b);
01382    mpiFree(&c);
01383    mpiFree(&q0);
01384    mpiFree(&r0);
01385    mpiFree(&t);
01386    mpiFree(&u);
01387    mpiFree(&v);
01388 
01389    //Return status code
01390    return error;
01391 }
01392 
01393 
01394 /**
01395  * @brief Modular exponentiation
01396  * @param[out] r Resulting integer R = A ^ E mod P
01397  * @param[in] a Pointer to a multiple precision integer
01398  * @param[in] e Exponent
01399  * @param[in] p Modulus
01400  * @return Error code
01401  **/
01402 
01403 error_t mpiExpMod(Mpi *r, const Mpi *a, const Mpi *e, const Mpi *p)
01404 {
01405    error_t error;
01406    int_t i;
01407    int_t j;
01408    int_t n;
01409    uint_t d;
01410    uint_t k;
01411    uint_t u;
01412    Mpi b;
01413    Mpi c2;
01414    Mpi t;
01415    Mpi s[8];
01416 
01417    //Initialize multiple precision integers
01418    mpiInit(&b);
01419    mpiInit(&c2);
01420    mpiInit(&t);
01421 
01422    //Initialize precomputed values
01423    for(i = 0; i < arraysize(s); i++)
01424       mpiInit(&s[i]);
01425 
01426    //Very small exponents are often selected with low Hamming weight.
01427    //The sliding window mechanism should be disabled in that case
01428    d = (mpiGetBitLength(e) <= 32) ? 1 : 4;
01429 
01430    //Even modulus?
01431    if(mpiIsEven(p))
01432    {
01433       //Let B = A^2
01434       MPI_CHECK(mpiMulMod(&b, a, a, p));
01435       //Let S[0] = A
01436       MPI_CHECK(mpiCopy(&s[0], a));
01437 
01438       //Precompute S[i] = A^(2 * i + 1)
01439       for(i = 1; i < (1 << (d - 1)); i++)
01440       {
01441          MPI_CHECK(mpiMulMod(&s[i], &s[i - 1], &b, p));
01442       }
01443 
01444       //Let R = 1
01445       MPI_CHECK(mpiSetValue(r, 1));
01446 
01447       //The exponent is processed in a right-to-left fashion
01448       i = mpiGetBitLength(e) - 1;
01449 
01450       //Perform sliding window exponentiation
01451       while(i >= 0)
01452       {
01453          //The sliding window exponentiation algorithm decomposes E
01454          //into zero and nonzero windows
01455          if(!mpiGetBitValue(e, i))
01456          {
01457             //Compute R = R^2
01458             MPI_CHECK(mpiMulMod(r, r, r, p));
01459             //Next bit to be processed
01460             i--;
01461          }
01462          else
01463          {
01464             //Find the longest window
01465             n = MAX(i - d + 1, 0);
01466 
01467             //The least significant bit of the window must be equal to 1
01468             while(!mpiGetBitValue(e, n)) n++;
01469 
01470             //The algorithm processes more than one bit per iteration
01471             for(u = 0, j = i; j >= n; j--)
01472             {
01473                //Compute R = R^2
01474                MPI_CHECK(mpiMulMod(r, r, r, p));
01475                //Compute the relevant index to be used in the precomputed table
01476                u = (u << 1) | mpiGetBitValue(e, j);
01477             }
01478 
01479             //Perform a single multiplication per iteration
01480             MPI_CHECK(mpiMulMod(r, r, &s[u >> 1], p));
01481             //Next bit to be processed
01482             i = n - 1;
01483          }
01484       }
01485    }
01486    else
01487    {
01488       //Compute the smaller C = (2^32)^k such as C > P
01489       k = mpiGetLength(p);
01490 
01491       //Compute C^2 mod P
01492       MPI_CHECK(mpiSetValue(&c2, 1));
01493       MPI_CHECK(mpiShiftLeft(&c2, 2 * k * (MPI_INT_SIZE * 8)));
01494       MPI_CHECK(mpiMod(&c2, &c2, p));
01495 
01496       //Let B = A * C mod P
01497       if(mpiComp(a, p) >= 0)
01498       {
01499          MPI_CHECK(mpiMod(&b, a, p));
01500          MPI_CHECK(mpiMontgomeryMul(&b, &b, &c2, k, p, &t));
01501       }
01502       else
01503       {
01504          MPI_CHECK(mpiMontgomeryMul(&b, a, &c2, k, p, &t));
01505       }
01506 
01507       //Let R = B^2 * C^-1 mod P
01508       MPI_CHECK(mpiMontgomeryMul(r, &b, &b, k, p, &t));
01509       //Let S[0] = B
01510       MPI_CHECK(mpiCopy(&s[0], &b));
01511 
01512       //Precompute S[i] = B^(2 * i + 1) * C^-1 mod P
01513       for(i = 1; i < (1 << (d - 1)); i++)
01514       {
01515          MPI_CHECK(mpiMontgomeryMul(&s[i], &s[i - 1], r, k, p, &t));
01516       }
01517 
01518       //Let R = C mod P
01519       MPI_CHECK(mpiCopy(r, &c2));
01520       MPI_CHECK(mpiMontgomeryRed(r, r, k, p, &t));
01521 
01522       //The exponent is processed in a right-to-left fashion
01523       i = mpiGetBitLength(e) - 1;
01524 
01525       //Perform sliding window exponentiation
01526       while(i >= 0)
01527       {
01528          //The sliding window exponentiation algorithm decomposes E
01529          //into zero and nonzero windows
01530          if(!mpiGetBitValue(e, i))
01531          {
01532             //Compute R = R^2 * C^-1 mod P
01533             MPI_CHECK(mpiMontgomeryMul(r, r, r, k, p, &t));
01534             //Next bit to be processed
01535             i--;
01536          }
01537          else
01538          {
01539             //Find the longest window
01540             n = MAX(i - d + 1, 0);
01541 
01542             //The least significant bit of the window must be equal to 1
01543             while(!mpiGetBitValue(e, n)) n++;
01544 
01545             //The algorithm processes more than one bit per iteration
01546             for(u = 0, j = i; j >= n; j--)
01547             {
01548                //Compute R = R^2 * C^-1 mod P
01549                MPI_CHECK(mpiMontgomeryMul(r, r, r, k, p, &t));
01550                //Compute the relevant index to be used in the precomputed table
01551                u = (u << 1) | mpiGetBitValue(e, j);
01552             }
01553 
01554             //Compute R = R * T[u/2] * C^-1 mod P
01555             MPI_CHECK(mpiMontgomeryMul(r, r, &s[u >> 1], k, p, &t));
01556             //Next bit to be processed
01557             i = n - 1;
01558          }
01559       }
01560 
01561       //Compute R = R * C^-1 mod P
01562       MPI_CHECK(mpiMontgomeryRed(r, r, k, p, &t));
01563    }
01564 
01565 end:
01566    //Release multiple precision integers
01567    mpiFree(&b);
01568    mpiFree(&c2);
01569    mpiFree(&t);
01570 
01571    //Release precomputed values
01572    for(i = 0; i < arraysize(s); i++)
01573       mpiFree(&s[i]);
01574 
01575    //Return status code
01576    return error;
01577 }
01578 
01579 
01580 /**
01581  * @brief Montgomery multiplication
01582  * @param[out] r Resulting integer R = A * B / 2^k mod P
01583  * @param[in] a An integer A such as 0 <= A < 2^k
01584  * @param[in] b An integer B such as 0 <= B < 2^k
01585  * @param[in] k An integer k such as P < 2^k
01586  * @param[in] p Modulus P
01587  * @param[in] t An preallocated integer T (for internal operation)
01588  * @return Error code
01589  **/
01590 
01591 error_t mpiMontgomeryMul(Mpi *r, const Mpi *a, const Mpi *b, uint_t k, const Mpi *p, Mpi *t)
01592 {
01593    error_t error;
01594    uint_t i;
01595    uint_t m;
01596    uint_t n;
01597    uint_t q;
01598 
01599    //Use Newton's method to compute the inverse of P[0] mod 2^32
01600    for(m = 2 - p->data[0], i = 0; i < 4; i++)
01601       m = m * (2 - m * p->data[0]);
01602 
01603    //Precompute -1/P[0] mod 2^32;
01604    m = ~m + 1;
01605 
01606    //We assume that B is always less than 2^k
01607    n = MIN(b->size, k);
01608 
01609    //Make sure T is large enough
01610    MPI_CHECK(mpiGrow(t, 2 * k + 1));
01611    //Let T = 0
01612    MPI_CHECK(mpiSetValue(t, 0));
01613 
01614    //Perform Montgomery multiplication
01615    for(i = 0; i < k; i++)
01616    {
01617       //Check current index
01618       if(i < a->size)
01619       {
01620          //Compute q = ((T[i] + A[i] * B[0]) * m) mod 2^32
01621          q = (t->data[i] + a->data[i] * b->data[0]) * m;
01622          //Compute T = T + A[i] * B
01623          mpiMulAccCore(t->data + i, b->data, n, a->data[i]);
01624       }
01625       else
01626       {
01627          //Compute q = (T[i] * m) mod 2^32
01628          q = t->data[i] * m;
01629       }
01630 
01631       //Compute T = T + q * P
01632       mpiMulAccCore(t->data + i, p->data, k, q);
01633    }
01634 
01635    //Compute R = T / 2^(32 * k)
01636    MPI_CHECK(mpiShiftRight(t, k * (MPI_INT_SIZE * 8)));
01637    MPI_CHECK(mpiCopy(r, t));
01638 
01639    //A final subtraction is required
01640    if(mpiComp(r, p) >= 0)
01641    {
01642       MPI_CHECK(mpiSub(r, r, p));
01643    }
01644 
01645 end:
01646    //Return status code
01647    return error;
01648 }
01649 
01650 
01651 /**
01652  * @brief Montgomery reduction
01653  * @param[out] r Resulting integer R = A / 2^k mod P
01654  * @param[in] a An integer A such as 0 <= A < 2^k
01655  * @param[in] k An integer k such as P < 2^k
01656  * @param[in] p Modulus P
01657  * @param[in] t An preallocated integer T (for internal operation)
01658  * @return Error code
01659  **/
01660 
01661 error_t mpiMontgomeryRed(Mpi *r, const Mpi *a, uint_t k, const Mpi *p, Mpi *t)
01662 {
01663    uint_t value;
01664    Mpi b;
01665 
01666    //Let B = 1
01667    value = 1;
01668    b.sign = 1;
01669    b.size = 1;
01670    b.data = &value;
01671 
01672    //Compute R = A / 2^k mod P
01673    return mpiMontgomeryMul(r, a, &b, k, p, t);
01674 }
01675 
01676 
01677 #if (MPI_ASM_SUPPORT == DISABLED)
01678 
01679 /**
01680  * @brief Multiply-accumulate operation
01681  * @param[out] r Resulting integer
01682  * @param[in] a First operand A
01683  * @param[in] m Size of A in words
01684  * @param[in] b Second operand B
01685  **/
01686 
01687 void mpiMulAccCore(uint_t *r, const uint_t *a, int_t m, const uint_t b)
01688 {
01689    int_t i;
01690    uint32_t c;
01691    uint32_t u;
01692    uint32_t v;
01693    uint64_t p;
01694 
01695    //Clear variables
01696    c = 0;
01697    u = 0;
01698    v = 0;
01699 
01700    //Perform multiplication
01701    for(i = 0; i < m; i++)
01702    {
01703       p = (uint64_t) a[i] * b;
01704       u = (uint32_t) p;
01705       v = (uint32_t) (p >> 32);
01706 
01707       u += c;
01708       if(u < c) v++;
01709 
01710       u += r[i];
01711       if(u < r[i]) v++;
01712 
01713       r[i] = u;
01714       c = v;
01715    }
01716 
01717    //Propagate carry
01718    for(; c != 0; i++)
01719    {
01720       r[i] += c;
01721       c = (r[i] < c);
01722    }
01723 }
01724 
01725 #endif
01726 
01727 
01728 /**
01729  * @brief Display the contents of a multiple precision integer
01730  * @param[in] stream Pointer to a FILE object that identifies an output stream
01731  * @param[in] prepend String to prepend to the left of each line
01732  * @param[in] a Pointer to a multiple precision integer
01733  **/
01734 
01735 void mpiDump(FILE *stream, const char_t *prepend, const Mpi *a)
01736 {
01737    uint_t i;
01738 
01739    //Process each word
01740    for(i = 0; i < a->size; i++)
01741    {
01742       //Beginning of a new line?
01743       if(i == 0 || ((a->size - i - 1) % 8) == 7)
01744          fprintf(stream, "%s", prepend);
01745 
01746       //Display current data
01747       fprintf(stream, "%08X ", a->data[a->size - 1 - i]);
01748 
01749       //End of current line?
01750       if(!((a->size - i - 1) % 8) || i == (a->size - 1))
01751          fprintf(stream, "\r\n");
01752    }
01753 }
01754 
01755 #endif
01756