overby.me
this repo has no description
1/*
2 * Copyright (c) 2002 Apple Computer, Inc. All rights reserved.
3 *
4 * @APPLE_LICENSE_HEADER_START@
5 *
6 * The contents of this file constitute Original Code as defined in and
7 * are subject to the Apple Public Source License Version 1.1 (the
8 * "License"). You may not use this file except in compliance with the
9 * License. Please obtain a copy of the License at
10 * http://www.apple.com/publicsource and read it before using this file.
11 *
12 * This Original Code and all software distributed under the License are
13 * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER
14 * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
15 * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT. Please see the
17 * License for the specific language governing rights and limitations
18 * under the License.
19 *
20 * @APPLE_LICENSE_HEADER_END@
21 */
22//
23// d3ops.c
24//
25// Double-double Function Library
26// Copyright: � 1995-96 by Apple Computer, Inc., all rights reserved
27//
28// Triple-double precision basic operations.
29//
30// (high, mid, low) <-- (a, b, c) <op> (d, e, f)
31// where <op> is add, multiply, or divide
32//
33// void _d3div(double a, double b, double c, double d, double e, double f,
34// double *high, double *mid, double *low);
35// void _d3mul(double a, double b, double c, double d, double e, double f,
36// double *high, double *mid, double *low);
37// void _d3add(double a, double b, double c, double d, double e, double f,
38// double *high, double *mid, double *low);
39//
40
41#include "math.h"
42#include "fp_private.h"
43#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
44#include "DD.h"
45
46static const hexdouble infinity = {{0x7ff00000, 0x00000000}};
47static const hexdouble tiny = {{0x00100000, 0x00000000}};
48
49
50void _d3div(double a, double b, double c, double d, double e, double f,
51 double *high, double *mid, double *low)
52{
53 // (high, mid, low) <-- (a, b, c) / (d, e, f)
54
55 double p, q, r, s, t, u, v, w, x, y, z;
56 double g, h, o, rmid, rlow;
57 // int lx, l1, l2;
58
59 // lx = (fabs(d) < tiny.d); // look out for overflow on 1/d
60 p = 1.0/d; // reciprocal of h.o. denominator term
61 if (fabs(d) < tiny.d) {
62 r=a/d; // if quotient is infinity,
63 if (fabs(r)==infinity.d) { // stop process, and return
64 *high = r;
65 *mid = 0.0e0;
66 *low = 0.0e0; // properly signed infinity
67 return;
68 }
69 }
70 else {
71 q = a*p;
72 r = q + (a - q*d)*p; // h.o. quotient term
73 }
74 s = a - r*d; // exact remainder
75 t = __FMUL( r, e); // second order remainder term
76 // l1 = (fabs(b) > fabs(t))
77 u = (r*e - t) + r*f; // third order remainder term
78 w = __FSUB( b, t ); // combine the
79 y = __FADD( s, w ); // second order terms
80 // l2 = (fabs(s) > fabs(w)); // and catch bits that may fall
81
82 if (fabs(b) > fabs(t)) // off the end to include in the
83 x = b - w - t; // third order term
84 else
85 x = b - (w + t);
86 if (fabs(s) > fabs(w))
87 z = s - y + w;
88 else
89 z = w - y + s;
90 g = c - u + x + z; // third order term
91 if (fabs(d) < tiny.d) { // caution against tiny quotients, again
92 o = y/d;
93 v = (y - __FMUL( o, d ) + g - (e*o + (f*o))) / d;
94 }
95 else {
96 h = y*p;
97 o = h + (y - h*d)*p; // second order term
98 v = (y - o*d + g - (e*o + (f*o)))*p;// third order quotient term
99 }
100 rmid = o + v; // start distillation of quotient
101 rlow = o - rmid + v;
102
103 *high = __FADD( r, rmid );
104 *mid = (r - *high + rmid) + rlow;
105 *low = (r - *high + rmid) - *mid + rlow;
106
107 return;
108}
109
110
111void _d3mul(double a, double b, double c, double d, double e, double f,
112 double *high, double *mid, double *low)
113{
114 // (high, mid, low) <-- (a, b, c) * (d, e, f)
115
116 double s, t, u, v, w, x, y, z, tlow, resbump;
117 // int l1, l2;
118
119 u = d*c + e*b + f*a; // sextuple precision multiplication
120 v = __FMUL( d, b ); // derived from _q_mp10.c
121 w = d*b - v;
122 x = __FMUL( e, a );
123 y = e*a - x;
124 // l1=(fabs(x) > fabs(v));
125 s = w + y + u;
126 z = __FMUL( d, a);
127 w = d*a - z;
128 u = __FADD( v, x );
129 // l2=(fabs(w) > fabs(u));
130 if (fabs(x) > fabs(v))
131 s = x - u + v + s;
132 else
133 s = v - u + x + s;
134 t = __FADD( u, w );
135 if (fabs(w) > fabs(u))
136 tlow = w - t + u + s;
137 else
138 tlow = u - t + w + s;
139 resbump = __FADD( t, tlow );
140 tlow = t - resbump + tlow;
141
142 *high = __FADD( z, resbump );
143 *mid = (z - *high + resbump) + tlow;
144 *low = (z - *high + resbump) - *mid + tlow;
145
146 return;
147}
148
149
150void _d3add(double a, double b, double c, double d, double e, double f,
151 double *high, double *mid, double *low)
152{
153 // (high, mid, low) <-- (a, b, c) + (d, e, f)
154
155 // int l1, l2, l3, l4, l5, l5, l6, l7, l8;
156 double p, q, r, s, t, u, v, w, x, y, z;
157 double res, carry, z1, z2, resmid, reslow;
158
159 // l1=(fabs(a) > fabs(d));
160 // l2=(fabs(b) > fabs(e));
161 // l3=(fabs(c) > fabs(f));
162
163 p = __FADD( a, d );
164 q = __FADD( b, e );
165 r = __FADD( c, f );
166
167 if (fabs(a) > fabs(d)) // a b c
168 s = a - p + d; // d e f
169 else // p s (replaces a and d
170 s = d - p + a;
171
172 if (fabs(b) > fabs(e))
173 t = b - q + e; // q t (replaces b and e
174 else
175 t = e - q + b;
176
177 // l4=(fabs(s) > fabs(q));
178 // l5=(fabs(t) > fabs(r));
179
180 if (fabs(c) > fabs(f))
181 u = c - r + f; // r u (replaces c and f
182 else
183 u = f - r + c;
184
185 v = __FADD( s, q );
186 x = __FADD( t, r );
187
188 // l6=(fabs(p) > fabs(v));
189
190 if (fabs(s) > fabs(q))
191 w = s - v + q; // v w (replaces s and q
192 else
193 w = q - v + s;
194
195 // l7=(fabs(x) > fabs(w));
196
197 if (fabs(t) > fabs(r))
198 y = t - x + r; // x y (replaces t and r
199 else
200 y = r - x + t;
201
202 z = __FADD( x, w );
203 res = __FADD( p, v );
204
205 if (fabs(p) > fabs(v))
206 carry = p - res + v; // res carry (replaces p and v
207 else
208 carry = v - res + p;
209
210 // l8=(fabs(carry) > fabs(z));
211
212 if (fabs(x) > fabs(w))
213 z1 = x - z + w; // z z1 (replaces x and w
214 else
215 z1 = w - z + x;
216
217 z2 = u + y + z1; // z2 (replaces u, y, and z1
218 resmid = __FADD( carry, z );
219
220 if (fabs(carry) > fabs(z))
221 reslow = carry - resmid + z; // resmid carry2 (replaces carry, z
222 else
223 reslow = z - resmid + carry;
224
225 *high = __FADD( res, resmid );
226 *mid = (res - *high + resmid) + reslow;
227 *low = (res - *high + resmid) - *mid + reslow + z2;
228
229 return;
230}
231#endif
232#if 0 // set to 1 to include code for self-contained sanity test
233
234#include <stdio.h>
235
236main()
237{
238 double a1,a2,a3, b1,b2,b3, c1,c2,c3, d1,d2,d3, e1,e2,e3;
239
240 _d3div( 2, 0, 0, 3, 0, 0, &a1, &a2, &a3); // a = 2/3
241 _d3div( 5, 0, 0, 7, 0, 0, &b1, &b2, &b3); // b = 5/7
242 _d3mul(a1, a2, a3, b1, b2, b3, &c1, &c2, &c3); // c = a*b = 10/21
243 _d3div(21, 0, 0, 10, 0, 0, &d1, &d2, &d3); // d = 21/10
244 _d3mul(c1, c2, c3, d1, d2, d3, &e1, &e2, &e3); // e = c*d (should be 1)
245
246 printf("e: %24.18e %24.18e %24.18e\n", e1, e2, e3);
247
248 _d3add(a1, a2, a3, b1, b2, b3, &c1, &c2, &c3); // c = a+b = 29/21
249 _d3div( 8, 0, 0, 21, 0, 0, &d1, &d2, &d3); // d = 8/21
250 _d3add(c1, c2, c3,-d1,-d2,-d3, &e1, &e2, &e3); // e = c-d (should be 1)
251
252 printf("e: %24.18e %24.18e %24.18e\n", e1, e2, e3);
253}
254
255// Correct test results:
256//
257// e: 1.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
258// e: 1.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
259
260#endif