1 /******************************************************************************
2 Copyright (c) 2007-2011, Intel Corp.
3 All rights reserved.
4
5 Redistribution and use in source and binary forms, with or without
6 modification, are permitted provided that the following conditions are met:
7
8 * Redistributions of source code must retain the above copyright notice,
9 this list of conditions and the following disclaimer.
10 * Redistributions in binary form must reproduce the above copyright
11 notice, this list of conditions and the following disclaimer in the
12 documentation and/or other materials provided with the distribution.
13 * Neither the name of Intel Corporation nor the names of its contributors
14 may be used to endorse or promote products derived from this software
15 without specific prior written permission.
16
17 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
21 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
27 THE POSSIBILITY OF SUCH DAMAGE.
28 ******************************************************************************/
29
30 #include "dpml_private.h"
31 #include "sqrt_macros.h"
32
33 #undef MAKE_ASINH
34 #undef MAKE_ACOSH
35
36 #if defined(ASINH)
37 # define MAKE_ASINH
38 # define BASE_NAME ASINH_BASE_NAME
39 # define _F_ENTRY_NAME F_ASINH_NAME
40 #elif defined(ACOSH)
41 # define BASE_NAME ACOSH_BASE_NAME
42 # define _F_ENTRY_NAME F_ACOSH_NAME
43 #else
44 # error "Must have one of ASINH, ACOSH defined"
45 #endif
46
47 #if !defined(F_ENTRY_NAME)
48 # define F_ENTRY_NAME _F_ENTRY_NAME
49 #endif
50
51 /*
52 Arcsinh & Arccosh
53 --------------------------------------
54
55
56 This source can be compiled into both Arcsine and Arccosine
57 routines. The definitions necessary to create the function follow.
58
59
60 Function Generation:
61
62 Along with any standard compile time definitions required by the dpml
63 the following items should be defined on the compilation command line to
64 create the indicated routine.
65
66 Arcsinh : ASINH
67
68 Arccosh : ACOSH
69
70
71 To create each routine's 'include' file an initial compilation should be
72 done using the following definition in addition to those above.
73
74 MAKE_INCLUDE
75
76
77 Selectable Build-time Parameters:
78
79 The definitions below define the minimum "overhang" limits for those
80 ranges of the routine with adjustable accuracy bounds. The numbers
81 specified in the definitions are the number of binary digits of
82 overhang. A complete discussion of these values and their use is
83 included in the individual routine documentation.
84 */
85
86 #define POLY_RANGE_OVERHANG 5
87 #define REDUCE_RANGE_OVERHANG 5
88 #define ASYM_RANGE_OVERHANG 7
89 #define LARGE_RANGE_OVERHANG 7
90
91
92
93
94
95 #if !defined(MAKE_INCLUDE)
96 #include STR(BUILD_FILE_NAME)
97 #endif
98
99 /*
100 Arcsinh
101 --------------------------
102
103 The Arcsinh designs described here are the result of
104 an effort to create a fast Arcsinh routine with error bounds near 1/2
105 lsb. The inherent conflict is that, to create fast routines we
106 generally need to give up some accuracy, and conversely, to increase
107 accuracy we often must give up speed. As a result, the design we're
108 presenting defines a user (builder) configureable routine. That is, it
109 is set up such that the builder of a routine may choose, through the
110 proper setting of parameters, the degree of accuracy of the generated
111 routine and hence, indirectly, its speed.
112
113
114 The Design:
115
116 The overall domain of the Arcsinh function has been divided up into
117 six regions or paths as follows:
118
119 (1) (2) (3) (4) (5) (6)
120 |--------|------------|-----------|-----------|-------|----------|
121 0 small polynomial reduction asymptotic large huge
122
123 (Note: Although the domain of Arcsinh extends from -infinite to
124 +infinite, the problem can be considered one of only positive arguments
125 through the application of the identity asinh(-x) = - asinh(x). )
126
127
128 Within each region a unique approximation to the Arcsinh function is
129 used. Each is chosen for its error characteristics, efficiency and the
130 range over which it can be applied.
131
132 1. Small region:
133
134 asinh(x) = x (x <= max_small)
135
136 Within the "small" region the Arcsinh function is approximated as
137 asinh(x) = x. This is a very quick approximation but it may only be
138 applied to small input values. There is effectively no associated
139 storage costs. By limiting the magnitude of x the error bound can
140 be limited to <= 1/2 lsb.
141
142 2. Polynomial region:
143
144 Within the "polynomial" region the function is approximated as
145
146 asinh(x) = x (1 + x^2 P(x)) (max_small_x <x <= max_poly_x)
147
148 where P(x) is a minimax polynomial approximation to (asinh(x)/x -1)/x^2,
149 given by Remes' algorithm and max_poly_x is the upper bound of the
150 polynomial region whose value satisfies:
151
152 (asinh(x)-x)/asinh(x) <= 2^(-POLY_RANGE_OVERHANG)
153
154
155 3. Reduction region:
156
157 In this region, asinh(x) is computed by
158
159 asinh(x) = asinh(x0) + asinh(x*sqrt(1+x0^2)-x0*sqrt(1+x^2))
160
161 max_poly_x < x <= max_reduce_x
162
163 i.e. asinh(x) is computed as the sum of two quantities: asinh(x0),
164 and a reduced value asinh(t), where
165
166 t = x * sqrt(1+x0^2) - x0*sqrt(1+x^2).
167
168 This approach incurs the cost of calculating t and also an
169 lookup table. The values x0, asinh(x0) and sqrt(1+xo^2) are
170 stored in the table to reduce the run time cost. Increased accuracy
171 and efficiency are gained by choosing the asinh(x0) table values such
172 that they have a predetermined number of trailing 0's or 1's beyond
173 the extent of the floating point precision. This reduces the error
174 bound and avoids the need to perform an extended addition between
175 the two quantities. The error bounds here can be established at a value
176 close to 1/2 lsb.
177
178
179 4. Asymptotic region: In this region, asinh(x) is computed as:
180
181 asinh(x) = ln(2x) + 1/4 x^-2 - 3/16 x^-4
182
183 + 5/96 x^-6 -...
184
185 where max_reduce_x < x <= max_asym_x
186
187 The upperbound of the reduction region, max_reduce_x, ( or the lowerbound
188 of the asymptotic region) is determined by finding a smallest x such that
189
190 ((asinh(x) - ln(2x)) / asinh(x)) < 2^-(ASYM_RANGE_OVERHANG)
191
192
193 5. Large region: In this region, asinh(x) is computed as:
194
195 asinh(x) = ln(2x), where max_asym_x < x <= HUGE/2
196
197 where HUGE is the largest floating number represented by the machine
198 and max_asym_x is determined by
199
200 (asinh(x) - ln(2x))/asinh(x) < 2^-(F_PRECISION+LARGE_RANGE_OVERHANG+1)
201
202
203 6. Huge region: In this region, to avoid overflow, asinh(x) is computed as:
204
205 asinh(x) = ln(2) + log(x), where HUGE/2 < x <= HUGE
206
207
208 Special cases:
209
210 Infinities and Nans passed as input result in an Infinities or
211 Nans being returned.
212
213
214
215 Configuring the implementation:
216
217 For polynomial, reduction, asymptotic, large and huge regions (paths 2, 3,
218 4, 5 and 6), the implementation has been set up so that the builder can
219 control the accuracy. This is accomplished by allowing the builder to
220 specify a lower limit to the floating-point alignment shift of the operation
221 which significantly affects round-off error in that range. By establishing
222 an alignment shift the builder determines the relative accuracy of the
223 approximation and thereby determines the effective rate at which the
224 routine will execute.
225
226 In the case of the approximations used for Arcsinh we have the following
227 situations:
228
229 Polynomial region:
230
231 The function is computed as
232
233 x + x^3 P(x)
234
235 where P(x) is approximately
236
237 -1/6 + 3/40 x^2 + ....
238
239 Thus overhang is
240
241 (asinh(x) - x) / asinh(x)
242
243 or approximately x^2 P(x).
244
245
246 Reduced region:
247
248 In this region the result are computed by
249
250 asinh(x) = asinh(xi) + asinh(ti),
251
252 where ti = x * sqrt(1+xi^2) - xi * sqrt(1+x^2). Obviously, we
253 have asinh(ti) = mx - mi, where mx = asinh(x) and mi = asinh(xi).
254 Thus the overhang is determined by
255
256 (mx-mi)/mi.
257
258 The minimum alignment shift can be properly controlled by adjusting
259 the size of interval (m1 - m0) and m0.
260
261 The builder of the Arcsinh routines can specify the overhang limits
262 for each of the above paths (The "small" regions have
263 error bounds pre-established to within 1/2 lsb). The larger the
264 overhang is, the more accurate of the results are. However, larger
265 overhang generates larger size of the lookup table. The values of these
266 overhangs are defined within the main source file header and can be set
267 there.
268
269
270
271 Design Specifics:
272
273 The following sections discuss the design and implementation details of
274 each path.
275
276
277 Note:
278
279
280 Interlaced along with the documentation is the source code necessary to
281 generate the boundary points, constants, etc. Understanding the design
282 should not require an understanding of the source.
283 */
284
285
286 #if defined(MAKE_INCLUDE)
287
288 @divert divertText
289
290 /*
291 The following command establishes the working precision required for
292 accurate computations in the current floating point precision.
293
294 Note: For most all our work, the required precision is bounded by
295 the accuracy needed to generate the "accurate" table values of the
296 'reduce' range. There we need an accuracy equivalent of the current
297 floating point precision plus the number of bits in the "accurate"
298 value overhang.
299
300 precision =
301 ceil( (2*floating_precision +
302 bits of overhang in the reduce range +
303 7-bit potential shift in a normalized MP number) / MP_radix)
304 */
305
306 working_prec = ceil( (2*F_PRECISION + REDUCE_RANGE_OVERHANG + MP_RADIX_BITS-1) / MP_RADIX_BITS ) + 1;
307
308 #define SET_MP_PREC(x) precision = x
309 #define RND_TO_FMT(x) bround(x, F_PRECISION)
310 #define PRINT_POLY_TERM(x) printf(" %#.4" STR(F_CHAR) ", \n", (x) )
311
312
313 /* INC_BIT & DEC_BIT increment and decrement x, respectively, by one
314 bit at position p.
315 */
316
317 #if defined(MAKE_ASINH)
318
319 #define BINARY_EXP(x) (bexp(x) - 1)
320 #define DEC_BIT(x, p) x -= 2^(BINARY_EXP(x) - p)
321 #define INC_BIT(x, p) x += 2^(BINARY_EXP(x) - p)
322 #define TRUNCATE(x, p) bchop(x, p)
323
324
325
326 /* "index" determines the lookup index of a floating point value.
327 It extracts the exponent and necessary fraction bits of a
328 floating point number and returns them as an integer.
329 */
330
index()331 function index()
332 {
333 v = $1;
334
335 /* get base 2 exponent of value */
336 n = BINARY_EXP(v);
337 /* get first K fraction bits of value (NOT
338 including the hidden bit) as an integer */
339 f = bextr(v, 2, K+1);
340
341 /* index = ((bias + n - norm)<< K) + f
342 as an integer */
343 return ( (F_EXP_BIAS + n - F_NORM) * 2^K) + f;
344 }
345
346
347 /* "make_accurate" determines a value, v, such that asinh(v) has
348 trailing 1's or 0's between the limit of the current precision
349 and the extent of the specified overhang.
350 lower_bound < v <= argument.
351 */
352
make_accurate()353 function make_accurate()
354 {
355 v = $1; lower_bound = $2; upper_bound = $3; overhang = $4;
356
357 ones_mask = (2^overhang) - 1;
358 v1 = v2 = v;
359 while (v1 > lower_bound || v2 < upper_bound) {
360 if (v1 > lower_bound) {
361 bits = bextr( asinh(v1), (F_PRECISION + 1), (F_PRECISION + overhang) );
362 if ( (bits == 0) || (bits == ones_mask) )
363 return v1;
364 DEC_BIT(v1, (F_PRECISION - 1));
365 };
366 if (v2 < upper_bound) {
367 INC_BIT(v2, (F_PRECISION - 1));
368 bits = bextr( asinh(v2), (F_PRECISION + 1), (F_PRECISION + overhang) );
369 if ( (bits == 0) || (bits == ones_mask) )
370 return v2;
371 }
372 }
373
374 printf("Couldn't find an ACCURATE value \n");
375 return -1;
376 }
377
378 SET_MP_PREC(working_prec);
379
380 /*
381 The "small" region:
382
383 Within the "small" region asinh(x) is approximated as asinh(x) = x.
384 The reasoning follows:
385
386 Given the Taylor Series expansion to asinh(x)
387
388 asinh(x) = x - x^3 (1/6) + x^5 (3/40) + ..., for x < 1.
389
390 we find that successive terms of the series decrease rapidly in
391 magnitude, and that as x goes to 0, the relative distance between
392 individual terms of the series becomes greater. It is meaningful to
393 ask, at what point does x^3/6, and hence all succeeding terms of the
394 series, become irrelevant with respect to the magnitude of x given a
395 fixed floating point precision? I.e. when is the ratio of x^3/6 to
396 x less than 1 / 2^(precision + 1)? Solving we find:
397
398 let,
399 p = current floating point precision.
400
401 (1/6) * (x^3)/x < 1/2^(p+1) ==>
402
403 x < sqrt( 6/2^(p+1) )
404
405 So, when x < sqrt( 6/2^(p+1) ), asinh(x) correctly rounded is
406 equivalent to x. The value sqrt( 6/2^(p+1) ) will thus be made the
407 upper bound of the "small" region.
408
409 error bounds:
410
411 Since x is equivalent to asinh(x) correctly rounded, the error bound
412 for this approximation is 1/2 lsb.
413 */
414
415 max_small_x = sqrt( 6 / 2^(F_PRECISION + 1) );
416
417
418 /*
419 The "polynomial" region:
420
421 Within the "polynomial" region, Arcsinh is approximated as
422 asinh(x) = x - x^3/6 + x^5 (3/40) + ..., or rather a truncated
423 polynomial approximation to this series. This is a reasonably quick
424 approximation and has storage costs limited to that of the
425 coefficients.
426
427 error bounds:
428
429 The error bound for this approach is roughly determined by two
430 things:
431
432 - The overhang between the first two terms of the series (or
433
434 more exactly, between x and the sum of the remaining terms of
435 the series = x - asinh(x)).
436
437 - The accuracy of the term largest in magnitude.
438
439 The overhang between the terms is a function of x and, in general,
440 will decrease as x gets larger. The largest term is x and, as an
441 input argument, is assumed to be exact. This implies that we can
442 enforce an error bound of our choosing (greater than 1/2 lsb) by
443 limiting the size of our input argument such that the leading two
444 terms maintain a chosen overhang.
445
446 If our desired overhang limit is V, we can compute a maximum X for
447 which that overhang is satisfied by determining when the following
448 is true:
449
450 (X - asinh(X))/asinh(X) < 2^-V
451
452 This point X will be upper bound of the "polynomial" range.
453
454 Note:
455 Since X can not be computed analytically it must be computed
456 numerically (e.g Newton's method).
457
458
459
460 The following code determines the polynomial range
461 */
462
463 rho = 2^-(POLY_RANGE_OVERHANG);
464 c = (1 - rho);
465 x = 0.5;
466 error = 1;
467 while (error > 2^-(2*F_PRECISION)) {
468 f = asinh(x) - x * c;
469 f1 = (1/sqrt(1+x*x)) - c;
470 next_x = x - f/f1;
471 error = abs(f);
472 x = next_x;
473 }
474 max_poly_x = x;
475
476
477
478 /* The following code determines the upper bound of the reduced range
479 (or the lower bound of the asymptotic region.) In the asymptotic region,
480 asinh(x) is determined by the following formula:
481
482 asinh(x) = ln(2x) + 1/4 x^-2 - 3/16 x^-4
483
484 + 5/96 x^-6 -...
485
486 The maximum of the reduction region( or the minimum of the asymptotic
487 region) is determined by finding a minimal x such that
488
489 ((asinh(x) - ln(2x)) / asinh(x)) < 2^-(ASYM_RANGE_OVERHANG);
490 */
491
492 rho = 2^-(ASYM_RANGE_OVERHANG);
493 c = (1 - rho);
494 x = 2.5;
495 error = 1;
496 while (error > 2^-(2*F_PRECISION)) {
497 f = c * asinh(x) - log(2*x);
498 f1 = c * (1/sqrt(1+x*x)) - 1/x ;
499 next_x = x - f/f1;
500 error = abs(f);
501 x = next_x;
502 }
503 max_reduce_x = x;
504
505
506
507 /* The following code determines the upperbound of the asymptotic
508
509 region. ( or the lower bound of the large region.) Within this region,
510 asinh(x) is approximated by
511
512 sign(x) asinh(x) = ln(2|x|) + 1/4 x^-2 - 3/16 x^-4
513
514 + 5/96 x^-6 -...
515
516 The upperbound of the asymptotic region( or the lower bound of the large
517 region) is determined by finding a minimal x such that
518
519 (asinh(x) - ln(2x))/asinh(x) < 2^-(F_PRECISION+LARGE_RANGE_OVERHANG+1)
520 */
521
522 SET_MP_PREC(2*working_prec+1);
523 rho = 2^-(F_PRECISION+LARGE_RANGE_OVERHANG+1);
524 c = (1 - rho);
525 x = 2.5;
526 error = 1;
527 while (error > 2^-(2*F_PRECISION)) {
528 f = c * asinh(x) - log(2*x);
529 f1 = c * (1/sqrt(1+x*x)) - 1/x ;
530 next_x = x - f/f1;
531 error = abs(f);
532 x = next_x;
533 }
534 max_asym_x = x;
535
536
537 /* The following code determines and prints the terms of the
538 asymptotic expansion in the asymptotic region.
539
540 asinh(x) = ln(2x) + 1/4 x^-2 - 3/16 x^-4
541
542 + 5/96 x^-6 -...
543 */
544
545
546 /* Approximation to the function: x^2(asinh(x) - ln(2x)) */
asinh_asym_func()547 function asinh_asym_func()
548 {
549 x = $1;
550 if (x == 0)
551 return 1/4;
552 else
553 return ((x*x)*(asinh(x) - ln(2*x)));
554 }
555
556
557 old_precision = precision;
558 precision = ceil(2*F_PRECISION/MP_RADIX_BITS) + 5;
559
560
561 remes(REMES_FIND_POLYNOMIAL+ REMES_RELATIVE_WEIGHT+ REMES_RECIP_SQUARE_ARG,
562 max_reduce_x, max_asym_x, asinh_asym_func, (F_PRECISION + 1),
563 &poly_deg, &poly_coef);
564
565 printf("#define NUM_ASYM_TERM %i\n", poly_deg + 1 );
566
567 printf("#define EVALUATE_ASYM_RANGE_POLYNOMIAL(x,c,y) \\\n");
568 printf(" POLY_%i(x,c,y)\n", poly_deg + 1);
569
570 printf("static const TABLE_UNION asym_range_coef[] = { \n");
571 for (i = 0; i <= poly_deg; i++)
572 printf("\t%#.4" STR(F_CHAR) ",\n", poly_coef[i] );
573 printf("}; \n\n");
574 precision = old_precision;
575
576
577
578 /* The "reduction" region: Within the "reduction" region asinh(x) is
579 approximated as
580
581 asinh(x) = asinh(m) + asinh(t),
582
583 where t = x*sqrt(1+m*m) - m*sqrt(1+x*x).
584
585 error bounds:
586
587 As in the "polynomial" regions the error bound here is
588 roughly a function of:
589
590 - The overhang between the final addition of asinh(m) and
591 asinh(t).
592
593 - The accuracy of the dominating (larger) value, asinh(m).
594
595 As for the overhang, given some value m, t will decrease
596 as x moves closer to m. This in turn implies asinh(t)
597 becomes smaller. Thus given some x in the "reduce" range we can
598 chose some m near x such that asinh(m) + asinh(t) has
599 at least the desired alignment shift.
600
601 As for the accuracy of asinh(m), we can chose the m above such that
602 asinh(m) has trailing 0's or 1's from the boundary of the current
603 precision to the extent of the chosen overhang. This will reduce
604 the overall error in our computation. These asinh(m) are known as
605 "accurate" values.
606
607 This implies that we can roughly enforce what error bound we choose
608 (greater than 1/2 lsb).
609
610
611
612 Determining the "m" and asinh(m) tables:
613
614 Ensuring a certain overhang, k, between asinh(m) and
615 asinh(t):
616
617 From the identity above, it is obvious that asinh(t) is the
618 difference between asinh(x) and asinh(m). Thus, if we choose a
619 net of equal length subintervals in [asinh(max_poly_x),
620 asihn(max_reduce_x)] such that the size of subinterval is sufficiently
621 small, we can ensure the alignment shift.
622
623 asinh(t)/asinh(m) = (ax - am)/am,
624
625 where ax = asinh(x) and am = asinh(m) and (ax-am) <= 1/2 of the
626 subinterval size.
627
628
629 The generation of these "m",s and x's divides the reduce range up
630 into subintervals like the following:
631
632
633 |--|--|---|---|----|----|-- ... --|
634 Xo Mo X1 M1 X2 M2 X3 Xn
635
636 For each input x such that Xi < x <= Xi+1 the asinh(x) is computed
637 using Mi and asinh(Mi). Note that these subdivisions are
638 non-uniform. As x moves to the right the relative size of these
639 intervals is increasing.
640
641 Indexing the "m" and asinh(m) table:
642
643 Given an x which lies within the 'reduce' range we need an efficient
644 way of determining which "m" and asinh(m) should be used in our our
645 calculations. We will accomplish this using a second table. For
646 each value of x such that Xi < x <= Xi+1 we will use the exponent
647 bits of x and a certain number of fraction bits (enough to uniquely
648 characterize which subinterval x lies within) to act as an index
649 into a second "indexing" table. The values stored within the
650 "indexing" table will in turn point to the appropriate value of "m"
651 (and asinh(m)) to use for the current x.
652
653
654 Mapping input arguments to the appropriate region:
655
656 As described throughout, each input argument x maps to one of the
657 three different approximations for asinh(x). An efficient way of
658 determining which approximation should be used for an input x is to
659 calculate the "index" of each argument x (as we do in the 'reduce'
660 region) and use it to make the choice of approximation.
661 This simply involves determining the index value for the boundary
662 points of each region and then comparing the index of the input
663 argument x to these values and branching accordingly.
664
665
666 The indicies:
667
668 Indices consist of exponents & fraction bits to uniquely characterize
669 an interval. The number of fraction bits indicates table size.
670
671 Calculation of the number of fraction bits needed for the index:
672
673 It is desirable to minimize the number of fraction bits used to
674 address the "indexing" table. Each additional bit we use will
675 increases the size of the "indexing" table by a power of 2. (This is
676 not completely true. Depending on what value is chosen for the
677 upper bound of the 'reduce' range it may not be necessary to store
678 indexes for values at the far end of the range. The increase in
679 table size, however, is still on the order of a power of 2.)
680
681 Since the size of the intervals Xi <--> Xi+1 decrease as x goes to
682 one, we find that the smallest interval for which we need to
683 uniquely specify each x is the interval X(n-1) <--> Xn. The minimum
684 number of fraction bits we need to characterize this interval is
685 given by the overhang difference between X(n-1) and Xn. Thus, the
686 minimum number of fraction bits, k, required to satisfy our index is
687 given by:
688
689 (X(n-1) - Xn)/Xn = 1/2^k
690
691 k = floor( log2 ((Xn - X(n-1))/Xn) )
692
693 as the number of fraction bits required for our index.
694
695
696 Mapping the index:
697
698 Realizing that an arguments exponent and fraction bits are going to
699 be looked at as an integer index to the "indexing" table leaves an
700 issue of addressing. Since we want Xo, the lower bound of the
701 'reduce' range, to map to the first element of the "indexing" table
702 it is necessary to map the generated index of Xo down to zero. This
703 is accomplished by predetermining the index of Xo and using it as a
704 offset (subtracting it off) from the generated index of Xo.
705 If we subtract this offset from all generated indexes we can
706 map the indexes of the 'reduce' range into a table look-up address
707 between 0 and tablesize-1.
708
709
710 The following code:
711 1. compute the minimum number of "index" bits required
712 to accurately determine the mapping of input values to table
713 values in the 'reduce' range.
714 2. compute the accurate table.
715 */
716
717
718 rho = 2^-REDUCE_RANGE_OVERHANG;
719 max_poly_ax = asinh(max_poly_x);
720 x0 = max_poly_x;
721 ax0 = max_poly_ax;
722 max_reduce_ax = asinh(max_reduce_x);
723 max_K = 1;
724 max_delta_ax = 0.0;
725 n = 0;
726 printf("static const TABLE_UNION asinh_tab[] = { \n");
727
728 while(ax0 < max_reduce_ax) {
729 delta_ax = rho * ax0;
730 ax1 = ax0 + delta_ax;
731 x1 = sinh(ax1);
732 if (ax1 > max_reduce_ax) {
733 ax1 = max_reduce_ax;
734 delta_ax = max_reduce_ax - ax0;
735 axm = ax0 + delta_ax/2;
736 x1 = max_reduce_x;
737 xm = sinh(axm);
738 xm = make_accurate( RND_TO_FMT(xm), x0, x1, REDUCE_RANGE_OVERHANG);
739 printf("\t%#.4" STR(F_CHAR) ",", xm);
740 printf("\t%#.4" STR(F_CHAR) ",", asinh(xm));
741 cosh_asinh_value = sqrt(1.0 + xm * xm);
742 printf("\t%#.4" STR(F_CHAR) ",\n", cosh_asinh_value);
743 n++;
744 if (max_delta_ax < delta_ax) max_delta_ax = delta_ax;
745 break;
746 }
747 axm = ax0 + delta_ax/2;
748 xm = sinh(axm);
749 xm = make_accurate( RND_TO_FMT(xm), x0, x1, REDUCE_RANGE_OVERHANG);
750 printf("\t%#.4" STR(F_CHAR) ",", xm);
751 printf("\t%#.4" STR(F_CHAR) ",", asinh(xm));
752 cosh_asinh_value = sqrt(1.0 + xm * xm);
753 printf("\t%#.4" STR(F_CHAR) ",\n", cosh_asinh_value);
754 delta_x = x1 - x0;
755 K = abs( floor ( log2 (delta_x/x0)));
756 if (max_K < K) max_K = K;
757 if (max_delta_ax < delta_ax) max_delta_ax = delta_ax;
758 ax0 = ax1;
759 x0 = x1;
760 n++;
761 }
762 count = n;
763 printf("}; \n\n");
764 printf("#define TABLE_ENTRY_SIZE %i \n", count);
765 max_delta_x = sinh(max_delta_ax);
766 K = max_K;
767 K = K - 1;
768 printf("#define K %i \n", K);
769
770 /* computation of the "index" represented by the maximum
771 value we will evaluate in the 'polynomial' range. Since the
772 'reduce' range begins beyond this value, we will use this
773 "index" to map all "indicies" so that those of the 'reduce'
774 range will take on values between 0 and the table size.
775 */
776
777 offset = index( TRUNCATE(max_poly_x, K+1) );
778 printf("#define OFFSET_IND %i \n", offset);
779
780 DEC_BIT(max_small_x, K);
781
782 printf("#define MAX_SMALL_INDEX %i \n", (index(max_small_x)-offset));
783
784 printf("#define MAX_POLY_INDEX %i \n", (index(max_poly_x)-offset));
785
786 printf("#define MAX_REDUCE_INDEX %i \n", (index(max_reduce_x)-offset));
787
788 printf("#define MAX_ASYM_INDEX %i \n", (index(max_asym_x)-offset));
789
790 half_huge_x = MPHOC_F_POS_HUGE/2;
791
792 printf("#define HALF_HUGE_INDEX %i \n", (index(half_huge_x)-offset));
793
794
795
796 /* computation of the "accurate" table "indecies" for values
797 within the 'reduce' range.
798 */
799
800 rho = 2^-REDUCE_RANGE_OVERHANG;
801
802 if (count <= 256)
803 printf("static const U_INT_8 asinh_index_table[] = { ");
804 else
805 printf("static const U_INT_16 asinh_index_table[] = { ");
806 bytes_used = 0;
807 i = 0;
808 j = 0;
809 m = 0;
810 x0 = max_poly_x;
811 ax0 = asinh(x0);
812
813 while(ax0 < max_reduce_ax) {
814 delta_ax = rho * ax0;
815 ax1 = ax0 + delta_ax;
816 if (ax1 > max_reduce_ax) {
817 ax1 = max_reduce_ax;
818 delta_ax = max_reduce_ax - ax0;
819 }
820 x1 = sinh(ax1);
821 b = TRUNCATE( x0, K+1);
822 next_b = TRUNCATE( x1, K+1);
823 cur_index = b;
824 while (cur_index < next_b) {
825 if ((j++ % 8) == 0) {
826 printf("\n\t");
827 bytes_used = 1;
828 }
829 else
830 bytes_used++;
831 printf("%i, ", i);
832 INC_BIT(cur_index, K);
833 }
834 i++;
835 ax0 = ax1;
836 x0 = x1;
837 }
838 i--;
839 /*
840 * Pad this table up to avoid alignment problem on HP machine.
841 */
842 bytes_to_pad = 8 - bytes_used;
843 for (k = 0; k < bytes_to_pad; k++)
844 printf("0, " );
845
846 printf("\n}; \n\n");
847
848
849
850
851 /* Generate the coefficients for the 'polynomial' range polynomial */
852
853
854 /* Approximation to the function: (asinh(x) - x) / x^3 */
asinh_func()855 function asinh_func()
856 {
857 if ($1 == 0)
858 return -1/6;
859 else
860 return (asinh($1) - ($1))/($1 * $1 * $1);
861 }
862
863
864 old_precision = precision;
865 precision = ceil(2*F_PRECISION/MP_RADIX_BITS) + 5;
866
867
868 remes(REMES_FIND_POLYNOMIAL + REMES_RELATIVE_WEIGHT + REMES_SQUARE_ARG,
869 0.0, max_poly_x, asinh_func, (F_PRECISION),
870 &poly_deg, &poly_coef);
871
872 printf("#define EVALUATE_POLY_RANGE_POLYNOMIAL(x,c,y) \\\n");
873 printf(" ODD_POLY_%i_U(x,c,y)\n", 2 * poly_deg + 3);
874
875 printf("static const TABLE_UNION poly_range_coef[] = { \n");
876 for (i = 0; i <= poly_deg; i++)
877 printf("\t%#.4" STR(F_CHAR) ",\n", poly_coef[i] );
878 printf("}; \n\n");
879
880
881
882 /* Generate the coefficients for the 'reduce' range polynomial */
883
884 remes(REMES_FIND_POLYNOMIAL + REMES_RELATIVE_WEIGHT + REMES_SQUARE_ARG,
885 0.0, max_delta_x, asinh_func, (F_PRECISION + 1),
886 &poly_deg, &poly_coef);
887
888 printf("#define EVALUATE_REDUCE_RANGE_POLYNOMIAL(x,c,y) \\\n");
889 printf(" ODD_POLY_%i_U(x,c,y)\n", 2 * poly_deg + 3);
890
891 printf("static const TABLE_UNION reduce_range_coef[] = { \n");
892 for (i = 0; i <= poly_deg; i++)
893 printf("\t%#.4" STR(F_CHAR) ",\n", poly_coef[i] );
894 printf("}; \n\n");
895
896
897 precision = old_precision;
898
899
900
901 #else /* ACOSH */
902
903 /* The computation of ACOSH is devided into four regions:
904
905 1. Direct region: In this region, ACOSH(x) is computed as
906
907 acosh(x) = log(x + sqrt((x-1)*(x+1)))
908
909 where 1 < x <= max_direct_x.
910
911
912 2. Asymptotic regions: In this region, ACOSH(x) is computed as
913
914 acosh(x) = ln(2x) - 1/4 x^-2 - 3/16 x^-4
915
916 - 5/96 x^-6 -...
917
918 where max_direct_x < x <= max_asym_x.
919
920 3. Large region: In this region, ACOSH(x) is computed as
921
922 acosh(x) = log(2x).
923
924 where max_asym_x < x <= HUGE/2.
925
926
927 4. Huge region: In this region, ACOSH(x) is computed as
928
929 acosh(x) = ln(2) + log(x).
930
931 where HUGE/2 < x <= HUGE.
932 */
933
934
935 /* The following code determines the upper bound of the direct region
936 (or the lower bound of the asymptotic region). In the asymptotic region,
937 asinh(x) is determined by the following formula:
938
939 acosh(x) = ln(2x) - 1/4 x^-2 - 3/16 x^-4
940
941 - 5/96 x^-6 -...
942
943 The upper bound of the direct region region is determined by finding
944 the smallest x such that
945
946 ((ln(2x)-acosh(x)) / acosh(x)) < 2^-(ASYM_RANGE_OVERHANG);
947 */
948
root_func()949 function root_func() {
950 acosh_x = acosh($1);
951 return (ln(2*($1)) - acosh_x)/acosh_x;
952 }
953
954 lo = binc(1,(F_PRECISION - 1));
955 hi = 10;
956 rho = 2^-(ASYM_RANGE_OVERHANG);
957 max_direct_x = find_root(MP_FIND_ROOT_NO_DERIV, lo, hi, rho, root_func);
958
959 printf("static const TABLE_UNION max_direct_x[] = { %#.4" STR(F_CHAR) " };\n\n", max_direct_x);
960
961
962
963 /* The following code determines the the upper bound of the asymptotic
964 region.(or the lower bound of the large region.) In the large region,
965 acosh(x) is approximated by
966
967 acosh(x) = ln(2x).
968
969 Thus, the lower bound of the large region (or max_asym_x) is determined
970 by finding the smallest x such that
971
972 -(ln(2x)-acosh(x))/acosh(x) < 2^(-F_PRECISION+LARGE_RANGE_OVER_HANG+1)
973 */
974
975 SET_MP_PREC(2 * working_prec+1);
976 rho = 2^-(F_PRECISION+LARGE_RANGE_OVERHANG+1);
977 c = (1 + rho);
978 x = 2.5;
979 error = 1;
980 while (error > 2^-(2*F_PRECISION)) {
981 f = log(2*x) - c * acosh(x);
982 f1 = 1/x - c * (1/sqrt(x*x-1));
983 next_x = x - f/f1;
984 error = abs(f);
985 x = next_x;
986 }
987 max_asym_x = x;
988
989 printf("static const TABLE_UNION max_asym_x[] = { %#.4" STR(F_CHAR) " };\n\n", max_asym_x);
990
991
992 /* The computation in asymptotic and large region is more efficient
993 than the computation in the direct region if the number of the
994 terms in the asymptotic expansion is kept reasonably small. The
995 following procedure computes and prints the terms of the
996 asymptotic series.
997 */
998
999 /* Approximation to the function: -[x^2(acosh(x) - ln(2x))] */
acosh_func()1000 function acosh_func()
1001 {
1002 x = $1;
1003 if (x == 0)
1004 return 1/4;
1005 else
1006 return -((x*x)*(acosh(x) - ln(2*x)));
1007 }
1008
1009
1010 old_precision = precision;
1011 precision = ceil(2*F_PRECISION/MP_RADIX_BITS) + 5;
1012
1013
1014 remes(REMES_FIND_POLYNOMIAL+ REMES_RELATIVE_WEIGHT+ REMES_RECIP_SQUARE_ARG,
1015 max_direct_x, max_asym_x, acosh_func, (F_PRECISION + 1),
1016 &poly_deg, &poly_coef);
1017
1018 printf("#define EVALUATE_ASYM_RANGE_POLYNOMIAL(x,c,y) \\\n");
1019 printf(" POLY_%i(x,c,y)\n", poly_deg + 1);
1020
1021 printf("static const TABLE_UNION asym_range_coef[] = { \n");
1022 for (i = 0; i <= poly_deg; i++)
1023 printf("\t%#.4" STR(F_CHAR) ",\n", poly_coef[i] );
1024 printf("}; \n\n");
1025 precision = old_precision;
1026
1027
1028 half_huge_x = MPHOC_F_POS_HUGE/2;
1029 printf("static const TABLE_UNION half_huge_x[] = { %#.4" STR(F_CHAR) " };\n\n", half_huge_x);
1030
1031 #endif /* MAKE_ASINH */
1032
1033 /* compute correctly rounded "high" and "low" parts of log(2) */
1034
1035 SET_MP_PREC(2 * working_prec);
1036 log_2 = log(2);
1037 printf("static const TABLE_UNION log_2[] = {\n");
1038 printf(" %#.4" STR(F_CHAR) " \n", log_2);
1039 printf("};\n\n");
1040
1041 SET_MP_PREC(working_prec);
1042
1043
1044 /* the following values are defined for use in performing
1045 automated testing with MTC.
1046 */
1047
1048 #ifdef MTC_DEFS
1049 huge_x = (MPHOC_F_POS_HUGE/10)*9;
1050 half_huge_x = MPHOC_F_POS_HUGE/2;
1051 #ifdef MAKE_ASINH
1052 printf("#define MAX_SMALL_PT m:%m\n", max_small_x );
1053 printf("#define MAX_POLY_PT m:%m\n", max_poly_x );
1054 printf("#define MAX_REDUCE_PT m:%m\n", max_reduce_x );
1055 printf("#define MAX_ASYM_PT m:%m\n", max_asym_x );
1056 printf("#define HUGE_PT m:%m\n", huge_x );
1057 printf("#define HALF_HUGE_PT m:%m\n", half_huge_x );
1058 #else
1059 printf("#define MAX_DIRECT_PT m:%m\n", max_direct_x );
1060 printf("#define MAX_ASYM_PT m:%m\n", max_asym_x );
1061 printf("#define HUGE_PT m:%m\n", huge_x );
1062 printf("#define HALF_HUGE_PT m:%m\n", half_huge_x );
1063 #endif
1064 #endif
1065
1066 @end_divert
1067 @eval my $outText = MphocEval( GetStream( "divertText" ) ); \
1068 my $headerText = GetHeaderText( STR(BUILD_FILE_NAME), \
1069 "Definitions and constants for " . \
1070 STR(F_ENTRY_NAME), __FILE__); \
1071 print "$headerText\n$outText";
1072
1073 #endif /* MAKE_INCLUDE */
1074
1075
1076
1077
1078 typedef struct { F_TYPE x_value, asinh_value, cosh_asinh_value;} TABLE_ENTRY;
1079
1080
1081
1082 /* macros specific to asinh */
1083
1084 #define MAX_DIRECT_X *(F_TYPE *)max_direct_x
1085 #define MAX_ASYM_X *(F_TYPE *)max_asym_x
1086 #define HALF_HUGE_X *(F_TYPE *)half_huge_x
1087
1088 #define POLY_RANGE_COEF (F_TYPE *)poly_range_coef
1089 #define REDUCE_RANGE_COEF (F_TYPE *)reduce_range_coef
1090 #define ASYM_RANGE_COEF (F_TYPE *)asym_range_coef
1091
1092 #define ASINH_TABLE(indx) *((TABLE_ENTRY *)asinh_tab + indx)
1093 #define LOG_2 *(F_TYPE *)log_2
1094
1095 #if VAX_FLOATING
1096 # if BITS_PER_WORD == 32
1097 # define INDEX_SHIFT (F_EXP_POS + (BITS_PER_WORD - 16) - K)
1098 # else
1099 # define INDEX_SHIFT (F_EXP_POS + (BITS_PER_F_TYPE - 16) - K)
1100 # endif
1101 # define IEEE_ONLY(x)
1102 #elif IEEE_FLOATING
1103 # define INDEX_SHIFT (F_EXP_POS - K)
1104 # define IEEE_ONLY(x) x
1105 #endif
1106
1107
1108 #define GET_INDEX(exp_word, indx) (indx) = (exp_word >> INDEX_SHIFT) - OFFSET_IND;
1109
1110 #define GET_MID_POINT(indx, m, sqrt_one_add_m2) \
1111 (indx) = asinh_index_table[(indx)]; \
1112 (m) = (ASINH_TABLE(indx)).x_value; \
1113 (sqrt_one_add_m2) = (ASINH_TABLE(indx)).cosh_asinh_value
1114
1115
1116
1117 #define ADD_ASINH_TABLE_VALUE(x, indx, y) \
1118 (y) = (F_TYPE)((ASINH_TABLE(indx)).asinh_value + (x))
1119
1120 #if defined(MAKE_ASINH)
1121
1122 F_F_PROTO( F_LN_NAME ) ;
1123
1124 F_TYPE
F_ENTRY_NAME(F_TYPE x)1125 F_ENTRY_NAME (F_TYPE x)
1126 {
1127 EXCEPTION_RECORD_DECLARATION
1128 F_TYPE z, m, t, u, sqrt_one_add_m2, sqrt_one_add_z2, one_add_z2;
1129 F_SIGN_TYPE sign;
1130 WORD indx;
1131 U_WORD exp_word_z;
1132
1133 F_SAVE_SIGN_AND_GET_ABS(x, sign, z);
1134 GET_EXP_WORD(z, exp_word_z);
1135 exp_word_z = PDP_SHUFFLE(exp_word_z);
1136 if (F_EXP_WORD_IS_INFINITE_OR_NAN(exp_word_z)) return(x);
1137
1138 GET_INDEX(exp_word_z, indx);
1139 if (indx >= MAX_REDUCE_INDEX) goto asym_or_ln2x;
1140
1141 /* Reduced region */
1142
1143 if (indx > MAX_POLY_INDEX) {
1144 GET_MID_POINT(indx, m, sqrt_one_add_m2);
1145 one_add_z2 = ((F_TYPE) 1.0 + z*z);
1146 t = (z - m)*(z + m);
1147 F_HW_OR_SW_SQRT(one_add_z2, sqrt_one_add_z2);
1148 t = t/(z * sqrt_one_add_m2 + m * sqrt_one_add_z2);
1149 EVALUATE_REDUCE_RANGE_POLYNOMIAL(t, REDUCE_RANGE_COEF, z);
1150 ADD_ASINH_TABLE_VALUE(z, indx, t);
1151 F_NEGATE_IF_SIGN_NEG(sign, t);
1152 return t;
1153 }
1154
1155 /* Polynomial region */
1156
1157 else if (indx > MAX_SMALL_INDEX) {
1158 EVALUATE_POLY_RANGE_POLYNOMIAL(z, POLY_RANGE_COEF, t);
1159 F_NEGATE_IF_SIGN_NEG(sign, t);
1160 return(t);
1161 }
1162
1163 /* Small region */
1164
1165 else return (x);
1166
1167
1168 asym_or_ln2x:
1169
1170 /* Asymptotic region */
1171
1172 if (indx < MAX_ASYM_INDEX) {
1173 u = 1/(z*z);
1174 EVALUATE_ASYM_RANGE_POLYNOMIAL(u, ASYM_RANGE_COEF, t);
1175 t += F_LN_NAME(2*z);
1176 F_NEGATE_IF_SIGN_NEG(sign, t);
1177 return(t);
1178 }
1179
1180 /* Large region */
1181
1182 else if ( indx < HALF_HUGE_INDEX) {
1183 t = F_LN_NAME(2*z);
1184 F_NEGATE_IF_SIGN_NEG(sign, t);
1185 return(t);
1186 }
1187
1188 /* Huge region */
1189
1190 else {
1191 t = LOG_2
1192 + F_LN_NAME(z);
1193 F_NEGATE_IF_SIGN_NEG(sign, t);
1194 return t;
1195 }
1196 }
1197
1198 #else /* ACOSH */
1199
1200 F_F_PROTO( F_LN_NAME ) ;
1201 F_F_PROTO( F_LOG1P_NAME ) ;
1202
1203 F_TYPE
F_ENTRY_NAME(F_TYPE x)1204 F_ENTRY_NAME (F_TYPE x)
1205 {
1206 EXCEPTION_RECORD_DECLARATION
1207 F_TYPE t, x2_sub_one, sqrt_x2_sub_one;
1208 F_TYPE u;
1209 U_WORD exp_word_x;
1210
1211
1212 GET_EXP_WORD(x, exp_word_x);
1213 IEEE_ONLY(if (F_EXP_WORD_IS_ABNORMAL(exp_word_x)) goto non_finite_x;)
1214
1215 if (x <= 1.0) goto out_of_range_or_one;
1216
1217 if (x > MAX_DIRECT_X) goto asym_or_ln2x;
1218 else {
1219 /* Direct region */
1220
1221 x2_sub_one = (x - (F_TYPE)1.0) * (x + (F_TYPE)1.0);
1222 F_HW_OR_SW_SQRT(x2_sub_one, sqrt_x2_sub_one);
1223 u = x - (F_TYPE)1.0;
1224 u += sqrt_x2_sub_one;
1225 t = F_LOG1P_NAME(u);
1226 return(t);
1227 }
1228
1229 asym_or_ln2x:
1230
1231 if (x < MAX_ASYM_X) {
1232
1233 /* Asymptotic region */
1234
1235 u = 1/(x*x);
1236 EVALUATE_ASYM_RANGE_POLYNOMIAL(u, ASYM_RANGE_COEF, t);
1237 t = F_LN_NAME(2*x) - t;
1238 return(t);
1239 }
1240 else if ( x < HALF_HUGE_X) {
1241
1242 /* Large region */
1243
1244 t = F_LN_NAME(2*x);
1245 return(t);
1246 }
1247 else {
1248 /* Huge region */
1249
1250 t = LOG_2
1251 + F_LN_NAME(x);
1252 return t;
1253 }
1254
1255 out_of_range_or_one:
1256
1257 if (x != 1.0) goto invalid_argument;
1258 return (F_TYPE) 0.0;
1259
1260 #if IEEE_FLOATING
1261
1262 non_finite_x:
1263
1264 F_CLASSIFY(x, exp_word_x);
1265 if ((exp_word_x == F_C_POS_INF) || (F_C_BASE_CLASS(exp_word_x) == F_C_NAN))
1266 return x;
1267
1268 #endif
1269
1270 invalid_argument:
1271
1272 GET_EXCEPTION_RESULT_1(ACOSH_ARG_LT_ONE, x, t);
1273 return t;
1274
1275 }
1276
1277 #endif /* MAKE_ASINH */
1278
1279
1280
1281 #ifdef MTC
1282
1283 #ifdef MAKE_ASINH
1284 @divert > dpml_asinh.mtc
1285 #else
1286 @divert > dpml_acosh.mtc
1287 #endif
1288 /*
1289 ** accuracy and key point tests for single
1290 ** and double precision asinh or acosh functions.
1291 */
1292
1293 #ifdef MAKE_ASINH
1294
1295 #if (F_PRECISION==24)
1296 build asinh_build = STR(PASTE(F_ENTRY_NAME , .o)) "my_asinh.o" "my_logf.o"
1297 "my_log.o" "dpml_globals.o" "dpml_exception.o" "dpml_sqrt_s_table.o"
1298 "dpml_sqrt_t_table.o";
1299 function asinh_func = F_CHAR F_ENTRY_NAME( F_CHAR.v.r );
1300 comparison_function asinh_backup = B_CHAR my_asinh( B_CHAR.v.r );
1301
1302 #else
1303 build asinh_build = STR(PASTE(F_ENTRY_NAME , .o)) "my_log.o"
1304 "dpml_globals.o" "dpml_exception.o" "dpml_sqrt_t_table.o";
1305 function asinh_func = F_CHAR F_ENTRY_NAME( F_CHAR.v.r );
1306 comparison_function asinh_backup = void mp_asinh(m.r.r, m.r.w);
1307 #endif
1308
1309 domain asinh_accuracy =
1310 { [ 0 , MAX_SMALL_PT ):uniform:1000 }
1311 { [ MAX_SMALL_PT, MAX_POLY_PT ):uniform:4000 }
1312 { [ MAX_POLY_PT, MAX_REDUCE_PT):uniform:4000 }
1313 { [ MAX_REDUCE_PT, MAX_ASYM_PT):uniform:4000 }
1314 { ( MAX_ASYM_PT, HALF_HUGE_PT):uniform:4000 }
1315 { ( HALF_HUGE_PT, HUGE_PT):uniform:1000 }
1316 ;
1317
1318 domain asinh_keypoint =
1319 lsb = 0.5;
1320 { -1.0 | der }
1321 { 0.0 | der }
1322 ;
1323
1324
1325 test asinh_acc =
1326 type = accuracy
1327 error = lsb;
1328 stats = max;
1329 points = 1024;
1330 ;
1331
1332 function = asinh_func;
1333
1334 comparison_function = asinh_backup;
1335
1336 domain = asinh_accuracy;
1337
1338 build = asinh_build;
1339
1340 output =
1341 file = STR(PASTE(F_ENTRY_NAME , _acc.out));
1342 ;
1343 ;
1344
1345
1346 test asinh_key =
1347 type = key_point;
1348
1349 function = asinh_func;
1350
1351 comparison_function = asinh_backup;
1352
1353 domain = asinh_keypoint;
1354
1355 build = asinh_build;
1356 output =
1357 file = STR(PASTE(F_ENTRY_NAME, _key.out));
1358 style = verbose;
1359 ;
1360 ;
1361 /* For testing acosh */
1362 #else
1363 build acosh_build = STR(PASTE(F_ENTRY_NAME, .o)) "my_log.o" "my_logf.o"
1364 "logp.o" "logpf.o" "dpml_globals.o" "dpml_exception.o" "dpml_sqrt_s_table.o" "dpml_sqrt_t_table.o";
1365
1366 function acosh_func = F_CHAR F_ENTRY_NAME( F_CHAR.v.r );
1367
1368 comparison_function acosh_backup = void mp_acosh(m.r.r, m.r.w);
1369
1370 domain acosh_accuracy =
1371 { [ 1.0, MAX_DIRECT_PT):uniform:4000 }
1372 { [ MAX_DIRECT_PT, MAX_ASYM_PT):uniform:4000 }
1373 { [ MAX_ASYM_PT, HALF_HUGE_PT):uniform:4000 }
1374 { [ HALF_HUGE_PT, HUGE_PT):uniform:4000 }
1375 ;
1376
1377 domain acosh_keypoint =
1378 lsb = 0.5;
1379 { 1.0 | der }
1380 { -1.0 | der }
1381 ;
1382
1383
1384 test acosh_acc =
1385 type = accuracy
1386 error = lsb;
1387 stats = max;
1388 points = 1024;
1389 ;
1390
1391 function = acosh_func;
1392
1393 comparison_function = acosh_backup;
1394
1395 domain = acosh_accuracy;
1396
1397 build = acosh_build;
1398
1399 output =
1400 file = STR(PASTE(F_ENTRY_NAME, _acc.out));
1401 ;
1402 ;
1403
1404
1405 test acosh_key =
1406 type = key_point;
1407
1408 function = acosh_func;
1409
1410 comparison_function = acosh_backup;
1411
1412 domain = acosh_keypoint;
1413
1414 build = acosh_build;
1415 output =
1416 file = STR(PASTE(F_ENTRY_NAME , _key.out));
1417 style = verbose;
1418 ;
1419 ;
1420
1421
1422 #endif /* MAKE_ASINH */
1423 @end_divert
1424
1425 #endif /* MTC */
1426
1427