summaryrefslogtreecommitdiffstats
path: root/pbn_div.c
blob: f7b3a565cf804043ddd7dbb1c644078a5cd0e768 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
/* ----------------------------------------------------------------------- *
 *
 *   Copyright 2007 H. Peter Anvin - All Rights Reserved
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU Lesser General Public License as
 *   published by the Free Software Foundation, Inc.,
 *   59 Temple Place Ste 330, Boston MA 02111-1307, USA; version 2.1,
 *   incorporated herein by reference.
 *
 * ----------------------------------------------------------------------- */

/*
 * pbn_div.c
 */

/* quotient, remainder, numerator, denominator

   The sign on the outputs = sign(numerator) XOR sign(denominator)

   Returns -1 on error (divide by zero) and frees the inputs.

   This is doing it the slow way, which is bitwise long division.
   If this turns out to be a bottleneck, then it should be replaced
   by smarter algorithms, but this algorithm should at least be
   correct, so it should be usable for now. */

#include "pbn.h"

int pbn_div(struct pbn **qp, struct pbn **rp, struct pbn *n, struct pbn *d)
{
    struct pbn *q, *r;
    int rv = d->bits == 0 ? -1 : 0;
    int cm;
    int shift, bits, len, i;
    pbn_atom_t *qn, qm;

    if (rv || (cm = pbn_abscmp(n, d)) < 0) {
	/* We either had division by zero, or |n| < |d| */
	int rv = d->bits == 0;

	if (rp) {
	    if (d->minus && n->bits) {
		n = pbn_cow(n, n->len);
		n->minus ^= 1;
	    }
	    *rp = n;
	} else {
	    pbn_free(n);
	}
	pbn_free(d);

	if (qp)
	    *qp = pbn_int(0);

	return rv;		/* Division by zero? */
	}

    if (cm == 0) {
	/* Not really a common case, but we got it anyway */
	int quot = (n->minus ^ d->minus) ? -1 : 1;
	pbn_free(d);
	pbn_free(n);

	if (rp)
	    *rp = pbn_int(0);
	if (qp)
	    *qp = pbn_int(quot);

	return 0;
    }

    shift = n->bits - d->bits;
    bits = shift+1;		/* Maximum possible number of quotient bits */
    len = (bits+PBN_ATOM_BITS-1)/PBN_ATOM_BITS;

    r = pbn_cow(n, n->len);
    q = pbn_new(len);
    d = pbn_shl(d, shift);
    q->minus = r->minus ^ d->minus;
    r->minus = d->minus;	/* Makes our lives easier for now... */

    qm = (pbn_atom_t)1 << (shift % PBN_ATOM_BITS);
    qn = &q->num[shift/PBN_ATOM_BITS];
    for (i = shift; i >= 0; i--) {
	if (pbn_abscmp(r, d) >= 0) {
	    r = pbn_sub(r, d);
	    *qn |= qm;
	}
	d = pbn_shr(d, 1);
	qm >>= 1;
	if (!qm) {
	    qm = (pbn_atom_t)1 << (PBN_ATOM_BITS-1);
	    qn--;
	}
    }

    r->minus = q->minus;

    if (rp)
	*rp = pbn_adjust_bits(r);
    else
	pbn_free(r);

    if (qp)
	*qp = pbn_adjust_bits(q);
    else
	pbn_free(q);

    return 0;
}