ghc/runtime/gmp/mpn_sub.c

   1 /* mpn_sub -- Subtract two low-level natural-number integers.
   2
   3 Copyright (C) 1991 Free Software Foundation, Inc.
   4
   5 This file is part of the GNU MP Library.
   6
   7 The GNU MP Library is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU General Public License as published by
   9 the Free Software Foundation; either version 2, or (at your option)
  10 any later version.
  11
  12 The GNU MP Library is distributed in the hope that it will be useful,
  13 but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 GNU General Public License for more details.
  16
  17 You should have received a copy of the GNU General Public License
  18 along with the GNU MP Library; see the file COPYING.  If not, write to
  19 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
  20
  21 #include "gmp.h"
  22 #include "gmp-impl.h"
  23
  24 /* Subtract SUB_PTR/SUB_SIZE from MIN_PTR/MIN_SIZE and store the
  25    result (MIN_SIZE words) at DIF_PTR.
  26
  27    Return 1 if min < sub (result is negative).  Otherwise, return the
  28    negative difference between the number of words in dif and min.
  29    (I.e.  return 0 if the result has MIN_SIZE words, -1 if it has
  30    MIN_SIZE - 1 words, etc.)
  31
  32    Argument constraint: MIN_SIZE >= SUB_SIZE.
  33
  34    The size of DIF can be calculated as MIN_SIZE + the return value.  */
  35
  36 mp_size
  37 #ifdef __STDC__
  38 mpn_sub (mp_ptr dif_ptr,
  39          mp_srcptr min_ptr, mp_size min_size,
  40          mp_srcptr sub_ptr, mp_size sub_size)
  41 #else
  42 mpn_sub (dif_ptr, min_ptr, min_size, sub_ptr, sub_size)
  43      mp_ptr dif_ptr;
  44      mp_srcptr min_ptr;
  45      mp_size min_size;
  46      mp_srcptr sub_ptr;
  47      mp_size sub_size;
  48 #endif
  49 {
  50   mp_limb m, s, dif;
  51   mp_size j;
  52
  53   /* The loop counter and index J goes from some negative value to zero.
  54      This way the loops are faster.  Need to offset the base pointers
  55      to take care of the negative indices.  */
  56
  57   j = -sub_size;
  58   if (j == 0)
  59     goto sub_finished;
  60
  61   min_ptr -= j;
  62   sub_ptr -= j;
  63   dif_ptr -= j;
  64
  65   /* There are two do-loops, marked NON-CARRY LOOP and CARRY LOOP that
  66      jump between each other.  The first loop is for when the previous
  67      subtraction didn't produce a carry-out; the second is for the
  68      complementary case.  */
  69
  70   /* NON-CARRY LOOP */
  71   do
  72     {
  73       m = min_ptr[j];
  74       s = sub_ptr[j];
  75       dif = m - s;
  76       dif_ptr[j] = dif;
  77       if (dif > m)
  78         goto cy_loop;
  79     ncy_loop:
  80       j++;
  81     }
  82   while (j < 0);
  83
  84   /* We have exhausted SUB, with no carry out.  Copy remaining part of
  85      MIN to DIF.  */
  86
  87  sub_finished:
  88   j = sub_size - min_size;
  89
  90   /* If there's no difference between the length of the operands, the
  91      last words might have become zero, and re-normalization is needed.  */
  92   if (j == 0)
  93     goto normalize;
  94
  95   min_ptr -= j;
  96   dif_ptr -= j;
  97
  98   goto copy;
  99
 100   /* CARRY LOOP */
 101   do
 102     {
 103       m = min_ptr[j];
 104       s = sub_ptr[j];
 105       dif = m - s - 1;
 106       dif_ptr[j] = dif;
 107       if (dif < m)
 108         goto ncy_loop;
 109     cy_loop:
 110       j++;
 111     }
 112   while (j < 0);
 113
 114   /* We have exhausted SUB, but need to propagate carry.  */
 115
 116   j = sub_size - min_size;
 117   if (j == 0)
 118     return 1;                   /* min < sub.  Flag it to the caller */
 119
 120   min_ptr -= j;
 121   dif_ptr -= j;
 122
 123   /* Propagate carry.  Sooner or later the carry will cancel with a
 124      non-zero word, because the minuend is normalized.  Considering this,
 125      there's no need to test the index J.  */
 126   for (;;)
 127     {
 128       m = min_ptr[j];
 129       dif = m - 1;
 130       dif_ptr[j] = dif;
 131       j++;
 132       if (dif < m)
 133         break;
 134     }
 135
 136   if (j == 0)
 137     goto normalize;
 138
 139  copy:
 140   /* Don't copy the remaining words of MIN to DIF if MIN_PTR and DIF_PTR
 141      are equal.  It would just be a no-op copying.  Return 0, as the length
 142      of the result equals that of the minuend.  */
 143   if (dif_ptr == min_ptr)
 144     return 0;
 145
 146   do
 147     {
 148       dif_ptr[j] = min_ptr[j];
 149       j++;
 150     }
 151   while (j < 0);
 152   return 0;
 153
 154  normalize:
 155   for (j = -1; j >= -min_size; j--)
 156     {
 157       if (dif_ptr[j] != 0)
 158         return j + 1;
 159     }
 160
 161   return -min_size;
 162 }