ghc/runtime/gmp/mpn_div.c

   1 /* mpn_div -- Divide natural numbers, producing both remainder and
   2    quotient.
   3
   4 Copyright (C) 1991 Free Software Foundation, Inc.
   5
   6 This file is part of the GNU MP Library.
   7
   8 The GNU MP Library 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, or (at your option)
  11 any later version.
  12
  13 The GNU MP Library 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 the GNU MP Library; see the file COPYING.  If not, write to
  20 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
  21
  22 #include "gmp.h"
  23 #include "gmp-impl.h"
  24 #include "longlong.h"
  25
  26 /* Divide num (NUM_PTR/NUM_SIZE) by den (DEN_PTR/DEN_SIZE) and write
  27    the quotient at QUOT_PTR and the remainder at NUM_PTR.
  28
  29    Return 0 or 1, depending on if the quotient size is (NSIZE - DSIZE)
  30    or (NSIZE - DSIZE + 1).
  31
  32    Argument constraints:
  33    1. The most significant bit of d must be set.
  34    2. QUOT_PTR != DEN_PTR and QUOT_PTR != NUM_PTR, i.e. the quotient storage
  35       area must be distinct from either input operands.
  36
  37    The exact sizes of the quotient and remainder must be determined
  38    by the caller, in spite of the return value.  The return value just
  39    informs the caller about if the highest digit is written or not, and
  40    it may very well be 0.  */
  41
  42 /* THIS WILL BE IMPROVED SOON.  MORE COMMENTS AND FASTER CODE.  */
  43
  44 mp_size
  45 #ifdef __STDC__
  46 mpn_div (mp_ptr quot_ptr,
  47          mp_ptr num_ptr, mp_size num_size,
  48          mp_srcptr den_ptr, mp_size den_size)
  49 #else
  50 mpn_div (quot_ptr, num_ptr, num_size, den_ptr, den_size)
  51      mp_ptr quot_ptr;
  52      mp_ptr num_ptr;
  53      mp_size num_size;
  54      mp_srcptr den_ptr;
  55      mp_size den_size;
  56 #endif
  57 {
  58   mp_size q_is_long = 0;
  59
  60   switch (den_size)
  61     {
  62     case 0:
  63       /* We are asked to divide by zero, so go ahead and do it!
  64          (To make the compiler not remove this statement, assign NUM_SIZE
  65          and fall through.)  */
  66       num_size = 1 / den_size;
  67
  68     case 1:
  69       {
  70         mp_size i;
  71         mp_limb n1, n0;
  72         mp_limb d;
  73
  74         d = den_ptr[0];
  75         i = num_size - 1;
  76         n1 = num_ptr[i];
  77         i--;
  78
  79         if (n1 >= d)
  80           {
  81             q_is_long = 1;
  82             n1 = 0;
  83             i++;
  84           }
  85
  86         for (; i >= 0; i--)
  87           {
  88             n0 = num_ptr[i];
  89             udiv_qrnnd (quot_ptr[i], n1, n1, n0, d);
  90           }
  91
  92         num_ptr[0] = n1;
  93       }
  94       break;
  95
  96     case 2:
  97       {
  98         mp_size i;
  99         mp_limb n0, n1, n2;
 100         mp_limb d0, d1;
 101
 102         num_ptr += num_size - 2;
 103         d0 = den_ptr[1];
 104         d1 = den_ptr[0];
 105         n0 = num_ptr[1];
 106         n1 = num_ptr[0];
 107
 108         if (n0 >= d0)
 109           {
 110             q_is_long = 1;
 111             n1 = n0;
 112             n0 = 0;
 113             num_ptr++;
 114             num_size++;
 115           }
 116
 117         for (i = num_size - den_size - 1; i >= 0; i--)
 118           {
 119             mp_limb q;
 120             mp_limb r;
 121
 122             num_ptr--;
 123             if (n0 == d0)
 124               {
 125                 /* Q should be either 111..111 or 111..110.  Need special
 126                    treatment of this rare case as normal division would
 127                    give overflow.  */
 128                 q = ~0;
 129
 130                 r = n1 + d0;
 131                 if (r < d0)     /* Carry in the addition? */
 132                   {
 133                     n2 = num_ptr[0];
 134
 135                     add_ssaaaa (n0, n1, r - d1, n2, 0, d1);
 136                     quot_ptr[i] = q;
 137                     continue;
 138                   }
 139                 n0 = d1 - (d1 != 0);
 140                 n1 = -d1;
 141               }
 142             else
 143               {
 144                 udiv_qrnnd (q, r, n0, n1, d0);
 145                 umul_ppmm (n0, n1, d1, q);
 146               }
 147
 148             n2 = num_ptr[0];
 149           q_test:
 150             if (n0 > r || (n0 == r && n1 > n2))
 151               {
 152                 /* The estimated Q was too large.  */
 153                 q--;
 154
 155                 sub_ddmmss (n0, n1, n0, n1, 0, d1);
 156                 r += d0;
 157                 if (r >= d0)    /* If not carry, test q again.  */
 158                   goto q_test;
 159               }
 160
 161             quot_ptr[i] = q;
 162             sub_ddmmss (n0, n1, r, n2, n0, n1);
 163           }
 164         num_ptr[1] = n0;
 165         num_ptr[0] = n1;
 166       }
 167       break;
 168
 169     default:
 170       {
 171         mp_size i;
 172         mp_limb d0 = den_ptr[den_size - 1];
 173         mp_limb d1 = den_ptr[den_size - 2];
 174         mp_limb n0 = num_ptr[num_size - 1];
 175         int ugly_hack_flag = 0;
 176
 177         if (n0 >= d0)
 178           {
 179
 180             /* There's a problem with this case, which shows up later in the
 181                code.  q becomes 1 (and sometimes 0) the first time when
 182                we've been here, and w_cy == 0 after the main do-loops below.
 183                But c = num_ptr[j] reads rubbish outside the num_ptr vector!
 184                Maybe I can solve this cleanly when I fix the early-end
 185                optimization here in the default case.  For now, I change the
 186                add_back entering condition, to kludge.  Leaving the stray
 187                memref behind!
 188
 189                HACK: Added ugly_hack_flag to make it work.  */
 190
 191             q_is_long = 1;
 192             n0 = 0;
 193             num_size++;
 194             ugly_hack_flag = 1;
 195           }
 196
 197         num_ptr += num_size;
 198         den_ptr += den_size;
 199         for (i = num_size - den_size - 1; i >= 0; i--)
 200           {
 201             mp_limb q;
 202             mp_limb n1;
 203             mp_limb w_cy;
 204             mp_limb d, c;
 205             mp_size j;
 206
 207             num_ptr--;
 208             if (n0 == d0)
 209               /* This might over-estimate q, but it's probably not worth
 210                  the extra code here to find out.  */
 211               q = ~0;
 212             else
 213               {
 214                 mp_limb r;
 215
 216                 udiv_qrnnd (q, r, n0, num_ptr[-1], d0);
 217                 umul_ppmm (n1, n0, d1, q);
 218
 219                 while (n1 > r || (n1 == r && n0 > num_ptr[-2]))
 220                   {
 221                     q--;
 222                     r += d0;
 223                     if (r < d0) /* I.e. "carry in previous addition?"  */
 224                       break;
 225                     n1 -= n0 < d1;
 226                     n0 -= d1;
 227                   }
 228               }
 229
 230             w_cy = 0;
 231             j = -den_size;
 232             do
 233               {
 234                 d = den_ptr[j];
 235                 c = num_ptr[j];
 236                 umul_ppmm (n1, n0, d, q);
 237                 n0 += w_cy;
 238                 w_cy = (n0 < w_cy) + n1;
 239                 n0 = c - n0;
 240                 num_ptr[j] = n0;
 241                 if (n0 > c)
 242                   goto cy_loop;
 243               ncy_loop:
 244                 j++;
 245               }
 246             while  (j < 0);
 247
 248             if (ugly_hack_flag)
 249               {
 250                 c = 0;
 251                 ugly_hack_flag = 0;
 252               }
 253             else
 254               c = num_ptr[j];
 255             if (c >= w_cy)
 256               goto store_q;
 257             goto add_back;
 258
 259             do
 260               {
 261                 d = den_ptr[j];
 262                 c = num_ptr[j];
 263                 umul_ppmm (n1, n0, d, q);
 264                 n0 += w_cy;
 265                 w_cy = (n0 < w_cy) + n1;
 266                 n0 = c - n0 - 1;
 267                 num_ptr[j] = n0;
 268                 if (n0 < c)
 269                   goto ncy_loop;
 270               cy_loop:
 271                 j++;
 272               }
 273             while  (j < 0);
 274
 275             if (ugly_hack_flag)
 276               {
 277                 c = 0;
 278                 ugly_hack_flag = 0;
 279               }
 280             else
 281               c = num_ptr[j];
 282             w_cy++;
 283             if (c >= w_cy)
 284               goto store_q;
 285
 286           add_back:
 287             j = -den_size;
 288             do
 289               {
 290                 d = den_ptr[j];
 291                 n0 = num_ptr[j] + d;
 292                 num_ptr[j] = n0;
 293                 if (n0 < d)
 294                   goto ab_cy_loop;
 295               ab_ncy_loop:
 296                 j++;
 297               }
 298             while  (j < 0);
 299             abort ();           /* We should always have a carry out! */
 300
 301             do
 302               {
 303                 d = den_ptr[j];
 304                 n0 = num_ptr[j] + d + 1;
 305                 num_ptr[j] = n0;
 306                 if (n0 > d)
 307                   goto ab_ncy_loop;
 308               ab_cy_loop:
 309                 j++;
 310               }
 311             while  (j < 0);
 312             q--;
 313
 314           store_q:
 315             quot_ptr[i] = q;
 316           }
 317       }
 318     }
 319
 320   return q_is_long;
 321 }