New implementation of mp_invmod for even moduli. 3x-500x faster than

xgcd for even moduli.

4 files changed, 623 insertions, 183 deletions

diff --git a/security/nss/lib/freebl/mpi/mpi-priv.h b/security/nss/lib/freebl/mpi/mpi-priv.h
index 69e2f3374..becc253dc 100644
--- a/security/nss/lib/freebl/mpi/mpi-priv.h
+++ b/security/nss/lib/freebl/mpi/mpi-priv.h
@@ -217,13 +217,16 @@ mp_err   s_mp_exptmod(const mp_int *a, const mp_int *b, const mp_int *m, mp_int
 mp_err   s_mp_2expt(mp_int *a, mp_digit k);    /* a = 2^k                 */
 int      s_mp_cmp(const mp_int *a, const mp_int *b); /* magnitude comparison */
 int      s_mp_cmp_d(const mp_int *a, mp_digit d); /* magnitude digit compare */
-int      s_mp_ispow2(mp_int *v);               /* is v a power of 2?      */
+int      s_mp_ispow2(const mp_int *v);         /* is v a power of 2?      */
 int      s_mp_ispow2d(mp_digit d);             /* is d a power of 2?      */
 
 int      s_mp_tovalue(char ch, int r);          /* convert ch to value    */
 char     s_mp_todigit(mp_digit val, int r, int low); /* convert val to digit */
 int      s_mp_outlen(int bits, int r);          /* output length in bytes */
 mp_digit s_mp_invmod_radix(mp_digit P);   /* returns (P ** -1) mod RADIX */
+mp_err   s_mp_invmod_odd_m( const mp_int *a, const mp_int *m, mp_int *c);
+mp_err   s_mp_invmod_2d(    const mp_int *a, mp_size k,       mp_int *c);
+mp_err   s_mp_invmod_even_m(const mp_int *a, const mp_int *m, mp_int *c);
 
 /* ------ mpv functions, operate on arrays of digits, not on mp_int's ------ */
 void     s_mpv_mul_d(const mp_digit *a, mp_size a_len, mp_digit b, mp_digit *c);
@@ -241,3 +244,4 @@ mp_err   s_mpv_div_2dx1d(mp_digit Nhi, mp_digit Nlo, mp_digit divisor,
 
 /* }}} */
 #endif
+
diff --git a/security/nss/lib/freebl/mpi/mpi-test.c b/security/nss/lib/freebl/mpi/mpi-test.c
index 16e020ffa..e0fce438c 100644
--- a/security/nss/lib/freebl/mpi/mpi-test.c
+++ b/security/nss/lib/freebl/mpi/mpi-test.c
@@ -65,23 +65,23 @@
   for the comparison tests accordingly.  Most of the other tests
   should be fine as long as you re-compute the solutions, though.
  */
-char   *mp1  = "639A868CDA0C569861B";
-char   *mp2  = "AAFC0A3FE45E5E09DBE2C29";
-char   *mp3  = "B55AA8DF8A7E83241F38AC7A9E479CAEF2E4D7C5";
-char   *mp4  = "-63DBC2265B88268DC801C10EA68476B7BDE0090F";
-char   *mp5  = "F595CB42";
-char   *mp5a = "-4B597E";
-char   *mp6  = "0";
-char   *mp7  = "EBFA7121CD838CE6439CC59DDB4CBEF3";
-char   *mp8  = "5";
-char   *mp9  = "F74A2876A1432698923B0767DA19DCF3D71795EE";
-char   *mp10 = "9184E72A000";
-char   *mp11 = "54D79A3557E8";
-char   *mp12 = "10000000000000000";
-char   *mp13 = 
+const char   *mp1  = "639A868CDA0C569861B";
+const char   *mp2  = "AAFC0A3FE45E5E09DBE2C29";
+const char   *mp3  = "B55AA8DF8A7E83241F38AC7A9E479CAEF2E4D7C5";
+const char   *mp4  = "-63DBC2265B88268DC801C10EA68476B7BDE0090F";
+const char   *mp5  = "F595CB42";
+const char   *mp5a = "-4B597E";
+const char   *mp6  = "0";
+const char   *mp7  = "EBFA7121CD838CE6439CC59DDB4CBEF3";
+const char   *mp8  = "5";
+const char   *mp9  = "F74A2876A1432698923B0767DA19DCF3D71795EE";
+const char   *mp10 = "9184E72A000";
+const char   *mp11 = "54D79A3557E8";
+const char   *mp12 = "10000000000000000";
+const char   *mp13 = 
 "34584F700C15A341E40BF7BFDD88A6630C8FF2B2067469372D391342BDAB6163963C"
 "D5A5C79F708BDE26E0CCF2DB66CD6D6089E29A877C45F2B050D226E6DA88";
-char   *mp14 =
+const char   *mp14 =
 "AC3FA0EABAAC45724814D798942A1E28E14C81E0DE8055CED630E7689DA648683645DB6E"
 "458D9F5338CC3D4E33A5D1C9BF42780133599E60DEE0049AFA8F9489501AE5C9AA2B8C13"
 "FD21285A538B2CA87A626BB56E0A654C8707535E637FF4E39174157402BDE3AA30C9F134"
@@ -97,7 +97,7 @@ char   *mp14 =
 "F224E6874926C8D24D34B457FD2C9A586C6B99582DC24F787A39E3942786CF1D494B6EB4"
 "A513498CDA0B217C4E80BCE7DA1C704C35E071AC21E0DA9F57C27C3533F46A8D20B04137"
 "C1B1384BE4B2EB46";
-char   *mp15 = 
+const char   *mp15 = 
 "39849CF7FD65AF2E3C4D87FE5526221103D90BA26A6642FFE3C3ECC0887BBBC57E011BF1"
 "05D822A841653509C68F79EBE51C0099B8CBB04DEF31F36F5954208A3209AC122F0E11D8"
 "4AE67A494D78336A2066D394D42E27EF6B03DDAF6D69F5112C93E714D27C94F82FC7EF77"
@@ -113,21 +113,23 @@ char   *mp15 =
 "434ADBED36D54ACDFDFF70A4EFB46E285131FE725F1C637D1C62115EDAD01C4189716327"
 "BFAA79618B1656CBFA22C2C965687D0381CC2FE0245913C4D8D96108213680BD8E93E821"
 "822AD9DDBFE4BD04";
-char   *mp16 = "4A724340668DB150339A70";
-char   *mp17 = "8ADB90F58";
-char   *mp18 = "C64C230AB20E5";
-char *mp19 = "F1C9DACDA287F2E3C88DCE2393B8F53DAAAC1196DC36510962B6B59454CFE64B";
-char *mp20 = "D445662C8B6FE394107B867797750C326E0F4A967E135FC430F6CD7207913AC7";
-
-mp_digit md1 = 0;
-mp_digit md2 = 0x1;
-mp_digit md3 = 0x80;
-mp_digit md4 = 0x9C97;
-mp_digit md5 = 0xF5BF;
-mp_digit md6 = 0x14A0;
-mp_digit md7 = 0x03E8;
-mp_digit md8 = 0x0101;
-mp_digit md9 = 0xA;
+const char   *mp16 = "4A724340668DB150339A70";
+const char   *mp17 = "8ADB90F58";
+const char   *mp18 = "C64C230AB20E5";
+const char *mp19 = 
+"F1C9DACDA287F2E3C88DCE2393B8F53DAAAC1196DC36510962B6B59454CFE64B";
+const char *mp20 = 
+"D445662C8B6FE394107B867797750C326E0F4A967E135FC430F6CD7207913AC7";
+
+const mp_digit md1 = 0;
+const mp_digit md2 = 0x1;
+const mp_digit md3 = 0x80;
+const mp_digit md4 = 0x9C97;
+const mp_digit md5 = 0xF5BF;
+const mp_digit md6 = 0x14A0;
+const mp_digit md7 = 0x03E8;
+const mp_digit md8 = 0x0101;
+const mp_digit md9 = 0xA;
 
 /* 
    Solutions of the form x_mpABC, where:
@@ -142,18 +144,18 @@ mp_digit md9 = 0xA;
    it is a constant; otherwise, it is a full integer.  
  */
 
-char *p_mp12   = "4286AD72E095C9FE009938750743174ADDD7FD1E53";
-char *p_mp34   = "-46BDBD66CA108C94A8CF46C325F7B6E2F2BA82D35"
-                 "A1BFD6934C441EE369B60CA29BADC26845E918B";
-char *p_mp57   = "E260C265A0A27C17AD5F4E59D6E0360217A2EBA6";
-char *p_mp22   = "7233B5C1097FFC77CCF55928FDC3A5D31B712FDE7A1E91";
-char *p_mp1d4  = "3CECEA2331F4220BEF68DED";
-char *p_mp8d6  = "6720";
-char *p_mp1113 =
+const char *p_mp12   = "4286AD72E095C9FE009938750743174ADDD7FD1E53";
+const char *p_mp34   = "-46BDBD66CA108C94A8CF46C325F7B6E2F2BA82D35"
+                       "A1BFD6934C441EE369B60CA29BADC26845E918B";
+const char *p_mp57   = "E260C265A0A27C17AD5F4E59D6E0360217A2EBA6";
+const char *p_mp22   = "7233B5C1097FFC77CCF55928FDC3A5D31B712FDE7A1E91";
+const char *p_mp1d4  = "3CECEA2331F4220BEF68DED";
+const char *p_mp8d6  = "6720";
+const char *p_mp1113 =
 "11590FC3831C8C3C51813142C88E566408DB04F9E27642F6471A1822E0100B12F7F1"
 "5699A127C0FA9D26DCBFF458522661F30C6ADA4A07C8C90F9116893F6DBFBF24C3A2"
 "4340";
-char *p_mp1415 = 
+const char *p_mp1415 = 
 "26B36540DE8B3586699CCEAE218A2842C7D5A01590E70C4A26E789107FBCDB06AA2C"
 "6DDC39E6FA18B16FCB2E934C9A5F844DAD60EE3B1EA82199EC5E9608F67F860FB965"
 "736055DF0E8F2540EB28D07F47E309B5F5D7C94FF190AB9C83A6970160CA700B1081"
@@ -186,15 +188,15 @@ char *p_mp1415 =
 "DD0C08D3E3EBDF0AF54203B43AFDFC40D8FC79C97A4B0A4E1BEB14D8FCEFDDED8758"
 "6ED65B18";
 
-char *mp_mp345 = "B9B6D3A3";
-char *mp_mp335 = "16609C2D";
+const char *mp_mp345 = "B9B6D3A3";
+const char *mp_mp335 = "16609C2D";
 
-char *s_mp13   = "B55AA8DF8A7E83241F38B2B446B06A4FB84E5DE0";
-char *s_mp34   = "517EE6B92EF65C965736EB6BF7C325F73504CEB6";
-char *s_mp46   = "-63DBC2265B88268DC801C10EA68476B7BDE0090F";
-char *s_mp5d4  = "F59667D9";
-char *s_mp2d5  = "AAFC0A3FE45E5E09DBF21E8";
-char *s_mp1415 = 
+const char *s_mp13   = "B55AA8DF8A7E83241F38B2B446B06A4FB84E5DE0";
+const char *s_mp34   = "517EE6B92EF65C965736EB6BF7C325F73504CEB6";
+const char *s_mp46   = "-63DBC2265B88268DC801C10EA68476B7BDE0090F";
+const char *s_mp5d4  = "F59667D9";
+const char *s_mp2d5  = "AAFC0A3FE45E5E09DBF21E8";
+const char *s_mp1415 = 
 "E5C43DE2B811F4A084625F96E9504039E5258D8348E698CEB9F4D4292622042DB446"
 "F75F4B65C1FB7A317257FA354BB5A45E789AEC254EAECE11F80A53E3B513822491DB"
 "D9399DEC4807A2A3A10360129AC93F4A42388D3BF20B310DD0E9E9F4BE07FC88D53A"
@@ -212,24 +214,24 @@ char *s_mp1415 =
 "48A37FB13F84ED4FB7ACA18C4639EE64309BDD3D552AEB4AAF44295943DC1229A497"
 "A84A";
 
-char *ms_mp345 = "1E71E292";
-
-char *d_mp12   = "-AAFBA6A55DD183FD854A60E";
-char *d_mp34   = "119366B05E606A9B1E73A6D8944CC1366B0C4E0D4";
-char *d_mp5d4  = "F5952EAB";
-char *d_mp6d2  = "-1";
-char *md_mp345 = "26596B86";
-
-char *q_mp42   = "-95825A1FFA1A155D5";
-char *r_mp42   = "-6312E99D7700A3DCB32ADF2";
-char *q_mp45a  = "15344CDA3D841F661D2B61B6EDF7828CE36";
-char *r_mp45a  = "-47C47B";
-char *q_mp7c2  = "75FD3890E6C1C67321CE62CEEDA65F79";
-char *q_mp3d6  = "8CAFD53C272BD6FE8B0847BDC3B539EFAB5C3";
-char *r_mp3d6  = "1E5";
-char *r_mp5d5  = "1257";
-char *r_mp47   = "B3A9018D970281A90FB729A181D95CB8";
-char *q_mp1404 = 
+const char *ms_mp345 = "1E71E292";
+
+const char *d_mp12   = "-AAFBA6A55DD183FD854A60E";
+const char *d_mp34   = "119366B05E606A9B1E73A6D8944CC1366B0C4E0D4";
+const char *d_mp5d4  = "F5952EAB";
+const char *d_mp6d2  = "-1";
+const char *md_mp345 = "26596B86";
+
+const char *q_mp42   = "-95825A1FFA1A155D5";
+const char *r_mp42   = "-6312E99D7700A3DCB32ADF2";
+const char *q_mp45a  = "15344CDA3D841F661D2B61B6EDF7828CE36";
+const char *r_mp45a  = "-47C47B";
+const char *q_mp7c2  = "75FD3890E6C1C67321CE62CEEDA65F79";
+const char *q_mp3d6  = "8CAFD53C272BD6FE8B0847BDC3B539EFAB5C3";
+const char *r_mp3d6  = "1E5";
+const char *r_mp5d5  = "1257";
+const char *r_mp47   = "B3A9018D970281A90FB729A181D95CB8";
+const char *q_mp1404 = 
 "-1B994D869142D3EF6123A3CBBC3C0114FA071CFCEEF4B7D231D65591D32501AD80F"
 "FF49AE4EC80514CC071EF6B42521C2508F4CB2FEAD69A2D2EF3934087DCAF88CC4C4"
 "659F1CA8A7F4D36817D802F778F1392337FE36302D6865BF0D4645625DF8BB044E19"
@@ -245,25 +247,26 @@ char *q_mp1404 =
 "422299D21899A22F853B0C93081CC9925E350132A0717A611DD932A68A0ACC6E4C7F"
 "7F685EF8C1F4910AEA5DC00BB5A36FCA07FFEAA490C547F6E14A08FE87041AB803E1"
 "BD9E23E4D367A2C35762F209073DFF48F3";
-char *r_mp1404 = "12FF98621ABF63144BFFC3207AC8FC10D8D1A09";
-
-char *q_mp13c  = "34584F700C15A341E40BF7BFDD88A6630C8FF2B2067469372D391342"
-                 "BDAB6163963CD5A5C79F708BDE26E0CCF2DB66CD6D6089E29A877C45";
-char *r_mp13c  = "F2B050D226E6DA88";
-char *q_mp9c16 = "F74A2876A1432698923B0767DA19DCF3D71795E";
-char *r_mp9c16 = "E";
-
-char *e_mp5d9 = "A8FD7145E727A20E52E73D22990D35D158090307A"
-                "13A5215AAC4E9AB1E96BD34E531209E03310400";
-char *e_mp78  = "AA5F72C737DFFD8CCD108008BFE7C79ADC01A819B"
-                "32B75FB82EC0FB8CA83311DA36D4063F1E57857A2"
-                "1AB226563D84A15BB63CE975FF1453BD6750C58D9"
-                "D113175764F5D0B3C89B262D4702F4D9640A3";
-char *me_mp817 = "E504493ACB02F7F802B327AB13BF25";
-char *me_mp5d47 = "1D45ED0D78F2778157992C951DD2734C";
-char *me_mp1512 = "FB5B2A28D902B9D9";
-char *me_mp161718 = "423C6AC6DBD74";
-char *me_mp5114 =
+const char *r_mp1404 = "12FF98621ABF63144BFFC3207AC8FC10D8D1A09";
+
+const char *q_mp13c  = 
+		"34584F700C15A341E40BF7BFDD88A6630C8FF2B2067469372D391342"
+		"BDAB6163963CD5A5C79F708BDE26E0CCF2DB66CD6D6089E29A877C45";
+const char *r_mp13c  = "F2B050D226E6DA88";
+const char *q_mp9c16 = "F74A2876A1432698923B0767DA19DCF3D71795E";
+const char *r_mp9c16 = "E";
+
+const char *e_mp5d9 = "A8FD7145E727A20E52E73D22990D35D158090307A"
+		      "13A5215AAC4E9AB1E96BD34E531209E03310400";
+const char *e_mp78  = "AA5F72C737DFFD8CCD108008BFE7C79ADC01A819B"
+		      "32B75FB82EC0FB8CA83311DA36D4063F1E57857A2"
+		      "1AB226563D84A15BB63CE975FF1453BD6750C58D9"
+		      "D113175764F5D0B3C89B262D4702F4D9640A3";
+const char *me_mp817 = "E504493ACB02F7F802B327AB13BF25";
+const char *me_mp5d47 = "1D45ED0D78F2778157992C951DD2734C";
+const char *me_mp1512 = "FB5B2A28D902B9D9";
+const char *me_mp161718 = "423C6AC6DBD74";
+const char *me_mp5114 =
 "64F0F72807993578BBA3C7C36FFB184028F9EB9A810C92079E1498D8A80FC848E1F0"
 "25F1DE43B7F6AC063F5CC29D8A7C2D7A66269D72BF5CDC327AF88AF8EF9E601DCB0A"
 "3F35BFF3525FB1B61CE3A25182F17C0A0633B4089EA15BDC47664A43FEF639748AAC"
@@ -281,16 +284,16 @@ char *me_mp5114 =
 "B9DDA0CF4DFF35BB8D31245912BF4497FD0BD95F0C604E26EA5A8EA4F5EAE870A5BD"
 "FE8C";
 
-char *e_mpc2d3 = "100000000000000000000000000000000";
+const char *e_mpc2d3 = "100000000000000000000000000000000";
 
-char *t_mp9    = "FB9B6E32FF0452A34746";
-char *i_mp27   = "B6AD8DCCDAF92B6FE57D062FFEE3A99";
-char *i_mp2019 = 
+const char *t_mp9    = "FB9B6E32FF0452A34746";
+const char *i_mp27   = "B6AD8DCCDAF92B6FE57D062FFEE3A99";
+const char *i_mp2019 = 
 "BDF3D88DC373A63EED92903115B03FC8501910AF68297B4C41870AED3EA9F839";
 /* "15E3FE09E8AE5523AABA197BD2D16318D3CA148EDF4AE1C1C52FC96AFAF5680B"; */
 
 
-char *t_mp15 =
+const char *t_mp15 =
 "795853094E59B0008093BCA8DECF68587C64BDCA2F3F7F8963DABC12F1CFFFA9B8C4"
 "365232FD4751870A0EF6CA619287C5D8B7F1747D95076AB19645EF309773E9EACEA0"
 "975FA4AE16251A8DA5865349C3A903E3B8A2C0DEA3C0720B6020C7FED69AFF62BB72"
@@ -300,14 +303,14 @@ char *t_mp15 =
 "2496882877B069E877B59740DC1226F18A5C0F66F64A5F59A9FAFC5E9FC45AEC0E7A"
 "BEE244F7DD3AC268CF512A0E52E4F5BE5B94";
 
-char *g_mp71   = "1";
-char *g_mp25   = "7";
-char *l_mp1011 = "C589E3D7D64A6942A000";
+const char *g_mp71   = "1";
+const char *g_mp25   = "7";
+const char *l_mp1011 = "C589E3D7D64A6942A000";
 
 /* mp9 in radices from 5 to 64 inclusive */
 #define LOW_RADIX   5
 #define HIGH_RADIX  64
-char *v_mp9[] = {
+const char *v_mp9[] = {
   "404041130042310320100141302000203430214122130002340212132414134210033",
   "44515230120451152500101352430105520150025145320010504454125502",
   "644641136612541136016610100564613624243140151310023515322",
@@ -370,7 +373,7 @@ char *v_mp9[] = {
   "FTAA7QXGoQOaZi7PzePtFFN5vNk"
 };
 
-unsigned char b_mp4[] = {
+const unsigned char b_mp4[] = {
   0x01, 
 #if MP_DIGIT_MAX > MP_32BIT_MAX
   0x00, 0x00, 0x00, 0x00,
@@ -390,6 +393,7 @@ void reason(char *fmt, ...);
 /*------------------------------------------------------------------------*/
 
 char g_intbuf[4096];  /* buffer for integer comparison   */
+char a_intbuf[4096];  /* buffer for integer comparison   */
 int  g_verbose = 1;   /* print out reasons for failure?  */
 int  res;
 
@@ -1342,9 +1346,18 @@ int test_exptmod_d(void)
 
 int test_invmod(void)
 {
-  mp_int  a, m;
+  mp_int  a, m, c;
+  mp_int  p1, p2, p3, p4, p5;
+  mp_int  t1, t2, t3, t4;
   mp_err  res;
 
+  /* 5 128-bit primes. */
+  static const char ivp1[] = { "AAD8A5A2A2BEF644BAEE7DB0CA643719" };
+  static const char ivp2[] = { "CB371AD2B79A90BCC88D0430663E40B9" };
+  static const char ivp3[] = { "C6C818D4DF2618406CA09280C0400099" };
+  static const char ivp4[] = { "CE949C04512E68918006B1F0D7E93F27" };
+  static const char ivp5[] = { "F8EE999B6416645040687440E0B89F51" };
+
   mp_init(&a); mp_init(&m);
   mp_read_radix(&a, mp2, 16); mp_read_radix(&m, mp7, 16);
 
@@ -1371,6 +1384,216 @@ int test_invmod(void)
     return 1;
   }
 
+/* Need the following test cases:
+  Odd modulus
+    - a is odd,      relatively prime to m
+    - a is odd,  not relatively prime to m
+    - a is even,     relatively prime to m
+    - a is even, not relatively prime to m
+  Even modulus
+    - a is even  (should fail)
+    - a is odd,  not relatively prime to m
+    - a is odd,      relatively prime to m,
+      m is not a power of 2
+	- m has factor 2**k, k < 32
+	- m has factor 2**k, k > 32
+      m is a power of 2, 2**k
+	- k < 32
+	- k > 32
+*/
+
+  mp_init(&a);  mp_init(&m);  mp_init(&c);  
+  mp_init(&p1); mp_init(&p2); mp_init(&p3); mp_init(&p4); mp_init(&p5); 
+  mp_init(&t1); mp_init(&t2); mp_init(&t3); mp_init(&t4); 
+
+  mp_read_radix(&p1, ivp1, 16);
+  mp_read_radix(&p2, ivp2, 16);
+  mp_read_radix(&p3, ivp3, 16);
+  mp_read_radix(&p4, ivp4, 16);
+  mp_read_radix(&p5, ivp5, 16);
+
+  IFOK( mp_2expt(&t2, 68) );	/* t2 = 2**68 */
+  IFOK( mp_2expt(&t3, 128) );	/* t3 = 2**128 */
+  IFOK( mp_2expt(&t4, 31) );	/* t4 = 2**31 */
+
+/* test 3: Odd modulus - a is odd, relatively prime to m */
+
+  IFOK( mp_mul(&p1, &p2, &a) );
+  IFOK( mp_mul(&p3, &p4, &m) );
+  IFOK( mp_invmod(&a, &m, &t1) );
+  IFOK( mp_invmod_xgcd(&a, &m, &c) );
+
+  if (mp_cmp(&t1, &c) != 0) {
+    mp_toradix(&t1, g_intbuf, 16);
+    mp_toradix(&c,  a_intbuf, 16);
+    reason("error: invmod test 3 computed %s, expected %s\n", 
+           g_intbuf, a_intbuf);
+    return 1;
+  }
+  mp_clear(&a);  mp_clear(&t1); mp_clear(&c);  
+  mp_init(&a);   mp_init(&t1);  mp_init(&c);   
+
+/* test 4: Odd modulus - a is odd, NOT relatively prime to m */
+
+  IFOK( mp_mul(&p1, &p3, &a) );
+  /* reuse same m as before */
+
+  res = mp_invmod_xgcd(&a, &m, &c);
+  if (res != MP_UNDEF) 
+    goto CLEANUP4;
+
+  res = mp_invmod(&a, &m, &t1); /* we expect this to fail. */
+  if (res != MP_UNDEF) {
+CLEANUP4:
+    reason("error: invmod test 4 succeeded, should have failed.\n");
+    return 1;
+  }
+  mp_clear(&a);  mp_clear(&t1); mp_clear(&c);  
+  mp_init(&a);   mp_init(&t1);  mp_init(&c);   
+
+/* test 5: Odd modulus - a is even, relatively prime to m */
+
+  IFOK( mp_mul(&p1, &t2, &a) );
+  /* reuse m */
+  IFOK( mp_invmod(&a, &m, &t1) );
+  IFOK( mp_invmod_xgcd(&a, &m, &c) );
+
+  if (mp_cmp(&t1, &c) != 0) {
+    mp_toradix(&t1, g_intbuf, 16);
+    mp_toradix(&c,  a_intbuf, 16);
+    reason("error: invmod test 5 computed %s, expected %s\n", 
+           g_intbuf, a_intbuf);
+    return 1;
+  }
+  mp_clear(&a);  mp_clear(&t1); mp_clear(&c);  
+  mp_init(&a);   mp_init(&t1);  mp_init(&c);   
+
+/* test 6: Odd modulus - a is odd, NOT relatively prime to m */
+
+  /* reuse t2 */
+  IFOK( mp_mul(&t2, &p3, &a) );
+  /* reuse same m as before */
+
+  res = mp_invmod_xgcd(&a, &m, &c);
+  if (res != MP_UNDEF) 
+    goto CLEANUP6;
+
+  res = mp_invmod(&a, &m, &t1); /* we expect this to fail. */
+  if (res != MP_UNDEF) {
+CLEANUP6:
+    reason("error: invmod test 6 succeeded, should have failed.\n");
+    return 1;
+  }
+  mp_clear(&a);  mp_clear(&m); mp_clear(&c);  mp_clear(&t1); 
+  mp_init(&a);   mp_init(&m);  mp_init(&c);   mp_init(&t1); 
+
+/* test 7: Even modulus, even a, should fail */
+
+  IFOK( mp_mul(&p3, &t3, &m) ); /* even m */
+  /* reuse t2 */
+  IFOK( mp_mul(&p1, &t2, &a) ); /* even a */
+
+  res = mp_invmod_xgcd(&a, &m, &c);
+  if (res != MP_UNDEF) 
+    goto CLEANUP7;
+
+  res = mp_invmod(&a, &m, &t1); /* we expect this to fail. */
+  if (res != MP_UNDEF) {
+CLEANUP7:
+    reason("error: invmod test 7 succeeded, should have failed.\n");
+    return 1;
+  }
+  mp_clear(&a);  mp_clear(&c);  mp_clear(&t1); 
+  mp_init(&a);   mp_init(&c);   mp_init(&t1); 
+
+/* test 8: Even modulus    - a is odd,  not relatively prime to m */
+
+  /* reuse m */
+  IFOK( mp_mul(&p3, &p1, &a) ); /* even a */
+
+  res = mp_invmod_xgcd(&a, &m, &c);
+  if (res != MP_UNDEF) 
+    goto CLEANUP8;
+
+  res = mp_invmod(&a, &m, &t1); /* we expect this to fail. */
+  if (res != MP_UNDEF) {
+CLEANUP8:
+    reason("error: invmod test 8 succeeded, should have failed.\n");
+    return 1;
+  }
+  mp_clear(&a);  mp_clear(&m); mp_clear(&c);  mp_clear(&t1); 
+  mp_init(&a);   mp_init(&m);  mp_init(&c);   mp_init(&t1); 
+
+/* test 9: Even modulus    - m has factor 2**k, k < 32
+ *	                   - a is odd, relatively prime to m,
+ */
+  IFOK( mp_mul(&p3, &t4, &m) ); /* even m */
+  IFOK( mp_mul(&p1, &p2, &a) );
+  IFOK( mp_invmod(&a, &m, &t1) );
+  IFOK( mp_invmod_xgcd(&a, &m, &c) );
+
+  if (mp_cmp(&t1, &c) != 0) {
+    mp_toradix(&t1, g_intbuf, 16);
+    mp_toradix(&c,  a_intbuf, 16);
+    reason("error: invmod test 9 computed %s, expected %s\n", 
+           g_intbuf, a_intbuf);
+    return 1;
+  }
+  mp_clear(&m);  mp_clear(&t1); mp_clear(&c);  
+  mp_init(&m);   mp_init(&t1);  mp_init(&c);   
+
+/* test 10: Even modulus    - m has factor 2**k, k > 32
+ *	                    - a is odd, relatively prime to m,
+ */
+  IFOK( mp_mul(&p3, &t3, &m) ); /* even m */
+  /* reuse a */
+  IFOK( mp_invmod(&a, &m, &t1) );
+  IFOK( mp_invmod_xgcd(&a, &m, &c) );
+
+  if (mp_cmp(&t1, &c) != 0) {
+    mp_toradix(&t1, g_intbuf, 16);
+    mp_toradix(&c,  a_intbuf, 16);
+    reason("error: invmod test 10 computed %s, expected %s\n", 
+           g_intbuf, a_intbuf);
+    return 1;
+  }
+  mp_clear(&t1); mp_clear(&c);  
+  mp_init(&t1);  mp_init(&c);   
+
+/* test 11: Even modulus    - m is a power of 2, 2**k | k < 32
+ *                          - a is odd, relatively prime to m,
+ */
+  IFOK( mp_invmod(&a, &t4, &t1) );
+  IFOK( mp_invmod_xgcd(&a, &t4, &c) );
+
+  if (mp_cmp(&t1, &c) != 0) {
+    mp_toradix(&t1, g_intbuf, 16);
+    mp_toradix(&c,  a_intbuf, 16);
+    reason("error: invmod test 11 computed %s, expected %s\n", 
+           g_intbuf, a_intbuf);
+    return 1;
+  }
+  mp_clear(&t1); mp_clear(&c);  
+  mp_init(&t1);  mp_init(&c);   
+
+/* test 12: Even modulus    - m is a power of 2, 2**k | k > 32
+ *                          - a is odd, relatively prime to m,
+ */
+  IFOK( mp_invmod(&a, &t3, &t1) );
+  IFOK( mp_invmod_xgcd(&a, &t3, &c) );
+
+  if (mp_cmp(&t1, &c) != 0) {
+    mp_toradix(&t1, g_intbuf, 16);
+    mp_toradix(&c,  a_intbuf, 16);
+    reason("error: invmod test 12 computed %s, expected %s\n", 
+           g_intbuf, a_intbuf);
+    return 1;
+  }
+
+  mp_clear(&a);  mp_clear(&m);  mp_clear(&c);  
+  mp_clear(&t1); mp_clear(&t2); mp_clear(&t3); mp_clear(&t4); 
+  mp_clear(&p1); mp_clear(&p2); mp_clear(&p3); mp_clear(&p4); mp_clear(&p5); 
+
   return 0;
 }
 
diff --git a/security/nss/lib/freebl/mpi/mpi.c b/security/nss/lib/freebl/mpi/mpi.c
index 0a62c5511..ff1becfa4 100644
--- a/security/nss/lib/freebl/mpi/mpi.c
+++ b/security/nss/lib/freebl/mpi/mpi.c
@@ -947,10 +947,7 @@ CLEANUP:
   Compute q = a / b and r = a mod b.  Input parameters may be re-used
   as output parameters.  If q or r is NULL, that portion of the
   computation will be discarded (although it will still be computed)
-
-  Pay no attention to the hacker behind the curtain.
  */
-
 mp_err mp_div(const mp_int *a, const mp_int *b, mp_int *q, mp_int *r)
 {
   mp_err   res;
@@ -1865,11 +1862,12 @@ mp_err mp_lcm(mp_int *a, mp_int *b, mp_int *c)
   mp_xgcd(a, b, g, x, y)
 
   Compute g = (a, b) and values x and y satisfying Bezout's identity
-  (that is, ax + by = g).  This uses the extended binary GCD algorithm
+  (that is, ax + by = g).  This uses the binary extended GCD algorithm
   based on the Stein algorithm used for mp_gcd()
+  See algorithm 14.61 in Handbook of Applied Cryptogrpahy.
  */
 
-mp_err mp_xgcd(mp_int *a, mp_int *b, mp_int *g, mp_int *x, mp_int *y)
+mp_err mp_xgcd(const mp_int *a, const mp_int *b, mp_int *g, mp_int *x, mp_int *y)
 {
   mp_int   gx, xc, yc, u, v, A, B, C, D;
   mp_int  *clean[9];
@@ -1880,35 +1878,37 @@ mp_err mp_xgcd(mp_int *a, mp_int *b, mp_int *g, mp_int *x, mp_int *y)
     return MP_RANGE;
 
   /* Initialize all these variables we need */
-  if((res = mp_init(&u)) != MP_OKAY) goto CLEANUP;
+  MP_CHECKOK( mp_init(&u) );
   clean[++last] = &u;
-  if((res = mp_init(&v)) != MP_OKAY) goto CLEANUP;
+  MP_CHECKOK( mp_init(&v) );
   clean[++last] = &v;
-  if((res = mp_init(&gx)) != MP_OKAY) goto CLEANUP;
+  MP_CHECKOK( mp_init(&gx) );
   clean[++last] = &gx;
-  if((res = mp_init(&A)) != MP_OKAY) goto CLEANUP;
+  MP_CHECKOK( mp_init(&A) );
   clean[++last] = &A;
-  if((res = mp_init(&B)) != MP_OKAY) goto CLEANUP;
+  MP_CHECKOK( mp_init(&B) );
   clean[++last] = &B;
-  if((res = mp_init(&C)) != MP_OKAY) goto CLEANUP;
+  MP_CHECKOK( mp_init(&C) );
   clean[++last] = &C;
-  if((res = mp_init(&D)) != MP_OKAY) goto CLEANUP;
+  MP_CHECKOK( mp_init(&D) );
   clean[++last] = &D;
-  if((res = mp_init_copy(&xc, a)) != MP_OKAY) goto CLEANUP;
+  MP_CHECKOK( mp_init_copy(&xc, a) );
   clean[++last] = &xc;
   mp_abs(&xc, &xc);
-  if((res = mp_init_copy(&yc, b)) != MP_OKAY) goto CLEANUP;
+  MP_CHECKOK( mp_init_copy(&yc, b) );
   clean[++last] = &yc;
   mp_abs(&yc, &yc);
 
   mp_set(&gx, 1);
 
-  /* Divide by two until at least one of them is even */
+  /* Divide by two until at least one of them is odd */
   while(mp_iseven(&xc) && mp_iseven(&yc)) {
-    s_mp_div_2(&xc);
-    s_mp_div_2(&yc);
-    if((res = s_mp_mul_2(&gx)) != MP_OKAY)
-      goto CLEANUP;
+    mp_size nx = mp_trailing_zeros(&xc);
+    mp_size ny = mp_trailing_zeros(&yc);
+    mp_size n  = MP_MIN(nx, ny);
+    s_mp_div_2d(&xc,n);
+    s_mp_div_2d(&yc,n);
+    MP_CHECKOK( s_mp_mul_2d(&gx,n) );
   }
 
   mp_copy(&xc, &u);
@@ -1916,16 +1916,16 @@ mp_err mp_xgcd(mp_int *a, mp_int *b, mp_int *g, mp_int *x, mp_int *y)
   mp_set(&A, 1); mp_set(&D, 1);
 
   /* Loop through binary GCD algorithm */
-  for(;;) {
+  do {
     while(mp_iseven(&u)) {
       s_mp_div_2(&u);
 
       if(mp_iseven(&A) && mp_iseven(&B)) {
 	s_mp_div_2(&A); s_mp_div_2(&B);
       } else {
-	if((res = mp_add(&A, &yc, &A)) != MP_OKAY) goto CLEANUP;
+	MP_CHECKOK( mp_add(&A, &yc, &A) );
 	s_mp_div_2(&A);
-	if((res = mp_sub(&B, &xc, &B)) != MP_OKAY) goto CLEANUP;
+	MP_CHECKOK( mp_sub(&B, &xc, &B) );
 	s_mp_div_2(&B);
       }
     }
@@ -1936,39 +1936,33 @@ mp_err mp_xgcd(mp_int *a, mp_int *b, mp_int *g, mp_int *x, mp_int *y)
       if(mp_iseven(&C) && mp_iseven(&D)) {
 	s_mp_div_2(&C); s_mp_div_2(&D);
       } else {
-	if((res = mp_add(&C, &yc, &C)) != MP_OKAY) goto CLEANUP;
+	MP_CHECKOK( mp_add(&C, &yc, &C) );
 	s_mp_div_2(&C);
-	if((res = mp_sub(&D, &xc, &D)) != MP_OKAY) goto CLEANUP;
+	MP_CHECKOK( mp_sub(&D, &xc, &D) );
 	s_mp_div_2(&D);
       }
     }
 
     if(mp_cmp(&u, &v) >= 0) {
-      if((res = mp_sub(&u, &v, &u)) != MP_OKAY) goto CLEANUP;
-      if((res = mp_sub(&A, &C, &A)) != MP_OKAY) goto CLEANUP;
-      if((res = mp_sub(&B, &D, &B)) != MP_OKAY) goto CLEANUP;
-
+      MP_CHECKOK( mp_sub(&u, &v, &u) );
+      MP_CHECKOK( mp_sub(&A, &C, &A) );
+      MP_CHECKOK( mp_sub(&B, &D, &B) );
     } else {
-      if((res = mp_sub(&v, &u, &v)) != MP_OKAY) goto CLEANUP;
-      if((res = mp_sub(&C, &A, &C)) != MP_OKAY) goto CLEANUP;
-      if((res = mp_sub(&D, &B, &D)) != MP_OKAY) goto CLEANUP;
-
+      MP_CHECKOK( mp_sub(&v, &u, &v) );
+      MP_CHECKOK( mp_sub(&C, &A, &C) );
+      MP_CHECKOK( mp_sub(&D, &B, &D) );
     }
+  } while (mp_cmp_z(&u) != 0);
 
-    /* If we're done, copy results to output */
-    if(mp_cmp_z(&u) == 0) {
-      if(x)
-	if((res = mp_copy(&C, x)) != MP_OKAY) goto CLEANUP;
+  /* copy results to output */
+  if(x)
+    MP_CHECKOK( mp_copy(&C, x) );
 
-      if(y)
-	if((res = mp_copy(&D, y)) != MP_OKAY) goto CLEANUP;
+  if(y)
+    MP_CHECKOK( mp_copy(&D, y) );
       
-      if(g)
-	if((res = mp_mul(&gx, &v, g)) != MP_OKAY) goto CLEANUP;
-
-      break;
-    }
-  }
+  if(g)
+    MP_CHECKOK( mp_mul(&gx, &v, g) );
 
  CLEANUP:
   while(last >= 0)
@@ -1980,6 +1974,51 @@ mp_err mp_xgcd(mp_int *a, mp_int *b, mp_int *g, mp_int *x, mp_int *y)
 
 /* }}} */
 
+mp_size mp_trailing_zeros(const mp_int *mp)
+{
+  mp_digit d;
+  mp_size  n = 0;
+  int      ix;
+
+  if (!mp || !MP_DIGITS(mp) || !mp_cmp_z(mp))
+    return n;
+
+  for (ix = 0; !(d = MP_DIGIT(mp,ix)) && (ix < MP_USED(mp)); ++ix)
+    n += MP_DIGIT_BIT;
+  if (!d)
+    return 0;	/* shouldn't happen, but ... */
+#if (MP_DIGIT_MAX > MP_32BIT_MAX) 
+  if (!(d & 0xffffffffU)) {
+    d >>= 32;
+    n  += 32;
+  }
+#endif
+  if (!(d & 0xffffU)) {
+    d >>= 16;
+    n  += 16;
+  }
+  if (!(d & 0xffU)) {
+    d >>= 8;
+    n  += 8;
+  }
+  if (!(d & 0xfU)) {
+    d >>= 4;
+    n  += 4;
+  }
+  if (!(d & 0x3U)) {
+    d >>= 2;
+    n  += 2;
+  }
+  if (!(d & 0x1U)) {
+    d >>= 1;
+    n  += 1;
+  }
+#if MP_ARGCHK == 2
+  assert(0 != (d & 1));
+#endif
+  return n;
+}
+
 /* Given a and prime p, computes c and k such that a*c == 2**k (mod p).
 ** Returns k (positive) or error (negative).
 ** This technique from the paper "Fast Modular Reciprocals" (unpublished)
@@ -2009,9 +2048,14 @@ mp_err s_mp_almost_inverse(const mp_int *a, const mp_int *p, mp_int *c)
   for (;;) {
     int diff_sign;
     while (mp_iseven(&f)) {
-      s_mp_div_2d(&f, 1);
-      MP_CHECKOK( s_mp_mul_2d(&d, 1) );
-      ++k;
+      mp_size n = mp_trailing_zeros(&f);
+      if (!n) {
+	res = MP_UNDEF;
+	goto CLEANUP;
+      }
+      s_mp_div_2d(&f, n);
+      MP_CHECKOK( s_mp_mul_2d(&d, n) );
+      k += n;
     }
     if (mp_cmp_d(&f, 1) == MP_EQ) {	/* f == 1 */
       res = k;
@@ -2047,7 +2091,7 @@ CLEANUP:
   return res;
 }
 
-/* Compute T = (P ** -1) mod (2 ** 32).  Also works for 16-bit mp_digits.
+/* Compute T = (P ** -1) mod MP_RADIX.  Also works for 16-bit mp_digits.
 ** This technique from the paper "Fast Modular Reciprocals" (unpublished)
 ** by Richard Schroeppel (a.k.a. Captain Nemo).
 */
@@ -2107,17 +2151,46 @@ CLEANUP:
   return res;
 }
 
-/* {{{ mp_invmod(a, m, c) */
+/* compute mod inverse using Schroeppel's method, only if m is odd */
+mp_err s_mp_invmod_odd_m(const mp_int *a, const mp_int *m, mp_int *c)
+{
+  int k;
+  mp_err  res;
+  mp_int  x;
 
-/*
-  mp_invmod(a, m, c)
+  ARGCHK(a && m && c, MP_BADARG);
 
-  Compute c = a^-1 (mod m), if there is an inverse for a (mod m).
-  This is equivalent to the question of whether (a, m) = 1.  If not,
-  MP_UNDEF is returned, and there is no inverse.
- */
+  if(mp_cmp_z(a) == 0 || mp_cmp_z(m) == 0)
+    return MP_RANGE;
+  if (mp_iseven(m))
+    return MP_UNDEF;
+
+  MP_DIGITS(&x) = 0;
+
+  if (a == c) {
+    if ((res = mp_init_copy(&x, a)) != MP_OKAY)
+      return res;
+    if (a == m) 
+      m = &x;
+    a = &x;
+  } else if (m == c) {
+    if ((res = mp_init_copy(&x, m)) != MP_OKAY)
+      return res;
+    m = &x;
+  } else {
+    MP_DIGITS(&x) = 0;
+  }
+
+  MP_CHECKOK( s_mp_almost_inverse(a, m, c) );
+  k = res;
+  MP_CHECKOK( s_mp_fixup_reciprocal(c, m, k, c) );
+CLEANUP:
+  mp_clear(&x);
+  return res;
+}
 
-mp_err mp_invmod(mp_int *a, mp_int *m, mp_int *c)
+/* Known good algorithm for computing modular inverse.  But slow. */
+mp_err mp_invmod_xgcd(const mp_int *a, const mp_int *m, mp_int *c)
 {
   mp_int  g, x;
   mp_err  res;
@@ -2129,36 +2202,12 @@ mp_err mp_invmod(mp_int *a, mp_int *m, mp_int *c)
 
   MP_DIGITS(&g) = 0;
   MP_DIGITS(&x) = 0;
-
-  if (mp_isodd(m)) {
-    int k;
-
-    if (a == c) {
-      if ((res = mp_init_copy(&x, a)) != MP_OKAY)
-	return res;
-      if (a == m) 
-	m = &x;
-      a = &x;
-    } else if (m == c) {
-      if ((res = mp_init_copy(&x, m)) != MP_OKAY)
-	return res;
-      m = &x;
-    } else {
-      MP_DIGITS(&x) = 0;
-    }
-
-    MP_CHECKOK( s_mp_almost_inverse(a, m, c) );
-    k = res;
-    MP_CHECKOK( s_mp_fixup_reciprocal(c, m, k, c) );
-    goto CLEANUP;
-  }
-
   MP_CHECKOK( mp_init(&x) );
   MP_CHECKOK( mp_init(&g) );
 
   MP_CHECKOK( mp_xgcd(a, m, &g, &x, NULL) );
 
-  if(mp_cmp_d(&g, 1) != MP_EQ) {
+  if (mp_cmp_d(&g, 1) != MP_EQ) {
     res = MP_UNDEF;
     goto CLEANUP;
   }
@@ -2171,6 +2220,166 @@ CLEANUP:
   mp_clear(&g);
 
   return res;
+}
+
+/* modular inverse where modulus is 2**k. */
+/* c = a**-1 mod 2**k */
+mp_err s_mp_invmod_2d(const mp_int *a, mp_size k, mp_int *c)
+{
+  mp_err res;
+  mp_size ix = k + 4;
+  mp_int t0, t1, val, tmp, two2k;
+
+  static const mp_digit d2 = 2;
+  static const mp_int two = { MP_ZPOS, 1, 1, (mp_digit *)&d2 };
+
+  if (mp_iseven(a))
+    return MP_UNDEF;
+  if (k <= MP_DIGIT_BIT) {
+    mp_digit i = s_mp_invmod_radix(MP_DIGIT(a,0));
+    if (k < MP_DIGIT_BIT)
+      i &= ((mp_digit)1 << k) - (mp_digit)1;
+    mp_set(c, i);
+    return MP_OKAY;
+  }
+  MP_DIGITS(&t0) = 0;
+  MP_DIGITS(&t1) = 0;
+  MP_DIGITS(&val) = 0;
+  MP_DIGITS(&tmp) = 0;
+  MP_DIGITS(&two2k) = 0;
+  MP_CHECKOK( mp_init_copy(&val, a) );
+  s_mp_mod_2d(&val, k);
+  MP_CHECKOK( mp_init_copy(&t0, &val) );
+  MP_CHECKOK( mp_init_copy(&t1, &t0)  );
+  MP_CHECKOK( mp_init(&tmp) );
+  MP_CHECKOK( mp_init(&two2k) );
+  MP_CHECKOK( s_mp_2expt(&two2k, k) );
+  do {
+    MP_CHECKOK( mp_mul(&val, &t1, &tmp)  );
+    MP_CHECKOK( mp_sub(&two, &tmp, &tmp) );
+    MP_CHECKOK( mp_mul(&t1, &tmp, &t1)   );
+    s_mp_mod_2d(&t1, k);
+    while (MP_SIGN(&t1) != MP_ZPOS) {
+      MP_CHECKOK( mp_add(&t1, &two2k, &t1) );
+    }
+    if (mp_cmp(&t1, &t0) == MP_EQ) 
+      break;
+    MP_CHECKOK( mp_copy(&t1, &t0) );
+  } while (--ix > 0);
+  if (!ix) {
+    res = MP_UNDEF;
+  } else {
+    mp_exch(c, &t1);
+  }
+
+CLEANUP:
+  mp_clear(&t0);
+  mp_clear(&t1);
+  mp_clear(&val);
+  mp_clear(&tmp);
+  mp_clear(&two2k);
+  return res;
+}
+
+mp_err s_mp_invmod_even_m(const mp_int *a, const mp_int *m, mp_int *c)
+{
+  mp_err res;
+  mp_size k;
+  mp_int oddFactor, evenFactor;	/* factors of the modulus */
+  mp_int oddPart, evenPart;	/* parts to combine via CRT. */
+  mp_int C2, tmp1, tmp2;
+
+  static const mp_digit d1 = 1;
+  static const mp_int one = { MP_ZPOS, 1, 1, (mp_digit *)&d1 };
+
+  if ((res = s_mp_ispow2(m)) >= 0) {
+    k = res;
+    return s_mp_invmod_2d(a, k, c);
+  }
+  MP_DIGITS(&oddFactor) = 0;
+  MP_DIGITS(&evenFactor) = 0;
+  MP_DIGITS(&oddPart) = 0;
+  MP_DIGITS(&evenPart) = 0;
+  MP_DIGITS(&C2)     = 0;
+  MP_DIGITS(&tmp1)   = 0;
+  MP_DIGITS(&tmp2)   = 0;
+
+  MP_CHECKOK( mp_init_copy(&oddFactor, m) );    /* oddFactor = m */
+  MP_CHECKOK( mp_init(&evenFactor) );
+  MP_CHECKOK( mp_init(&oddPart) );
+  MP_CHECKOK( mp_init(&evenPart) );
+  MP_CHECKOK( mp_init(&C2)     );
+  MP_CHECKOK( mp_init(&tmp1)   );
+  MP_CHECKOK( mp_init(&tmp2)   );
+
+  k = mp_trailing_zeros(m);
+  s_mp_div_2d(&oddFactor, k);
+  MP_CHECKOK( s_mp_2expt(&evenFactor, k) );
+
+  /* compute a**-1 mod oddFactor. */
+  MP_CHECKOK( s_mp_invmod_odd_m(a, &oddFactor, &oddPart) );
+  /* compute a**-1 mod evenFactor, where evenFactor == 2**k. */
+  MP_CHECKOK( s_mp_invmod_2d(   a,       k,    &evenPart) );
+
+  /* Use Chinese Remainer theorem to compute a**-1 mod m. */
+  /* let m1 = oddFactor,  v1 = oddPart, 
+   * let m2 = evenFactor, v2 = evenPart.
+   */
+
+  /* Compute C2 = m1**-1 mod m2. */
+  MP_CHECKOK( s_mp_invmod_2d(&oddFactor, k,    &C2) );
+
+  /* compute u = (v2 - v1)*C2 mod m2 */
+  MP_CHECKOK( mp_sub(&evenPart, &oddPart,   &tmp1) );
+  MP_CHECKOK( mp_mul(&tmp1,     &C2,        &tmp2) );
+  s_mp_mod_2d(&tmp2, k);
+  while (MP_SIGN(&tmp2) != MP_ZPOS) {
+    MP_CHECKOK( mp_add(&tmp2, &evenFactor, &tmp2) );
+  }
+
+  /* compute answer = v1 + u*m1 */
+  MP_CHECKOK( mp_mul(&tmp2,     &oddFactor, c) );
+  MP_CHECKOK( mp_add(&oddPart,  c,          c) );
+  /* not sure this is necessary, but it's low cost if not. */
+  MP_CHECKOK( mp_mod(c,         m,          c) );
+
+CLEANUP:
+  mp_clear(&oddFactor);
+  mp_clear(&evenFactor);
+  mp_clear(&oddPart);
+  mp_clear(&evenPart);
+  mp_clear(&C2);
+  mp_clear(&tmp1);
+  mp_clear(&tmp2);
+  return res;
+}
+
+
+/* {{{ mp_invmod(a, m, c) */
+
+/*
+  mp_invmod(a, m, c)
+
+  Compute c = a^-1 (mod m), if there is an inverse for a (mod m).
+  This is equivalent to the question of whether (a, m) = 1.  If not,
+  MP_UNDEF is returned, and there is no inverse.
+ */
+
+mp_err mp_invmod(const mp_int *a, const mp_int *m, mp_int *c)
+{
+
+  ARGCHK(a && m && c, MP_BADARG);
+
+  if(mp_cmp_z(a) == 0 || mp_cmp_z(m) == 0)
+    return MP_RANGE;
+
+  if (mp_isodd(m)) {
+    return s_mp_invmod_odd_m(a, m, c);
+  }
+  if (mp_iseven(a))
+    return MP_UNDEF;	/* not invertable */
+
+  return s_mp_invmod_even_m(a, m, c);
 
 } /* end mp_invmod() */
 
@@ -2304,7 +2513,7 @@ mp_err mp_toraw(mp_int *mp, char *str)
   character or the end of the string.
  */
 
-mp_err  mp_read_radix(mp_int *mp, char *str, int radix)
+mp_err  mp_read_radix(mp_int *mp, const char *str, int radix)
 {
   int     ix = 0, val = 0;
   mp_err  res;
@@ -4167,7 +4376,7 @@ int      s_mp_cmp_d(const mp_int *a, mp_digit d)
   Returns -1 if the value is not a power of two; otherwise, it returns
   k such that v = 2^k, i.e. lg(v).
  */
-int      s_mp_ispow2(mp_int *v)
+int      s_mp_ispow2(const mp_int *v)
 {
   mp_digit d;
   int      extra = 0, ix;
diff --git a/security/nss/lib/freebl/mpi/mpi.h b/security/nss/lib/freebl/mpi/mpi.h
index 768c71f67..6f5e87128 100644
--- a/security/nss/lib/freebl/mpi/mpi.h
+++ b/security/nss/lib/freebl/mpi/mpi.h
@@ -246,8 +246,9 @@ int    mp_iseven(const mp_int *a);
 #if MP_NUMTH
 mp_err mp_gcd(mp_int *a, mp_int *b, mp_int *c);
 mp_err mp_lcm(mp_int *a, mp_int *b, mp_int *c);
-mp_err mp_xgcd(mp_int *a, mp_int *b, mp_int *g, mp_int *x, mp_int *y);
-mp_err mp_invmod(mp_int *a, mp_int *m, mp_int *c);
+mp_err mp_xgcd(const mp_int *a, const mp_int *b, mp_int *g, mp_int *x, mp_int *y);
+mp_err mp_invmod(const mp_int *a, const mp_int *m, mp_int *c);
+mp_err mp_invmod_xgcd(const mp_int *a, const mp_int *m, mp_int *c);
 #endif /* end MP_NUMTH */
 
 /* Input and output        */
@@ -259,7 +260,7 @@ void   mp_print(mp_int *mp, FILE *ofp);
 mp_err mp_read_raw(mp_int *mp, char *str, int len);
 int    mp_raw_size(mp_int *mp);
 mp_err mp_toraw(mp_int *mp, char *str);
-mp_err mp_read_radix(mp_int *mp, char *str, int radix);
+mp_err mp_read_radix(mp_int *mp, const char *str, int radix);
 int    mp_radix_size(mp_int *mp, int radix);
 mp_err mp_toradix(mp_int *mp, char *str, int radix);
 int    mp_tovalue(char ch, int r);
@@ -279,6 +280,9 @@ mp_err mp_to_unsigned_octets(const mp_int *mp, unsigned char *str, mp_size maxle
 mp_err mp_to_signed_octets(const mp_int *mp, unsigned char *str, mp_size maxlen);
 mp_err mp_to_fixlen_octets(const mp_int *mp, unsigned char *str, mp_size len);
 
+/* Miscellaneous */
+mp_size mp_trailing_zeros(const mp_int *mp);
+
 #define MP_CHECKOK(x)  if (MP_OKAY > (res = (x))) goto CLEANUP
 #define MP_CHECKERR(x) if (MP_OKAY > (res = (x))) goto CLEANUP