usr/src/lib/libmp/common/mdiv.c

   1 /*      Copyright (c) 1984, 1986, 1987, 1988, 1989 AT&T */
   2 /*        All Rights Reserved   */
   3
   4
   5 /*
   6  * Copyright (c) 1980 Regents of the University of California.
   7  * All rights reserved.  The Berkeley software License Agreement
   8  * specifies the terms and conditions for redistribution.
   9  */
  10 /*      Portions Copyright(c) 1988, Sun Microsystems Inc.       */
  11 /*      All Rights Reserved                                     */
  12
  13 /*
  14  * Copyright (c) 1997, by Sun Microsystems, Inc.
  15  * All rights reserved.
  16  */
  17
  18 #ident  "%Z%%M% %I%     %E% SMI"        /* SVr4.0 1.1   */
  19
  20 /* LINTLIBRARY */
  21
  22 #include <mp.h>
  23 #include <stdio.h>
  24 #include <stdlib.h>
  25 #include <sys/types.h>
  26 #include "libmp.h"
  27
  28 static void m_div(MINT *, MINT *, MINT *, MINT *);
  29
  30 void
  31 mp_mdiv(MINT *a, MINT *b, MINT *q, MINT *r)
  32 {
  33         MINT x, y;
  34         int sign;
  35
  36         sign = 1;
  37         x.len = y.len = 0;
  38         _mp_move(a, &x);
  39         _mp_move(b, &y);
  40         if (x.len < 0) {
  41                 sign = -1;
  42                 x.len = -x.len;
  43         }
  44         if (y.len < 0) {
  45                 sign = -sign;
  46                 y.len = -y.len;
  47         }
  48         _mp_xfree(q);
  49         _mp_xfree(r);
  50         m_div(&x, &y, q, r);
  51         if (sign == -1) {
  52                 q->len = -q->len;
  53                 r->len = -r->len;
  54         }
  55         _mp_xfree(&x);
  56         _mp_xfree(&y);
  57 }
  58
  59 static int
  60 m_dsb(int qx, int n, short *a, short *b)
  61 {
  62         int borrow;
  63         int s3b2shit;
  64         int j;
  65         short fifteen = 15;
  66         short *aptr, *bptr;
  67 #ifdef DEBUGDSB
  68         (void) printf("m_dsb %d %d %d %d\n", qx, n, *a, *b);
  69 #endif
  70
  71         borrow = 0;
  72         aptr = a;
  73         bptr = b;
  74         for (j = n; j > 0; j--) {
  75 #ifdef DEBUGDSB
  76                 (void) printf("1 borrow=%x %d %d %d\n", borrow, (*aptr * qx),
  77                     *bptr, *aptr);
  78 #endif
  79                 borrow -= (*aptr++) * qx - *bptr;
  80 #ifdef DEBUGDSB
  81                 (void) printf("2 borrow=%x %d %d %d\n", borrow, (*aptr * qx),
  82                     *bptr, *aptr);
  83 #endif
  84                 *bptr++ = (short)(borrow & 077777);
  85 #ifdef DEBUGDSB
  86                 (void) printf("3 borrow=%x %d %d %d\n", borrow, (*aptr * qx),
  87                     *bptr, *aptr);
  88 #endif
  89                 if (borrow >= 0) borrow >>= fifteen; /* 3b2 */
  90                 else borrow = 0xfffe0000 | (borrow >> fifteen);
  91 #ifdef DEBUGDSB
  92                 (void) printf("4 borrow=%x %d %d %d\n", borrow, (*aptr * qx),
  93                     *bptr, *aptr);
  94 #endif
  95         }
  96         borrow += *bptr;
  97         *bptr = (short)(borrow & 077777);
  98                 if (borrow >= 0) s3b2shit = borrow >> fifteen; /* 3b2 */
  99                 else s3b2shit = 0xfffe0000 | (borrow >> fifteen);
 100         if (s3b2shit == 0) {
 101 #ifdef DEBUGDSB
 102         (void) printf("mdsb 0\n");
 103 #endif
 104                 return (0);
 105         }
 106         borrow = 0;
 107         aptr = a;
 108         bptr = b;
 109         for (j = n; j > 0; j--) {
 110                 borrow += *aptr++ + *bptr;
 111                 *bptr++ = (short)(borrow & 077777);
 112                 if (borrow >= 0) borrow >>= fifteen; /* 3b2 */
 113                 else borrow = 0xfffe0000 | (borrow >>fifteen);
 114         }
 115 #ifdef DEBUGDSB
 116         (void) printf("mdsb 1\n");
 117 #endif
 118         return (1);
 119 }
 120
 121 static int
 122 m_trq(short v1, short v2, short u1, short u2, short u3)
 123 {
 124         short d;
 125         int x1;
 126         int c1;
 127
 128         c1 = u1 * 0100000 + u2;
 129         if (u1 == v1) {
 130                 d = 077777;
 131         } else {
 132                 d = (short)(c1 / v1);
 133         }
 134         do {
 135                 x1 = c1 - v1 * d;
 136                 x1 = x1 * 0100000 + u3 - v2 * d;
 137                 --d;
 138         } while (x1 < 0);
 139 #ifdef DEBUGMTRQ
 140         (void) printf("mtrq %d %d %d %d %d %d\n", v1, v2, u1, u2, u3, (d+1));
 141 #endif
 142         return ((int)d + 1);
 143 }
 144
 145 static void
 146 m_div(MINT *a, MINT *b, MINT *q, MINT *r)
 147 {
 148         MINT u, v, x, w;
 149         short d;
 150         short *qval;
 151         short *uval;
 152         int j;
 153         int qq;
 154         int n;
 155         short v1;
 156         short v2;
 157
 158         u.len = v.len = x.len = w.len = 0;
 159         if (b->len == 0) {
 160                 _mp_fatal("mdiv divide by zero");
 161                 return;
 162         }
 163         if (b->len == 1) {
 164                 r->val = _mp_xalloc(1, "m_div1");
 165                 mp_sdiv(a, b->val[0], q, r->val);
 166                 if (r->val[0] == 0) {
 167                         free(r->val);
 168                         r->len = 0;
 169                 } else {
 170                         r->len = 1;
 171                 }
 172                 return;
 173         }
 174         if (a -> len < b -> len) {
 175                 q->len = 0;
 176                 r->len = a->len;
 177                 r->val = _mp_xalloc(r->len, "m_div2");
 178                 for (qq = 0; qq < r->len; qq++) {
 179                         r->val[qq] = a->val[qq];
 180                 }
 181                 return;
 182         }
 183         x.len = 1;
 184         x.val = &d;
 185         n = b->len;
 186         d = 0100000 / (b->val[n - 1] + 1);
 187         mp_mult(a, &x, &u); /* subtle: relies on mult allocing extra space */
 188         mp_mult(b, &x, &v);
 189 #ifdef DEBUG_MDIV
 190         (void) printf("  u=%s\n", mtox(&u));
 191         (void) printf("  v=%s\n", mtox(&v));
 192 #endif
 193         v1 = v.val[n - 1];
 194         v2 = v.val[n - 2];
 195         qval = _mp_xalloc(a -> len - n + 1, "m_div3");
 196         uval = u.val;
 197         for (j = a->len - n; j >= 0; j--) {
 198                 qq = m_trq(v1, v2, uval[j + n], uval[j + n - 1],
 199                                                         uval[j + n - 2]);
 200                 if (m_dsb(qq, n, v.val, uval + j))
 201                         qq -= 1;
 202                 qval[j] = (short)qq;
 203         }
 204         x.len = n;
 205         x.val = u.val;
 206         _mp_mcan(&x);
 207 #ifdef DEBUG_MDIV
 208         (void) printf("  x=%s\n", mtox(&x));
 209         (void) printf("  d(in)=%d\n", (d));
 210 #endif
 211         mp_sdiv(&x, d, &w, &d);
 212 #ifdef DEBUG_MDIV
 213         (void) printf("  w=%s\n", mtox(&w));
 214         (void) printf("  d(out)=%d\n", (d));
 215 #endif
 216         r->len = w.len;
 217         r->val = w.val;
 218         q->val = qval;
 219         qq = a->len - n + 1;
 220         if (qq > 0 && qval[qq - 1] == 0)
 221                 qq -= 1;
 222         q->len = qq;
 223         if (qq == 0)
 224                 free(qval);
 225         if (x.len != 0)
 226                 _mp_xfree(&u);
 227         _mp_xfree(&v);
 228 }