aboutsummaryrefslogtreecommitdiff
path: root/src/benchmarks/cfrac/pdivmod.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/benchmarks/cfrac/pdivmod.c')
-rw-r--r--src/benchmarks/cfrac/pdivmod.c315
1 files changed, 315 insertions, 0 deletions
diff --git a/src/benchmarks/cfrac/pdivmod.c b/src/benchmarks/cfrac/pdivmod.c
new file mode 100644
index 0000000..ce4a24b
--- /dev/null
+++ b/src/benchmarks/cfrac/pdivmod.c
@@ -0,0 +1,315 @@
+#include "pdefs.h"
+#include "precision.h"
+
+#ifdef DEBUG
+#include <stdio.h>
+#endif
+
+#ifdef ASM_16BIT
+#include "asm16bit.h"
+#endif
+
+/*
+ * Divide u (dividend) by v (divisor); If non-null, qp and rp are set to
+ * quotient and remainder. The result returned will be *qp, unless qp is
+ * NULL, then *rp will be returned if non-null, otherwise pUndef is returned.
+ *
+ * Produce:
+ *
+ * q (quotient) = u div v (v != 0)
+ * truncation is toward zero
+ *
+ * r (remainder) = u mod v
+ * = u - u div v * v (v != 0)
+ * = u (v == 0)
+ * ( e.g. u == q*v + r )
+ * remainder has same sign and dividend
+ *
+ * Note: this has opposite convention than the C standard div fuction,
+ * but the same convention of the typical C "/" operator
+ * It is also inconvienient for the mod function.
+ */
+/*
+ * This algorithm is taken almost verbatum from Knuth Vol 2.
+ * Please note the following trivial(?) array index
+ * transformations (since MSD to LSD order is reversed):
+ *
+ * q[0..m] to Q[0..m] thus q[i] == Q[m-i]
+ * r[1..n] R[0..n-1] r[i] == R[n+1-i]
+ * u[0..m+n] w[0..m+n] u[i] == w[m+n-i]
+ * v[1..n] x[0..n-1] v[i] == x[n-i]
+ *
+ * let N == n - 1 so that n == N + 1 thus:
+ *
+ * q[0..m] to Q[0..m] thus q[i] == Q[m-i]
+ * r[1..n] R[0..N] r[i] == R[N+2-i]
+ * u[0..m+n] w[0..m+N+1] u[i] == w[m+N+1-i]
+ * v[1..n] x[0..N] v[i] == x[N+1-i]
+ */
+
+/*
+ * Note: Be very observent of the usage of uPtr, and vPtr.
+ * They are used to point to u, v, w, q or r as necessary.
+ */
+precision pdivmod(u, v, qp, rp)
+ precision u, v, *qp, *rp;
+{
+ register digitPtr uPtr, vPtr, qPtr, LastPtr;
+
+ register accumulator temp; /* 0 <= temp < base^2 */
+ register digit carry; /* 0 <= carry < 2 */
+ register digit hi; /* 0 <= hi < base */
+
+ register posit n, m;
+ digit d; /* 0 <= d < base */
+ digit qd; /* 0 <= qd < base */
+#ifdef DEBUG
+ int i;
+#endif
+
+ precision q, r, w; /* quotient, remainder, temporary */
+
+ n = v->size; /* size of v and r */
+
+ (void) pparm(u);
+ (void) pparm(v);
+ if (u->size < n) {
+ q = pUndef;
+ r = pUndef;
+ pset(&q, pzero);
+ pset(&r, u);
+ goto done;
+ }
+
+ m = u->size - n;
+
+ uPtr = u->value + m + n;
+ vPtr = v->value + n;
+
+ q = palloc(m + 1);
+ if (q == pUndef) return q;
+
+ q->sign = (u->sign != v->sign); /* can generate -0 */
+
+ r = palloc(n);
+ if (r == pUndef) {
+ pdestroy(q);
+ return r;
+ }
+ r->sign = u->sign;
+/*
+ * watch out! does this function return: q=floor(a/b) or trunc(a/b)?
+ * it's currently the latter, but every mathmaticion I have talked to
+ * prefers the former so that a % b returns between 0 to b-1. The
+ * problem is that this is slower and disagrees with C common practice.
+ */
+ qPtr = q->value + m + 1;
+
+ if (n == 1) {
+ d = *--vPtr; /* d is only digit of v */
+ if (d == 0) { /* divide by zero? */
+ q = pnew(errorp(PDOMAIN, "pdivmod", "divide by zero"));
+ } else { /* single digit divide */
+#ifndef ASM_16BIT
+ vPtr = r->value + n;
+ hi = 0; /* hi is current remainder */
+ do {
+ temp = mulBase(hi); /* 0 <= temp <= (base-1)^2 */
+ temp += *--uPtr; /* 0 <= temp <= base(base-1) */
+ hi = uModDiv(temp, d, --qPtr); /* 0 <= hi < base */
+ } while (uPtr > u->value);
+ *--vPtr = hi;
+#else
+ qPtr -= m + 1;
+ *(r->value) = memdivw1(qPtr, u->value, m + 1, d);
+#endif
+ }
+ } else { /* muti digit divide */
+ /*
+ * normalize: multiply u and v by d so hi digit of v > b/2
+ */
+ d = BASE / (*--vPtr+1); /* high digit of v */
+
+ w = palloc(n); /* size of v */
+ if (w == pUndef) return w;
+
+#ifndef ASM_16BIT
+ vPtr = v->value;
+ uPtr = w->value; /* very confusing. just a temp */
+ LastPtr = vPtr + n;
+ hi = 0;
+ do { /* single digit multiply */
+ temp = uMul(*vPtr++, d); /* 0<= temp <= base(base-1)/2 */
+ temp += hi; /* 0 <= temp <= (base^2-1)/2 */
+ hi = divBase(temp); /* 0 <= hi < base / 2 */
+ *uPtr++ = modBase(temp); /* 0 <= hi < base / 2 */
+ } while (vPtr < LastPtr); /* on exit hi == 0 */
+#else
+ hi = memmulw1(w->value, v->value, n, d);
+#endif
+
+ pset(&v, w);
+ pdestroy(w);
+
+ w = palloc(m + n + 1);
+ if (w == pUndef) return w;
+
+#ifndef ASM_16BIT
+ uPtr = u->value;
+ vPtr = w->value; /* very confusing. just a temp */
+ LastPtr = uPtr + m + n;
+ do { /* single digit multiply */
+ temp = uMul(*uPtr++, d);
+ temp += hi;
+ hi = divBase(temp);
+ *vPtr++ = modBase(temp);
+ } while (uPtr < LastPtr);
+ *vPtr = hi; /* note extra digit */
+#else
+ hi = memmulw1(w->value, u->value, m + n, d);
+ w->value[m + n] = hi;
+#endif
+
+ pset(&u, w);
+ pdestroy(w);
+
+#ifdef DEBUG
+ printf("m = %d n = %d\nd = %d\n", m, n, d);
+ printf("norm u = "); pshow(u);
+ printf("norm v = "); pshow(v);
+#endif
+
+ uPtr = u->value + m + 1; /* current least significant digit */
+ do {
+ --uPtr;
+#ifdef DEBUG
+ printf(" u = ");
+ for (i = n; i >= 0; --i) printf("%.*x ", sizeof(digit) * 2, uPtr[i]);
+ putchar('\n');
+ printf(" v = ");
+ for (i = 1; i < 3; i++) printf("%.*x ", sizeof(digit) * 2,
+ v->value[n-i]);
+ putchar('\n');
+#endif
+#ifndef ASM_16BIT
+ vPtr = v->value + n;
+ LastPtr = uPtr + n;
+ if (*LastPtr == *--vPtr) { /* guess next digit */
+ qd = BASE - 1;
+ } else {
+ temp = mulBase(*LastPtr);
+ temp += *--LastPtr; /* 0 <= temp< base^2 */
+ temp = uModDiv(temp, *vPtr, &qd);
+ --vPtr;
+ --LastPtr;
+ while (uMul(*vPtr, qd) > mulBase(temp) + *LastPtr) {
+ --qd;
+ temp += vPtr[1];
+ if (temp >= BASE) break; /* if so, vPtr*qd <= temp*base */
+ }
+ LastPtr += 2;
+ }
+ /*
+ * Single digit Multiply then Subtract
+ */
+ vPtr = v->value;
+ carry = 1; /* noborrow bit */
+ hi = 0; /* hi digit of multiply */
+ do {
+ /* multiply */
+ temp = uMul(qd, *vPtr++); /* 0 <= temp <= (base-1)^2 */
+ temp += hi; /* 0 <= temp <= base(base-1) */
+ hi = divBase(temp);
+ temp = modBase(temp);
+ /* subtract */
+ temp = (BASE-1) - temp; /* 0 <= temp < base */
+ temp += *uPtr + carry; /* 0 <= temp < 2*base */
+ carry = divBase(temp);
+ *uPtr++ = modBase(temp); /* 0 <= carry < 2 */
+ } while (uPtr < LastPtr);
+ temp = (BASE-1) - hi;
+ temp += *uPtr + carry;
+ carry = divBase(temp);
+ *uPtr = modBase(temp);
+ uPtr -= n;
+#else
+#if 0
+ carry = !memmulsubw(uPtr, v->value, n, qd); /* 1 if noborrow */
+#endif
+ carry = !memdivw(uPtr, v->value, n, &qd); /* 1 if noborrow */
+#endif
+#ifdef DEBUG
+ printf(" qhat = %.*x\n", sizeof(digit) * 2, qd);
+ printf(" new u = ");
+ for (i = n; i >= 0; --i) printf("%.*x ", sizeof(digit) * 2, uPtr[i]);
+ putchar('\n');
+#endif
+ if (carry == 0) { /* Test remainder, add back */
+ vPtr = v->value;
+ LastPtr = uPtr + n;
+ do {
+ temp = *uPtr + *vPtr++;
+ temp += carry;
+ carry = divBase(temp);
+ *uPtr++ = modBase(temp);
+ } while (uPtr < LastPtr);
+ *uPtr += carry - BASE; /* real strange but works */
+ uPtr -= n;
+ --qd;
+#ifdef DEBUG
+ printf(" decrementing q...adding back\n");
+ printf(" fixed u = ");
+ for (i = n; i >= 0; --i) printf("%.*x ", sizeof(digit) * 2, uPtr[i]);
+ putchar('\n');
+ printf(" newq = %.*x\n", sizeof(digit) * 2, qd);
+#endif
+ }
+ *--qPtr = qd; /* one leading zero possible */
+#ifdef DEBUG
+ putchar('\n');
+#endif
+ } while (uPtr > u->value);
+
+ /*
+ * Un-normalize to get remainder
+ */
+#ifndef ASM_16BIT
+ uPtr = u->value + n; /* skip hi digit (it's zero) */
+ vPtr = r->value + n;
+ hi = 0; /* hi is current remainder */
+ do { /* single digit divide */
+ temp = mulBase(hi); /* 0<=temp < base^2-(base-1) */
+ temp += *--uPtr; /* 0 <= temp < base^2 */
+ hi = uModDiv(temp, d, --vPtr);
+ } while (uPtr > u->value); /* carry will be zero */
+#else
+ carry = memdivw1(r->value, u->value, n, d); /* always 0 */
+#endif
+ pnorm(r); /* remainder may have many leading 0's */
+ }
+
+ if (m > 0 && qPtr[m] == 0) {
+ --(q->size); /* normalize */
+ }
+ if (q->size == 1 && *qPtr == 0) q->sign = false;
+
+done:
+
+ pdestroy(u);
+ pdestroy(v);
+
+ if (rp == (precision *) -1) {
+ if (qp != pNull) pset(qp, q);
+ pdestroy(q);
+ return presult(r);
+ } else if (qp == (precision *) -1) {
+ if (rp != pNull) pset(rp, r);
+ pdestroy(r);
+ return presult(q);
+ }
+ if (qp != pNull) pset(qp, q);
+ if (rp != pNull) pset(rp, r);
+ pdestroy(q);
+ pdestroy(r);
+ return pUndef;
+}
generated by cgit v1.2.3 (git 2.25.1) at 2026-04-25 12:44:18 +0000