bn_mp_karatsuba_sqr.c 3.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121
  1. #include <tommath.h>
  2. #ifdef BN_MP_KARATSUBA_SQR_C
  3. /* LibTomMath, multiple-precision integer library -- Tom St Denis
  4. *
  5. * LibTomMath is a library that provides multiple-precision
  6. * integer arithmetic as well as number theoretic functionality.
  7. *
  8. * The library was designed directly after the MPI library by
  9. * Michael Fromberger but has been written from scratch with
  10. * additional optimizations in place.
  11. *
  12. * The library is free for all purposes without any express
  13. * guarantee it works.
  14. *
  15. * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com
  16. */
  17. /* Karatsuba squaring, computes b = a*a using three
  18. * half size squarings
  19. *
  20. * See comments of karatsuba_mul for details. It
  21. * is essentially the same algorithm but merely
  22. * tuned to perform recursive squarings.
  23. */
  24. int mp_karatsuba_sqr (mp_int * a, mp_int * b)
  25. {
  26. mp_int x0, x1, t1, t2, x0x0, x1x1;
  27. int B, err;
  28. err = MP_MEM;
  29. /* min # of digits */
  30. B = a->used;
  31. /* now divide in two */
  32. B = B >> 1;
  33. /* init copy all the temps */
  34. if (mp_init_size (&x0, B) != MP_OKAY)
  35. goto ERR;
  36. if (mp_init_size (&x1, a->used - B) != MP_OKAY)
  37. goto X0;
  38. /* init temps */
  39. if (mp_init_size (&t1, a->used * 2) != MP_OKAY)
  40. goto X1;
  41. if (mp_init_size (&t2, a->used * 2) != MP_OKAY)
  42. goto T1;
  43. if (mp_init_size (&x0x0, B * 2) != MP_OKAY)
  44. goto T2;
  45. if (mp_init_size (&x1x1, (a->used - B) * 2) != MP_OKAY)
  46. goto X0X0;
  47. {
  48. register int x;
  49. register mp_digit *dst, *src;
  50. src = a->dp;
  51. /* now shift the digits */
  52. dst = x0.dp;
  53. for (x = 0; x < B; x++) {
  54. *dst++ = *src++;
  55. }
  56. dst = x1.dp;
  57. for (x = B; x < a->used; x++) {
  58. *dst++ = *src++;
  59. }
  60. }
  61. x0.used = B;
  62. x1.used = a->used - B;
  63. mp_clamp (&x0);
  64. /* now calc the products x0*x0 and x1*x1 */
  65. if (mp_sqr (&x0, &x0x0) != MP_OKAY)
  66. goto X1X1; /* x0x0 = x0*x0 */
  67. if (mp_sqr (&x1, &x1x1) != MP_OKAY)
  68. goto X1X1; /* x1x1 = x1*x1 */
  69. /* now calc (x1+x0)**2 */
  70. if (s_mp_add (&x1, &x0, &t1) != MP_OKAY)
  71. goto X1X1; /* t1 = x1 - x0 */
  72. if (mp_sqr (&t1, &t1) != MP_OKAY)
  73. goto X1X1; /* t1 = (x1 - x0) * (x1 - x0) */
  74. /* add x0y0 */
  75. if (s_mp_add (&x0x0, &x1x1, &t2) != MP_OKAY)
  76. goto X1X1; /* t2 = x0x0 + x1x1 */
  77. if (s_mp_sub (&t1, &t2, &t1) != MP_OKAY)
  78. goto X1X1; /* t1 = (x1+x0)**2 - (x0x0 + x1x1) */
  79. /* shift by B */
  80. if (mp_lshd (&t1, B) != MP_OKAY)
  81. goto X1X1; /* t1 = (x0x0 + x1x1 - (x1-x0)*(x1-x0))<<B */
  82. if (mp_lshd (&x1x1, B * 2) != MP_OKAY)
  83. goto X1X1; /* x1x1 = x1x1 << 2*B */
  84. if (mp_add (&x0x0, &t1, &t1) != MP_OKAY)
  85. goto X1X1; /* t1 = x0x0 + t1 */
  86. if (mp_add (&t1, &x1x1, b) != MP_OKAY)
  87. goto X1X1; /* t1 = x0x0 + t1 + x1x1 */
  88. err = MP_OKAY;
  89. X1X1:mp_clear (&x1x1);
  90. X0X0:mp_clear (&x0x0);
  91. T2:mp_clear (&t2);
  92. T1:mp_clear (&t1);
  93. X1:mp_clear (&x1);
  94. X0:mp_clear (&x0);
  95. ERR:
  96. return err;
  97. }
  98. #endif
  99. /* $Source: /cvs/libtom/libtommath/bn_mp_karatsuba_sqr.c,v $ */
  100. /* $Revision: 1.5 $ */
  101. /* $Date: 2006/03/31 14:18:44 $ */