/******************************************************************************//*                                                                            *//* Functions for arithmetic and number theory with large integers in C        *//* Software supplement to the book "Cryptography in C and C++"                *//* by Michael Welschenbach, published by Apress Berkeley CA, 2005             *//*                                                                            *//* Module kmul.c           Revision: 19.12.2000                               *//*                                                                            *//*  Copyright (C) 1998-2005 by Michael Welschenbach                           *//*  Copyright (C) 2001-2005 by Springer-Verlag Berlin, Heidelberg             *//*  Copyright (C) 2001-2005 by Apress L.P., Berkeley, CA                      *//*  Copyright (C) 2002-2005 by Wydawnictwa MIKOM, Poland                      *//*  Copyright (C) 2002-2005 by PHEI, P.R.China                                *//*  Copyright (C) 2002-2005 by InfoBook, Korea                                *//*  Copyright (C) 2002-2005 by Triumph Publishing, Russia                     *//*                                                                            *//*  All Rights Reserved                                                       *//*                                                                            *//*  The software may be used for noncommercial purposes and may be altered,   *//*  as long as the following conditions are accepted without any              *//*  qualification:                                                            *//*                                                                            *//*  (1) All changes to the sources must be identified in such a way that the  *//*      changed software cannot be misinterpreted as the original software.   *//*                                                                            *//*  (2) The statements of copyright may not be removed or altered.            *//*                                                                            *//*  (3) The following DISCLAIMER is accepted:                                 *//*                                                                            *//*  DISCLAIMER:                                                               *//*                                                                            *//*  There is no warranty for the software contained in this distribution, to  *//*  the extent permitted by applicable law. The copyright holders provide the *//*  software `as is' without warranty of any kind, either expressed or        *//*  implied, including, but not limited to, the implied warranty of fitness   *//*  for a particular purpose. The entire risk as to the quality and           *//*  performance of the program is with you.                                   *//*                                                                            *//*  In no event unless required by applicable law or agreed to in writing     *//*  will the copyright holders, or any of the individual authors named in     *//*  the source files, be liable to you for damages, including any general,    *//*  special, incidental or consequential damages arising out of any use of    *//*  the software or out of inability to use the software (including but not   *//*  limited to any financial losses, loss of data or data being rendered      *//*  inaccurate or losses sustained by you or by third parties as a result of  *//*  a failure of the software to operate with any other programs), even if    *//*  such holder or other party has been advised of the possibility of such    *//*  damages.                                                                  *//*                                                                            *//******************************************************************************/#include <memory.h>#include "flint.h"#include "kmul.h"#define NO_ASSERTS#ifdef FLINT_DEBUG#undef NO_ASSERTS#define ASSERT_LOG_AND_QUIT#include "_assert.h"#include "_alloc.h"#ifdef COVERAGE#include "utclog.h"#endif#endif#ifdef NO_ASSERTS#define Assert(a) (void)0#endifstatic void kmul (clint *, clint *, int, int, CLINT);static void ksqr (clint *, int, CLINT);static void shiftadd (CLINT a_l, CLINT b_l, int l, CLINT s_l);static void addkar (clint *a_l, clint *b_l, int len_arg, CLINT s_l);#define MUL_THRESHOLD 40#define SQR_THRESHOLD 40CLINTD tmp_l;/******************************************************************************//*                                                                            *//*  Function:  Interface to Karatsuba multiplication                          *//*  Syntax:    void kmul_l (CLINT a_l, CLINT b_l, CLINT c_l);                 *//*  Input:     aa_l, bb_l (Factors)                                           *//*  Output:    p_l (Product)                                                  *//*  Returns:   E_CLINT_OK if everything is O.K.                               *//*             E_CLINT_OFL in case of Overflow                                *//*                                                                            *//******************************************************************************/int __FLINT_APIkmul_l (CLINT a_l, CLINT b_l, CLINT p_l){  CLINT aa_l, bb_l;  CLINTD pp_l;  int OFL = E_CLINT_OK;  cpy_l (aa_l, a_l);  cpy_l (bb_l, b_l);  kmul (LSDPTR_L(aa_l), LSDPTR_L(bb_l), DIGITS_L (aa_l), DIGITS_L (bb_l), pp_l);  if (DIGITS_L (pp_l) > (USHORT) CLINTMAXDIGIT)                /* Overflow ? */    {      ANDMAX_L (pp_l);                            /* Reduction modulo Nmax+1 */      OFL = E_CLINT_OFL;    }  cpy_l (p_l, pp_l);  ZEROCLINT_L (aa_l);  ZEROCLINT_L (bb_l);  ZEROCLINTD_L (pp_l);  ZEROCLINTD_L (tmp_l);  return OFL;}/******************************************************************************//*                                                                            *//*  Function:  Interface to Karatsuba squaring                                *//*  Syntax:    void ksqr_l (CLINT a_l, CLINT c_l);                            *//*  Input:     aa_l (Factor)                                                  *//*  Output:    p_l (Square)                                                   *//*  Returns:   E_CLINT_OK if everything is O.K.                               *//*             E_CLINT_OFL in case of overflow                                *//*                                                                            *//******************************************************************************/int __FLINT_APIksqr_l (CLINT a_l, CLINT p_l){  CLINT aa_l;  CLINTD pp_l;  int OFL = E_CLINT_OK;  cpy_l (aa_l, a_l);  ksqr (LSDPTR_L(aa_l), DIGITS_L (aa_l), pp_l);  if (DIGITS_L (pp_l) > (USHORT) CLINTMAXDIGIT)                /* Overflow ? */    {      ANDMAX_L (pp_l);                            /* Reduction modulo Nmax+1 */      OFL = E_CLINT_OFL;    }  cpy_l (p_l, pp_l);  ZEROCLINT_L (aa_l);  ZEROCLINTD_L (pp_l);  ZEROCLINTD_L (tmp_l);  return OFL;}/******************************************************************************//*                                                                            *//*  Function:  Karatsuba multiplication of two factors a_l and b_l            *//*             with 2k digits to base B                                       *//*  Syntax:    void kmul (clint *aptr_l, clint *bptr_l,                       *//*                        int len_a, int len_b, CLINT p_l);                   *//*  Input:     aptr_l (Pointer to least significant digit of a_l)             *//*             bptr_l (Pointer to least significant digit of b_l)             *//*             len_a (Number of digits of a_l)                                *//*             len_b (Number of digits of b_l)                                *//*  Output:    p_l (Product)                                                  *//*  Returns:   -                                                              *//*                                                                            *//******************************************************************************/static voidkmul (clint *aptr_l, clint *bptr_l, int len_a, int len_b, CLINT p_l){  CLINT c01_l, c10_l;  clint c0_l[CLINTMAXSHORT + 2], c1_l[CLINTMAXSHORT + 2], c2_l[CLINTMAXSHORT + 2];  clint *a1ptr_l, *b1ptr_l;  int l2;  if ((len_a == len_b) && (len_a >= MUL_THRESHOLD) && (0 == (len_a & 1)))    {      l2 = len_a/2;      a1ptr_l = aptr_l + l2;      b1ptr_l = bptr_l + l2;      kmul (aptr_l, bptr_l, l2, l2, c0_l);      kmul (a1ptr_l, b1ptr_l, l2, l2, c1_l);      addkar (a1ptr_l, aptr_l, l2, c01_l);      addkar (b1ptr_l, bptr_l, l2, c10_l);      kmul (LSDPTR_L (c01_l), LSDPTR_L (c10_l), DIGITS_L (c01_l), DIGITS_L (c10_l), c2_l);      sub (c2_l, c1_l, tmp_l);      sub (tmp_l, c0_l, c2_l);      shiftadd (c1_l, c2_l, l2, tmp_l);      shiftadd (tmp_l, c0_l, l2, p_l);    }  else /* Fallback to nonrecursive multiplication */    {      memcpy (LSDPTR_L (c1_l), aptr_l, len_a * sizeof (clint));      memcpy (LSDPTR_L (c2_l), bptr_l, len_b * sizeof (clint));      SETDIGITS_L (c1_l, len_a);      SETDIGITS_L (c2_l, len_b);      mult (c1_l, c2_l, p_l);    }  ZEROCLINT_L (c01_l);  ZEROCLINT_L (c10_l);  ZEROCLINT_L (c0_l);  ZEROCLINT_L (c1_l);  ZEROCLINT_L (c2_l);}/******************************************************************************//*                                                                            *//*  Function:  Karatsuba squaring of a factor a_l                             *//*             with 2k digits to base B                                       *//*  Syntax:    void kmul (clint *aptr_l, int len_a, CLINT p_l);               *//*  Input:     aptr_l (Pointer to least significant digit of a_l)             *//*             len_a (Number of digits of a_l)                                *//*  Output:    p_l (Square)                                                   *//*  Returns:   -                                                              *//*                                                                            *//******************************************************************************/static voidksqr (clint *aptr_l, int len_a, CLINT p_l){  CLINT c01_l;  clint c0_l[CLINTMAXSHORT + 2], c1_l[CLINTMAXSHORT + 2], c2_l[CLINTMAXSHORT + 2];  clint *a1ptr_l;  int l2;  if ((len_a >= SQR_THRESHOLD) && (0 == (len_a & 1)))    {      l2 = len_a/2;      a1ptr_l = aptr_l + l2;      ksqr (aptr_l, l2, c0_l);      ksqr (a1ptr_l, l2, c1_l);      addkar (a1ptr_l, aptr_l, l2, c01_l);      ksqr (LSDPTR_L (c01_l), DIGITS_L (c01_l), c2_l);      sub (c2_l, c1_l, tmp_l);      sub (tmp_l, c0_l, c2_l);      shiftadd (c1_l, c2_l, l2, tmp_l);      shiftadd (tmp_l, c0_l, l2, p_l);    }  else /* Fallback to nonrecursive squaring */    {      memcpy (LSDPTR_L (c1_l), aptr_l, len_a * sizeof (clint));      SETDIGITS_L (c1_l, len_a);      sqr (c1_l, p_l);    }  ZEROCLINT_L (c01_l);  ZEROCLINT_L (c0_l);  ZEROCLINT_L (c1_l);  ZEROCLINT_L (c2_l);}/******************************************************************************//*                                                                            *//*  Function:  Addition of two arguments                                      *//*             Parameters are pointers to LSD of arguments                    *//*  Syntax:    void addkar (clint *a_l, clint *b_l, int dgts, CLINT s_l);     *//*  Input:     a_l, b_l (Arguments with equal number of digits)               *//*  Output:    s_l (Sum)                                                      *//*  Returns:   -                                                              *//*                                                                            *//******************************************************************************/static voidaddkar (clint *a_l, clint *b_l, int dgts, CLINT s_l){  clint *msdptra_l;  clint *aptr_l, *bptr_l, *sptr_l = LSDPTR_L (s_l);  ULONG carry = 0L;  aptr_l = a_l;  bptr_l = b_l;  msdptra_l = a_l + dgts - 1;  SETDIGITS_L (s_l, dgts);  while (aptr_l <= msdptra_l)    {      *sptr_l++ = (USHORT) (carry = (ULONG) * aptr_l++ + (ULONG) * bptr_l++                            + (ULONG) (USHORT) (carry >> BITPERDGT));    }  if (carry & BASE)    {      *sptr_l = 1;      INCDIGITS_L (s_l);    }}/******************************************************************************//*                                                                            *//*  Function:  Addition of two arguments combined with shift of first         *//*             argument, w/o check for leading zeros                          *//*  Syntax:    void shiftadd (CLINT a_l, CLINT b_l, int dgts, CLINT s_l);     *//*  Input:     a_l, b_l (arguments)                                           *//*             dgts (Exponent to base 2)                                      *//*             Assumption: DIGITS_L (a_l) + dgts >= DIGITS_L (b_l))           *//*  Output:    s_l (Sum = (2^dgts)*a_l + b_l)                                 *//*  Returns:   -                                                              *//*                                                                            *//******************************************************************************/static voidshiftadd (CLINT a_l, CLINT b_l, int dgts, CLINT s_l){  clint *msdptra_l, *msdptrb_l;  clint *aptr_l, *bptr_l, *sptr_l = LSDPTR_L (s_l);  ULONG carry = 0L;  aptr_l = LSDPTR_L (a_l);  bptr_l = LSDPTR_L (b_l);  msdptra_l = MSDPTR_L (a_l);  msdptrb_l = MSDPTR_L (b_l);  SETDIGITS_L (s_l, DIGITS_L (a_l) + dgts);  while (bptr_l <= MIN (msdptrb_l, b_l + dgts))    {      *sptr_l++ = *bptr_l++;    }  while (sptr_l <= s_l + dgts)    {      *sptr_l++ = 0;    }  while (bptr_l <= msdptrb_l)    {      *sptr_l++ = (USHORT) (carry = (ULONG) * aptr_l++ + (ULONG) * bptr_l++                            + (ULONG) (USHORT) (carry >> BITPERDGT));    }  while (aptr_l <= msdptra_l)    {      *sptr_l++ = (USHORT) (carry = (ULONG) * aptr_l++                            + (ULONG) (USHORT) (carry >> BITPERDGT));    }  if (carry & BASE)    {      *sptr_l = 1;      INCDIGITS_L (s_l);    }  RMLDZRS_L (s_l);}