d0p1/ack
1
0
Fork 0
Code Issues Pull requests Packages Projects Releases Wiki Activity

Simplify flt_arith now that mantissa uses uint32_t.

Browse source 
It seems that someone wanted to build flt_arith with a compiler that
had long but not unsigned long.  This required extra code to
accomplish unsigned right shift, unsigned division, and unsigned
comparison using the signed operations.  Now that we use uint32_t, we
can simply use the unsigned operations and remove the ucmp() function.
We have similar code in mach/proto/fp/ and in
lang/cem/libcc.ansi/stdlib/ext_comp.c where we use the unsigned
operations.

Some long variables become uint32_t, and some masks with 0xFFFFFFFF
disappear because uint32_t has only 32 bits.

Update flt_arith.3 to show that mantissa uses uint32_t.

Provide a target to install modules/src/flt_arith/test.c as flt_test
so I can run the tests.
This commit is contained in:
George Koehler 2016-11-05 21:29:03 -04:00
parent 3bb41d3910
commit daeeb5aca3
11 changed files with 64 additions and 82 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

10
modules/src/flt_arith/b64_add.c
View file

@ -5,6 +5,7 @@
/* $Id$ */
#include <stdint.h>
#include "flt_misc.h"
int
@ -15,16 +16,15 @@ flt_b64_add(e1,e2)
int carry;
/* add higher pair of 32 bits */
overflow = ucmp((long)0xFFFFFFFF - e1->flt_h_32, e2->flt_h_32) < 0;
overflow = (0xFFFFFFFFUL - e1->flt_h_32 < e2->flt_h_32);
e1->flt_h_32 += e2->flt_h_32;
/* add lower pair of 32 bits */
carry = ucmp((long)0xFFFFFFFF - e1->flt_l_32, e2->flt_l_32) < 0;
carry = (0xFFFFFFFFUL - e1->flt_l_32 < e2->flt_l_32);
e1->flt_l_32 += e2->flt_l_32;
if ((carry) && ((++e1->flt_h_32 &~0xFFFFFFFF) || e1->flt_h_32 == 0)) {
e1->flt_h_32 = 0;
if ((carry) && (++e1->flt_h_32 == 0))
return(1); /* had a 64 bit overflow */
}
return(overflow); /* return status from higher add */
}

15
modules/src/flt_arith/b64_sft.c
View file

@ -23,12 +23,10 @@ flt_b64_sft(e,n)
n -= 32;
}
if (n > 0) {
e->flt_l_32 = (e->flt_l_32 >> 1) & 0x7FFFFFFF;
e->flt_l_32 >>= (n - 1);
e->flt_l_32 >>= n;
if (e->flt_h_32 != 0) {
e->flt_l_32 |= (e->flt_h_32 << (32 - n)) & 0xFFFFFFFF;
e->flt_h_32 = (e->flt_h_32 >> 1) & 0x7FFFFFFF;
e->flt_h_32 >>= (n - 1);
e->flt_l_32 |= (e->flt_h_32 << (32 - n));
e->flt_h_32 >>= n;
}
}
n = -n;
@ -38,11 +36,10 @@ flt_b64_sft(e,n)
n -= 32;
}
if (n > 0) {
e->flt_h_32 = (e->flt_h_32 << n) & 0xFFFFFFFF;
e->flt_h_32 <<= n;
if (e->flt_l_32 != 0) {
long l = (e->flt_l_32 >> 1) & 0x7FFFFFFF;
e->flt_h_32 |= (l >> (31 - n));
e->flt_l_32 = (e->flt_l_32 << n) & 0xFFFFFFFF;
e->flt_h_32 |= (e->flt_l_32 >> (32 - n));
e->flt_l_32 <<= n;
}
}
}

15
modules/src/flt_arith/build.lua
View file

@ -15,7 +15,6 @@ clibrary {
"./flt_umin.c",
"./flt_chk.c",
"./split.c",
"./ucmp.c",
},
hdrs = { "./flt_arith.h" },
deps = {
@ -24,4 +23,18 @@ clibrary {
}
}
-- The test program isn't built by default. Here is a target
-- modules/src/flt_arith+pkg to install it as flt_test
cprogram {
name = "test",
srcs = { "./test.c" },
deps = { "+lib" },
}
installable {
name = "pkg",
map = {
["$(INSDIR)/bin/flt_test"] = "+test",
},
}

9
modules/src/flt_arith/flt_add.c
View file

@ -44,11 +44,10 @@ flt_add(e1,e2,e3)
}
if (e1->flt_sign != e2->flt_sign) {
/* e2 + e1 = e2 - (-e1) */
int tmp = ucmp(e1->m1, e2->m1);
int tmp2 = ucmp(e1->m2, e2->m2);
if (tmp > 0 || (tmp == 0 && tmp2 > 0)) {
if (e1->m1 > e2->m1 ||
(e1->m1 == e2->m1 && e1->m2 > e2->m2)) {
/* abs(e1) > abs(e2) */
if (tmp2 < 0) {
if (e1->m2 < e2->m2) {
e1->m1 -= 1; /* carry in */
}
e1->m1 -= e2->m1;
@ -56,7 +55,7 @@ flt_add(e1,e2,e3)
*e3 = *e1;
}
else {
if (tmp2 > 0)
if (e1->m2 > e2->m2)
e2->m1 -= 1; /* carry in */
e2->m1 -= e1->m1;
e2->m2 -= e1->m2;

4
modules/src/flt_arith/flt_arith.3
View file

@ -9,8 +9,8 @@ flt_arith \- high precision floating point arithmetic
.if t .ta 3m 13m 22m
.if n .ta 5m 25m 40m
struct flt_mantissa {
long flt_h_32; /* high order 32 bits of mantissa */
long flt_l_32; /* low order 32 bits of mantissa */
uint32_t flt_h_32; /* high order 32 bits of mantissa */
uint32_t flt_l_32; /* low order 32 bits of mantissa */
};
typedef struct {

47
modules/src/flt_arith/flt_div.c
View file

@ -5,14 +5,15 @@
/* $Id$ */
#include <stdint.h>
#include "flt_misc.h"
void
flt_div(e1,e2,e3)
register flt_arith *e1,*e2,*e3;
{
long result[2];
register long *lp;
uint32_t result[2];
register uint32_t *rp;
unsigned short u[9], v[5];
register int j;
register unsigned short *u_p = u;
@ -44,7 +45,7 @@ flt_div(e1,e2,e3)
while (! v[maxv]) maxv--;
result[0] = 0;
result[1] = 0;
lp = result;
rp = result;
/*
* Use an algorithm of Knuth (The art of programming, Seminumerical
@ -52,35 +53,25 @@ flt_div(e1,e2,e3)
* with base 65536.
*/
for (j = 0; j <= 3; j++, u_p++) {
long q_est, temp;
long v1 = v[1];
uint32_t q_est, temp;
if (j == 2) lp++;
if (j == 2) rp++;
if (u_p[0] == 0 && u_p[1] < v[1]) continue;
temp = ((long)u_p[0] << 16) + u_p[1];
temp = ((uint32_t)u_p[0] << 16) + u_p[1];
if (u_p[0] >= v[1]) {
q_est = 0x0000FFFFL;
q_est = 0x0000FFFFUL;
}
else if (v[1] == 1) {
q_est = temp;
}
else if (temp >= 0) {
q_est = temp / v1;
}
else {
long rem;
q_est = (0x7FFFFFFF/v1)+((temp&0x7FFFFFFF)/v1);
rem = (0x7FFFFFFF%v1)+((temp&0x7FFFFFFF)%v1)+1;
while (rem >= v1) {
q_est++;
rem -= v1;
}
q_est = temp / v[1];
}
temp -= q_est * v1;
temp -= q_est * v[1];
while (!(temp&0xFFFF0000) &&
ucmp((long)v[2]*q_est,(temp<<16)+(long)u_p[2]) > 0) {
v[2]*q_est > (temp<<16)+u_p[2]) {
q_est--;
temp += v1;
temp += v[1];
}
/* Now, according to Knuth, we have an estimate of the
quotient, that is either correct or one too big, but
@ -88,11 +79,11 @@ flt_div(e1,e2,e3)
*/
if (q_est != 0) {
int i;
long k = 0;
uint32_t k = 0;
int borrow = 0;
for (i = maxv; i > 0; i--) {
long tmp = q_est * (long)v[i] + k + borrow;
uint32_t tmp = q_est * v[i] + k + borrow;
unsigned short md = tmp & 0xFFFF;
borrow = (md > u_p[i]);
@ -100,7 +91,7 @@ flt_div(e1,e2,e3)
k = (tmp >> 16) & 0xFFFF;
}
k += borrow;
borrow = (long)u_p[0] < k;
borrow = u_p[0] < k;
u_p[0] -= k;
if (borrow) {
@ -110,15 +101,15 @@ flt_div(e1,e2,e3)
q_est--;
borrow = 0;
for (i = maxv; i > 0; i--) {
long tmp
= v[i]+(long)u_p[i]+borrow;
uint32_t tmp
= v[i]+(uint32_t)u_p[i]+borrow;
u_p[i] = tmp & 0xFFFF;
borrow = (tmp >> 16) & 0xFFFF;
}
u_p[0] += borrow;
}
*lp |= (j & 1) ? q_est : (q_est<<16);
*rp |= (j & 1) ? q_est : (q_est<<16);
}
}
e3->m1 = result[0];

2
modules/src/flt_arith/flt_misc.h
View file

@ -15,13 +15,11 @@
#define EXT_MIN (-16384) /* min exponent */
/* hiding of names: */
#define ucmp _flt_ucmp
#define flt_nrm _flt_nrm
#define flt_chk _flt_chk
#define flt_b64_add _flt_64add
#define flt_split _flt_split
int ucmp(long, long);
void flt_nrm(flt_arith *);
void flt_chk(flt_arith *);
int flt_b64_add(struct flt_mantissa *, struct flt_mantissa *);

13
modules/src/flt_arith/flt_mul.c
View file

@ -5,6 +5,7 @@
/* $Id$ */
#include <stdint.h>
#include "flt_misc.h"
void
@ -43,9 +44,9 @@ flt_mul(e1,e2,e3)
*/
for(i=4, pres = &result[4];i--;pres--) if (mp[i]) {
unsigned short k = 0;
long mpi = mp[i];
uint32_t mpi = mp[i];
for(j=4;j--;) {
long tmp = (long)pres[j] + k;
long tmp = (uint32_t)pres[j] + k;
if (mc[j]) tmp += mpi * mc[j];
pres[j] = tmp & 0xFFFF;
k = (tmp >> 16) & 0xFFFF;
@ -64,12 +65,12 @@ flt_mul(e1,e2,e3)
/*
* combine the registers to a total
*/
e3->m1 = ((long)result[0] << 16) + result[1];
e3->m2 = ((long)result[2] << 16) + result[3];
e3->m1 = ((uint32_t)result[0] << 16) + result[1];
e3->m2 = ((uint32_t)result[2] << 16) + result[3];
if (result[4] & 0x8000) {
if (++e3->m2 == 0 || (e3->m2 & ~ 0xFFFFFFFF)) {
if (++e3->m2 == 0) {
e3->m2 = 0;
if (++e3->m1 == 0 || (e3->m1 & ~ 0xFFFFFFFF)) {
if (++e3->m1 == 0) {
e3->m1 = 0x80000000;
e3->flt_exp++;
}

5
modules/src/flt_arith/flt_nrm.c
View file

@ -5,6 +5,7 @@
/* $Id$ */
#include <stdint.h>
#include "flt_misc.h"
void
@ -24,8 +25,8 @@ flt_nrm(e)
e->m2 = 0L;
e->flt_exp -= 32;
}
if ((e->m1 & 0x80000000) == 0) {
long l = 0x40000000;
if ((e->m1 & 0x80000000UL) == 0) {
uint32_t l = 0x40000000UL;
int cnt = -1;
while (! (l & e->m1)) {

5
modules/src/flt_arith/test.c
View file

@ -1,3 +1,5 @@
#include <stdio.h>
#include <string.h>
#include "flt_arith.h"
struct tests {
@ -22,6 +24,7 @@ struct tests {
{ 0, 0, 0, 0}
};
int
main()
{
register struct tests *p = tests;
@ -31,7 +34,7 @@ main()
if (! dotest(p)) exit_status = 1;
p++;
}
exit(exit_status);
return exit_status;
}
int

21
modules/src/flt_arith/ucmp.c
View file

@ -1,21 +0,0 @@
/*
(c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands.
See the copyright notice in the ACK home directory, in the file "Copyright".
*/
/* $Id$ */
#include "flt_misc.h"
int
ucmp(l1,l2)
long l1,l2;
{
if (l1 == l2) return 0;
if (l2 >= 0) {
if (l1 > l2 || l1 < 0) return 1;
return -1;
}
if (l1 >= 0 || l1 < l2) return -1;
return 1;
}