bigdiv.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <string.h>
   4
   5 #define MAX 100
   6 #define BASE 10
   7
   8 int first(int *X)
   9 {
  10         /* Find the leftmost non-zero digit of X */
  11         int n = MAX-1;
  12         while (n >= 0 && !X[n])
  13                 n--;
  14         return n;
  15 }
  16
  17 int mul2(int *X, int n)
  18 {
  19         /* Multiply X by 2 and return its new length */
  20         int c = 0;
  21         for (int i=0; i<=n; i++) {
  22                 int x = 2*X[i] + c;
  23                 X[i] = x % BASE;
  24                 c = x / BASE;
  25         }
  26         if (c)
  27                 X[++n] = c;
  28         return n;
  29 }
  30
  31 int submul(int *X, int xn, int *Y, int yn, int z, int *T)
  32 {
  33         /* T = X-z*Y, returns -1 if it's <0; T and X can be the same array */
  34         if (!z)
  35                 return 0;
  36         if (yn > xn)
  37                 return -1;
  38         int c = 1, d = 0;
  39         for (int i=0; i<=xn; i++) {
  40                 int y = z*((i <= yn) ? Y[i] : 0) + d;
  41                 d = y / BASE;
  42                 y %= BASE;
  43                 int x = X[i] + (BASE-1-y) + c;
  44                 T[i] = x % BASE;
  45                 c = x / BASE;
  46         }
  47         return (c ? 0 : -1);
  48 }
  49
  50 void divide(int *X0, int *Y0, int *Z)
  51 {
  52         /* A copy of X and Y will be needed to avoid overwriting inputs */
  53         int X[MAX], Y[MAX];
  54         memcpy(X, X0, MAX * sizeof(*X));
  55         memcpy(Y, Y0, MAX * sizeof(*Y));
  56
  57         /* Find the leftmost digit of X and Y */
  58         int xn = first(X), yn = first(Y);
  59         if (yn < 0)             // Division by 0
  60                 exit(1);
  61
  62         /* While Y starts with <BASE/2, multiply both X and Y by 2 */
  63         while (Y[yn] < BASE/2) {
  64                 xn = mul2(X, xn);
  65                 yn = mul2(Y, yn);
  66         }
  67
  68         /* Find the size of result and fill Z with zeroes */
  69         for (int i=0; i<MAX; i++)
  70                 Z[i] = 0;
  71         int zn = xn - yn;
  72         if (zn < 0)             /* The result will obviously round to 0 */
  73                 return;
  74
  75         /* The main division loop */
  76         while (zn >= 0) {
  77                 int T[MAX];
  78                 while (xn > 0 && !X[xn])        /* Fix the estimate of size of X */
  79                         xn--;
  80                 int z = X[xn];
  81                 if (xn-zn > yn)                 /* First 1 or 2 digits of X */
  82                         z = 10*z + X[xn-1];
  83                 Z[zn] = z / Y[yn];
  84                 if (submul(X+zn, xn-zn, Y, yn, Z[zn], T) < 0)
  85                         Z[zn]--;
  86                 submul(X+zn, xn-zn, Y, yn, Z[zn], X+zn);
  87                 zn--;
  88         }
  89 }
  90
  91 int main(void)
  92 {
  93         int A[MAX] = { 6, 7, 5, 8, 4, 0, 1 };
  94         int B[MAX] = { 4, 6 };
  95         int C[MAX];
  96         divide(A, B, C);
  97         for (int i=MAX-1; i>=0; i--)
  98                 printf("%d", C[i]);
  99         putchar('\n');
 100         return 0;
 101 }