summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/regress/lib/libc/cephes/Makefile7
-rw-r--r--src/regress/lib/libc/cephes/drand.c174
-rw-r--r--src/regress/lib/libc/cephes/econst.c114
-rw-r--r--src/regress/lib/libc/cephes/eexp.c86
-rw-r--r--src/regress/lib/libc/cephes/ehead.h59
-rw-r--r--src/regress/lib/libc/cephes/elog.c110
-rw-r--r--src/regress/lib/libc/cephes/epow.c187
-rw-r--r--src/regress/lib/libc/cephes/etanh.c70
-rw-r--r--src/regress/lib/libc/cephes/etodec.c199
-rw-r--r--src/regress/lib/libc/cephes/ieee.c4153
-rw-r--r--src/regress/lib/libc/cephes/ieetst.c875
-rw-r--r--src/regress/lib/libc/cephes/mconf.h187
-rw-r--r--src/regress/lib/libc/cephes/mtherr.c114
13 files changed, 6335 insertions, 0 deletions
diff --git a/src/regress/lib/libc/cephes/Makefile b/src/regress/lib/libc/cephes/Makefile
new file mode 100644
index 0000000000..75cc85f4a8
--- /dev/null
+++ b/src/regress/lib/libc/cephes/Makefile
@@ -0,0 +1,7 @@
1# $OpenBSD: Makefile,v 1.1 2011/07/02 18:11:01 martynas Exp $
2
3PROG = ieetst
4SRCS = drand.c econst.c eexp.c elog.c epow.c etanh.c etodec.c ieee.c \
5 ieetst.c mtherr.c
6
7.include <bsd.regress.mk>
diff --git a/src/regress/lib/libc/cephes/drand.c b/src/regress/lib/libc/cephes/drand.c
new file mode 100644
index 0000000000..7f7000b4e8
--- /dev/null
+++ b/src/regress/lib/libc/cephes/drand.c
@@ -0,0 +1,174 @@
1/* $OpenBSD: drand.c,v 1.1 2011/07/02 18:11:01 martynas Exp $ */
2
3/*
4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19/* drand.c
20 *
21 * Pseudorandom number generator
22 *
23 *
24 *
25 * SYNOPSIS:
26 *
27 * double y, drand();
28 *
29 * drand( &y );
30 *
31 *
32 *
33 * DESCRIPTION:
34 *
35 * Yields a random number 1.0 <= y < 2.0.
36 *
37 * The three-generator congruential algorithm by Brian
38 * Wichmann and David Hill (BYTE magazine, March, 1987,
39 * pp 127-8) is used. The period, given by them, is
40 * 6953607871644.
41 *
42 * Versions invoked by the different arithmetic compile
43 * time options DEC, IBMPC, and MIEEE, produce
44 * approximately the same sequences, differing only in the
45 * least significant bits of the numbers. The UNK option
46 * implements the algorithm as recommended in the BYTE
47 * article. It may be used on all computers. However,
48 * the low order bits of a double precision number may
49 * not be adequately random, and may vary due to arithmetic
50 * implementation details on different computers.
51 *
52 * The other compile options generate an additional random
53 * integer that overwrites the low order bits of the double
54 * precision number. This reduces the period by a factor of
55 * two but tends to overcome the problems mentioned.
56 *
57 */
58
59#include "mconf.h"
60
61
62/* Three-generator random number algorithm
63 * of Brian Wichmann and David Hill
64 * BYTE magazine, March, 1987 pp 127-8
65 *
66 * The period, given by them, is (p-1)(q-1)(r-1)/4 = 6.95e12.
67 */
68
69static int sx = 1;
70static int sy = 10000;
71static int sz = 3000;
72
73static union {
74 double d;
75 unsigned short s[4];
76} unkans;
77
78/* This function implements the three
79 * congruential generators.
80 */
81
82static int ranwh()
83{
84int r, s;
85
86/* sx = sx * 171 mod 30269 */
87r = sx/177;
88s = sx - 177 * r;
89sx = 171 * s - 2 * r;
90if( sx < 0 )
91 sx += 30269;
92
93
94/* sy = sy * 172 mod 30307 */
95r = sy/176;
96s = sy - 176 * r;
97sy = 172 * s - 35 * r;
98if( sy < 0 )
99 sy += 30307;
100
101/* sz = 170 * sz mod 30323 */
102r = sz/178;
103s = sz - 178 * r;
104sz = 170 * s - 63 * r;
105if( sz < 0 )
106 sz += 30323;
107/* The results are in static sx, sy, sz. */
108return 0;
109}
110
111/* drand.c
112 *
113 * Random double precision floating point number between 1 and 2.
114 *
115 * C callable:
116 * drand( &x );
117 */
118
119int drand( a )
120double *a;
121{
122unsigned short r;
123#ifdef DEC
124unsigned short s, t;
125#endif
126
127/* This algorithm of Wichmann and Hill computes a floating point
128 * result:
129 */
130ranwh();
131unkans.d = sx/30269.0 + sy/30307.0 + sz/30323.0;
132r = unkans.d;
133unkans.d -= r;
134unkans.d += 1.0;
135
136/* if UNK option, do nothing further.
137 * Otherwise, make a random 16 bit integer
138 * to overwrite the least significant word
139 * of unkans.
140 */
141#ifdef UNK
142/* do nothing */
143#else
144ranwh();
145r = sx * sy + sz;
146#endif
147
148#ifdef DEC
149/* To make the numbers as similar as possible
150 * in all arithmetics, the random integer has
151 * to be inserted 3 bits higher up in a DEC number.
152 * An alternative would be put it 3 bits lower down
153 * in all the other number types.
154 */
155s = unkans.s[2];
156t = s & 07; /* save these bits to put in at the bottom */
157s &= 0177770;
158s |= (r >> 13) & 07;
159unkans.s[2] = s;
160t |= r << 3;
161unkans.s[3] = t;
162#endif
163
164#ifdef IBMPC
165unkans.s[0] = r;
166#endif
167
168#ifdef MIEEE
169unkans.s[3] = r;
170#endif
171
172*a = unkans.d;
173return 0;
174}
diff --git a/src/regress/lib/libc/cephes/econst.c b/src/regress/lib/libc/cephes/econst.c
new file mode 100644
index 0000000000..4232059e4c
--- /dev/null
+++ b/src/regress/lib/libc/cephes/econst.c
@@ -0,0 +1,114 @@
1/* $OpenBSD: econst.c,v 1.1 2011/07/02 18:11:01 martynas Exp $ */
2
3/*
4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19/* econst.c */
20/* e type constants used by high precision check routines */
21
22#include "ehead.h"
23
24
25#if NE == 10
26/* 0.0 */
27unsigned short ezero[NE] =
28 {0x0000, 0x0000, 0x0000, 0x0000,
29 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
30
31/* 5.0E-1 */
32unsigned short ehalf[NE] =
33 {0x0000, 0x0000, 0x0000, 0x0000,
34 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
35
36/* 1.0E0 */
37unsigned short eone[NE] =
38 {0x0000, 0x0000, 0x0000, 0x0000,
39 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
40
41/* 2.0E0 */
42unsigned short etwo[NE] =
43 {0x0000, 0x0000, 0x0000, 0x0000,
44 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
45
46/* 3.2E1 */
47unsigned short e32[NE] =
48 {0x0000, 0x0000, 0x0000, 0x0000,
49 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
50
51/* 6.93147180559945309417232121458176568075500134360255E-1 */
52unsigned short elog2[NE] =
53 {0x40f3, 0xf6af, 0x03f2, 0xb398,
54 0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
55
56/* 1.41421356237309504880168872420969807856967187537695E0 */
57unsigned short esqrt2[NE] =
58 {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
59 0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
60
61/* 3.14159265358979323846264338327950288419716939937511E0 */
62unsigned short epi[NE] =
63 {0x2902, 0x1cd1, 0x80dc, 0x628b,
64 0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
65
66/* 5.7721566490153286060651209008240243104215933593992E-1 */
67unsigned short eeul[NE] = {
680xd1be,0xc7a4,0076660,0063743,0111704,0x3ffe,};
69
70#else
71
72/* 0.0 */
73unsigned short ezero[NE] = {
740, 0000000,0000000,0000000,0000000,0000000,};
75/* 5.0E-1 */
76unsigned short ehalf[NE] = {
770, 0000000,0000000,0000000,0100000,0x3ffe,};
78/* 1.0E0 */
79unsigned short eone[NE] = {
800, 0000000,0000000,0000000,0100000,0x3fff,};
81/* 2.0E0 */
82unsigned short etwo[NE] = {
830, 0000000,0000000,0000000,0100000,0040000,};
84/* 3.2E1 */
85unsigned short e32[NE] = {
860, 0000000,0000000,0000000,0100000,0040004,};
87/* 6.93147180559945309417232121458176568075500134360255E-1 */
88unsigned short elog2[NE] = {
890xc9e4,0x79ab,0150717,0013767,0130562,0x3ffe,};
90/* 1.41421356237309504880168872420969807856967187537695E0 */
91unsigned short esqrt2[NE] = {
920x597e,0x6484,0174736,0171463,0132404,0x3fff,};
93/* 2/sqrt(PI) =
94 * 1.12837916709551257389615890312154517168810125865800E0 */
95unsigned short eoneopi[NE] = {
960x71d5,0x688d,0012333,0135202,0110156,0x3fff,};
97/* 3.14159265358979323846264338327950288419716939937511E0 */
98unsigned short epi[NE] = {
990xc4c6,0xc234,0020550,0155242,0144417,0040000,};
100/* 5.7721566490153286060651209008240243104215933593992E-1 */
101unsigned short eeul[NE] = {
1020xd1be,0xc7a4,0076660,0063743,0111704,0x3ffe,};
103#endif
104extern unsigned short ezero[];
105extern unsigned short ehalf[];
106extern unsigned short eone[];
107extern unsigned short etwo[];
108extern unsigned short e32[];
109extern unsigned short elog2[];
110extern unsigned short esqrt2[];
111extern unsigned short eoneopi[];
112extern unsigned short epi[];
113extern unsigned short eeul[];
114
diff --git a/src/regress/lib/libc/cephes/eexp.c b/src/regress/lib/libc/cephes/eexp.c
new file mode 100644
index 0000000000..74f0d6adb3
--- /dev/null
+++ b/src/regress/lib/libc/cephes/eexp.c
@@ -0,0 +1,86 @@
1/* $OpenBSD: eexp.c,v 1.1 2011/07/02 18:11:01 martynas Exp $ */
2
3/*
4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19/* xexp.c */
20/* exponential function check routine */
21/* by Stephen L. Moshier. */
22
23
24#include "ehead.h"
25
26void eexp( x, y )
27unsigned short *x, *y;
28{
29unsigned short num[NE], den[NE], x2[NE];
30long i;
31unsigned short sign, expchk;
32
33/* range reduction theory: x = i + f, 0<=f<1;
34 * e**x = e**i * e**f
35 * e**i = 2**(i/log 2).
36 * Let i/log2 = i1 + f1, 0<=f1<1.
37 * Then e**i = 2**i1 * 2**f1, so
38 * e**x = 2**i1 * e**(log 2 * f1) * e**f.
39 */
40if( ecmp(x, ezero) == 0 )
41 {
42 emov( eone, y );
43 return;
44 }
45emov(x, x2);
46expchk = x2[NE-1];
47sign = expchk & 0x8000;
48x2[NE-1] &= 0x7fff;
49
50/* Test for excessively large argument */
51expchk &= 0x7fff;
52if( expchk > (EXONE + 15) )
53 {
54 eclear( y );
55 if( sign == 0 )
56 einfin( y );
57 return;
58 }
59
60eifrac( x2, &i, num ); /* x = i + f */
61
62if( i != 0 )
63 {
64 ltoe( &i, den ); /* floating point i */
65 ediv( elog2, den, den ); /* i/log 2 */
66 eifrac( den, &i, den ); /* i/log 2 = i1 + f1 */
67 emul( elog2, den, den ); /* log 2 * f1 */
68 eadd( den, num, x2 ); /* log 2 * f1 + f */
69 }
70
71/*x2[NE-1] -= 1;*/
72eldexp( x2, -1L, x2 ); /* divide by 2 */
73etanh( x2, x2 ); /* tanh( x/2 ) */
74eadd( x2, eone, num ); /* 1 + tanh */
75eneg( x2 );
76eadd( x2, eone, den ); /* 1 - tanh */
77ediv( den, num, y ); /* (1 + tanh)/(1 - tanh) */
78
79/*y[NE-1] += i;*/
80if( sign )
81 {
82 ediv( y, eone, y );
83 i = -i;
84 }
85eldexp( y, i, y ); /* multiply by 2**i */
86}
diff --git a/src/regress/lib/libc/cephes/ehead.h b/src/regress/lib/libc/cephes/ehead.h
new file mode 100644
index 0000000000..009bcf89cc
--- /dev/null
+++ b/src/regress/lib/libc/cephes/ehead.h
@@ -0,0 +1,59 @@
1/* $OpenBSD: ehead.h,v 1.1 2011/07/02 18:11:01 martynas Exp $ */
2
3/*
4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19/* Include file for extended precision arithmetic programs.
20 */
21
22/* Number of 16 bit words in external x type format */
23#define NE 10
24
25/* Number of 16 bit words in internal format */
26#define NI (NE+3)
27
28/* Array offset to exponent */
29#define E 1
30
31/* Array offset to high guard word */
32#define M 2
33
34/* Number of bits of precision */
35#define NBITS ((NI-4)*16)
36
37/* Maximum number of decimal digits in ASCII conversion
38 * = NBITS*log10(2)
39 */
40#define NDEC (NBITS*8/27)
41
42/* The exponent of 1.0 */
43#define EXONE (0x3fff)
44
45void eadd(), esub(), emul(), ediv();
46int ecmp(), enormlz(), eshift();
47void eshup1(), eshup8(), eshup6(), eshdn1(), eshdn8(), eshdn6();
48void eabs(), eneg(), emov(), eclear(), einfin(), efloor();
49void eldexp(), efrexp(), eifrac(), ltoe();
50void esqrt(), elog(), eexp(), etanh(), epow();
51void asctoe(), asctoe24(), asctoe53(), asctoe64();
52void etoasc(), e24toasc(), e53toasc(), e64toasc();
53void etoe64(), etoe53(), etoe24(), e64toe(), e53toe(), e24toe();
54int mtherr();
55extern unsigned short ezero[], ehalf[], eone[], etwo[];
56extern unsigned short elog2[], esqrt2[];
57
58
59/* by Stephen L. Moshier. */
diff --git a/src/regress/lib/libc/cephes/elog.c b/src/regress/lib/libc/cephes/elog.c
new file mode 100644
index 0000000000..079cc754f4
--- /dev/null
+++ b/src/regress/lib/libc/cephes/elog.c
@@ -0,0 +1,110 @@
1/* $OpenBSD: elog.c,v 1.1 2011/07/02 18:11:01 martynas Exp $ */
2
3/*
4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19/* xlog.c */
20/* natural logarithm */
21/* by Stephen L. Moshier. */
22
23#include "mconf.h"
24#include "ehead.h"
25
26
27
28void elog( x, y )
29unsigned short *x, *y;
30{
31unsigned short xx[NE], z[NE], a[NE], b[NE], t[NE], qj[NE];
32long ex;
33int fex;
34
35
36if( x[NE-1] & (unsigned short )0x8000 )
37 {
38 eclear(y);
39 mtherr( "elog", DOMAIN );
40 return;
41 }
42if( ecmp( x, ezero ) == 0 )
43 {
44 einfin( y );
45 eneg(y);
46 mtherr( "elog", SING );
47 return;
48 }
49if( ecmp( x, eone ) == 0 )
50 {
51 eclear( y );
52 return;
53 }
54
55/* range reduction: log x = log( 2**ex * m ) = ex * log2 + log m */
56efrexp( x, &fex, xx );
57/*
58emov(x, xx );
59ex = xx[NX-1] & 0x7fff;
60ex -= 0x3ffe;
61xx[NX-1] = 0x3ffe;
62*/
63
64/* Adjust range to 1/sqrt(2), sqrt(2) */
65esqrt2[NE-1] -= 1;
66if( ecmp( xx, esqrt2 ) < 0 )
67 {
68 fex -= 1;
69 emul( xx, etwo, xx );
70 }
71esqrt2[NE-1] += 1;
72
73esub( eone, xx, a );
74if( a[NE-1] == 0 )
75 {
76 eclear( y );
77 goto logdon;
78 }
79eadd( eone, xx, b );
80ediv( b, a, y ); /* store (x-1)/(x+1) in y */
81
82emul( y, y, z );
83
84emov( eone, a );
85emov( eone, b );
86emov( eone, qj );
87do
88 {
89 eadd( etwo, qj, qj ); /* 2 * i + 1 */
90 emul( z, a, a );
91 ediv( qj, a, t );
92 eadd( t, b, b );
93 }
94while( ((b[NE-1] & 0x7fff) - (t[NE-1] & 0x7fff)) < NBITS );
95
96
97emul( b, y, y );
98emul( y, etwo, y );
99
100logdon:
101
102/* now add log of 2**ex */
103if( fex != 0 )
104 {
105 ex = fex;
106 ltoe( &ex, b );
107 emul( elog2, b, b );
108 eadd( b, y, y );
109 }
110}
diff --git a/src/regress/lib/libc/cephes/epow.c b/src/regress/lib/libc/cephes/epow.c
new file mode 100644
index 0000000000..646268fce7
--- /dev/null
+++ b/src/regress/lib/libc/cephes/epow.c
@@ -0,0 +1,187 @@
1/* $OpenBSD: epow.c,v 1.1 2011/07/02 18:11:01 martynas Exp $ */
2
3/*
4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19/* epow.c */
20/* power function: z = x**y */
21/* by Stephen L. Moshier. */
22
23
24#include "ehead.h"
25
26extern int rndprc;
27void epowi();
28
29void epow( x, y, z )
30unsigned short *x, *y, *z;
31{
32unsigned short w[NE];
33int rndsav;
34long li;
35
36efloor( y, w );
37if( ecmp(y,w) == 0 )
38 {
39 eifrac( y, &li, w );
40 if( li < 0 )
41 li = -li;
42 if( (li < 0x7fffffff) && (li != 0x80000000) )
43 {
44 epowi( x, y, z );
45 return;
46 }
47 }
48/* z = exp( y * log(x) ) */
49rndsav = rndprc;
50rndprc = NBITS;
51elog( x, w );
52emul( y, w, w );
53eexp( w, z );
54rndprc = rndsav;
55emul( eone, z, z );
56}
57
58
59/* y is integer valued. */
60
61void epowi( x, y, z )
62unsigned short x[], y[], z[];
63{
64unsigned short w[NE];
65long li, lx;
66unsigned long lu;
67int rndsav;
68unsigned short signx;
69/* unsigned short signy; */
70
71rndsav = rndprc;
72eifrac( y, &li, w );
73if( li < 0 )
74 lx = -li;
75else
76 lx = li;
77
78if( (lx == 0x7fffffff) || (lx == 0x80000000) )
79 {
80 epow( x, y, z );
81 goto done;
82 }
83
84if( (x[NE-1] & (unsigned short )0x7fff) == 0 )
85 {
86 if( li == 0 )
87 {
88 emov( eone, z );
89 return;
90 }
91 else if( li < 0 )
92 {
93 einfin( z );
94 return;
95 }
96 else
97 {
98 eclear( z );
99 return;
100 }
101 }
102
103if( li == 0L )
104 {
105 emov( eone, z );
106 return;
107 }
108
109emov( x, w );
110signx = w[NE-1] & (unsigned short )0x8000;
111w[NE-1] &= (unsigned short )0x7fff;
112
113/* Overflow detection */
114/*
115lx = li * (w[NE-1] - 0x3fff);
116if( lx > 16385L )
117 {
118 einfin( z );
119 mtherr( "epowi", OVERFLOW );
120 goto done;
121 }
122if( lx < -16450L )
123 {
124 eclear( z );
125 return;
126 }
127*/
128rndprc = NBITS;
129
130if( li < 0 )
131 {
132 lu = (unsigned int )( -li );
133/* signy = 0xffff;*/
134 ediv( w, eone, w );
135 }
136else
137 {
138 lu = (unsigned int )li;
139/* signy = 0;*/
140 }
141
142/* First bit of the power */
143if( lu & 1 )
144 {
145 emov( w, z );
146 }
147else
148 {
149 emov( eone, z );
150 signx = 0;
151 }
152
153
154lu >>= 1;
155while( lu != 0L )
156 {
157 emul( w, w, w ); /* arg to the 2-to-the-kth power */
158 if( lu & 1L ) /* if that bit is set, then include in product */
159 emul( w, z, z );
160 lu >>= 1;
161 }
162
163
164done:
165
166if( signx )
167 eneg( z ); /* odd power of negative number */
168
169/*
170if( signy )
171 {
172 if( ecmp( z, ezero ) != 0 )
173 {
174 ediv( z, eone, z );
175 }
176 else
177 {
178 einfin( z );
179 printf( "epowi OVERFLOW\n" );
180 }
181 }
182*/
183rndprc = rndsav;
184emul( eone, z, z );
185}
186
187
diff --git a/src/regress/lib/libc/cephes/etanh.c b/src/regress/lib/libc/cephes/etanh.c
new file mode 100644
index 0000000000..4ac5ff1c21
--- /dev/null
+++ b/src/regress/lib/libc/cephes/etanh.c
@@ -0,0 +1,70 @@
1/* $OpenBSD: etanh.c,v 1.1 2011/07/02 18:11:01 martynas Exp $ */
2
3/*
4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19/* xtanh.c */
20/* hyperbolic tangent check routine */
21/* this subroutine is used by the exponential function routine */
22/* by Stephen L. Moshier. */
23
24
25
26#include "ehead.h"
27
28
29void etanh( x, y )
30unsigned short *x, *y;
31{
32unsigned short e[NE], r[NE], j[NE], xx[NE], m2[NE];
33short i, n;
34long lj;
35
36emov( x, r );
37r[NE-1] &= (unsigned short )0x7fff;
38if( ecmp(r, eone) >= 0 )
39 {
40/* tanh(x) = (exp(x) - exp(-x)) / (exp(x) + exp(-x))
41 * Note eexp() calls xtanh, but with an argument less than (1 + log 2)/2.
42 */
43 eexp( r, e );
44 ediv( e, eone, r );
45 esub( r, e, xx );
46 eadd( r, e, j );
47 ediv( j, xx, y );
48 return;
49 }
50
51emov( etwo, m2 );
52eneg( m2 );
53
54n = NBITS/8; /* Number of terms to do in the continued fraction */
55lj = 2 * n + 1;
56ltoe( &lj, j );
57
58emov( j, e );
59emul( x, x, xx );
60
61/* continued fraction */
62for( i=0; i<n; i++)
63 {
64 ediv( e, xx, r );
65 eadd( m2, j, j );
66 eadd( r, j, e );
67 }
68
69ediv( e, x, y );
70}
diff --git a/src/regress/lib/libc/cephes/etodec.c b/src/regress/lib/libc/cephes/etodec.c
new file mode 100644
index 0000000000..a15845efb6
--- /dev/null
+++ b/src/regress/lib/libc/cephes/etodec.c
@@ -0,0 +1,199 @@
1/* $OpenBSD: etodec.c,v 1.1 2011/07/02 18:11:01 martynas Exp $ */
2
3/*
4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19#include "ehead.h"
20void emovi(), emovo(), ecleaz(), eshdn8(), emdnorm();
21void todec();
22/*
23; convert DEC double precision to e type
24; double d;
25; short e[NE];
26; dectoe( &d, e );
27*/
28void dectoe( d, e )
29unsigned short *d;
30unsigned short *e;
31{
32unsigned short y[NI];
33register unsigned short r, *p;
34
35ecleaz(y); /* start with a zero */
36p = y; /* point to our number */
37r = *d; /* get DEC exponent word */
38if( *d & (unsigned int )0x8000 )
39 *p = 0xffff; /* fill in our sign */
40++p; /* bump pointer to our exponent word */
41r &= 0x7fff; /* strip the sign bit */
42if( r == 0 ) /* answer = 0 if high order DEC word = 0 */
43 goto done;
44
45
46r >>= 7; /* shift exponent word down 7 bits */
47r += EXONE - 0201; /* subtract DEC exponent offset */
48 /* add our e type exponent offset */
49*p++ = r; /* to form our exponent */
50
51r = *d++; /* now do the high order mantissa */
52r &= 0177; /* strip off the DEC exponent and sign bits */
53r |= 0200; /* the DEC understood high order mantissa bit */
54*p++ = r; /* put result in our high guard word */
55
56*p++ = *d++; /* fill in the rest of our mantissa */
57*p++ = *d++;
58*p = *d;
59
60eshdn8(y); /* shift our mantissa down 8 bits */
61done:
62emovo( y, e );
63}
64
65
66
67/*
68; convert e type to DEC double precision
69; double d;
70; short e[NE];
71; etodec( e, &d );
72*/
73#if 0
74static unsigned short decbit[NI] = {0,0,0,0,0,0,0200,0};
75void etodec( x, d )
76unsigned short *x, *d;
77{
78unsigned short xi[NI];
79register unsigned short r;
80int i, j;
81
82emovi( x, xi );
83*d = 0;
84if( xi[0] != 0 )
85 *d = 0100000;
86r = xi[E];
87if( r < (EXONE - 128) )
88 goto zout;
89i = xi[M+4];
90if( (i & 0200) != 0 )
91 {
92 if( (i & 0377) == 0200 )
93 {
94 if( (i & 0400) != 0 )
95 {
96 /* check all less significant bits */
97 for( j=M+5; j<NI; j++ )
98 {
99 if( xi[j] != 0 )
100 goto yesrnd;
101 }
102 }
103 goto nornd;
104 }
105yesrnd:
106 eaddm( decbit, xi );
107 r -= enormlz(xi);
108 }
109
110nornd:
111
112r -= EXONE;
113r += 0201;
114if( r < 0 )
115 {
116zout:
117 *d++ = 0;
118 *d++ = 0;
119 *d++ = 0;
120 *d++ = 0;
121 return;
122 }
123if( r >= 0377 )
124 {
125 *d++ = 077777;
126 *d++ = -1;
127 *d++ = -1;
128 *d++ = -1;
129 return;
130 }
131r &= 0377;
132r <<= 7;
133eshup8( xi );
134xi[M] &= 0177;
135r |= xi[M];
136*d++ |= r;
137*d++ = xi[M+1];
138*d++ = xi[M+2];
139*d++ = xi[M+3];
140}
141#else
142
143extern int rndprc;
144
145void etodec( x, d )
146unsigned short *x, *d;
147{
148unsigned short xi[NI];
149long exp;
150int rndsav;
151
152emovi( x, xi );
153exp = (long )xi[E] - (EXONE - 0201); /* adjust exponent for offsets */
154/* round off to nearest or even */
155rndsav = rndprc;
156rndprc = 56;
157emdnorm( xi, 0, 0, exp, 64 );
158rndprc = rndsav;
159todec( xi, d );
160}
161
162void todec( x, y )
163unsigned short *x, *y;
164{
165unsigned short i;
166unsigned short *p;
167
168p = x;
169*y = 0;
170if( *p++ )
171 *y = 0100000;
172i = *p++;
173if( i == 0 )
174 {
175 *y++ = 0;
176 *y++ = 0;
177 *y++ = 0;
178 *y++ = 0;
179 return;
180 }
181if( i > 0377 )
182 {
183 *y++ |= 077777;
184 *y++ = 0xffff;
185 *y++ = 0xffff;
186 *y++ = 0xffff;
187 return;
188 }
189i &= 0377;
190i <<= 7;
191eshup8( x );
192x[M] &= 0177;
193i |= x[M];
194*y++ |= i;
195*y++ = x[M+1];
196*y++ = x[M+2];
197*y++ = x[M+3];
198}
199#endif
diff --git a/src/regress/lib/libc/cephes/ieee.c b/src/regress/lib/libc/cephes/ieee.c
new file mode 100644
index 0000000000..e2b8aa7b99
--- /dev/null
+++ b/src/regress/lib/libc/cephes/ieee.c
@@ -0,0 +1,4153 @@
1/* $OpenBSD: ieee.c,v 1.1 2011/07/02 18:11:01 martynas Exp $ */
2
3/*
4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19/* ieee.c
20 *
21 * Extended precision IEEE binary floating point arithmetic routines
22 *
23 * Numbers are stored in C language as arrays of 16-bit unsigned
24 * short integers. The arguments of the routines are pointers to
25 * the arrays.
26 *
27 *
28 * External e type data structure, simulates Intel 8087 chip
29 * temporary real format but possibly with a larger significand:
30 *
31 * NE-1 significand words (least significant word first,
32 * most significant bit is normally set)
33 * exponent (value = EXONE for 1.0,
34 * top bit is the sign)
35 *
36 *
37 * Internal data structure of a number (a "word" is 16 bits):
38 *
39 * ei[0] sign word (0 for positive, 0xffff for negative)
40 * ei[1] biased exponent (value = EXONE for the number 1.0)
41 * ei[2] high guard word (always zero after normalization)
42 * ei[3]
43 * to ei[NI-2] significand (NI-4 significand words,
44 * most significant word first,
45 * most significant bit is set)
46 * ei[NI-1] low guard word (0x8000 bit is rounding place)
47 *
48 *
49 *
50 * Routines for external format numbers
51 *
52 * asctoe( string, e ) ASCII string to extended double e type
53 * asctoe64( string, &d ) ASCII string to long double
54 * asctoe53( string, &d ) ASCII string to double
55 * asctoe24( string, &f ) ASCII string to single
56 * asctoeg( string, e, prec ) ASCII string to specified precision
57 * e24toe( &f, e ) IEEE single precision to e type
58 * e53toe( &d, e ) IEEE double precision to e type
59 * e64toe( &d, e ) IEEE long double precision to e type
60 * eabs(e) absolute value
61 * eadd( a, b, c ) c = b + a
62 * eclear(e) e = 0
63 * ecmp (a, b) Returns 1 if a > b, 0 if a == b,
64 * -1 if a < b, -2 if either a or b is a NaN.
65 * ediv( a, b, c ) c = b / a
66 * efloor( a, b ) truncate to integer, toward -infinity
67 * efrexp( a, exp, s ) extract exponent and significand
68 * eifrac( e, &l, frac ) e to long integer and e type fraction
69 * euifrac( e, &l, frac ) e to unsigned long integer and e type fraction
70 * einfin( e ) set e to infinity, leaving its sign alone
71 * eldexp( a, n, b ) multiply by 2**n
72 * emov( a, b ) b = a
73 * emul( a, b, c ) c = b * a
74 * eneg(e) e = -e
75 * eround( a, b ) b = nearest integer value to a
76 * esub( a, b, c ) c = b - a
77 * e24toasc( &f, str, n ) single to ASCII string, n digits after decimal
78 * e53toasc( &d, str, n ) double to ASCII string, n digits after decimal
79 * e64toasc( &d, str, n ) long double to ASCII string
80 * etoasc( e, str, n ) e to ASCII string, n digits after decimal
81 * etoe24( e, &f ) convert e type to IEEE single precision
82 * etoe53( e, &d ) convert e type to IEEE double precision
83 * etoe64( e, &d ) convert e type to IEEE long double precision
84 * ltoe( &l, e ) long (32 bit) integer to e type
85 * ultoe( &l, e ) unsigned long (32 bit) integer to e type
86 * eisneg( e ) 1 if sign bit of e != 0, else 0
87 * eisinf( e ) 1 if e has maximum exponent (non-IEEE)
88 * or is infinite (IEEE)
89 * eisnan( e ) 1 if e is a NaN
90 * esqrt( a, b ) b = square root of a
91 *
92 *
93 * Routines for internal format numbers
94 *
95 * eaddm( ai, bi ) add significands, bi = bi + ai
96 * ecleaz(ei) ei = 0
97 * ecleazs(ei) set ei = 0 but leave its sign alone
98 * ecmpm( ai, bi ) compare significands, return 1, 0, or -1
99 * edivm( ai, bi ) divide significands, bi = bi / ai
100 * emdnorm(ai,l,s,exp) normalize and round off
101 * emovi( a, ai ) convert external a to internal ai
102 * emovo( ai, a ) convert internal ai to external a
103 * emovz( ai, bi ) bi = ai, low guard word of bi = 0
104 * emulm( ai, bi ) multiply significands, bi = bi * ai
105 * enormlz(ei) left-justify the significand
106 * eshdn1( ai ) shift significand and guards down 1 bit
107 * eshdn8( ai ) shift down 8 bits
108 * eshdn6( ai ) shift down 16 bits
109 * eshift( ai, n ) shift ai n bits up (or down if n < 0)
110 * eshup1( ai ) shift significand and guards up 1 bit
111 * eshup8( ai ) shift up 8 bits
112 * eshup6( ai ) shift up 16 bits
113 * esubm( ai, bi ) subtract significands, bi = bi - ai
114 *
115 *
116 * The result is always normalized and rounded to NI-4 word precision
117 * after each arithmetic operation.
118 *
119 * Exception flags are NOT fully supported.
120 *
121 * Define INFINITY in mconf.h for support of infinity; otherwise a
122 * saturation arithmetic is implemented.
123 *
124 * Define NANS for support of Not-a-Number items; otherwise the
125 * arithmetic will never produce a NaN output, and might be confused
126 * by a NaN input.
127 * If NaN's are supported, the output of ecmp(a,b) is -2 if
128 * either a or b is a NaN. This means asking if(ecmp(a,b) < 0)
129 * may not be legitimate. Use if(ecmp(a,b) == -1) for less-than
130 * if in doubt.
131 * Signaling NaN's are NOT supported; they are treated the same
132 * as quiet NaN's.
133 *
134 * Denormals are always supported here where appropriate (e.g., not
135 * for conversion to DEC numbers).
136 */
137
138/*
139 * Revision history:
140 *
141 * 5 Jan 84 PDP-11 assembly language version
142 * 2 Mar 86 fixed bug in asctoq()
143 * 6 Dec 86 C language version
144 * 30 Aug 88 100 digit version, improved rounding
145 * 15 May 92 80-bit long double support
146 *
147 * Author: S. L. Moshier.
148 */
149
150#include <stdio.h>
151#include "mconf.h"
152#include "ehead.h"
153
154/* Change UNK into something else. */
155#ifdef UNK
156#undef UNK
157#if BIGENDIAN
158#define MIEEE 1
159#else
160#define IBMPC 1
161#endif
162#endif
163
164/* NaN's require infinity support. */
165#ifdef NANS
166#ifndef INFINITY
167#define INFINITY
168#endif
169#endif
170
171/* This handles 64-bit long ints. */
172#define LONGBITS (8 * sizeof(long))
173
174/* Control register for rounding precision.
175 * This can be set to 80 (if NE=6), 64, 56, 53, or 24 bits.
176 */
177int rndprc = NBITS;
178extern int rndprc;
179
180void eaddm(), esubm(), emdnorm(), asctoeg(), enan();
181static void toe24(), toe53(), toe64(), toe113();
182void eremain(), einit(), eiremain();
183int ecmpm(), edivm(), emulm(), eisneg(), eisinf();
184void emovi(), emovo(), emovz(), ecleaz(), eadd1();
185void etodec(), todec(), dectoe();
186int eisnan(), eiisnan();
187
188
189
190void einit()
191{
192}
193
194/*
195; Clear out entire external format number.
196;
197; unsigned short x[];
198; eclear( x );
199*/
200
201void eclear( x )
202register unsigned short *x;
203{
204register int i;
205
206for( i=0; i<NE; i++ )
207 *x++ = 0;
208}
209
210
211
212/* Move external format number from a to b.
213 *
214 * emov( a, b );
215 */
216
217void emov( a, b )
218register unsigned short *a, *b;
219{
220register int i;
221
222for( i=0; i<NE; i++ )
223 *b++ = *a++;
224}
225
226
227/*
228; Absolute value of external format number
229;
230; short x[NE];
231; eabs( x );
232*/
233
234void eabs(x)
235unsigned short x[]; /* x is the memory address of a short */
236{
237
238x[NE-1] &= 0x7fff; /* sign is top bit of last word of external format */
239}
240
241
242
243
244/*
245; Negate external format number
246;
247; unsigned short x[NE];
248; eneg( x );
249*/
250
251void eneg(x)
252unsigned short x[];
253{
254
255#ifdef NANS
256if( eisnan(x) )
257 return;
258#endif
259x[NE-1] ^= 0x8000; /* Toggle the sign bit */
260}
261
262
263
264/* Return 1 if external format number is negative,
265 * else return zero.
266 */
267int eisneg(x)
268unsigned short x[];
269{
270
271#ifdef NANS
272if( eisnan(x) )
273 return( 0 );
274#endif
275if( x[NE-1] & 0x8000 )
276 return( 1 );
277else
278 return( 0 );
279}
280
281
282/* Return 1 if external format number has maximum possible exponent,
283 * else return zero.
284 */
285int eisinf(x)
286unsigned short x[];
287{
288
289if( (x[NE-1] & 0x7fff) == 0x7fff )
290 {
291#ifdef NANS
292 if( eisnan(x) )
293 return( 0 );
294#endif
295 return( 1 );
296 }
297else
298 return( 0 );
299}
300
301/* Check if e-type number is not a number.
302 */
303int eisnan(x)
304unsigned short x[];
305{
306
307#ifdef NANS
308int i;
309/* NaN has maximum exponent */
310if( (x[NE-1] & 0x7fff) != 0x7fff )
311 return (0);
312/* ... and non-zero significand field. */
313for( i=0; i<NE-1; i++ )
314 {
315 if( *x++ != 0 )
316 return (1);
317 }
318#endif
319return (0);
320}
321
322/*
323; Fill entire number, including exponent and significand, with
324; largest possible number. These programs implement a saturation
325; value that is an ordinary, legal number. A special value
326; "infinity" may also be implemented; this would require tests
327; for that value and implementation of special rules for arithmetic
328; operations involving inifinity.
329*/
330
331void einfin(x)
332register unsigned short *x;
333{
334register int i;
335
336#ifdef INFINITY
337for( i=0; i<NE-1; i++ )
338 *x++ = 0;
339*x |= 32767;
340#else
341for( i=0; i<NE-1; i++ )
342 *x++ = 0xffff;
343*x |= 32766;
344if( rndprc < NBITS )
345 {
346 if (rndprc == 113)
347 {
348 *(x - 9) = 0;
349 *(x - 8) = 0;
350 }
351 if( rndprc == 64 )
352 {
353 *(x-5) = 0;
354 }
355 if( rndprc == 53 )
356 {
357 *(x-4) = 0xf800;
358 }
359 else
360 {
361 *(x-4) = 0;
362 *(x-3) = 0;
363 *(x-2) = 0xff00;
364 }
365 }
366#endif
367}
368
369
370
371/* Move in external format number,
372 * converting it to internal format.
373 */
374void emovi( a, b )
375unsigned short *a, *b;
376{
377register unsigned short *p, *q;
378int i;
379
380q = b;
381p = a + (NE-1); /* point to last word of external number */
382/* get the sign bit */
383if( *p & 0x8000 )
384 *q++ = 0xffff;
385else
386 *q++ = 0;
387/* get the exponent */
388*q = *p--;
389*q++ &= 0x7fff; /* delete the sign bit */
390#ifdef INFINITY
391if( (*(q-1) & 0x7fff) == 0x7fff )
392 {
393#ifdef NANS
394 if( eisnan(a) )
395 {
396 *q++ = 0;
397 for( i=3; i<NI; i++ )
398 *q++ = *p--;
399 return;
400 }
401#endif
402 for( i=2; i<NI; i++ )
403 *q++ = 0;
404 return;
405 }
406#endif
407/* clear high guard word */
408*q++ = 0;
409/* move in the significand */
410for( i=0; i<NE-1; i++ )
411 *q++ = *p--;
412/* clear low guard word */
413*q = 0;
414}
415
416
417/* Move internal format number out,
418 * converting it to external format.
419 */
420void emovo( a, b )
421unsigned short *a, *b;
422{
423register unsigned short *p, *q;
424unsigned short i;
425
426p = a;
427q = b + (NE-1); /* point to output exponent */
428/* combine sign and exponent */
429i = *p++;
430if( i )
431 *q-- = *p++ | 0x8000;
432else
433 *q-- = *p++;
434#ifdef INFINITY
435if( *(p-1) == 0x7fff )
436 {
437#ifdef NANS
438 if( eiisnan(a) )
439 {
440 enan( b, NBITS );
441 return;
442 }
443#endif
444 einfin(b);
445 return;
446 }
447#endif
448/* skip over guard word */
449++p;
450/* move the significand */
451for( i=0; i<NE-1; i++ )
452 *q-- = *p++;
453}
454
455
456
457
458/* Clear out internal format number.
459 */
460
461void ecleaz( xi )
462register unsigned short *xi;
463{
464register int i;
465
466for( i=0; i<NI; i++ )
467 *xi++ = 0;
468}
469
470/* same, but don't touch the sign. */
471
472void ecleazs( xi )
473register unsigned short *xi;
474{
475register int i;
476
477++xi;
478for(i=0; i<NI-1; i++)
479 *xi++ = 0;
480}
481
482
483
484
485/* Move internal format number from a to b.
486 */
487void emovz( a, b )
488register unsigned short *a, *b;
489{
490register int i;
491
492for( i=0; i<NI-1; i++ )
493 *b++ = *a++;
494/* clear low guard word */
495*b = 0;
496}
497
498/* Return nonzero if internal format number is a NaN.
499 */
500
501int eiisnan (x)
502unsigned short x[];
503{
504int i;
505
506if( (x[E] & 0x7fff) == 0x7fff )
507 {
508 for( i=M+1; i<NI; i++ )
509 {
510 if( x[i] != 0 )
511 return(1);
512 }
513 }
514return(0);
515}
516
517#ifdef INFINITY
518/* Return nonzero if internal format number is infinite. */
519
520static int
521eiisinf (x)
522 unsigned short x[];
523{
524
525#ifdef NANS
526 if (eiisnan (x))
527 return (0);
528#endif
529 if ((x[E] & 0x7fff) == 0x7fff)
530 return (1);
531 return (0);
532}
533#endif
534
535/*
536; Compare significands of numbers in internal format.
537; Guard words are included in the comparison.
538;
539; unsigned short a[NI], b[NI];
540; cmpm( a, b );
541;
542; for the significands:
543; returns +1 if a > b
544; 0 if a == b
545; -1 if a < b
546*/
547int ecmpm( a, b )
548register unsigned short *a, *b;
549{
550int i;
551
552a += M; /* skip up to significand area */
553b += M;
554for( i=M; i<NI; i++ )
555 {
556 if( *a++ != *b++ )
557 goto difrnt;
558 }
559return(0);
560
561difrnt:
562if( *(--a) > *(--b) )
563 return(1);
564else
565 return(-1);
566}
567
568
569/*
570; Shift significand down by 1 bit
571*/
572
573void eshdn1(x)
574register unsigned short *x;
575{
576register unsigned short bits;
577int i;
578
579x += M; /* point to significand area */
580
581bits = 0;
582for( i=M; i<NI; i++ )
583 {
584 if( *x & 1 )
585 bits |= 1;
586 *x >>= 1;
587 if( bits & 2 )
588 *x |= 0x8000;
589 bits <<= 1;
590 ++x;
591 }
592}
593
594
595
596/*
597; Shift significand up by 1 bit
598*/
599
600void eshup1(x)
601register unsigned short *x;
602{
603register unsigned short bits;
604int i;
605
606x += NI-1;
607bits = 0;
608
609for( i=M; i<NI; i++ )
610 {
611 if( *x & 0x8000 )
612 bits |= 1;
613 *x <<= 1;
614 if( bits & 2 )
615 *x |= 1;
616 bits <<= 1;
617 --x;
618 }
619}
620
621
622
623/*
624; Shift significand down by 8 bits
625*/
626
627void eshdn8(x)
628register unsigned short *x;
629{
630register unsigned short newbyt, oldbyt;
631int i;
632
633x += M;
634oldbyt = 0;
635for( i=M; i<NI; i++ )
636 {
637 newbyt = *x << 8;
638 *x >>= 8;
639 *x |= oldbyt;
640 oldbyt = newbyt;
641 ++x;
642 }
643}
644
645/*
646; Shift significand up by 8 bits
647*/
648
649void eshup8(x)
650register unsigned short *x;
651{
652int i;
653register unsigned short newbyt, oldbyt;
654
655x += NI-1;
656oldbyt = 0;
657
658for( i=M; i<NI; i++ )
659 {
660 newbyt = *x >> 8;
661 *x <<= 8;
662 *x |= oldbyt;
663 oldbyt = newbyt;
664 --x;
665 }
666}
667
668/*
669; Shift significand up by 16 bits
670*/
671
672void eshup6(x)
673register unsigned short *x;
674{
675int i;
676register unsigned short *p;
677
678p = x + M;
679x += M + 1;
680
681for( i=M; i<NI-1; i++ )
682 *p++ = *x++;
683
684*p = 0;
685}
686
687/*
688; Shift significand down by 16 bits
689*/
690
691void eshdn6(x)
692register unsigned short *x;
693{
694int i;
695register unsigned short *p;
696
697x += NI-1;
698p = x + 1;
699
700for( i=M; i<NI-1; i++ )
701 *(--p) = *(--x);
702
703*(--p) = 0;
704}
705
706/*
707; Add significands
708; x + y replaces y
709*/
710
711void eaddm( x, y )
712unsigned short *x, *y;
713{
714register unsigned long a;
715int i;
716unsigned int carry;
717
718x += NI-1;
719y += NI-1;
720carry = 0;
721for( i=M; i<NI; i++ )
722 {
723 a = (unsigned long )(*x) + (unsigned long )(*y) + carry;
724 if( a & 0x10000 )
725 carry = 1;
726 else
727 carry = 0;
728 *y = (unsigned short )a;
729 --x;
730 --y;
731 }
732}
733
734/*
735; Subtract significands
736; y - x replaces y
737*/
738
739void esubm( x, y )
740unsigned short *x, *y;
741{
742unsigned long a;
743int i;
744unsigned int carry;
745
746x += NI-1;
747y += NI-1;
748carry = 0;
749for( i=M; i<NI; i++ )
750 {
751 a = (unsigned long )(*y) - (unsigned long )(*x) - carry;
752 if( a & 0x10000 )
753 carry = 1;
754 else
755 carry = 0;
756 *y = (unsigned short )a;
757 --x;
758 --y;
759 }
760}
761
762
763/* Divide significands */
764
765static unsigned short equot[NI] = {0}; /* was static */
766
767#if 0
768int edivm( den, num )
769unsigned short den[], num[];
770{
771int i;
772register unsigned short *p, *q;
773unsigned short j;
774
775p = &equot[0];
776*p++ = num[0];
777*p++ = num[1];
778
779for( i=M; i<NI; i++ )
780 {
781 *p++ = 0;
782 }
783
784/* Use faster compare and subtraction if denominator
785 * has only 15 bits of significance.
786 */
787p = &den[M+2];
788if( *p++ == 0 )
789 {
790 for( i=M+3; i<NI; i++ )
791 {
792 if( *p++ != 0 )
793 goto fulldiv;
794 }
795 if( (den[M+1] & 1) != 0 )
796 goto fulldiv;
797 eshdn1(num);
798 eshdn1(den);
799
800 p = &den[M+1];
801 q = &num[M+1];
802
803 for( i=0; i<NBITS+2; i++ )
804 {
805 if( *p <= *q )
806 {
807 *q -= *p;
808 j = 1;
809 }
810 else
811 {
812 j = 0;
813 }
814 eshup1(equot);
815 equot[NI-2] |= j;
816 eshup1(num);
817 }
818 goto divdon;
819 }
820
821/* The number of quotient bits to calculate is
822 * NBITS + 1 scaling guard bit + 1 roundoff bit.
823 */
824fulldiv:
825
826p = &equot[NI-2];
827for( i=0; i<NBITS+2; i++ )
828 {
829 if( ecmpm(den,num) <= 0 )
830 {
831 esubm(den, num);
832 j = 1; /* quotient bit = 1 */
833 }
834 else
835 j = 0;
836 eshup1(equot);
837 *p |= j;
838 eshup1(num);
839 }
840
841divdon:
842
843eshdn1( equot );
844eshdn1( equot );
845
846/* test for nonzero remainder after roundoff bit */
847p = &num[M];
848j = 0;
849for( i=M; i<NI; i++ )
850 {
851 j |= *p++;
852 }
853if( j )
854 j = 1;
855
856
857for( i=0; i<NI; i++ )
858 num[i] = equot[i];
859return( (int )j );
860}
861
862/* Multiply significands */
863int emulm( a, b )
864unsigned short a[], b[];
865{
866unsigned short *p, *q;
867int i, j, k;
868
869equot[0] = b[0];
870equot[1] = b[1];
871for( i=M; i<NI; i++ )
872 equot[i] = 0;
873
874p = &a[NI-2];
875k = NBITS;
876while( *p == 0 ) /* significand is not supposed to be all zero */
877 {
878 eshdn6(a);
879 k -= 16;
880 }
881if( (*p & 0xff) == 0 )
882 {
883 eshdn8(a);
884 k -= 8;
885 }
886
887q = &equot[NI-1];
888j = 0;
889for( i=0; i<k; i++ )
890 {
891 if( *p & 1 )
892 eaddm(b, equot);
893/* remember if there were any nonzero bits shifted out */
894 if( *q & 1 )
895 j |= 1;
896 eshdn1(a);
897 eshdn1(equot);
898 }
899
900for( i=0; i<NI; i++ )
901 b[i] = equot[i];
902
903/* return flag for lost nonzero bits */
904return(j);
905}
906
907#else
908
909/* Multiply significand of e-type number b
910by 16-bit quantity a, e-type result to c. */
911
912void m16m( a, b, c )
913unsigned short a;
914unsigned short b[], c[];
915{
916register unsigned short *pp;
917register unsigned long carry;
918unsigned short *ps;
919unsigned short p[NI];
920unsigned long aa, m;
921int i;
922
923aa = a;
924pp = &p[NI-2];
925*pp++ = 0;
926*pp = 0;
927ps = &b[NI-1];
928
929for( i=M+1; i<NI; i++ )
930 {
931 if( *ps == 0 )
932 {
933 --ps;
934 --pp;
935 *(pp-1) = 0;
936 }
937 else
938 {
939 m = (unsigned long) aa * *ps--;
940 carry = (m & 0xffff) + *pp;
941 *pp-- = (unsigned short )carry;
942 carry = (carry >> 16) + (m >> 16) + *pp;
943 *pp = (unsigned short )carry;
944 *(pp-1) = carry >> 16;
945 }
946 }
947for( i=M; i<NI; i++ )
948 c[i] = p[i];
949}
950
951
952/* Divide significands. Neither the numerator nor the denominator
953is permitted to have its high guard word nonzero. */
954
955
956int edivm( den, num )
957unsigned short den[], num[];
958{
959int i;
960register unsigned short *p;
961unsigned long tnum;
962unsigned short j, tdenm, tquot;
963unsigned short tprod[NI+1];
964
965p = &equot[0];
966*p++ = num[0];
967*p++ = num[1];
968
969for( i=M; i<NI; i++ )
970 {
971 *p++ = 0;
972 }
973eshdn1( num );
974tdenm = den[M+1];
975for( i=M; i<NI; i++ )
976 {
977 /* Find trial quotient digit (the radix is 65536). */
978 tnum = (((unsigned long) num[M]) << 16) + num[M+1];
979
980 /* Do not execute the divide instruction if it will overflow. */
981 if( (tdenm * 0xffffL) < tnum )
982 tquot = 0xffff;
983 else
984 tquot = tnum / tdenm;
985
986 /* Prove that the divide worked. */
987/*
988 tcheck = (unsigned long )tquot * tdenm;
989 if( tnum - tcheck > tdenm )
990 tquot = 0xffff;
991*/
992 /* Multiply denominator by trial quotient digit. */
993 m16m( tquot, den, tprod );
994 /* The quotient digit may have been overestimated. */
995 if( ecmpm( tprod, num ) > 0 )
996 {
997 tquot -= 1;
998 esubm( den, tprod );
999 if( ecmpm( tprod, num ) > 0 )
1000 {
1001 tquot -= 1;
1002 esubm( den, tprod );
1003 }
1004 }
1005/*
1006 if( ecmpm( tprod, num ) > 0 )
1007 {
1008 eshow( "tprod", tprod );
1009 eshow( "num ", num );
1010 printf( "tnum = %08lx, tden = %04x, tquot = %04x\n",
1011 tnum, den[M+1], tquot );
1012 }
1013*/
1014 esubm( tprod, num );
1015/*
1016 if( ecmpm( num, den ) >= 0 )
1017 {
1018 eshow( "num ", num );
1019 eshow( "den ", den );
1020 printf( "tnum = %08lx, tden = %04x, tquot = %04x\n",
1021 tnum, den[M+1], tquot );
1022 }
1023*/
1024 equot[i] = tquot;
1025 eshup6(num);
1026 }
1027/* test for nonzero remainder after roundoff bit */
1028p = &num[M];
1029j = 0;
1030for( i=M; i<NI; i++ )
1031 {
1032 j |= *p++;
1033 }
1034if( j )
1035 j = 1;
1036
1037for( i=0; i<NI; i++ )
1038 num[i] = equot[i];
1039
1040return( (int )j );
1041}
1042
1043
1044
1045/* Multiply significands */
1046int emulm( a, b )
1047unsigned short a[], b[];
1048{
1049unsigned short *p, *q;
1050unsigned short pprod[NI];
1051unsigned short j;
1052int i;
1053
1054equot[0] = b[0];
1055equot[1] = b[1];
1056for( i=M; i<NI; i++ )
1057 equot[i] = 0;
1058
1059j = 0;
1060p = &a[NI-1];
1061q = &equot[NI-1];
1062for( i=M+1; i<NI; i++ )
1063 {
1064 if( *p == 0 )
1065 {
1066 --p;
1067 }
1068 else
1069 {
1070 m16m( *p--, b, pprod );
1071 eaddm(pprod, equot);
1072 }
1073 j |= *q;
1074 eshdn6(equot);
1075 }
1076
1077for( i=0; i<NI; i++ )
1078 b[i] = equot[i];
1079
1080/* return flag for lost nonzero bits */
1081return( (int)j );
1082}
1083
1084
1085/*
1086eshow(str, x)
1087char *str;
1088unsigned short *x;
1089{
1090int i;
1091
1092printf( "%s ", str );
1093for( i=0; i<NI; i++ )
1094 printf( "%04x ", *x++ );
1095printf( "\n" );
1096}
1097*/
1098#endif
1099
1100
1101
1102/*
1103 * Normalize and round off.
1104 *
1105 * The internal format number to be rounded is "s".
1106 * Input "lost" indicates whether the number is exact.
1107 * This is the so-called sticky bit.
1108 *
1109 * Input "subflg" indicates whether the number was obtained
1110 * by a subtraction operation. In that case if lost is nonzero
1111 * then the number is slightly smaller than indicated.
1112 *
1113 * Input "exp" is the biased exponent, which may be negative.
1114 * the exponent field of "s" is ignored but is replaced by
1115 * "exp" as adjusted by normalization and rounding.
1116 *
1117 * Input "rcntrl" is the rounding control.
1118 */
1119
1120static int rlast = -1;
1121static int rw = 0;
1122static unsigned short rmsk = 0;
1123static unsigned short rmbit = 0;
1124static unsigned short rebit = 0;
1125static int re = 0;
1126static unsigned short rbit[NI] = {0,0,0,0,0,0,0,0};
1127
1128void emdnorm( s, lost, subflg, exp, rcntrl )
1129unsigned short s[];
1130int lost;
1131int subflg;
1132long exp;
1133int rcntrl;
1134{
1135int i, j;
1136unsigned short r;
1137
1138/* Normalize */
1139j = enormlz( s );
1140
1141/* a blank significand could mean either zero or infinity. */
1142#ifndef INFINITY
1143if( j > NBITS )
1144 {
1145 ecleazs( s );
1146 return;
1147 }
1148#endif
1149exp -= j;
1150#ifndef INFINITY
1151if( exp >= 32767L )
1152 goto overf;
1153#else
1154if( (j > NBITS) && (exp < 32767L) )
1155 {
1156 ecleazs( s );
1157 return;
1158 }
1159#endif
1160if( exp < 0L )
1161 {
1162 if( exp > (long )(-NBITS-1) )
1163 {
1164 j = (int )exp;
1165 i = eshift( s, j );
1166 if( i )
1167 lost = 1;
1168 }
1169 else
1170 {
1171 ecleazs( s );
1172 return;
1173 }
1174 }
1175/* Round off, unless told not to by rcntrl. */
1176if( rcntrl == 0 )
1177 goto mdfin;
1178/* Set up rounding parameters if the control register changed. */
1179if( rndprc != rlast )
1180 {
1181 ecleaz( rbit );
1182 switch( rndprc )
1183 {
1184 default:
1185 case NBITS:
1186 rw = NI-1; /* low guard word */
1187 rmsk = 0xffff;
1188 rmbit = 0x8000;
1189 rebit = 1;
1190 re = rw - 1;
1191 break;
1192 case 113:
1193 rw = 10;
1194 rmsk = 0x7fff;
1195 rmbit = 0x4000;
1196 rebit = 0x8000;
1197 re = rw;
1198 break;
1199 case 64:
1200 rw = 7;
1201 rmsk = 0xffff;
1202 rmbit = 0x8000;
1203 rebit = 1;
1204 re = rw-1;
1205 break;
1206/* For DEC arithmetic */
1207 case 56:
1208 rw = 6;
1209 rmsk = 0xff;
1210 rmbit = 0x80;
1211 rebit = 0x100;
1212 re = rw;
1213 break;
1214 case 53:
1215 rw = 6;
1216 rmsk = 0x7ff;
1217 rmbit = 0x0400;
1218 rebit = 0x800;
1219 re = rw;
1220 break;
1221 case 24:
1222 rw = 4;
1223 rmsk = 0xff;
1224 rmbit = 0x80;
1225 rebit = 0x100;
1226 re = rw;
1227 break;
1228 }
1229 rbit[re] = rebit;
1230 rlast = rndprc;
1231 }
1232
1233/* Shift down 1 temporarily if the data structure has an implied
1234 * most significant bit and the number is denormal.
1235 * For rndprc = 64 or NBITS, there is no implied bit.
1236 * But Intel long double denormals lose one bit of significance even so.
1237 */
1238#ifdef IBMPC
1239if( (exp <= 0) && (rndprc != NBITS) )
1240#else
1241if( (exp <= 0) && (rndprc != 64) && (rndprc != NBITS) )
1242#endif
1243 {
1244 lost |= s[NI-1] & 1;
1245 eshdn1(s);
1246 }
1247/* Clear out all bits below the rounding bit,
1248 * remembering in r if any were nonzero.
1249 */
1250r = s[rw] & rmsk;
1251if( rndprc < NBITS )
1252 {
1253 i = rw + 1;
1254 while( i < NI )
1255 {
1256 if( s[i] )
1257 r |= 1;
1258 s[i] = 0;
1259 ++i;
1260 }
1261 }
1262s[rw] &= ~rmsk;
1263if( (r & rmbit) != 0 )
1264 {
1265 if( r == rmbit )
1266 {
1267 if( lost == 0 )
1268 { /* round to even */
1269 if( (s[re] & rebit) == 0 )
1270 goto mddone;
1271 }
1272 else
1273 {
1274 if( subflg != 0 )
1275 goto mddone;
1276 }
1277 }
1278 eaddm( rbit, s );
1279 }
1280mddone:
1281#ifdef IBMPC
1282if( (exp <= 0) && (rndprc != NBITS) )
1283#else
1284if( (exp <= 0) && (rndprc != 64) && (rndprc != NBITS) )
1285#endif
1286 {
1287 eshup1(s);
1288 }
1289if( s[2] != 0 )
1290 { /* overflow on roundoff */
1291 eshdn1(s);
1292 exp += 1;
1293 }
1294mdfin:
1295s[NI-1] = 0;
1296if( exp >= 32767L )
1297 {
1298#ifndef INFINITY
1299overf:
1300#endif
1301#ifdef INFINITY
1302 s[1] = 32767;
1303 for( i=2; i<NI-1; i++ )
1304 s[i] = 0;
1305#else
1306 s[1] = 32766;
1307 s[2] = 0;
1308 for( i=M+1; i<NI-1; i++ )
1309 s[i] = 0xffff;
1310 s[NI-1] = 0;
1311 if( (rndprc < 64) || (rndprc == 113) )
1312 {
1313 s[rw] &= ~rmsk;
1314 if( rndprc == 24 )
1315 {
1316 s[5] = 0;
1317 s[6] = 0;
1318 }
1319 }
1320#endif
1321 return;
1322 }
1323if( exp < 0 )
1324 s[1] = 0;
1325else
1326 s[1] = (unsigned short )exp;
1327}
1328
1329
1330
1331/*
1332; Subtract external format numbers.
1333;
1334; unsigned short a[NE], b[NE], c[NE];
1335; esub( a, b, c ); c = b - a
1336*/
1337
1338static int subflg = 0;
1339
1340void esub( a, b, c )
1341unsigned short *a, *b, *c;
1342{
1343
1344#ifdef NANS
1345if( eisnan(a) )
1346 {
1347 emov (a, c);
1348 return;
1349 }
1350if( eisnan(b) )
1351 {
1352 emov(b,c);
1353 return;
1354 }
1355/* Infinity minus infinity is a NaN.
1356 * Test for subtracting infinities of the same sign.
1357 */
1358if( eisinf(a) && eisinf(b) && ((eisneg (a) ^ eisneg (b)) == 0))
1359 {
1360 mtherr( "esub", DOMAIN );
1361 enan( c, NBITS );
1362 return;
1363 }
1364#endif
1365subflg = 1;
1366eadd1( a, b, c );
1367}
1368
1369
1370/*
1371; Add.
1372;
1373; unsigned short a[NE], b[NE], c[NE];
1374; eadd( a, b, c ); c = b + a
1375*/
1376void eadd( a, b, c )
1377unsigned short *a, *b, *c;
1378{
1379
1380#ifdef NANS
1381/* NaN plus anything is a NaN. */
1382if( eisnan(a) )
1383 {
1384 emov(a,c);
1385 return;
1386 }
1387if( eisnan(b) )
1388 {
1389 emov(b,c);
1390 return;
1391 }
1392/* Infinity minus infinity is a NaN.
1393 * Test for adding infinities of opposite signs.
1394 */
1395if( eisinf(a) && eisinf(b)
1396 && ((eisneg(a) ^ eisneg(b)) != 0) )
1397 {
1398 mtherr( "eadd", DOMAIN );
1399 enan( c, NBITS );
1400 return;
1401 }
1402#endif
1403subflg = 0;
1404eadd1( a, b, c );
1405}
1406
1407void eadd1( a, b, c )
1408unsigned short *a, *b, *c;
1409{
1410unsigned short ai[NI], bi[NI], ci[NI];
1411int i, lost, j, k;
1412long lt, lta, ltb;
1413
1414#ifdef INFINITY
1415if( eisinf(a) )
1416 {
1417 emov(a,c);
1418 if( subflg )
1419 eneg(c);
1420 return;
1421 }
1422if( eisinf(b) )
1423 {
1424 emov(b,c);
1425 return;
1426 }
1427#endif
1428emovi( a, ai );
1429emovi( b, bi );
1430if( subflg )
1431 ai[0] = ~ai[0];
1432
1433/* compare exponents */
1434lta = ai[E];
1435ltb = bi[E];
1436lt = lta - ltb;
1437if( lt > 0L )
1438 { /* put the larger number in bi */
1439 emovz( bi, ci );
1440 emovz( ai, bi );
1441 emovz( ci, ai );
1442 ltb = bi[E];
1443 lt = -lt;
1444 }
1445lost = 0;
1446if( lt != 0L )
1447 {
1448 if( lt < (long )(-NBITS-1) )
1449 goto done; /* answer same as larger addend */
1450 k = (int )lt;
1451 lost = eshift( ai, k ); /* shift the smaller number down */
1452 }
1453else
1454 {
1455/* exponents were the same, so must compare significands */
1456 i = ecmpm( ai, bi );
1457 if( i == 0 )
1458 { /* the numbers are identical in magnitude */
1459 /* if different signs, result is zero */
1460 if( ai[0] != bi[0] )
1461 {
1462 eclear(c);
1463 return;
1464 }
1465 /* if same sign, result is double */
1466 /* double denomalized tiny number */
1467 if( (bi[E] == 0) && ((bi[3] & 0x8000) == 0) )
1468 {
1469 eshup1( bi );
1470 goto done;
1471 }
1472 /* add 1 to exponent unless both are zero! */
1473 for( j=1; j<NI-1; j++ )
1474 {
1475 if( bi[j] != 0 )
1476 {
1477 ltb += 1;
1478 if( ltb >= 0x7fff )
1479 {
1480 eclear(c);
1481 einfin(c);
1482 if( ai[0] != 0 )
1483 eneg(c);
1484 return;
1485 }
1486 break;
1487 }
1488 }
1489 bi[E] = (unsigned short )ltb;
1490 goto done;
1491 }
1492 if( i > 0 )
1493 { /* put the larger number in bi */
1494 emovz( bi, ci );
1495 emovz( ai, bi );
1496 emovz( ci, ai );
1497 }
1498 }
1499if( ai[0] == bi[0] )
1500 {
1501 eaddm( ai, bi );
1502 subflg = 0;
1503 }
1504else
1505 {
1506 esubm( ai, bi );
1507 subflg = 1;
1508 }
1509emdnorm( bi, lost, subflg, ltb, 64 );
1510
1511done:
1512emovo( bi, c );
1513}
1514
1515
1516
1517/*
1518; Divide.
1519;
1520; unsigned short a[NE], b[NE], c[NE];
1521; ediv( a, b, c ); c = b / a
1522*/
1523void ediv( a, b, c )
1524unsigned short *a, *b, *c;
1525{
1526unsigned short ai[NI], bi[NI];
1527int i, sign;
1528long lt, lta, ltb;
1529
1530/* IEEE says if result is not a NaN, the sign is "-" if and only if
1531 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
1532sign = eisneg(a) ^ eisneg(b);
1533
1534#ifdef NANS
1535/* Return any NaN input. */
1536if( eisnan(a) )
1537 {
1538 emov(a,c);
1539 return;
1540 }
1541if( eisnan(b) )
1542 {
1543 emov(b,c);
1544 return;
1545 }
1546/* Zero over zero, or infinity over infinity, is a NaN. */
1547if( ((ecmp(a,ezero) == 0) && (ecmp(b,ezero) == 0))
1548 || (eisinf (a) && eisinf (b)) )
1549 {
1550 mtherr( "ediv", DOMAIN );
1551 enan( c, NBITS );
1552 return;
1553 }
1554#endif
1555/* Infinity over anything else is infinity. */
1556#ifdef INFINITY
1557if( eisinf(b) )
1558 {
1559 einfin(c);
1560 goto divsign;
1561 }
1562if( eisinf(a) )
1563 {
1564 eclear(c);
1565 goto divsign;
1566 }
1567#endif
1568emovi( a, ai );
1569emovi( b, bi );
1570lta = ai[E];
1571ltb = bi[E];
1572if( bi[E] == 0 )
1573 { /* See if numerator is zero. */
1574 for( i=1; i<NI-1; i++ )
1575 {
1576 if( bi[i] != 0 )
1577 {
1578 ltb -= enormlz( bi );
1579 goto dnzro1;
1580 }
1581 }
1582 eclear(c);
1583 goto divsign;
1584 }
1585dnzro1:
1586
1587if( ai[E] == 0 )
1588 { /* possible divide by zero */
1589 for( i=1; i<NI-1; i++ )
1590 {
1591 if( ai[i] != 0 )
1592 {
1593 lta -= enormlz( ai );
1594 goto dnzro2;
1595 }
1596 }
1597 einfin(c);
1598 mtherr( "ediv", SING );
1599 goto divsign;
1600 }
1601dnzro2:
1602
1603i = edivm( ai, bi );
1604/* calculate exponent */
1605lt = ltb - lta + EXONE;
1606emdnorm( bi, i, 0, lt, 64 );
1607emovo( bi, c );
1608
1609divsign:
1610
1611if( sign )
1612 *(c+(NE-1)) |= 0x8000;
1613else
1614 *(c+(NE-1)) &= ~0x8000;
1615}
1616
1617
1618
1619/*
1620; Multiply.
1621;
1622; unsigned short a[NE], b[NE], c[NE];
1623; emul( a, b, c ); c = b * a
1624*/
1625void emul( a, b, c )
1626unsigned short *a, *b, *c;
1627{
1628unsigned short ai[NI], bi[NI];
1629int i, j, sign;
1630long lt, lta, ltb;
1631
1632/* IEEE says if result is not a NaN, the sign is "-" if and only if
1633 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
1634sign = eisneg(a) ^ eisneg(b);
1635
1636#ifdef NANS
1637/* NaN times anything is the same NaN. */
1638if( eisnan(a) )
1639 {
1640 emov(a,c);
1641 return;
1642 }
1643if( eisnan(b) )
1644 {
1645 emov(b,c);
1646 return;
1647 }
1648/* Zero times infinity is a NaN. */
1649if( (eisinf(a) && (ecmp(b,ezero) == 0))
1650 || (eisinf(b) && (ecmp(a,ezero) == 0)) )
1651 {
1652 mtherr( "emul", DOMAIN );
1653 enan( c, NBITS );
1654 return;
1655 }
1656#endif
1657/* Infinity times anything else is infinity. */
1658#ifdef INFINITY
1659if( eisinf(a) || eisinf(b) )
1660 {
1661 einfin(c);
1662 goto mulsign;
1663 }
1664#endif
1665emovi( a, ai );
1666emovi( b, bi );
1667lta = ai[E];
1668ltb = bi[E];
1669if( ai[E] == 0 )
1670 {
1671 for( i=1; i<NI-1; i++ )
1672 {
1673 if( ai[i] != 0 )
1674 {
1675 lta -= enormlz( ai );
1676 goto mnzer1;
1677 }
1678 }
1679 eclear(c);
1680 goto mulsign;
1681 }
1682mnzer1:
1683
1684if( bi[E] == 0 )
1685 {
1686 for( i=1; i<NI-1; i++ )
1687 {
1688 if( bi[i] != 0 )
1689 {
1690 ltb -= enormlz( bi );
1691 goto mnzer2;
1692 }
1693 }
1694 eclear(c);
1695 goto mulsign;
1696 }
1697mnzer2:
1698
1699/* Multiply significands */
1700j = emulm( ai, bi );
1701/* calculate exponent */
1702lt = lta + ltb - (EXONE - 1);
1703emdnorm( bi, j, 0, lt, 64 );
1704emovo( bi, c );
1705/* IEEE says sign is "-" if and only if operands have opposite signs. */
1706mulsign:
1707if( sign )
1708 *(c+(NE-1)) |= 0x8000;
1709else
1710 *(c+(NE-1)) &= ~0x8000;
1711}
1712
1713
1714
1715
1716/*
1717; Convert IEEE double precision to e type
1718; double d;
1719; unsigned short x[N+2];
1720; e53toe( &d, x );
1721*/
1722void e53toe( pe, y )
1723unsigned short *pe, *y;
1724{
1725#ifdef DEC
1726
1727dectoe( pe, y ); /* see etodec.c */
1728
1729#else
1730
1731register unsigned short r;
1732register unsigned short *p, *e;
1733unsigned short yy[NI];
1734int denorm, k;
1735
1736e = pe;
1737denorm = 0; /* flag if denormalized number */
1738ecleaz(yy);
1739#ifdef IBMPC
1740e += 3;
1741#endif
1742r = *e;
1743yy[0] = 0;
1744if( r & 0x8000 )
1745 yy[0] = 0xffff;
1746yy[M] = (r & 0x0f) | 0x10;
1747r &= ~0x800f; /* strip sign and 4 significand bits */
1748#ifdef INFINITY
1749if( r == 0x7ff0 )
1750 {
1751#ifdef NANS
1752#ifdef IBMPC
1753 if( ((pe[3] & 0xf) != 0) || (pe[2] != 0)
1754 || (pe[1] != 0) || (pe[0] != 0) )
1755 {
1756 enan( y, NBITS );
1757 return;
1758 }
1759#else
1760 if( ((pe[0] & 0xf) != 0) || (pe[1] != 0)
1761 || (pe[2] != 0) || (pe[3] != 0) )
1762 {
1763 enan( y, NBITS );
1764 return;
1765 }
1766#endif
1767#endif /* NANS */
1768 eclear( y );
1769 einfin( y );
1770 if( yy[0] )
1771 eneg(y);
1772 return;
1773 }
1774#endif
1775r >>= 4;
1776/* If zero exponent, then the significand is denormalized.
1777 * So, take back the understood high significand bit. */
1778if( r == 0 )
1779 {
1780 denorm = 1;
1781 yy[M] &= ~0x10;
1782 }
1783r += EXONE - 01777;
1784yy[E] = r;
1785p = &yy[M+1];
1786#ifdef IBMPC
1787*p++ = *(--e);
1788*p++ = *(--e);
1789*p++ = *(--e);
1790#endif
1791#ifdef MIEEE
1792++e;
1793*p++ = *e++;
1794*p++ = *e++;
1795*p++ = *e++;
1796#endif
1797(void )eshift( yy, -5 );
1798if( denorm )
1799 { /* if zero exponent, then normalize the significand */
1800 if( (k = enormlz(yy)) > NBITS )
1801 ecleazs(yy);
1802 else
1803 yy[E] -= (unsigned short )(k-1);
1804 }
1805emovo( yy, y );
1806#endif /* not DEC */
1807}
1808
1809void e64toe( pe, y )
1810unsigned short *pe, *y;
1811{
1812unsigned short yy[NI];
1813unsigned short *p, *q, *e;
1814int i;
1815
1816e = pe;
1817p = yy;
1818for( i=0; i<NE-5; i++ )
1819 *p++ = 0;
1820#ifdef IBMPC
1821for( i=0; i<5; i++ )
1822 *p++ = *e++;
1823#endif
1824#ifdef DEC
1825for( i=0; i<5; i++ )
1826 *p++ = *e++;
1827#endif
1828#ifdef MIEEE
1829p = &yy[0] + (NE-1);
1830*p-- = *e++;
1831++e;
1832for( i=0; i<4; i++ )
1833 *p-- = *e++;
1834#endif
1835
1836#ifdef IBMPC
1837/* For Intel long double, shift denormal significand up 1
1838 -- but only if the top significand bit is zero. */
1839if((yy[NE-1] & 0x7fff) == 0 && (yy[NE-2] & 0x8000) == 0)
1840 {
1841 unsigned short temp[NI+1];
1842 emovi(yy, temp);
1843 eshup1(temp);
1844 emovo(temp,y);
1845 return;
1846 }
1847#endif
1848#ifdef INFINITY
1849/* Point to the exponent field. */
1850p = &yy[NE-1];
1851if ((*p & 0x7fff) == 0x7fff)
1852 {
1853#ifdef NANS
1854#ifdef IBMPC
1855 for( i=0; i<4; i++ )
1856 {
1857 if((i != 3 && pe[i] != 0)
1858 /* Check for Intel long double infinity pattern. */
1859 || (i == 3 && pe[i] != 0x8000))
1860 {
1861 enan( y, NBITS );
1862 return;
1863 }
1864 }
1865#else
1866 /* In Motorola extended precision format, the most significant
1867 bit of an infinity mantissa could be either 1 or 0. It is
1868 the lower order bits that tell whether the value is a NaN. */
1869 if ((pe[2] & 0x7fff) != 0)
1870 goto bigend_nan;
1871
1872 for( i=3; i<=5; i++ )
1873 {
1874 if( pe[i] != 0 )
1875 {
1876bigend_nan:
1877 enan( y, NBITS );
1878 return;
1879 }
1880 }
1881#endif
1882#endif /* NANS */
1883 eclear( y );
1884 einfin( y );
1885 if( *p & 0x8000 )
1886 eneg(y);
1887 return;
1888 }
1889#endif
1890p = yy;
1891q = y;
1892for( i=0; i<NE; i++ )
1893 *q++ = *p++;
1894}
1895
1896void e113toe(pe,y)
1897unsigned short *pe, *y;
1898{
1899register unsigned short r;
1900unsigned short *e, *p;
1901unsigned short yy[NI];
1902int denorm, i;
1903
1904e = pe;
1905denorm = 0;
1906ecleaz(yy);
1907#ifdef IBMPC
1908e += 7;
1909#endif
1910r = *e;
1911yy[0] = 0;
1912if( r & 0x8000 )
1913 yy[0] = 0xffff;
1914r &= 0x7fff;
1915#ifdef INFINITY
1916if( r == 0x7fff )
1917 {
1918#ifdef NANS
1919#ifdef IBMPC
1920 for( i=0; i<7; i++ )
1921 {
1922 if( pe[i] != 0 )
1923 {
1924 enan( y, NBITS );
1925 return;
1926 }
1927 }
1928#else
1929 for( i=1; i<8; i++ )
1930 {
1931 if( pe[i] != 0 )
1932 {
1933 enan( y, NBITS );
1934 return;
1935 }
1936 }
1937#endif
1938#endif /* NANS */
1939 eclear( y );
1940 einfin( y );
1941 if( *e & 0x8000 )
1942 eneg(y);
1943 return;
1944 }
1945#endif /* INFINITY */
1946yy[E] = r;
1947p = &yy[M + 1];
1948#ifdef IBMPC
1949for( i=0; i<7; i++ )
1950 *p++ = *(--e);
1951#endif
1952#ifdef MIEEE
1953++e;
1954for( i=0; i<7; i++ )
1955 *p++ = *e++;
1956#endif
1957/* If denormal, remove the implied bit; else shift down 1. */
1958if( r == 0 )
1959 {
1960 yy[M] = 0;
1961 }
1962else
1963 {
1964 yy[M] = 1;
1965 eshift( yy, -1 );
1966 }
1967emovo(yy,y);
1968}
1969
1970
1971/*
1972; Convert IEEE single precision to e type
1973; float d;
1974; unsigned short x[N+2];
1975; dtox( &d, x );
1976*/
1977void e24toe( pe, y )
1978unsigned short *pe, *y;
1979{
1980register unsigned short r;
1981register unsigned short *p, *e;
1982unsigned short yy[NI];
1983int denorm, k;
1984
1985e = pe;
1986denorm = 0; /* flag if denormalized number */
1987ecleaz(yy);
1988#ifdef IBMPC
1989e += 1;
1990#endif
1991#ifdef DEC
1992e += 1;
1993#endif
1994r = *e;
1995yy[0] = 0;
1996if( r & 0x8000 )
1997 yy[0] = 0xffff;
1998yy[M] = (r & 0x7f) | 0200;
1999r &= ~0x807f; /* strip sign and 7 significand bits */
2000#ifdef INFINITY
2001if( r == 0x7f80 )
2002 {
2003#ifdef NANS
2004#ifdef MIEEE
2005 if( ((pe[0] & 0x7f) != 0) || (pe[1] != 0) )
2006 {
2007 enan( y, NBITS );
2008 return;
2009 }
2010#else
2011 if( ((pe[1] & 0x7f) != 0) || (pe[0] != 0) )
2012 {
2013 enan( y, NBITS );
2014 return;
2015 }
2016#endif
2017#endif /* NANS */
2018 eclear( y );
2019 einfin( y );
2020 if( yy[0] )
2021 eneg(y);
2022 return;
2023 }
2024#endif
2025r >>= 7;
2026/* If zero exponent, then the significand is denormalized.
2027 * So, take back the understood high significand bit. */
2028if( r == 0 )
2029 {
2030 denorm = 1;
2031 yy[M] &= ~0200;
2032 }
2033r += EXONE - 0177;
2034yy[E] = r;
2035p = &yy[M+1];
2036#ifdef IBMPC
2037*p++ = *(--e);
2038#endif
2039#ifdef DEC
2040*p++ = *(--e);
2041#endif
2042#ifdef MIEEE
2043++e;
2044*p++ = *e++;
2045#endif
2046(void )eshift( yy, -8 );
2047if( denorm )
2048 { /* if zero exponent, then normalize the significand */
2049 if( (k = enormlz(yy)) > NBITS )
2050 ecleazs(yy);
2051 else
2052 yy[E] -= (unsigned short )(k-1);
2053 }
2054emovo( yy, y );
2055}
2056
2057void etoe113(x,e)
2058unsigned short *x, *e;
2059{
2060unsigned short xi[NI];
2061long exp;
2062int rndsav;
2063
2064#ifdef NANS
2065if( eisnan(x) )
2066 {
2067 enan( e, 113 );
2068 return;
2069 }
2070#endif
2071emovi( x, xi );
2072exp = (long )xi[E];
2073#ifdef INFINITY
2074if( eisinf(x) )
2075 goto nonorm;
2076#endif
2077/* round off to nearest or even */
2078rndsav = rndprc;
2079rndprc = 113;
2080emdnorm( xi, 0, 0, exp, 64 );
2081rndprc = rndsav;
2082nonorm:
2083toe113 (xi, e);
2084}
2085
2086/* move out internal format to ieee long double */
2087static void toe113(a,b)
2088unsigned short *a, *b;
2089{
2090register unsigned short *p, *q;
2091unsigned short i;
2092
2093#ifdef NANS
2094if( eiisnan(a) )
2095 {
2096 enan( b, 113 );
2097 return;
2098 }
2099#endif
2100p = a;
2101#ifdef MIEEE
2102q = b;
2103#else
2104q = b + 7; /* point to output exponent */
2105#endif
2106
2107/* If not denormal, delete the implied bit. */
2108if( a[E] != 0 )
2109 {
2110 eshup1 (a);
2111 }
2112/* combine sign and exponent */
2113i = *p++;
2114#ifdef MIEEE
2115if( i )
2116 *q++ = *p++ | 0x8000;
2117else
2118 *q++ = *p++;
2119#else
2120if( i )
2121 *q-- = *p++ | 0x8000;
2122else
2123 *q-- = *p++;
2124#endif
2125/* skip over guard word */
2126++p;
2127/* move the significand */
2128#ifdef MIEEE
2129for (i = 0; i < 7; i++)
2130 *q++ = *p++;
2131#else
2132for (i = 0; i < 7; i++)
2133 *q-- = *p++;
2134#endif
2135}
2136
2137
2138void etoe64( x, e )
2139unsigned short *x, *e;
2140{
2141unsigned short xi[NI];
2142long exp;
2143int rndsav;
2144
2145#ifdef NANS
2146if( eisnan(x) )
2147 {
2148 enan( e, 64 );
2149 return;
2150 }
2151#endif
2152emovi( x, xi );
2153exp = (long )xi[E]; /* adjust exponent for offset */
2154#ifdef INFINITY
2155if( eisinf(x) )
2156 goto nonorm;
2157#endif
2158/* round off to nearest or even */
2159rndsav = rndprc;
2160rndprc = 64;
2161emdnorm( xi, 0, 0, exp, 64 );
2162rndprc = rndsav;
2163nonorm:
2164toe64( xi, e );
2165}
2166
2167/* move out internal format to ieee long double */
2168static void toe64( a, b )
2169unsigned short *a, *b;
2170{
2171register unsigned short *p, *q;
2172unsigned short i;
2173
2174#ifdef NANS
2175if( eiisnan(a) )
2176 {
2177 enan( b, 64 );
2178 return;
2179 }
2180#endif
2181#ifdef IBMPC
2182/* Shift Intel denormal significand down 1. */
2183if( a[E] == 0 )
2184 eshdn1(a);
2185#endif
2186p = a;
2187#ifdef MIEEE
2188q = b;
2189#else
2190q = b + 4; /* point to output exponent */
2191#if 1
2192/* NOTE: if data type is 96 bits wide, clear the last word here. */
2193*(q+1)= 0;
2194#endif
2195#endif
2196
2197/* combine sign and exponent */
2198i = *p++;
2199#ifdef MIEEE
2200if( i )
2201 *q++ = *p++ | 0x8000;
2202else
2203 *q++ = *p++;
2204*q++ = 0;
2205#else
2206if( i )
2207 *q-- = *p++ | 0x8000;
2208else
2209 *q-- = *p++;
2210#endif
2211/* skip over guard word */
2212++p;
2213/* move the significand */
2214#ifdef MIEEE
2215for( i=0; i<4; i++ )
2216 *q++ = *p++;
2217#else
2218#ifdef INFINITY
2219if (eiisinf (a))
2220 {
2221 /* Intel long double infinity. */
2222 *q-- = 0x8000;
2223 *q-- = 0;
2224 *q-- = 0;
2225 *q = 0;
2226 return;
2227 }
2228#endif
2229for( i=0; i<4; i++ )
2230 *q-- = *p++;
2231#endif
2232}
2233
2234
2235/*
2236; e type to IEEE double precision
2237; double d;
2238; unsigned short x[NE];
2239; etoe53( x, &d );
2240*/
2241
2242#ifdef DEC
2243
2244void etoe53( x, e )
2245unsigned short *x, *e;
2246{
2247etodec( x, e ); /* see etodec.c */
2248}
2249
2250static void toe53( x, y )
2251unsigned short *x, *y;
2252{
2253todec( x, y );
2254}
2255
2256#else
2257
2258void etoe53( x, e )
2259unsigned short *x, *e;
2260{
2261unsigned short xi[NI];
2262long exp;
2263int rndsav;
2264
2265#ifdef NANS
2266if( eisnan(x) )
2267 {
2268 enan( e, 53 );
2269 return;
2270 }
2271#endif
2272emovi( x, xi );
2273exp = (long )xi[E] - (EXONE - 0x3ff); /* adjust exponent for offsets */
2274#ifdef INFINITY
2275if( eisinf(x) )
2276 goto nonorm;
2277#endif
2278/* round off to nearest or even */
2279rndsav = rndprc;
2280rndprc = 53;
2281emdnorm( xi, 0, 0, exp, 64 );
2282rndprc = rndsav;
2283nonorm:
2284toe53( xi, e );
2285}
2286
2287
2288static void toe53( x, y )
2289unsigned short *x, *y;
2290{
2291unsigned short i;
2292unsigned short *p;
2293
2294
2295#ifdef NANS
2296if( eiisnan(x) )
2297 {
2298 enan( y, 53 );
2299 return;
2300 }
2301#endif
2302p = &x[0];
2303#ifdef IBMPC
2304y += 3;
2305#endif
2306*y = 0; /* output high order */
2307if( *p++ )
2308 *y = 0x8000; /* output sign bit */
2309
2310i = *p++;
2311if( i >= (unsigned int )2047 )
2312 { /* Saturate at largest number less than infinity. */
2313#ifdef INFINITY
2314 *y |= 0x7ff0;
2315#ifdef IBMPC
2316 *(--y) = 0;
2317 *(--y) = 0;
2318 *(--y) = 0;
2319#endif
2320#ifdef MIEEE
2321 ++y;
2322 *y++ = 0;
2323 *y++ = 0;
2324 *y++ = 0;
2325#endif
2326#else
2327 *y |= (unsigned short )0x7fef;
2328#ifdef IBMPC
2329 *(--y) = 0xffff;
2330 *(--y) = 0xffff;
2331 *(--y) = 0xffff;
2332#endif
2333#ifdef MIEEE
2334 ++y;
2335 *y++ = 0xffff;
2336 *y++ = 0xffff;
2337 *y++ = 0xffff;
2338#endif
2339#endif
2340 return;
2341 }
2342if( i == 0 )
2343 {
2344 (void )eshift( x, 4 );
2345 }
2346else
2347 {
2348 i <<= 4;
2349 (void )eshift( x, 5 );
2350 }
2351i |= *p++ & (unsigned short )0x0f; /* *p = xi[M] */
2352*y |= (unsigned short )i; /* high order output already has sign bit set */
2353#ifdef IBMPC
2354*(--y) = *p++;
2355*(--y) = *p++;
2356*(--y) = *p;
2357#endif
2358#ifdef MIEEE
2359++y;
2360*y++ = *p++;
2361*y++ = *p++;
2362*y++ = *p++;
2363#endif
2364}
2365
2366#endif /* not DEC */
2367
2368
2369
2370/*
2371; e type to IEEE single precision
2372; float d;
2373; unsigned short x[N+2];
2374; xtod( x, &d );
2375*/
2376void etoe24( x, e )
2377unsigned short *x, *e;
2378{
2379long exp;
2380unsigned short xi[NI];
2381int rndsav;
2382
2383#ifdef NANS
2384if( eisnan(x) )
2385 {
2386 enan( e, 24 );
2387 return;
2388 }
2389#endif
2390emovi( x, xi );
2391exp = (long )xi[E] - (EXONE - 0177); /* adjust exponent for offsets */
2392#ifdef INFINITY
2393if( eisinf(x) )
2394 goto nonorm;
2395#endif
2396/* round off to nearest or even */
2397rndsav = rndprc;
2398rndprc = 24;
2399emdnorm( xi, 0, 0, exp, 64 );
2400rndprc = rndsav;
2401nonorm:
2402toe24( xi, e );
2403}
2404
2405static void toe24( x, y )
2406unsigned short *x, *y;
2407{
2408unsigned short i;
2409unsigned short *p;
2410
2411#ifdef NANS
2412if( eiisnan(x) )
2413 {
2414 enan( y, 24 );
2415 return;
2416 }
2417#endif
2418p = &x[0];
2419#ifdef IBMPC
2420y += 1;
2421#endif
2422#ifdef DEC
2423y += 1;
2424#endif
2425*y = 0; /* output high order */
2426if( *p++ )
2427 *y = 0x8000; /* output sign bit */
2428
2429i = *p++;
2430if( i >= 255 )
2431 { /* Saturate at largest number less than infinity. */
2432#ifdef INFINITY
2433 *y |= (unsigned short )0x7f80;
2434#ifdef IBMPC
2435 *(--y) = 0;
2436#endif
2437#ifdef DEC
2438 *(--y) = 0;
2439#endif
2440#ifdef MIEEE
2441 ++y;
2442 *y = 0;
2443#endif
2444#else
2445 *y |= (unsigned short )0x7f7f;
2446#ifdef IBMPC
2447 *(--y) = 0xffff;
2448#endif
2449#ifdef DEC
2450 *(--y) = 0xffff;
2451#endif
2452#ifdef MIEEE
2453 ++y;
2454 *y = 0xffff;
2455#endif
2456#endif
2457 return;
2458 }
2459if( i == 0 )
2460 {
2461 (void )eshift( x, 7 );
2462 }
2463else
2464 {
2465 i <<= 7;
2466 (void )eshift( x, 8 );
2467 }
2468i |= *p++ & (unsigned short )0x7f; /* *p = xi[M] */
2469*y |= i; /* high order output already has sign bit set */
2470#ifdef IBMPC
2471*(--y) = *p;
2472#endif
2473#ifdef DEC
2474*(--y) = *p;
2475#endif
2476#ifdef MIEEE
2477++y;
2478*y = *p;
2479#endif
2480}
2481
2482
2483/* Compare two e type numbers.
2484 *
2485 * unsigned short a[NE], b[NE];
2486 * ecmp( a, b );
2487 *
2488 * returns +1 if a > b
2489 * 0 if a == b
2490 * -1 if a < b
2491 * -2 if either a or b is a NaN.
2492 */
2493int ecmp( a, b )
2494unsigned short *a, *b;
2495{
2496unsigned short ai[NI], bi[NI];
2497register unsigned short *p, *q;
2498register int i;
2499int msign;
2500
2501#ifdef NANS
2502if (eisnan (a) || eisnan (b))
2503 return( -2 );
2504#endif
2505emovi( a, ai );
2506p = ai;
2507emovi( b, bi );
2508q = bi;
2509
2510if( *p != *q )
2511 { /* the signs are different */
2512/* -0 equals + 0 */
2513 for( i=1; i<NI-1; i++ )
2514 {
2515 if( ai[i] != 0 )
2516 goto nzro;
2517 if( bi[i] != 0 )
2518 goto nzro;
2519 }
2520 return(0);
2521nzro:
2522 if( *p == 0 )
2523 return( 1 );
2524 else
2525 return( -1 );
2526 }
2527/* both are the same sign */
2528if( *p == 0 )
2529 msign = 1;
2530else
2531 msign = -1;
2532i = NI-1;
2533do
2534 {
2535 if( *p++ != *q++ )
2536 {
2537 goto diff;
2538 }
2539 }
2540while( --i > 0 );
2541
2542return(0); /* equality */
2543
2544
2545
2546diff:
2547
2548if( *(--p) > *(--q) )
2549 return( msign ); /* p is bigger */
2550else
2551 return( -msign ); /* p is littler */
2552}
2553
2554
2555
2556
2557/* Find nearest integer to x = floor( x + 0.5 )
2558 *
2559 * unsigned short x[NE], y[NE]
2560 * eround( x, y );
2561 */
2562void eround( x, y )
2563unsigned short *x, *y;
2564{
2565
2566eadd( ehalf, x, y );
2567efloor( y, y );
2568}
2569
2570
2571
2572
2573/*
2574; convert long (32-bit) integer to e type
2575;
2576; long l;
2577; unsigned short x[NE];
2578; ltoe( &l, x );
2579; note &l is the memory address of l
2580*/
2581void ltoe( lp, y )
2582long *lp; /* lp is the memory address of a long integer */
2583unsigned short *y; /* y is the address of a short */
2584{
2585unsigned short yi[NI];
2586unsigned long ll;
2587int k;
2588
2589ecleaz( yi );
2590if( *lp < 0 )
2591 {
2592 ll = (unsigned long )( -(*lp) ); /* make it positive */
2593 yi[0] = 0xffff; /* put correct sign in the e type number */
2594 }
2595else
2596 {
2597 ll = (unsigned long )( *lp );
2598 }
2599/* move the long integer to yi significand area */
2600if( sizeof(long) == 8 )
2601 {
2602 yi[M] = (unsigned short) (ll >> (LONGBITS - 16));
2603 yi[M + 1] = (unsigned short) (ll >> (LONGBITS - 32));
2604 yi[M + 2] = (unsigned short) (ll >> 16);
2605 yi[M + 3] = (unsigned short) ll;
2606 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
2607 }
2608else
2609 {
2610 yi[M] = (unsigned short )(ll >> 16);
2611 yi[M+1] = (unsigned short )ll;
2612 yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
2613 }
2614if( (k = enormlz( yi )) > NBITS ) /* normalize the significand */
2615 ecleaz( yi ); /* it was zero */
2616else
2617 yi[E] -= (unsigned short )k; /* subtract shift count from exponent */
2618emovo( yi, y ); /* output the answer */
2619}
2620
2621/*
2622; convert unsigned long (32-bit) integer to e type
2623;
2624; unsigned long l;
2625; unsigned short x[NE];
2626; ltox( &l, x );
2627; note &l is the memory address of l
2628*/
2629void ultoe( lp, y )
2630unsigned long *lp; /* lp is the memory address of a long integer */
2631unsigned short *y; /* y is the address of a short */
2632{
2633unsigned short yi[NI];
2634unsigned long ll;
2635int k;
2636
2637ecleaz( yi );
2638ll = *lp;
2639
2640/* move the long integer to ayi significand area */
2641if( sizeof(long) == 8 )
2642 {
2643 yi[M] = (unsigned short) (ll >> (LONGBITS - 16));
2644 yi[M + 1] = (unsigned short) (ll >> (LONGBITS - 32));
2645 yi[M + 2] = (unsigned short) (ll >> 16);
2646 yi[M + 3] = (unsigned short) ll;
2647 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
2648 }
2649else
2650 {
2651 yi[M] = (unsigned short )(ll >> 16);
2652 yi[M+1] = (unsigned short )ll;
2653 yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
2654 }
2655if( (k = enormlz( yi )) > NBITS ) /* normalize the significand */
2656 ecleaz( yi ); /* it was zero */
2657else
2658 yi[E] -= (unsigned short )k; /* subtract shift count from exponent */
2659emovo( yi, y ); /* output the answer */
2660}
2661
2662
2663/*
2664; Find long integer and fractional parts
2665
2666; long i;
2667; unsigned short x[NE], frac[NE];
2668; xifrac( x, &i, frac );
2669
2670 The integer output has the sign of the input. The fraction is
2671 the positive fractional part of abs(x).
2672*/
2673void eifrac( x, i, frac )
2674unsigned short *x;
2675long *i;
2676unsigned short *frac;
2677{
2678unsigned short xi[NI];
2679int j, k;
2680unsigned long ll;
2681
2682emovi( x, xi );
2683k = (int )xi[E] - (EXONE - 1);
2684if( k <= 0 )
2685 {
2686/* if exponent <= 0, integer = 0 and real output is fraction */
2687 *i = 0L;
2688 emovo( xi, frac );
2689 return;
2690 }
2691if( k > (8 * sizeof(long) - 1) )
2692 {
2693/*
2694; long integer overflow: output large integer
2695; and correct fraction
2696*/
2697 j = 8 * sizeof(long) - 1;
2698 if( xi[0] )
2699 *i = (long) ((unsigned long) 1) << j;
2700 else
2701 *i = (long) (((unsigned long) (~(0L))) >> 1);
2702 (void )eshift( xi, k );
2703 }
2704if( k > 16 )
2705 {
2706/*
2707 Shift more than 16 bits: shift up k-16 mod 16
2708 then shift by 16's.
2709*/
2710 j = k - ((k >> 4) << 4);
2711 eshift (xi, j);
2712 ll = xi[M];
2713 k -= j;
2714 do
2715 {
2716 eshup6 (xi);
2717 ll = (ll << 16) | xi[M];
2718 }
2719 while ((k -= 16) > 0);
2720 *i = ll;
2721 if (xi[0])
2722 *i = -(*i);
2723 }
2724else
2725 {
2726/* shift not more than 16 bits */
2727 eshift( xi, k );
2728 *i = (long )xi[M] & 0xffff;
2729 if( xi[0] )
2730 *i = -(*i);
2731 }
2732xi[0] = 0;
2733xi[E] = EXONE - 1;
2734xi[M] = 0;
2735if( (k = enormlz( xi )) > NBITS )
2736 ecleaz( xi );
2737else
2738 xi[E] -= (unsigned short )k;
2739
2740emovo( xi, frac );
2741}
2742
2743
2744/*
2745; Find unsigned long integer and fractional parts
2746
2747; unsigned long i;
2748; unsigned short x[NE], frac[NE];
2749; xifrac( x, &i, frac );
2750
2751 A negative e type input yields integer output = 0
2752 but correct fraction.
2753*/
2754void euifrac( x, i, frac )
2755unsigned short *x;
2756unsigned long *i;
2757unsigned short *frac;
2758{
2759unsigned short xi[NI];
2760int j, k;
2761unsigned long ll;
2762
2763emovi( x, xi );
2764k = (int )xi[E] - (EXONE - 1);
2765if( k <= 0 )
2766 {
2767/* if exponent <= 0, integer = 0 and argument is fraction */
2768 *i = 0L;
2769 emovo( xi, frac );
2770 return;
2771 }
2772if( k > (8 * sizeof(long)) )
2773 {
2774/*
2775; long integer overflow: output large integer
2776; and correct fraction
2777*/
2778 *i = ~(0L);
2779 (void )eshift( xi, k );
2780 }
2781else if( k > 16 )
2782 {
2783/*
2784 Shift more than 16 bits: shift up k-16 mod 16
2785 then shift up by 16's.
2786*/
2787 j = k - ((k >> 4) << 4);
2788 eshift (xi, j);
2789 ll = xi[M];
2790 k -= j;
2791 do
2792 {
2793 eshup6 (xi);
2794 ll = (ll << 16) | xi[M];
2795 }
2796 while ((k -= 16) > 0);
2797 *i = ll;
2798 }
2799else
2800 {
2801/* shift not more than 16 bits */
2802 eshift( xi, k );
2803 *i = (long )xi[M] & 0xffff;
2804 }
2805
2806if( xi[0] ) /* A negative value yields unsigned integer 0. */
2807 *i = 0L;
2808
2809xi[0] = 0;
2810xi[E] = EXONE - 1;
2811xi[M] = 0;
2812if( (k = enormlz( xi )) > NBITS )
2813 ecleaz( xi );
2814else
2815 xi[E] -= (unsigned short )k;
2816
2817emovo( xi, frac );
2818}
2819
2820
2821
2822/*
2823; Shift significand
2824;
2825; Shifts significand area up or down by the number of bits
2826; given by the variable sc.
2827*/
2828int eshift( x, sc )
2829unsigned short *x;
2830int sc;
2831{
2832unsigned short lost;
2833unsigned short *p;
2834
2835if( sc == 0 )
2836 return( 0 );
2837
2838lost = 0;
2839p = x + NI-1;
2840
2841if( sc < 0 )
2842 {
2843 sc = -sc;
2844 while( sc >= 16 )
2845 {
2846 lost |= *p; /* remember lost bits */
2847 eshdn6(x);
2848 sc -= 16;
2849 }
2850
2851 while( sc >= 8 )
2852 {
2853 lost |= *p & 0xff;
2854 eshdn8(x);
2855 sc -= 8;
2856 }
2857
2858 while( sc > 0 )
2859 {
2860 lost |= *p & 1;
2861 eshdn1(x);
2862 sc -= 1;
2863 }
2864 }
2865else
2866 {
2867 while( sc >= 16 )
2868 {
2869 eshup6(x);
2870 sc -= 16;
2871 }
2872
2873 while( sc >= 8 )
2874 {
2875 eshup8(x);
2876 sc -= 8;
2877 }
2878
2879 while( sc > 0 )
2880 {
2881 eshup1(x);
2882 sc -= 1;
2883 }
2884 }
2885if( lost )
2886 lost = 1;
2887return( (int )lost );
2888}
2889
2890
2891
2892/*
2893; normalize
2894;
2895; Shift normalizes the significand area pointed to by argument
2896; shift count (up = positive) is returned.
2897*/
2898int enormlz(x)
2899unsigned short x[];
2900{
2901register unsigned short *p;
2902int sc;
2903
2904sc = 0;
2905p = &x[M];
2906if( *p != 0 )
2907 goto normdn;
2908++p;
2909if( *p & 0x8000 )
2910 return( 0 ); /* already normalized */
2911while( *p == 0 )
2912 {
2913 eshup6(x);
2914 sc += 16;
2915/* With guard word, there are NBITS+16 bits available.
2916 * return true if all are zero.
2917 */
2918 if( sc > NBITS )
2919 return( sc );
2920 }
2921/* see if high byte is zero */
2922while( (*p & 0xff00) == 0 )
2923 {
2924 eshup8(x);
2925 sc += 8;
2926 }
2927/* now shift 1 bit at a time */
2928while( (*p & 0x8000) == 0)
2929 {
2930 eshup1(x);
2931 sc += 1;
2932 if( sc > (NBITS+16) )
2933 {
2934 mtherr( "enormlz", UNDERFLOW );
2935 return( sc );
2936 }
2937 }
2938return( sc );
2939
2940/* Normalize by shifting down out of the high guard word
2941 of the significand */
2942normdn:
2943
2944if( *p & 0xff00 )
2945 {
2946 eshdn8(x);
2947 sc -= 8;
2948 }
2949while( *p != 0 )
2950 {
2951 eshdn1(x);
2952 sc -= 1;
2953
2954 if( sc < -NBITS )
2955 {
2956 mtherr( "enormlz", OVERFLOW );
2957 return( sc );
2958 }
2959 }
2960return( sc );
2961}
2962
2963
2964
2965
2966/* Convert e type number to decimal format ASCII string.
2967 * The constants are for 64 bit precision.
2968 */
2969
2970#define NTEN 12
2971#define MAXP 4096
2972
2973#if NE == 10
2974static unsigned short etens[NTEN + 1][NE] =
2975{
2976 {0x6576, 0x4a92, 0x804a, 0x153f,
2977 0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
2978 {0x6a32, 0xce52, 0x329a, 0x28ce,
2979 0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
2980 {0x526c, 0x50ce, 0xf18b, 0x3d28,
2981 0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
2982 {0x9c66, 0x58f8, 0xbc50, 0x5c54,
2983 0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
2984 {0x851e, 0xeab7, 0x98fe, 0x901b,
2985 0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
2986 {0x0235, 0x0137, 0x36b1, 0x336c,
2987 0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
2988 {0x50f8, 0x25fb, 0xc76b, 0x6b71,
2989 0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
2990 {0x0000, 0x0000, 0x0000, 0x0000,
2991 0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
2992 {0x0000, 0x0000, 0x0000, 0x0000,
2993 0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
2994 {0x0000, 0x0000, 0x0000, 0x0000,
2995 0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
2996 {0x0000, 0x0000, 0x0000, 0x0000,
2997 0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
2998 {0x0000, 0x0000, 0x0000, 0x0000,
2999 0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
3000 {0x0000, 0x0000, 0x0000, 0x0000,
3001 0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
3002};
3003
3004static unsigned short emtens[NTEN + 1][NE] =
3005{
3006 {0x2030, 0xcffc, 0xa1c3, 0x8123,
3007 0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
3008 {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
3009 0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
3010 {0xf53f, 0xf698, 0x6bd3, 0x0158,
3011 0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
3012 {0xe731, 0x04d4, 0xe3f2, 0xd332,
3013 0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
3014 {0xa23e, 0x5308, 0xfefb, 0x1155,
3015 0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
3016 {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
3017 0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
3018 {0x2a20, 0x6224, 0x47b3, 0x98d7,
3019 0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
3020 {0x0b5b, 0x4af2, 0xa581, 0x18ed,
3021 0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
3022 {0xbf71, 0xa9b3, 0x7989, 0xbe68,
3023 0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
3024 {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
3025 0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
3026 {0xc155, 0xa4a8, 0x404e, 0x6113,
3027 0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
3028 {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
3029 0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
3030 {0xcccd, 0xcccc, 0xcccc, 0xcccc,
3031 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
3032};
3033#else
3034static unsigned short etens[NTEN+1][NE] = {
3035{0xc94c,0x979a,0x8a20,0x5202,0xc460,0x7525,},/* 10**4096 */
3036{0xa74d,0x5de4,0xc53d,0x3b5d,0x9e8b,0x5a92,},/* 10**2048 */
3037{0x650d,0x0c17,0x8175,0x7586,0xc976,0x4d48,},
3038{0xcc65,0x91c6,0xa60e,0xa0ae,0xe319,0x46a3,},
3039{0xddbc,0xde8d,0x9df9,0xebfb,0xaa7e,0x4351,},
3040{0xc66f,0x8cdf,0x80e9,0x47c9,0x93ba,0x41a8,},
3041{0x3cbf,0xa6d5,0xffcf,0x1f49,0xc278,0x40d3,},
3042{0xf020,0xb59d,0x2b70,0xada8,0x9dc5,0x4069,},
3043{0x0000,0x0000,0x0400,0xc9bf,0x8e1b,0x4034,},
3044{0x0000,0x0000,0x0000,0x2000,0xbebc,0x4019,},
3045{0x0000,0x0000,0x0000,0x0000,0x9c40,0x400c,},
3046{0x0000,0x0000,0x0000,0x0000,0xc800,0x4005,},
3047{0x0000,0x0000,0x0000,0x0000,0xa000,0x4002,}, /* 10**1 */
3048};
3049
3050static unsigned short emtens[NTEN+1][NE] = {
3051{0x2de4,0x9fde,0xd2ce,0x04c8,0xa6dd,0x0ad8,}, /* 10**-4096 */
3052{0x4925,0x2de4,0x3436,0x534f,0xceae,0x256b,}, /* 10**-2048 */
3053{0x87a6,0xc0bd,0xda57,0x82a5,0xa2a6,0x32b5,},
3054{0x7133,0xd21c,0xdb23,0xee32,0x9049,0x395a,},
3055{0xfa91,0x1939,0x637a,0x4325,0xc031,0x3cac,},
3056{0xac7d,0xe4a0,0x64bc,0x467c,0xddd0,0x3e55,},
3057{0x3f24,0xe9a5,0xa539,0xea27,0xa87f,0x3f2a,},
3058{0x67de,0x94ba,0x4539,0x1ead,0xcfb1,0x3f94,},
3059{0x4c2f,0xe15b,0xc44d,0x94be,0xe695,0x3fc9,},
3060{0xfdc2,0xcefc,0x8461,0x7711,0xabcc,0x3fe4,},
3061{0xd3c3,0x652b,0xe219,0x1758,0xd1b7,0x3ff1,},
3062{0x3d71,0xd70a,0x70a3,0x0a3d,0xa3d7,0x3ff8,},
3063{0xcccd,0xcccc,0xcccc,0xcccc,0xcccc,0x3ffb,}, /* 10**-1 */
3064};
3065#endif
3066
3067void e24toasc( x, string, ndigs )
3068unsigned short x[];
3069char *string;
3070int ndigs;
3071{
3072unsigned short w[NI];
3073
3074e24toe( x, w );
3075etoasc( w, string, ndigs );
3076}
3077
3078
3079void e53toasc( x, string, ndigs )
3080unsigned short x[];
3081char *string;
3082int ndigs;
3083{
3084unsigned short w[NI];
3085
3086e53toe( x, w );
3087etoasc( w, string, ndigs );
3088}
3089
3090
3091void e64toasc( x, string, ndigs )
3092unsigned short x[];
3093char *string;
3094int ndigs;
3095{
3096unsigned short w[NI];
3097
3098e64toe( x, w );
3099etoasc( w, string, ndigs );
3100}
3101
3102void e113toasc (x, string, ndigs)
3103unsigned short x[];
3104char *string;
3105int ndigs;
3106{
3107unsigned short w[NI];
3108
3109e113toe (x, w);
3110etoasc (w, string, ndigs);
3111}
3112
3113
3114void etoasc( x, string, ndigs )
3115unsigned short x[];
3116char *string;
3117int ndigs;
3118{
3119long digit;
3120unsigned short y[NI], t[NI], u[NI], w[NI];
3121unsigned short *p, *r, *ten;
3122unsigned short sign;
3123int i, j, k, expon, rndsav;
3124char *s, *ss;
3125unsigned short m;
3126
3127rndsav = rndprc;
3128#ifdef NANS
3129if( eisnan(x) )
3130 {
3131 sprintf( string, " NaN " );
3132 goto bxit;
3133 }
3134#endif
3135rndprc = NBITS; /* set to full precision */
3136emov( x, y ); /* retain external format */
3137if( y[NE-1] & 0x8000 )
3138 {
3139 sign = 0xffff;
3140 y[NE-1] &= 0x7fff;
3141 }
3142else
3143 {
3144 sign = 0;
3145 }
3146expon = 0;
3147ten = &etens[NTEN][0];
3148emov( eone, t );
3149/* Test for zero exponent */
3150if( y[NE-1] == 0 )
3151 {
3152 for( k=0; k<NE-1; k++ )
3153 {
3154 if( y[k] != 0 )
3155 goto tnzro; /* denormalized number */
3156 }
3157 goto isone; /* legal all zeros */
3158 }
3159tnzro:
3160
3161/* Test for infinity.
3162 */
3163if( y[NE-1] == 0x7fff )
3164 {
3165 if( sign )
3166 sprintf( string, " -Infinity " );
3167 else
3168 sprintf( string, " Infinity " );
3169 goto bxit;
3170 }
3171
3172/* Test for exponent nonzero but significand denormalized.
3173 * This is an error condition.
3174 */
3175if( (y[NE-1] != 0) && ((y[NE-2] & 0x8000) == 0) )
3176 {
3177 mtherr( "etoasc", DOMAIN );
3178 sprintf( string, "NaN" );
3179 goto bxit;
3180 }
3181
3182/* Compare to 1.0 */
3183i = ecmp( eone, y );
3184if( i == 0 )
3185 goto isone;
3186
3187if( i < 0 )
3188 { /* Number is greater than 1 */
3189/* Convert significand to an integer and strip trailing decimal zeros. */
3190 emov( y, u );
3191 u[NE-1] = EXONE + NBITS - 1;
3192
3193 p = &etens[NTEN-4][0];
3194 m = 16;
3195do
3196 {
3197 ediv( p, u, t );
3198 efloor( t, w );
3199 for( j=0; j<NE-1; j++ )
3200 {
3201 if( t[j] != w[j] )
3202 goto noint;
3203 }
3204 emov( t, u );
3205 expon += (int )m;
3206noint:
3207 p += NE;
3208 m >>= 1;
3209 }
3210while( m != 0 );
3211
3212/* Rescale from integer significand */
3213 u[NE-1] += y[NE-1] - (unsigned int )(EXONE + NBITS - 1);
3214 emov( u, y );
3215/* Find power of 10 */
3216 emov( eone, t );
3217 m = MAXP;
3218 p = &etens[0][0];
3219 while( ecmp( ten, u ) <= 0 )
3220 {
3221 if( ecmp( p, u ) <= 0 )
3222 {
3223 ediv( p, u, u );
3224 emul( p, t, t );
3225 expon += (int )m;
3226 }
3227 m >>= 1;
3228 if( m == 0 )
3229 break;
3230 p += NE;
3231 }
3232 }
3233else
3234 { /* Number is less than 1.0 */
3235/* Pad significand with trailing decimal zeros. */
3236 if( y[NE-1] == 0 )
3237 {
3238 while( (y[NE-2] & 0x8000) == 0 )
3239 {
3240 emul( ten, y, y );
3241 expon -= 1;
3242 }
3243 }
3244 else
3245 {
3246 emovi( y, w );
3247 for( i=0; i<NDEC+1; i++ )
3248 {
3249 if( (w[NI-1] & 0x7) != 0 )
3250 break;
3251/* multiply by 10 */
3252 emovz( w, u );
3253 eshdn1( u );
3254 eshdn1( u );
3255 eaddm( w, u );
3256 u[1] += 3;
3257 while( u[2] != 0 )
3258 {
3259 eshdn1(u);
3260 u[1] += 1;
3261 }
3262 if( u[NI-1] != 0 )
3263 break;
3264 if( eone[NE-1] <= u[1] )
3265 break;
3266 emovz( u, w );
3267 expon -= 1;
3268 }
3269 emovo( w, y );
3270 }
3271 k = -MAXP;
3272 p = &emtens[0][0];
3273 r = &etens[0][0];
3274 emov( y, w );
3275 emov( eone, t );
3276 while( ecmp( eone, w ) > 0 )
3277 {
3278 if( ecmp( p, w ) >= 0 )
3279 {
3280 emul( r, w, w );
3281 emul( r, t, t );
3282 expon += k;
3283 }
3284 k /= 2;
3285 if( k == 0 )
3286 break;
3287 p += NE;
3288 r += NE;
3289 }
3290 ediv( t, eone, t );
3291 }
3292isone:
3293/* Find the first (leading) digit. */
3294emovi( t, w );
3295emovz( w, t );
3296emovi( y, w );
3297emovz( w, y );
3298eiremain( t, y );
3299digit = equot[NI-1];
3300while( (digit == 0) && (ecmp(y,ezero) != 0) )
3301 {
3302 eshup1( y );
3303 emovz( y, u );
3304 eshup1( u );
3305 eshup1( u );
3306 eaddm( u, y );
3307 eiremain( t, y );
3308 digit = equot[NI-1];
3309 expon -= 1;
3310 }
3311s = string;
3312if( sign )
3313 *s++ = '-';
3314else
3315 *s++ = ' ';
3316/* Examine number of digits requested by caller. */
3317if( ndigs < 0 )
3318 ndigs = 0;
3319if( ndigs > NDEC )
3320 ndigs = NDEC;
3321if( digit == 10 )
3322 {
3323 *s++ = '1';
3324 *s++ = '.';
3325 if( ndigs > 0 )
3326 {
3327 *s++ = '0';
3328 ndigs -= 1;
3329 }
3330 expon += 1;
3331 }
3332else
3333 {
3334 *s++ = (char )digit + '0';
3335 *s++ = '.';
3336 }
3337/* Generate digits after the decimal point. */
3338for( k=0; k<=ndigs; k++ )
3339 {
3340/* multiply current number by 10, without normalizing */
3341 eshup1( y );
3342 emovz( y, u );
3343 eshup1( u );
3344 eshup1( u );
3345 eaddm( u, y );
3346 eiremain( t, y );
3347 *s++ = (char )equot[NI-1] + '0';
3348 }
3349digit = equot[NI-1];
3350--s;
3351ss = s;
3352/* round off the ASCII string */
3353if( digit > 4 )
3354 {
3355/* Test for critical rounding case in ASCII output. */
3356 if( digit == 5 )
3357 {
3358 emovo( y, t );
3359 if( ecmp(t,ezero) != 0 )
3360 goto roun; /* round to nearest */
3361 if( (*(s-1) & 1) == 0 )
3362 goto doexp; /* round to even */
3363 }
3364/* Round up and propagate carry-outs */
3365roun:
3366 --s;
3367 k = *s & 0x7f;
3368/* Carry out to most significant digit? */
3369 if( k == '.' )
3370 {
3371 --s;
3372 k = *s;
3373 k += 1;
3374 *s = (char )k;
3375/* Most significant digit carries to 10? */
3376 if( k > '9' )
3377 {
3378 expon += 1;
3379 *s = '1';
3380 }
3381 goto doexp;
3382 }
3383/* Round up and carry out from less significant digits */
3384 k += 1;
3385 *s = (char )k;
3386 if( k > '9' )
3387 {
3388 *s = '0';
3389 goto roun;
3390 }
3391 }
3392doexp:
3393/*
3394if( expon >= 0 )
3395 sprintf( ss, "e+%d", expon );
3396else
3397 sprintf( ss, "e%d", expon );
3398*/
3399 sprintf( ss, "E%d", expon );
3400bxit:
3401rndprc = rndsav;
3402}
3403
3404
3405
3406
3407/*
3408; ASCTOQ
3409; ASCTOQ.MAC LATEST REV: 11 JAN 84
3410; SLM, 3 JAN 78
3411;
3412; Convert ASCII string to quadruple precision floating point
3413;
3414; Numeric input is free field decimal number
3415; with max of 15 digits with or without
3416; decimal point entered as ASCII from teletype.
3417; Entering E after the number followed by a second
3418; number causes the second number to be interpreted
3419; as a power of 10 to be multiplied by the first number
3420; (i.e., "scientific" notation).
3421;
3422; Usage:
3423; asctoq( string, q );
3424*/
3425
3426/* ASCII to single */
3427void asctoe24( s, y )
3428char *s;
3429unsigned short *y;
3430{
3431asctoeg( s, y, 24 );
3432}
3433
3434
3435/* ASCII to double */
3436void asctoe53( s, y )
3437char *s;
3438unsigned short *y;
3439{
3440#ifdef DEC
3441asctoeg( s, y, 56 );
3442#else
3443asctoeg( s, y, 53 );
3444#endif
3445}
3446
3447
3448/* ASCII to long double */
3449void asctoe64( s, y )
3450char *s;
3451unsigned short *y;
3452{
3453asctoeg( s, y, 64 );
3454}
3455
3456/* ASCII to 128-bit long double */
3457void asctoe113 (s, y)
3458char *s;
3459unsigned short *y;
3460{
3461asctoeg( s, y, 113 );
3462}
3463
3464/* ASCII to super double */
3465void asctoe( s, y )
3466char *s;
3467unsigned short *y;
3468{
3469asctoeg( s, y, NBITS );
3470}
3471
3472/* Space to make a copy of the input string: */
3473static char lstr[82] = {0};
3474
3475void asctoeg( ss, y, oprec )
3476char *ss;
3477unsigned short *y;
3478int oprec;
3479{
3480unsigned short yy[NI], xt[NI], tt[NI];
3481int esign, decflg, sgnflg, nexp, exp, prec, lost;
3482int k, trail, c, rndsav;
3483long lexp;
3484unsigned short nsign, *p;
3485char *sp, *s;
3486
3487/* Copy the input string. */
3488s = ss;
3489while( *s == ' ' ) /* skip leading spaces */
3490 ++s;
3491sp = lstr;
3492for( k=0; k<79; k++ )
3493 {
3494 if( (*sp++ = *s++) == '\0' )
3495 break;
3496 }
3497*sp = '\0';
3498s = lstr;
3499
3500rndsav = rndprc;
3501rndprc = NBITS; /* Set to full precision */
3502lost = 0;
3503nsign = 0;
3504decflg = 0;
3505sgnflg = 0;
3506nexp = 0;
3507exp = 0;
3508prec = 0;
3509ecleaz( yy );
3510trail = 0;
3511
3512nxtcom:
3513k = *s - '0';
3514if( (k >= 0) && (k <= 9) )
3515 {
3516/* Ignore leading zeros */
3517 if( (prec == 0) && (decflg == 0) && (k == 0) )
3518 goto donchr;
3519/* Identify and strip trailing zeros after the decimal point. */
3520 if( (trail == 0) && (decflg != 0) )
3521 {
3522 sp = s;
3523 while( (*sp >= '0') && (*sp <= '9') )
3524 ++sp;
3525/* Check for syntax error */
3526 c = *sp & 0x7f;
3527 if( (c != 'e') && (c != 'E') && (c != '\0')
3528 && (c != '\n') && (c != '\r') && (c != ' ')
3529 && (c != ',') )
3530 goto error;
3531 --sp;
3532 while( *sp == '0' )
3533 *sp-- = 'z';
3534 trail = 1;
3535 if( *s == 'z' )
3536 goto donchr;
3537 }
3538/* If enough digits were given to more than fill up the yy register,
3539 * continuing until overflow into the high guard word yy[2]
3540 * guarantees that there will be a roundoff bit at the top
3541 * of the low guard word after normalization.
3542 */
3543 if( yy[2] == 0 )
3544 {
3545 if( decflg )
3546 nexp += 1; /* count digits after decimal point */
3547 eshup1( yy ); /* multiply current number by 10 */
3548 emovz( yy, xt );
3549 eshup1( xt );
3550 eshup1( xt );
3551 eaddm( xt, yy );
3552 ecleaz( xt );
3553 xt[NI-2] = (unsigned short )k;
3554 eaddm( xt, yy );
3555 }
3556 else
3557 {
3558 /* Mark any lost non-zero digit. */
3559 lost |= k;
3560 /* Count lost digits before the decimal point. */
3561 if (decflg == 0)
3562 nexp -= 1;
3563 }
3564 prec += 1;
3565 goto donchr;
3566 }
3567
3568switch( *s )
3569 {
3570 case 'z':
3571 break;
3572 case 'E':
3573 case 'e':
3574 goto expnt;
3575 case '.': /* decimal point */
3576 if( decflg )
3577 goto error;
3578 ++decflg;
3579 break;
3580 case '-':
3581 nsign = 0xffff;
3582 if( sgnflg )
3583 goto error;
3584 ++sgnflg;
3585 break;
3586 case '+':
3587 if( sgnflg )
3588 goto error;
3589 ++sgnflg;
3590 break;
3591 case ',':
3592 case ' ':
3593 case '\0':
3594 case '\n':
3595 case '\r':
3596 goto daldone;
3597 case 'i':
3598 case 'I':
3599 goto infinite;
3600 default:
3601 error:
3602#ifdef NANS
3603 enan( yy, NI*16 );
3604#else
3605 mtherr( "asctoe", DOMAIN );
3606 ecleaz(yy);
3607#endif
3608 goto aexit;
3609 }
3610donchr:
3611++s;
3612goto nxtcom;
3613
3614/* Exponent interpretation */
3615expnt:
3616
3617esign = 1;
3618exp = 0;
3619++s;
3620/* check for + or - */
3621if( *s == '-' )
3622 {
3623 esign = -1;
3624 ++s;
3625 }
3626if( *s == '+' )
3627 ++s;
3628while( (*s >= '0') && (*s <= '9') )
3629 {
3630 exp *= 10;
3631 exp += *s++ - '0';
3632 if (exp > 4977)
3633 {
3634 if (esign < 0)
3635 goto zero;
3636 else
3637 goto infinite;
3638 }
3639 }
3640if( esign < 0 )
3641 exp = -exp;
3642if( exp > 4932 )
3643 {
3644infinite:
3645 ecleaz(yy);
3646 yy[E] = 0x7fff; /* infinity */
3647 goto aexit;
3648 }
3649if( exp < -4977 )
3650 {
3651zero:
3652 ecleaz(yy);
3653 goto aexit;
3654 }
3655
3656daldone:
3657nexp = exp - nexp;
3658/* Pad trailing zeros to minimize power of 10, per IEEE spec. */
3659while( (nexp > 0) && (yy[2] == 0) )
3660 {
3661 emovz( yy, xt );
3662 eshup1( xt );
3663 eshup1( xt );
3664 eaddm( yy, xt );
3665 eshup1( xt );
3666 if( xt[2] != 0 )
3667 break;
3668 nexp -= 1;
3669 emovz( xt, yy );
3670 }
3671if( (k = enormlz(yy)) > NBITS )
3672 {
3673 ecleaz(yy);
3674 goto aexit;
3675 }
3676lexp = (EXONE - 1 + NBITS) - k;
3677emdnorm( yy, lost, 0, lexp, 64 );
3678/* convert to external format */
3679
3680
3681/* Multiply by 10**nexp. If precision is 64 bits,
3682 * the maximum relative error incurred in forming 10**n
3683 * for 0 <= n <= 324 is 8.2e-20, at 10**180.
3684 * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
3685 * For 0 >= n >= -999, it is -1.55e-19 at 10**-435.
3686 */
3687lexp = yy[E];
3688if( nexp == 0 )
3689 {
3690 k = 0;
3691 goto expdon;
3692 }
3693esign = 1;
3694if( nexp < 0 )
3695 {
3696 nexp = -nexp;
3697 esign = -1;
3698 if( nexp > 4096 )
3699 { /* Punt. Can't handle this without 2 divides. */
3700 emovi( etens[0], tt );
3701 lexp -= tt[E];
3702 k = edivm( tt, yy );
3703 lexp += EXONE;
3704 nexp -= 4096;
3705 }
3706 }
3707p = &etens[NTEN][0];
3708emov( eone, xt );
3709exp = 1;
3710do
3711 {
3712 if( exp & nexp )
3713 emul( p, xt, xt );
3714 p -= NE;
3715 exp = exp + exp;
3716 }
3717while( exp <= MAXP );
3718
3719emovi( xt, tt );
3720if( esign < 0 )
3721 {
3722 lexp -= tt[E];
3723 k = edivm( tt, yy );
3724 lexp += EXONE;
3725 }
3726else
3727 {
3728 lexp += tt[E];
3729 k = emulm( tt, yy );
3730 lexp -= EXONE - 1;
3731 }
3732
3733expdon:
3734
3735/* Round and convert directly to the destination type */
3736if( oprec == 53 )
3737 lexp -= EXONE - 0x3ff;
3738else if( oprec == 24 )
3739 lexp -= EXONE - 0177;
3740#ifdef DEC
3741else if( oprec == 56 )
3742 lexp -= EXONE - 0201;
3743#endif
3744rndprc = oprec;
3745emdnorm( yy, k, 0, lexp, 64 );
3746
3747aexit:
3748
3749rndprc = rndsav;
3750yy[0] = nsign;
3751switch( oprec )
3752 {
3753#ifdef DEC
3754 case 56:
3755 todec( yy, y ); /* see etodec.c */
3756 break;
3757#endif
3758 case 53:
3759 toe53( yy, y );
3760 break;
3761 case 24:
3762 toe24( yy, y );
3763 break;
3764 case 64:
3765 toe64( yy, y );
3766 break;
3767 case 113:
3768 toe113( yy, y );
3769 break;
3770 case NBITS:
3771 emovo( yy, y );
3772 break;
3773 }
3774}
3775
3776
3777
3778/* y = largest integer not greater than x
3779 * (truncated toward minus infinity)
3780 *
3781 * unsigned short x[NE], y[NE]
3782 *
3783 * efloor( x, y );
3784 */
3785static unsigned short bmask[] = {
37860xffff,
37870xfffe,
37880xfffc,
37890xfff8,
37900xfff0,
37910xffe0,
37920xffc0,
37930xff80,
37940xff00,
37950xfe00,
37960xfc00,
37970xf800,
37980xf000,
37990xe000,
38000xc000,
38010x8000,
38020x0000,
3803};
3804
3805void efloor( x, y )
3806unsigned short x[], y[];
3807{
3808register unsigned short *p;
3809int e, expon, i;
3810unsigned short f[NE];
3811
3812emov( x, f ); /* leave in external format */
3813expon = (int )f[NE-1];
3814e = (expon & 0x7fff) - (EXONE - 1);
3815if( e <= 0 )
3816 {
3817 eclear(y);
3818 goto isitneg;
3819 }
3820/* number of bits to clear out */
3821e = NBITS - e;
3822emov( f, y );
3823if( e <= 0 )
3824 return;
3825
3826p = &y[0];
3827while( e >= 16 )
3828 {
3829 *p++ = 0;
3830 e -= 16;
3831 }
3832/* clear the remaining bits */
3833*p &= bmask[e];
3834/* truncate negatives toward minus infinity */
3835isitneg:
3836
3837if( (unsigned short )expon & (unsigned short )0x8000 )
3838 {
3839 for( i=0; i<NE-1; i++ )
3840 {
3841 if( f[i] != y[i] )
3842 {
3843 esub( eone, y, y );
3844 break;
3845 }
3846 }
3847 }
3848}
3849
3850
3851/* unsigned short x[], s[];
3852 * long *exp;
3853 *
3854 * efrexp( x, exp, s );
3855 *
3856 * Returns s and exp such that s * 2**exp = x and .5 <= s < 1.
3857 * For example, 1.1 = 0.55 * 2**1
3858 * Handles denormalized numbers properly using long integer exp.
3859 */
3860void efrexp( x, exp, s )
3861unsigned short x[];
3862long *exp;
3863unsigned short s[];
3864{
3865unsigned short xi[NI];
3866long li;
3867
3868emovi( x, xi );
3869li = (long )((short )xi[1]);
3870
3871if( li == 0 )
3872 {
3873 li -= enormlz( xi );
3874 }
3875xi[1] = 0x3ffe;
3876emovo( xi, s );
3877*exp = li - 0x3ffe;
3878}
3879
3880
3881
3882/* unsigned short x[], y[];
3883 * long pwr2;
3884 *
3885 * eldexp( x, pwr2, y );
3886 *
3887 * Returns y = x * 2**pwr2.
3888 */
3889void eldexp( x, pwr2, y )
3890unsigned short x[];
3891long pwr2;
3892unsigned short y[];
3893{
3894unsigned short xi[NI];
3895long li;
3896int i;
3897
3898emovi( x, xi );
3899li = xi[1];
3900li += pwr2;
3901i = 0;
3902emdnorm( xi, i, i, li, 64 );
3903emovo( xi, y );
3904}
3905
3906
3907/* c = remainder after dividing b by a
3908 * Least significant integer quotient bits left in equot[].
3909 */
3910void eremain( a, b, c )
3911unsigned short a[], b[], c[];
3912{
3913unsigned short den[NI], num[NI];
3914
3915#ifdef NANS
3916if( eisinf(b) || (ecmp(a,ezero) == 0) || eisnan(a) || eisnan(b))
3917 {
3918 enan( c, NBITS );
3919 return;
3920 }
3921#endif
3922if( ecmp(a,ezero) == 0 )
3923 {
3924 mtherr( "eremain", SING );
3925 eclear( c );
3926 return;
3927 }
3928emovi( a, den );
3929emovi( b, num );
3930eiremain( den, num );
3931/* Sign of remainder = sign of quotient */
3932if( a[0] == b[0] )
3933 num[0] = 0;
3934else
3935 num[0] = 0xffff;
3936emovo( num, c );
3937}
3938
3939
3940void eiremain( den, num )
3941unsigned short den[], num[];
3942{
3943long ld, ln;
3944unsigned short j;
3945
3946ld = den[E];
3947ld -= enormlz( den );
3948ln = num[E];
3949ln -= enormlz( num );
3950ecleaz( equot );
3951while( ln >= ld )
3952 {
3953 if( ecmpm(den,num) <= 0 )
3954 {
3955 esubm(den, num);
3956 j = 1;
3957 }
3958 else
3959 {
3960 j = 0;
3961 }
3962 eshup1(equot);
3963 equot[NI-1] |= j;
3964 eshup1(num);
3965 ln -= 1;
3966 }
3967emdnorm( num, 0, 0, ln, 0 );
3968}
3969
3970/* NaN bit patterns
3971 */
3972#ifdef MIEEE
3973unsigned short nan113[8] = {
3974 0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
3975unsigned short nan64[6] = {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
3976unsigned short nan53[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
3977unsigned short nan24[2] = {0x7fff, 0xffff};
3978#endif
3979
3980#ifdef IBMPC
3981unsigned short nan113[8] = {0, 0, 0, 0, 0, 0, 0xc000, 0xffff};
3982unsigned short nan64[6] = {0, 0, 0, 0xc000, 0xffff, 0};
3983unsigned short nan53[4] = {0, 0, 0, 0xfff8};
3984unsigned short nan24[2] = {0, 0xffc0};
3985#endif
3986
3987
3988void enan (nan, size)
3989unsigned short *nan;
3990int size;
3991{
3992int i, n;
3993unsigned short *p;
3994
3995switch( size )
3996 {
3997#ifndef DEC
3998 case 113:
3999 n = 8;
4000 p = nan113;
4001 break;
4002
4003 case 64:
4004 n = 6;
4005 p = nan64;
4006 break;
4007
4008 case 53:
4009 n = 4;
4010 p = nan53;
4011 break;
4012
4013 case 24:
4014 n = 2;
4015 p = nan24;
4016 break;
4017
4018 case NBITS:
4019 for( i=0; i<NE-2; i++ )
4020 *nan++ = 0;
4021 *nan++ = 0xc000;
4022 *nan++ = 0x7fff;
4023 return;
4024
4025 case NI*16:
4026 *nan++ = 0;
4027 *nan++ = 0x7fff;
4028 *nan++ = 0;
4029 *nan++ = 0xc000;
4030 for( i=4; i<NI; i++ )
4031 *nan++ = 0;
4032 return;
4033#endif
4034 default:
4035 mtherr( "enan", DOMAIN );
4036 return;
4037 }
4038for (i=0; i < n; i++)
4039 *nan++ = *p++;
4040}
4041
4042
4043
4044/* Longhand square root. */
4045
4046static int esqinited = 0;
4047static unsigned short sqrndbit[NI];
4048
4049void esqrt( x, y )
4050short *x, *y;
4051{
4052unsigned short temp[NI], num[NI], sq[NI], xx[NI];
4053int i, j, k, n, nlups;
4054long m, exp;
4055
4056if( esqinited == 0 )
4057 {
4058 ecleaz( sqrndbit );
4059 sqrndbit[NI-2] = 1;
4060 esqinited = 1;
4061 }
4062/* Check for arg <= 0 */
4063i = ecmp( x, ezero );
4064if( i <= 0 )
4065 {
4066#ifdef NANS
4067 if (i == -2)
4068 {
4069 enan (y, NBITS);
4070 return;
4071 }
4072#endif
4073 eclear(y);
4074 if( i < 0 )
4075 mtherr( "esqrt", DOMAIN );
4076 return;
4077 }
4078
4079#ifdef INFINITY
4080if( eisinf(x) )
4081 {
4082 eclear(y);
4083 einfin(y);
4084 return;
4085 }
4086#endif
4087/* Bring in the arg and renormalize if it is denormal. */
4088emovi( x, xx );
4089m = (long )xx[1]; /* local long word exponent */
4090if( m == 0 )
4091 m -= enormlz( xx );
4092
4093/* Divide exponent by 2 */
4094m -= 0x3ffe;
4095exp = (unsigned short )( (m / 2) + 0x3ffe );
4096
4097/* Adjust if exponent odd */
4098if( (m & 1) != 0 )
4099 {
4100 if( m > 0 )
4101 exp += 1;
4102 eshdn1( xx );
4103 }
4104
4105ecleaz( sq );
4106ecleaz( num );
4107n = 8; /* get 8 bits of result per inner loop */
4108nlups = rndprc;
4109j = 0;
4110
4111while( nlups > 0 )
4112 {
4113/* bring in next word of arg */
4114 if( j < NE )
4115 num[NI-1] = xx[j+3];
4116/* Do additional bit on last outer loop, for roundoff. */
4117 if( nlups <= 8 )
4118 n = nlups + 1;
4119 for( i=0; i<n; i++ )
4120 {
4121/* Next 2 bits of arg */
4122 eshup1( num );
4123 eshup1( num );
4124/* Shift up answer */
4125 eshup1( sq );
4126/* Make trial divisor */
4127 for( k=0; k<NI; k++ )
4128 temp[k] = sq[k];
4129 eshup1( temp );
4130 eaddm( sqrndbit, temp );
4131/* Subtract and insert answer bit if it goes in */
4132 if( ecmpm( temp, num ) <= 0 )
4133 {
4134 esubm( temp, num );
4135 sq[NI-2] |= 1;
4136 }
4137 }
4138 nlups -= n;
4139 j += 1;
4140 }
4141
4142/* Adjust for extra, roundoff loop done. */
4143exp += (NBITS - 1) - rndprc;
4144
4145/* Sticky bit = 1 if the remainder is nonzero. */
4146k = 0;
4147for( i=3; i<NI; i++ )
4148 k |= (int )num[i];
4149
4150/* Renormalize and round off. */
4151emdnorm( sq, k, 0, exp, 64 );
4152emovo( sq, y );
4153}
diff --git a/src/regress/lib/libc/cephes/ieetst.c b/src/regress/lib/libc/cephes/ieetst.c
new file mode 100644
index 0000000000..fdce864f40
--- /dev/null
+++ b/src/regress/lib/libc/cephes/ieetst.c
@@ -0,0 +1,875 @@
1/* $OpenBSD: ieetst.c,v 1.1 2011/07/02 18:11:01 martynas Exp $ */
2
3/*
4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19/* Floating point to ASCII input and output string test program.
20 *
21 * Numbers in the native machine data structure are converted
22 * to e type, then to and from decimal ASCII strings. Native
23 * printf() and scanf() functions are also used to produce
24 * and read strings. The resulting e type binary values
25 * are compared, with diagnostic printouts of any discrepancies.
26 *
27 * Steve Moshier, 16 Dec 88
28 * last revision: 16 May 92
29 */
30
31#include <stdio.h>
32
33#include "mconf.h"
34#include "ehead.h"
35
36/* Include tests of 80-bit long double precision: */
37#define LDOUBLE 0
38/* Abort subtest after getting this many errors: */
39#define MAXERR 5
40/* Number of random arguments to try (set as large as you have
41 * patience for): */
42#define NRAND 100
43/* Perform internal consistency test: */
44#define CHKINTERNAL 0
45
46static unsigned short fullp[NE], rounded[NE];
47float prec24, sprec24, ssprec24;
48double prec53, sprec53, ssprec53;
49#if LDOUBLE
50long double prec64, sprec64, ssprec64;
51#endif
52
53static unsigned short rprint[NE], rscan[NE];
54static unsigned short q1[NE], q2[NE], q5[NE];
55static unsigned short e1[NE], e2[NE], e3[NE];
56static double d1, d2;
57static int errprint = 0;
58static int errscan = 0;
59static int identerr = 0;
60static int errtot = 0;
61static int count = 0;
62static char str0[80], str1[80], str2[80], str3[80];
63static unsigned short eten[NE], maxm[NE];
64
65int m, n, k2, mprec, SPREC;
66
67char *Ten = "10.0";
68char tformat[10];
69char *format24 = "%.8e";
70#ifdef DEC
71char *format53 = "%.17e";
72#else
73char *format53 = "%.16e";
74#endif
75char *fformat24 = "%e";
76char *fformat53 = "%le";
77char *pct = "%";
78char *quo = "\042";
79#if LDOUBLE
80char *format64 = "%.20Le";
81char *fformat64 = "%Le";
82#endif
83char *format;
84char *fformat;
85char *toomany = "Too many errors; aborting this test.\n";
86
87static int mnrflag;
88static int etrflag;
89void chkit(), printerr(), mnrand(), etrand(), shownoncrit();
90void chkid(), pvec();
91
92main()
93{
94int i, iprec, retval = 0;
95
96printf( "Steve Moshier's printf/scanf tester, version 0.2.\n\n" );
97#ifdef DEC
98 /* DEC PDP-11/VAX single precision not yet implemented */
99for( iprec = 1; iprec<2; iprec++ )
100#else
101for( iprec = 0; iprec<3; iprec++ )
102/*for( iprec = 2; iprec<3; iprec++ )*/
103#endif
104 {
105 errscan = 0;
106 identerr = 0;
107 errprint = 0;
108 eclear( rprint );
109 eclear( rscan );
110
111switch( iprec )
112 {
113 case 0:
114 SPREC = 8; /* # digits after the decimal point */
115 mprec = 24; /* # bits in the significand */
116 m = 9; /* max # decimal digits for correct rounding */
117 n = 13; /* max power of ten for correct rounding */
118 k2 = -125; /* underflow beyond 2^-k2 */
119 format = format24; /* printf format string */
120 fformat = fformat24; /* scanf format string */
121 mnrflag = 1; /* sets interval for random numbers */
122 etrflag = 1;
123 printf( "Testing FLOAT precision.\n" );
124 break;
125
126 case 1:
127#ifdef DEC
128 SPREC = 17;
129 mprec = 56;
130 m = 17;
131 n = 27;
132 k2 = -125;
133 format = format53;
134 fformat = fformat53;
135 mnrflag = 2;
136 etrflag = 1;
137 printf( "Testing DEC DOUBLE precision.\n" );
138 break;
139#else
140 SPREC = 16;
141 mprec = 53;
142 m = 17;
143 n = 27;
144 k2 = -1021;
145 format = format53;
146 fformat = fformat53;
147 mnrflag = 2;
148 etrflag = 2;
149 printf( "Testing DOUBLE precision.\n" );
150 break;
151#endif
152 case 2:
153#if LDOUBLE
154 SPREC = 20;
155 mprec = 64;
156 m = 20;
157 n = 34;
158 k2 = -16382;
159 format = format64;
160 fformat = fformat64;
161 mnrflag = 3;
162 etrflag = 3;
163 printf( "Testing LONG DOUBLE precision.\n" );
164 break;
165#else
166 goto nodenorm;
167#endif
168 }
169
170 asctoe( Ten, eten );
171/* 10^m - 1 */
172 d2 = m;
173 e53toe( &d2, e1 );
174 epow( eten, e1, maxm );
175 esub( eone, maxm, maxm );
176
177/* test 1 */
178 printf( "1. Checking 10^n - 1 for n = %d to %d.\n", -m, m );
179 emov( eone, q5 );
180 for( count=0; count<=m; count++ )
181 {
182 esub( eone, q5, fullp );
183 chkit( 1 );
184 ediv( q5, eone, q2 );
185 esub( eone, q2, fullp );
186 chkit( 1 );
187 emul( eten, q5, q5 );
188 if( errtot >= MAXERR )
189 {
190 printf( "%s", toomany );
191 goto end1;
192 }
193 }
194end1:
195 printerr();
196
197
198/* test 2 */
199 printf( "2. Checking powers of 10 from 10^-%d to 10^%d.\n", n, n );
200 emov( eone, q5 );
201 for( count=0; count<=n; count++ )
202 {
203 emov( q5, fullp );
204 chkit( 2 );
205 ediv( q5, eone, fullp );
206 chkit( 2 );
207 emul( eten, q5, q5 );
208 if( errtot >= MAXERR )
209 {
210 printf( "%s", toomany );
211 goto end2;
212 }
213 }
214end2:
215 printerr();
216
217/* test 3 */
218 printf( "3. Checking (10^%d-1)*10^n from n = -%d to %d.\n", m, n, n );
219 emov( eone, q5 );
220 for( count= -n; count<=n; count++ )
221 {
222 emul( maxm, q5, fullp );
223 chkit( 3 );
224 emov( q5, fullp );
225 ediv( fullp, eone, fullp );
226 emul( maxm, fullp, fullp );
227 chkit( 3 );
228 emul( eten, q5, q5 );
229 if( errtot >= MAXERR )
230 {
231 printf( "%s", toomany );
232 goto end3;
233 }
234 }
235end3:
236 printerr();
237
238
239
240/* test 4 */
241 printf( "4. Checking powers of 2 from 2^-24 to 2^+56.\n" );
242 d1 = -24.0;
243 e53toe( &d1, q1 );
244 epow( etwo, q1, q5 );
245
246 for( count = -24; count <= 56; count++ )
247 {
248 emov( q5, fullp );
249 chkit( 4 );
250 emul( etwo, q5, q5 );
251 if( errtot >= MAXERR )
252 {
253 printf( "%s", toomany );
254 goto end4;
255 }
256 }
257end4:
258 printerr();
259
260
261/* test 5 */
262 printf( "5. Checking 2^n - 1 for n = 0 to %d.\n", mprec );
263 emov( eone, q5 );
264 for( count=0; count<=mprec; count++ )
265 {
266 esub( eone, q5, fullp );
267 chkit( 5 );
268 emul( etwo, q5, q5 );
269 if( errtot >= MAXERR )
270 {
271 printf( "%s", toomany );
272 goto end5;
273 }
274 }
275end5:
276 printerr();
277
278/* test 6 */
279 printf( "6. Checking 2^n + 1 for n = 0 to %d.\n", mprec );
280 emov( eone, q5 );
281 for( count=0; count<=mprec; count++ )
282 {
283 eadd( eone, q5, fullp );
284 chkit( 6 );
285 emul( etwo, q5, q5 );
286 if( errtot >= MAXERR )
287 {
288 printf( "%s", toomany );
289 goto end6;
290 }
291 }
292end6:
293 printerr();
294
295/* test 7 */
296 printf(
297 "7. Checking %d values M * 10^N with random integer M and N,\n",
298 NRAND );
299 printf(" 1 <= M <= 10^%d - 1 and -%d <= N <= +%d.\n", m, n, n );
300 for( i=0; i<NRAND; i++ )
301 {
302 mnrand( fullp );
303 chkit( 7 );
304 if( errtot >= MAXERR )
305 {
306 printf( "%s", toomany );
307 goto end7;
308 }
309 }
310end7:
311 printerr();
312
313/* test 8 */
314 printf("8. Checking critical rounding cases.\n" );
315 for( i=0; i<20; i++ )
316 {
317 mnrand( fullp );
318 eabs( fullp );
319 if( ecmp( fullp, eone ) < 0 )
320 ediv( fullp, eone, fullp );
321 efloor( fullp, fullp );
322 eadd( ehalf, fullp, fullp );
323 chkit( 8 );
324 if( errtot >= MAXERR )
325 {
326 printf( "%s", toomany );
327 goto end8;
328 }
329 }
330end8:
331 printerr();
332
333
334
335/* test 9 */
336 printf("9. Testing on %d random non-denormal values.\n", NRAND );
337 for( i=0; i<NRAND; i++ )
338 {
339 etrand( fullp );
340 chkit( 9 );
341 }
342 printerr();
343 shownoncrit();
344
345/* test 10 */
346#if 0
347 printf(
348 "Do you want to check denormal numbers in this precision ? (y/n) " );
349 gets( str0 );
350 if( str0[0] != 'y' )
351 goto nodenorm;
352#endif
353
354 printf( "10. Checking denormal numbers.\n" );
355
356/* Form 2^-starting power */
357 d1 = k2;
358 e53toe( &d1, q1 );
359 epow( etwo, q1, e1 );
360
361/* Find 2^-mprec less than starting power */
362 d1 = -mprec + 4;
363 e53toe( &d1, q1 );
364 epow( etwo, q1, e3 );
365 emul( e1, e3, e3 );
366 emov( e3, e2 );
367 ediv( etwo, e2, e2 );
368
369 while( ecmp(e1,e2) != 0 )
370 {
371 eadd( e1, e2, fullp );
372 switch( mprec )
373 {
374#if LDOUBLE
375 case 64:
376 etoe64( e1, &sprec64 );
377 e64toe( &sprec64, q1 );
378 etoe64( fullp, &prec64 );
379 e64toe( &prec64, q2 );
380 break;
381#endif
382#ifdef DEC
383 case 56:
384#endif
385 case 53:
386 etoe53( e1, &sprec53 );
387 e53toe( &sprec53, q1 );
388 etoe53( fullp, &prec53 );
389 e53toe( &prec53, q2 );
390 break;
391
392 case 24:
393 etoe24( e1, &sprec24 );
394 e24toe( &sprec24, q1 );
395 etoe24( fullp, &prec24 );
396 e24toe( &prec24, q2 );
397 break;
398 }
399 if( ecmp( q2, ezero ) == 0 )
400 goto maxden;
401 chkit(10);
402 if( ecmp(q1,q2) == 0 )
403 {
404 ediv( etwo, e1, e1 );
405 emov( e3, e2 );
406 }
407 if( errtot >= MAXERR )
408 {
409 printf( "%s", toomany );
410 goto maxden;
411 }
412 ediv( etwo, e2, e2 );
413 }
414maxden:
415 printerr();
416nodenorm:
417 printf( "\n" );
418 retval |= errscan | identerr | errprint;
419 } /* loop on precision */
420printf( "End of test.\n" );
421return (retval);
422}
423
424#if CHKINTERNAL
425long double xprec64;
426double xprec53;
427float xprec24;
428
429/* Check binary -> printf -> scanf -> binary identity
430 * of internal routines
431 */
432void chkinternal( ref, tst, string )
433unsigned short ref[], tst[];
434char *string;
435{
436
437if( ecmp(ref,tst) != 0 )
438 {
439 printf( "internal identity compare error!\n" );
440 chkid( ref, tst, string );
441 }
442}
443#endif
444
445
446/* Check binary -> printf -> scanf -> binary identity
447 */
448void chkid( print, scan, string )
449unsigned short print[], scan[];
450char *string;
451{
452/* Test printf-scanf identity */
453if( ecmp( print, scan ) != 0 )
454 {
455 pvec( print, NE );
456 printf( " ->printf-> %s ->scanf->\n", string );
457 pvec( scan, NE );
458 printf( " is not an identity.\n" );
459 ++identerr;
460 }
461}
462
463
464/* Check scanf result
465 */
466void chkscan( ref, tst, string )
467unsigned short ref[], tst[];
468char *string;
469{
470/* Test scanf() */
471if( ecmp( ref, tst ) != 0 )
472 {
473 printf( "scanf(%s) -> ", string );
474 pvec( tst, NE );
475 printf( "\n should be " );
476 pvec( ref, NE );
477 printf( ".\n" );
478 ++errscan;
479 ++errtot;
480 }
481}
482
483
484/* Test printf() result
485 */
486void chkprint( ref, tst, string )
487unsigned short ref[], tst[];
488char *string;
489{
490if( ecmp(ref, tst) != 0 )
491 {
492 printf( "printf( ");
493 pvec( ref, NE );
494 printf( ") -> %s\n", string );
495 printf( " = " );
496 pvec( tst, NE );
497 printf( ".\n" );
498 ++errprint;
499 ++errtot;
500 }
501}
502
503
504/* Print array of n 16-bit shorts
505 */
506void pvec( x, n )
507unsigned short x[];
508int n;
509{
510int i;
511
512for( i=0; i<n; i++ )
513 {
514 printf( "%04x ", x[i] );
515 }
516}
517
518/* Measure worst case printf rounding error
519 */
520void cmpprint( ref, tst )
521unsigned short ref[], tst[];
522{
523unsigned short e[NE];
524
525if( ecmp( ref, ezero ) != 0 )
526 {
527 esub( ref, tst, e );
528 ediv( ref, e, e );
529 eabs( e );
530 if( ecmp( e, rprint ) > 0 )
531 emov( e, rprint );
532 }
533}
534
535/* Measure worst case scanf rounding error
536 */
537void cmpscan( ref, tst )
538unsigned short ref[], tst[];
539{
540unsigned short er[NE];
541
542if( ecmp( ref, ezero ) != 0 )
543 {
544 esub( ref, tst, er );
545 ediv( ref, er, er );
546 eabs( er );
547 if( ecmp( er, rscan ) > 0 )
548 emov( er, rscan );
549 if( ecmp( er, ehalf ) > 0 )
550 {
551 etoasc( tst, str1, 21 );
552 printf( "Bad error: scanf(%s) = %s !\n", str0, str1 );
553 }
554 }
555}
556
557/* Check rounded-down decimal string output of printf
558 */
559void cmptrunc( ref, tst )
560unsigned short ref[], tst[];
561{
562if( ecmp( ref, tst ) != 0 )
563 {
564 printf( "printf(%s%s%s, %s) -> %s\n", quo, tformat, quo, str1, str2 );
565 printf( "should be %s .\n", str3 );
566 errprint += 1;
567 }
568}
569
570
571void shownoncrit()
572{
573
574etoasc( rprint, str0, 3 );
575printf( "Maximum relative printf error found = %s .\n", str0 );
576etoasc( rscan, str0, 3 );
577printf( "Maximum relative scanf error found = %s .\n", str0 );
578}
579
580
581
582/* Produce arguments and call comparison subroutines.
583 */
584void chkit( testno )
585int testno;
586{
587unsigned short t[NE], u[NE], v[NE];
588int j;
589
590switch( mprec )
591 {
592#if LDOUBLE
593 case 64:
594 etoe64( fullp, &prec64 );
595 e64toe( &prec64, rounded );
596#if CHKINTERNAL
597 e64toasc( &prec64, str1, SPREC );
598 asctoe64( str1, &xprec64 );
599 e64toe( &xprec64, t );
600 chkinternal( rounded, t, str1 );
601#endif
602/* check printf and scanf */
603 sprintf( str2, format, prec64 );
604 sscanf( str2, fformat, &sprec64 );
605 e64toe( &sprec64, u );
606 chkid( rounded, u, str2 );
607 asctoe64( str2, &ssprec64 );
608 e64toe( &ssprec64, v );
609 chkscan( v, u, str2 );
610 chkprint( rounded, v, str2 );
611 if( testno < 8 )
612 break;
613/* rounding error measurement */
614 etoasc( fullp, str0, 24 );
615 etoe64( fullp, &ssprec64 );
616 e64toe( &ssprec64, u );
617 sprintf( str2, format, ssprec64 );
618 asctoe( str2, t );
619 cmpprint( u, t );
620 sscanf( str0, fformat, &sprec64 );
621 e64toe( &sprec64, t );
622 cmpscan( fullp, t );
623 if( testno < 8 )
624 break;
625/* strings rounded to less than maximum precision */
626 e64toasc( &ssprec64, str1, 24 );
627 for( j=SPREC-1; j>0; j-- )
628 {
629 e64toasc( &ssprec64, str3, j );
630 asctoe( str3, v );
631 sprintf( tformat, "%s.%dLe", pct, j );
632 sprintf( str2, tformat, ssprec64 );
633 asctoe( str2, t );
634 cmptrunc( v, t );
635 }
636 break;
637#endif
638#ifdef DEC
639 case 56:
640#endif
641 case 53:
642 etoe53( fullp, &prec53 );
643 e53toe( &prec53, rounded );
644#if CHKINTERNAL
645 e53toasc( &prec53, str1, SPREC );
646 asctoe53( str1, &xprec53 );
647 e53toe( &xprec53, t );
648 chkinternal( rounded, t, str1 );
649#endif
650 sprintf( str2, format, prec53 );
651 sscanf( str2, fformat, &sprec53 );
652 e53toe( &sprec53, u );
653 chkid( rounded, u, str2 );
654 asctoe53( str2, &ssprec53 );
655 e53toe( &ssprec53, v );
656 chkscan( v, u, str2 );
657 chkprint( rounded, v, str2 );
658 if( testno < 8 )
659 break;
660/* rounding error measurement */
661 etoasc( fullp, str0, 24 );
662 etoe53( fullp, &ssprec53 );
663 e53toe( &ssprec53, u );
664 sprintf( str2, format, ssprec53 );
665 asctoe( str2, t );
666 cmpprint( u, t );
667 sscanf( str0, fformat, &sprec53 );
668 e53toe( &sprec53, t );
669 cmpscan( fullp, t );
670 if( testno < 8 )
671 break;
672 e53toasc( &ssprec53, str1, 24 );
673 for( j=SPREC-1; j>0; j-- )
674 {
675 e53toasc( &ssprec53, str3, j );
676 asctoe( str3, v );
677 sprintf( tformat, "%s.%de", pct, j );
678 sprintf( str2, tformat, ssprec53 );
679 asctoe( str2, t );
680 cmptrunc( v, t );
681 }
682 break;
683
684 case 24:
685 etoe24( fullp, &prec24 );
686 e24toe( &prec24, rounded );
687#if CHKINTERNAL
688 e24toasc( &prec24, str1, SPREC );
689 asctoe24( str1, &xprec24 );
690 e24toe( &xprec24, t );
691 chkinternal( rounded, t, str1 );
692#endif
693 sprintf( str2, format, prec24 );
694 sscanf( str2, fformat, &sprec24 );
695 e24toe( &sprec24, u );
696 chkid( rounded, u, str2 );
697 asctoe24( str2, &ssprec24 );
698 e24toe( &ssprec24, v );
699 chkscan( v, u, str2 );
700 chkprint( rounded, v, str2 );
701 if( testno < 8 )
702 break;
703/* rounding error measurement */
704 etoasc( fullp, str0, 24 );
705 etoe24( fullp, &ssprec24 );
706 e24toe( &ssprec24, u );
707 sprintf( str2, format, ssprec24 );
708 asctoe( str2, t );
709 cmpprint( u, t );
710 sscanf( str0, fformat, &sprec24 );
711 e24toe( &sprec24, t );
712 cmpscan( fullp, t );
713/*
714 if( testno < 8 )
715 break;
716*/
717 e24toasc( &ssprec24, str1, 24 );
718 for( j=SPREC-1; j>0; j-- )
719 {
720 e24toasc( &ssprec24, str3, j );
721 asctoe( str3, v );
722 sprintf( tformat, "%s.%de", pct, j );
723 sprintf( str2, tformat, ssprec24 );
724 asctoe( str2, t );
725 cmptrunc( v, t );
726 }
727 break;
728 }
729}
730
731
732void printerr()
733{
734if( (errscan == 0) && (identerr == 0) && (errprint == 0) )
735 printf( "No errors found.\n" );
736else
737 {
738 printf( "%d binary -> decimal errors found.\n", errprint );
739 printf( "%d decimal -> binary errors found.\n", errscan );
740 }
741errscan = 0; /* reset for next test */
742identerr = 0;
743errprint = 0;
744errtot = 0;
745}
746
747
748/* Random number generator
749 * in the range M * 10^N, where 1 <= M <= 10^17 - 1
750 * and -27 <= N <= +27. Test values of M are logarithmically distributed
751 * random integers; test values of N are uniformly distributed random integers.
752 */
753
754static char *fwidth = "1.036163291797320557783096e1"; /* log(sqrt(10^9-1)) */
755static char *dwidth = "1.957197329044938830915E1"; /* log(sqrt(10^17-1)) */
756static char *ldwidth = "2.302585092994045684017491e1"; /* log(sqrt(10^20-1)) */
757
758static char *a13 = "13.0";
759static char *a27 = "27.0";
760static char *a34 = "34.0";
761static char *a10m13 = "1.0e-13";
762static unsigned short LOW[ NE ], WIDTH[NE], e27[NE], e10m13[NE];
763
764
765void mnrand( erand )
766unsigned short erand[];
767{
768unsigned short ea[NE], em[NE], en[NE], ex[NE];
769double x, a;
770
771if( mnrflag )
772 {
773 if( mnrflag == 3 )
774 {
775 asctoe( ldwidth, WIDTH );
776 asctoe( a34, e27 );
777 }
778 if( mnrflag == 2 )
779 {
780 asctoe( dwidth, WIDTH );
781 asctoe( a27, e27 );
782 }
783 if( mnrflag == 1 )
784 {
785 asctoe( fwidth, WIDTH );
786 asctoe( a13, e27 );
787 }
788 asctoe( a10m13, e10m13 );
789 mnrflag = 0;
790 }
791drand( &x );
792e53toe( &x, ex ); /* x = WIDTH * ( x - 1.0 ) + LOW; */
793esub( eone, ex, ex );
794emul( WIDTH, ex, ex );
795eexp( ex, ex ); /* x = exp(x); */
796
797drand( &a );
798e53toe( &a, ea );
799emul( ea, ex, ea ); /* a = 1.0e-13 * x * a; */
800emul( e10m13, ea, ea );
801eabs( ea );
802eadd( ea, ex, ex ); /* add fuzz */
803emul( ex, ex, ex ); /* square it, to get range to 10^17 - 1 */
804efloor( ex, em ); /* this is M */
805
806/* Random power of 10 */
807drand( &a );
808e53toe( &a, ex );
809esub( eone, ex, ex ); /* y3 = 54.0 * ( y3 - 1.0 ) + 0.5; */
810emul( e27, ex, ex );
811eadd( ex, ex, ex );
812eadd( ehalf, ex, ex );
813efloor( ex, ex ); /* y3 = floor( y3 ) - 27.0; */
814esub( e27, ex, en ); /* this is N */
815epow( eten, en, ex );
816emul( ex, em, erand );
817}
818
819/* -ln 2^16382 */
820char *ldemin = "-1.1355137111933024058873097E4";
821char *ldewid = "2.2710274223866048117746193E4";
822/* -ln 2^1022 */
823char *demin = "-7.0839641853226410622441123E2";
824char *dewid = "1.4167928370645282124488225E3";
825/* -ln 2^125 */
826char *femin = "-8.6643397569993163677154015E1";
827char *fewid = "1.7328679513998632735430803E2";
828
829void etrand( erand )
830unsigned short erand[];
831{
832unsigned short ea[NE], ex[NE];
833double x, a;
834
835if( etrflag )
836 {
837 if( etrflag == 3 )
838 {
839 asctoe( ldemin, LOW );
840 asctoe( ldewid, WIDTH );
841 asctoe( a34, e27 );
842 }
843 if( etrflag == 2 )
844 {
845 asctoe( demin, LOW );
846 asctoe( dewid, WIDTH );
847 asctoe( a27, e27 );
848 }
849 if( etrflag == 1 )
850 {
851 asctoe( femin, LOW );
852 asctoe( fewid, WIDTH );
853 asctoe( a13, e27 );
854 }
855 asctoe( a10m13, e10m13 );
856 etrflag = 0;
857 }
858drand( &x );
859e53toe( &x, ex ); /* x = WIDTH * ( x - 1.0 ) + LOW; */
860esub( eone, ex, ex );
861emul( WIDTH, ex, ex );
862eadd( LOW, ex, ex );
863eexp( ex, ex ); /* x = exp(x); */
864
865/* add fuzz
866 */
867drand( &a );
868e53toe( &a, ea );
869emul( ea, ex, ea ); /* a = 1.0e-13 * x * a; */
870emul( e10m13, ea, ea );
871if( ecmp( ex, ezero ) > 0 )
872 eneg( ea );
873eadd( ea, ex, erand );
874}
875
diff --git a/src/regress/lib/libc/cephes/mconf.h b/src/regress/lib/libc/cephes/mconf.h
new file mode 100644
index 0000000000..a92bd3ab64
--- /dev/null
+++ b/src/regress/lib/libc/cephes/mconf.h
@@ -0,0 +1,187 @@
1/* $OpenBSD: mconf.h,v 1.1 2011/07/02 18:11:01 martynas Exp $ */
2
3/*
4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19/* mconf.h
20 *
21 * Common include file for math routines
22 *
23 *
24 *
25 * SYNOPSIS:
26 *
27 * #include "mconf.h"
28 *
29 *
30 *
31 * DESCRIPTION:
32 *
33 * This file contains definitions for error codes that are
34 * passed to the common error handling routine mtherr()
35 * (which see).
36 *
37 * The file also includes a conditional assembly definition
38 * for the type of computer arithmetic (IEEE, DEC, Motorola
39 * IEEE, or UNKnown).
40 *
41 * For Digital Equipment PDP-11 and VAX computers, certain
42 * IBM systems, and others that use numbers with a 56-bit
43 * significand, the symbol DEC should be defined. In this
44 * mode, most floating point constants are given as arrays
45 * of octal integers to eliminate decimal to binary conversion
46 * errors that might be introduced by the compiler.
47 *
48 * For little-endian computers, such as IBM PC, that follow the
49 * IEEE Standard for Binary Floating Point Arithmetic (ANSI/IEEE
50 * Std 754-1985), the symbol IBMPC should be defined. These
51 * numbers have 53-bit significands. In this mode, constants
52 * are provided as arrays of hexadecimal 16 bit integers.
53 *
54 * Big-endian IEEE format is denoted MIEEE. On some RISC
55 * systems such as Sun SPARC, double precision constants
56 * must be stored on 8-byte address boundaries. Since integer
57 * arrays may be aligned differently, the MIEEE configuration
58 * may fail on such machines.
59 *
60 * To accommodate other types of computer arithmetic, all
61 * constants are also provided in a normal decimal radix
62 * which one can hope are correctly converted to a suitable
63 * format by the available C language compiler. To invoke
64 * this mode, define the symbol UNK.
65 *
66 * An important difference among these modes is a predefined
67 * set of machine arithmetic constants for each. The numbers
68 * MACHEP (the machine roundoff error), MAXNUM (largest number
69 * represented), and several other parameters are preset by
70 * the configuration symbol. Check the file const.c to
71 * ensure that these values are correct for your computer.
72 *
73 * Configurations NANS, INFINITIES, MINUSZERO, and DENORMAL
74 * may fail on many systems. Verify that they are supposed
75 * to work on your computer.
76 */
77
78#include <sys/types.h>
79#include <sys/endian.h>
80
81/* Constant definitions for math error conditions
82 */
83
84#define DOMAIN 1 /* argument domain error */
85#define SING 2 /* argument singularity */
86#define OVERFLOW 3 /* overflow range error */
87#define UNDERFLOW 4 /* underflow range error */
88#define TLOSS 5 /* total loss of precision */
89#define PLOSS 6 /* partial loss of precision */
90
91#define EDOM 33
92#define ERANGE 34
93
94/* Complex numeral. */
95typedef struct
96 {
97 double r;
98 double i;
99 } cmplx;
100
101/* Long double complex numeral. */
102typedef struct
103 {
104 double r;
105 double i;
106 } cmplxl;
107
108/* Type of computer arithmetic */
109
110/* PDP-11, Pro350, VAX:
111 */
112#ifdef __vax__
113#define DEC 1
114#endif /* __vax__ */
115
116/* Intel IEEE, low order words come first:
117 */
118/* #define IBMPC 1 */
119
120/* Motorola IEEE, high order words come first
121 * (Sun 680x0 workstation):
122 */
123/* #define MIEEE 1 */
124
125/* UNKnown arithmetic, invokes coefficients given in
126 * normal decimal format. Beware of range boundary
127 * problems (MACHEP, MAXLOG, etc. in const.c) and
128 * roundoff problems in pow.c:
129 * (Sun SPARCstation)
130 */
131#ifndef __vax__
132#define UNK 1
133#endif /* !__vax__ */
134
135/* If you define UNK, then be sure to set BIGENDIAN properly. */
136#if BYTE_ORDER == BIG_ENDIAN
137#define BIGENDIAN 1
138#endif /* BYTE_ORDER == BIG_ENDIAN */
139
140/* Define this `volatile' if your compiler thinks
141 * that floating point arithmetic obeys the associative
142 * and distributive laws. It will defeat some optimizations
143 * (but probably not enough of them).
144 *
145 * #define VOLATILE volatile
146 */
147#define VOLATILE
148
149/* For 12-byte long doubles on an i386, pad a 16-bit short 0
150 * to the end of real constants initialized by integer arrays.
151 *
152 * #define XPD 0,
153 *
154 * Otherwise, the type is 10 bytes long and XPD should be
155 * defined blank (e.g., Microsoft C).
156 *
157 * #define XPD
158 */
159#define XPD 0,
160
161/* Define to support tiny denormal numbers, else undefine. */
162#ifndef __vax__
163#define DENORMAL 1
164#endif /* !__vax__ */
165
166/* Define to ask for infinity support, else undefine. */
167#ifndef __vax__
168#define INFINITIES 1
169#endif /* !__vax__ */
170
171/* Define to ask for support of numbers that are Not-a-Number,
172 else undefine. This may automatically define INFINITIES in some files. */
173#ifndef __vax__
174#define NANS 1
175#endif /* !__vax__ */
176
177/* Define to distinguish between -0.0 and +0.0. */
178#define MINUSZERO 1
179
180/* Define 1 for ANSI C atan2() function
181 See atan.c and clog.c. */
182#define ANSIC 1
183
184int mtherr();
185
186/* Variable for error reporting. See mtherr.c. */
187extern int merror;
diff --git a/src/regress/lib/libc/cephes/mtherr.c b/src/regress/lib/libc/cephes/mtherr.c
new file mode 100644
index 0000000000..9a47a198bd
--- /dev/null
+++ b/src/regress/lib/libc/cephes/mtherr.c
@@ -0,0 +1,114 @@
1/* $OpenBSD: mtherr.c,v 1.1 2011/07/02 18:11:01 martynas Exp $ */
2
3/*
4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19/* mtherr.c
20 *
21 * Library common error handling routine
22 *
23 *
24 *
25 * SYNOPSIS:
26 *
27 * char *fctnam;
28 * int code;
29 * int mtherr();
30 *
31 * mtherr( fctnam, code );
32 *
33 *
34 *
35 * DESCRIPTION:
36 *
37 * This routine may be called to report one of the following
38 * error conditions (in the include file mconf.h).
39 *
40 * Mnemonic Value Significance
41 *
42 * DOMAIN 1 argument domain error
43 * SING 2 function singularity
44 * OVERFLOW 3 overflow range error
45 * UNDERFLOW 4 underflow range error
46 * TLOSS 5 total loss of precision
47 * PLOSS 6 partial loss of precision
48 * EDOM 33 Unix domain error code
49 * ERANGE 34 Unix range error code
50 *
51 * The default version of the file prints the function name,
52 * passed to it by the pointer fctnam, followed by the
53 * error condition. The display is directed to the standard
54 * output device. The routine then returns to the calling
55 * program. Users may wish to modify the program to abort by
56 * calling exit() under severe error conditions such as domain
57 * errors.
58 *
59 * Since all error conditions pass control to this function,
60 * the display may be easily changed, eliminated, or directed
61 * to an error logging device.
62 *
63 * SEE ALSO:
64 *
65 * mconf.h
66 *
67 */
68
69#include <stdio.h>
70#include "mconf.h"
71
72int merror = 0;
73
74/* Notice: the order of appearance of the following
75 * messages is bound to the error codes defined
76 * in mconf.h.
77 */
78static char *ermsg[7] = {
79"unknown", /* error code 0 */
80"domain", /* error code 1 */
81"singularity", /* et seq. */
82"overflow",
83"underflow",
84"total loss of precision",
85"partial loss of precision"
86};
87
88
89int mtherr( name, code )
90char *name;
91int code;
92{
93
94/* Display string passed by calling program,
95 * which is supposed to be the name of the
96 * function in which the error occurred:
97 */
98printf( "\n%s ", name );
99
100/* Set global error message word */
101merror = code;
102
103/* Display error message defined
104 * by the code argument.
105 */
106if( (code <= 0) || (code >= 7) )
107 code = 0;
108printf( "%s error\n", ermsg[code] );
109
110/* Return to calling
111 * program
112 */
113return( 0 );
114}
generated by cgit v1.2.3-55-g6feb (git 2.39.1) at 2025-03-21 15:35:43 +0000