crypto/mpi/generic_udiv-w-sdiv.c

   1 /* mpihelp_udiv_w_sdiv -- implement udiv_qrnnd on machines with only signed
   2  *                        division.
   3  * Copyright (C) 1992, 1994, 1996, 1998 Free Software Foundation, Inc.
   4  * Contributed by Peter L. Montgomery.
   5  *
   6  * This file is part of GnuPG.
   7  *
   8  * GnuPG is free software; you can redistribute it and/or modify
   9  * it under the terms of the GNU General Public License as published by
  10  * the Free Software Foundation; either version 2 of the License, or
  11  * (at your option) any later version.
  12  *
  13  * GnuPG is distributed in the hope that it will be useful,
  14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  * GNU General Public License for more details.
  17  *
  18  * You should have received a copy of the GNU General Public License
  19  * along with this program; if not, write to the Free Software
  20  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA
  21  */
  22
  23 #include "mpi-internal.h"
  24 #include "longlong.h"
  25
  26
  27 #if 0  /* not yet ported to MPI */
  28
  29 mpi_limb_t
  30 mpihelp_udiv_w_sdiv( mpi_limp_t *rp,
  31                      mpi_limp_t *a1,
  32                      mpi_limp_t *a0,
  33                      mpi_limp_t *d   )
  34 {
  35   mp_limb_t q, r;
  36   mp_limb_t c0, c1, b1;
  37
  38   if ((mpi_limb_signed_t) d >= 0)
  39     {
  40       if (a1 < d - a1 - (a0 >> (BITS_PER_MP_LIMB - 1)))
  41         {
  42           /* dividend, divisor, and quotient are nonnegative */
  43           sdiv_qrnnd (q, r, a1, a0, d);
  44         }
  45       else
  46         {
  47           /* Compute c1*2^32 + c0 = a1*2^32 + a0 - 2^31*d */
  48           sub_ddmmss (c1, c0, a1, a0, d >> 1, d << (BITS_PER_MP_LIMB - 1));
  49           /* Divide (c1*2^32 + c0) by d */
  50           sdiv_qrnnd (q, r, c1, c0, d);
  51           /* Add 2^31 to quotient */
  52           q += (mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1);
  53         }
  54     }
  55   else
  56     {
  57       b1 = d >> 1;                      /* d/2, between 2^30 and 2^31 - 1 */
  58       c1 = a1 >> 1;                     /* A/2 */
  59       c0 = (a1 << (BITS_PER_MP_LIMB - 1)) + (a0 >> 1);
  60
  61       if (a1 < b1)                      /* A < 2^32*b1, so A/2 < 2^31*b1 */
  62         {
  63           sdiv_qrnnd (q, r, c1, c0, b1); /* (A/2) / (d/2) */
  64
  65           r = 2*r + (a0 & 1);           /* Remainder from A/(2*b1) */
  66           if ((d & 1) != 0)
  67             {
  68               if (r >= q)
  69                 r = r - q;
  70               else if (q - r <= d)
  71                 {
  72                   r = r - q + d;
  73                   q--;
  74                 }
  75               else
  76                 {
  77                   r = r - q + 2*d;
  78                   q -= 2;
  79                 }
  80             }
  81         }
  82       else if (c1 < b1)                 /* So 2^31 <= (A/2)/b1 < 2^32 */
  83         {
  84           c1 = (b1 - 1) - c1;
  85           c0 = ~c0;                     /* logical NOT */
  86
  87           sdiv_qrnnd (q, r, c1, c0, b1); /* (A/2) / (d/2) */
  88
  89           q = ~q;                       /* (A/2)/b1 */
  90           r = (b1 - 1) - r;
  91
  92           r = 2*r + (a0 & 1);           /* A/(2*b1) */
  93
  94           if ((d & 1) != 0)
  95             {
  96               if (r >= q)
  97                 r = r - q;
  98               else if (q - r <= d)
  99                 {
 100                   r = r - q + d;
 101                   q--;
 102                 }
 103               else
 104                 {
 105                   r = r - q + 2*d;
 106                   q -= 2;
 107                 }
 108             }
 109         }
 110       else                              /* Implies c1 = b1 */
 111         {                               /* Hence a1 = d - 1 = 2*b1 - 1 */
 112           if (a0 >= -d)
 113             {
 114               q = -1;
 115               r = a0 + d;
 116             }
 117           else
 118             {
 119               q = -2;
 120               r = a0 + 2*d;
 121             }
 122         }
 123     }
 124
 125   *rp = r;
 126   return q;
 127 }
 128
 129 #endif
 130