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 #define BASE_NAME exp
31 #include "dpml_ux.h"
32
33 #if !defined(MAKE_INCLUDE)
34 # include STR(BUILD_FILE_NAME)
35 #endif
36
37 extern _X_FLOAT PACKED_CONSTANT_TABLE[ LAST_CONS_INDEX ];
38
39 /*
40 ** UX_EXP_REDUCE performs argument reduction for the exponential family of
41 ** functions. Given and input argument, x, UX_EXP_REDUCE computes the reduced
42 ** argument, z, and the scale factor, s, as:
43 **
44 ** lnb*x = s*ln2 + z, |z| < ln2/2
45 **
46 ** where b is equal to e or 10. If |x| > 2^16, UX_EXP_REDUCE returns a value of
47 ** s and z that will force underflow or overflow in the pack routine.
48 */
49
50 #if !defined(UX_EXP_REDUCE)
51 # define UX_EXP_REDUCE __INTERNAL_NAME(ux_exp_reduce__)
52 #endif
53
54 static WORD
UX_EXP_REDUCE(UX_FLOAT * orig_argument,UX_FLOAT * reduced_argument,UX_FRACTION_DIGIT_TYPE * constants)55 UX_EXP_REDUCE(UX_FLOAT * orig_argument, UX_FLOAT * reduced_argument,
56 UX_FRACTION_DIGIT_TYPE * constants )
57 {
58 U_WORD shift, reduce_constant_exp;
59 UX_SIGN_TYPE sign;
60 UX_EXPONENT_TYPE exponent, scale_exponent;
61 UX_FRACTION_DIGIT_TYPE scale, msd, lsd;
62 UX_FLOAT ux_scale, tmp;
63
64 exponent = G_UX_EXPONENT(orig_argument);
65 sign = G_UX_SIGN(orig_argument);
66
67 reduce_constant_exp = constants[2];
68 if ( (UX_UNSIGNED_EXPONENT_TYPE) (exponent + 1 - reduce_constant_exp) > 18)
69 { /* Either no reduction is necessary, or exponent > 17 */
70
71 scale = 0;
72 UX_COPY(orig_argument, reduced_argument);
73 if (exponent > 0)
74 { /* exponent > 17, force underflow or overflow */
75 P_UX_EXPONENT(reduced_argument, -128);
76 scale = sign ? UX_UNDERFLOW_EXPONENT : UX_OVERFLOW_EXPONENT;
77 }
78 return scale;
79 }
80
81 /*
82 ** Given an input argument of the form x = 2^n*f, we want to compute
83 ** lnb*x = scale*ln2 + z, |z| <= ln2/2. Or equivalently, scale =
84 ** nint(x*lnb/ln2) and z = scale*ln2. Suppose, the number of bits in a
85 ** fraction digit is k, and we define K = 2^k. Further suppose that F is
86 ** the high k-1 bits of f and L is the high k bits of lnb/ln2. Then
87 **
88 ** scale = nint(x*lnb/ln2)
89 ** = nint[ 2^n*f*(lnb/ln2) ]
90 ** ~ nint{ 2^n*[F/(K/2)]*[L/(K/4)] }
91 ** = nint{ 2^(n+3)*(F*L)/K^2 }
92 ** = nint{ 2^(n+3)*[ Hi(F*L)*K + Lo(F*L) ]/K^2 }
93 ** ~ nint{ 2^(n+3)*H(F*L)/K }
94 ** = nint{ H(F*L)/2^(k - n - 3) }
95 **
96 ** so that we can compute scale by computing the high k bits of F*L and
97 ** "shifting" right k-n-3 bits. Since we want to multiply by scale,
98 ** we actually mask out the low order bits after rounding. Note that
99 ** since we took only the high k-1 bits of f, there is no possibility
100 ** of a carry out on the round.
101 */
102
103 msd = G_UX_MSD(orig_argument) >> 1;
104 UMULH( msd, constants[0], scale);
105 shift = (BITS_PER_UX_FRACTION_DIGIT_TYPE - 3) - exponent;
106 scale += SET_BIT(shift - 1);
107 scale &= -SET_BIT(shift);
108
109 /*
110 ** Now compute (x - scale*high_bits_of_ln2) - scale*low_bits_of_ln2
111 ** Begin by make sure scale is normalized. It could have at most two
112 ** leading zeros
113 */
114
115 while ((UX_SIGNED_FRACTION_DIGIT_TYPE) scale > 0)
116 {
117 scale += scale;
118 shift++;
119 }
120
121
122 /*
123 ** Get scale*high_bits_of_ln2 and subtract from x. Theres a small
124 ** complication that needs to be dealt with here: When computing
125 ** scale*high_bits_of_ln2, it may be unnormalized by one bit. Which
126 ** causes x to be right shifted one bit on the subtraction, there by
127 ** losing the last bit of x. Most of the time, this is unimportant.
128 ** However, for very large arguments with a non-zero lsb, this results
129 ** in very large error in the final answer, so we need to normalize
130 ** scale*high_bits_of_ln2 before subtracting
131 */
132
133 scale_exponent = BITS_PER_UX_FRACTION_DIGIT_TYPE - shift;
134 EXTENDED_DIGIT_MULTIPLY(scale, constants[1], msd, lsd);
135 exponent = scale_exponent;
136 if (((UX_SIGNED_FRACTION_DIGIT_TYPE) msd) > 0)
137 {
138 exponent--;
139 msd = (msd + msd) + (lsd >> (BITS_PER_UX_FRACTION_DIGIT_TYPE - 1));
140 lsd += lsd;
141 }
142
143 /* adjust the product exponent by the exponent of the constant */
144 UX_SET_SIGN_EXP_MSD(&tmp, sign, exponent + reduce_constant_exp, msd);
145 P_UX_FRACTION_DIGIT(&tmp, 1, lsd);
146 ADDSUB(orig_argument, &tmp, SUB, &tmp);
147
148 /* scale*low_bits_of_ln2 and subtract from x - scale*high_bits_of_ln2 */
149
150 UX_SET_SIGN_EXP_MSD(&ux_scale, sign, scale_exponent, scale);
151 MULTIPLY(&ux_scale, (UX_FLOAT *)&constants[3], reduced_argument);
152 ADDSUB(&tmp, reduced_argument, SUB | NO_NORMALIZATION, reduced_argument);
153
154 scale >>= shift;
155 scale = (sign) ? -scale : scale;
156 return scale;
157 }
158
159 /*
160 ** UX_EXP_COMMON is the unpacked interface to routine that will compute b^x for
161 ** b = e or 10. It calls UX_EXP_REDUCE to get the exponent and reduced
162 ** argument and then evaluates the exp polynomial
163 */
164
165 #if !defined(UX_EXP_COMMON)
166 # define UX_EXP_COMMON __INTERNAL_NAME(ux_exp_common__)
167 #endif
168
169 void
UX_EXP_COMMON(UX_FLOAT * unpacked_argument,UX_FLOAT * unpacked_result,UX_FRACTION_DIGIT_TYPE * constant_table)170 UX_EXP_COMMON( UX_FLOAT * unpacked_argument, UX_FLOAT * unpacked_result,
171 UX_FRACTION_DIGIT_TYPE * constant_table)
172 {
173 UX_EXPONENT_TYPE scale;
174 UX_FLOAT reduced_argument;
175
176 /* Get reduced argument */
177 scale = UX_EXP_REDUCE(unpacked_argument, &reduced_argument, constant_table);
178
179 /* Compute e^reduced_argument */
180
181 EVALUATE_RATIONAL(
182 &reduced_argument,
183 (FIXED_128 *) &constant_table[EXP_COEF_INDEX],
184 constant_table[EXP_DEGREE_INDEX],
185 NUMERATOR_FLAGS(STANDARD),
186 unpacked_result);
187
188 /* Scale e^reduced_argument */
189
190 UX_INCR_EXPONENT(unpacked_result, scale);
191 }
192
193 /*
194 ** UX_EXP is the unpacked interface to the exponential routine. It calls
195 ** UX_EXP_COMMONN routine to compute its result.
196 */
197
198 #if !defined(UX_EXP)
199 # define UX_EXP __INTERNAL_NAME(ux_exp__)
200 #endif
201
202 void
UX_EXP(UX_FLOAT * unpacked_argument,UX_FLOAT * unpacked_result)203 UX_EXP( UX_FLOAT * unpacked_argument, UX_FLOAT * unpacked_result)
204 {
205 UX_EXP_COMMON(unpacked_argument, unpacked_result,
206 EXP_CONSTANT_TABLE_ADDRESS);
207 }
208
209
210 /*
211 ** F_EXP_NAME is the user level packed x-float exp routine
212 */
213
214 #undef F_ENTRY_NAME
215 #define F_ENTRY_NAME F_EXP_NAME
216
X_X_PROTO(F_ENTRY_NAME,packed_result,packed_argument)217 X_X_PROTO(F_ENTRY_NAME, packed_result, packed_argument)
218 {
219 WORD fp_class;
220 UX_FLOAT unpacked_argument, unpacked_result;
221 EXCEPTION_INFO_DECL
222 DECLARE_X_FLOAT(packed_result)
223
224 INIT_EXCEPTION_INFO;
225 fp_class = UNPACK(
226 PASS_ARG_X_FLOAT(packed_argument),
227 & unpacked_argument,
228 EXP_CLASS_TO_ACTION_MAP,
229 PASS_RET_X_FLOAT(packed_result)
230 OPT_EXCEPTION_INFO );
231
232 if (0 > fp_class)
233 RETURN_X_FLOAT(packed_result);
234
235 UX_EXP( &unpacked_argument, &unpacked_result);
236
237 PACK(
238 &unpacked_result,
239 PASS_RET_X_FLOAT(packed_result),
240 EXP_UNDERFLOW,
241 EXP_OVERFLOW
242 OPT_EXCEPTION_INFO );
243
244 RETURN_X_FLOAT(packed_result);
245
246 }
247
248
249 /*
250 ** F_EXPM1_NAME is the packed x-float expm1 function. F_EXPM1_NAME exam the
251 ** size of the reduced argument. If it is small enough, a direct polynomial
252 ** evaluation is perform. Otherwise, UX_EXP computes expm1(x) = exp(x) - 1
253 */
254
255 #undef F_ENTRY_NAME
256 #define F_ENTRY_NAME F_EXPM1_NAME
257
X_X_PROTO(F_ENTRY_NAME,packed_result,packed_argument)258 X_X_PROTO(F_ENTRY_NAME, packed_result, packed_argument)
259 {
260 WORD fp_class;
261 UX_EXPONENT_TYPE scale;
262 UX_FLOAT unpacked_argument, unpacked_result, reduced_argument, one;
263 UX_FRACTION_DIGIT_TYPE * constants;
264 EXCEPTION_INFO_DECL
265 DECLARE_X_FLOAT(packed_result)
266
267 INIT_EXCEPTION_INFO;
268 fp_class = UNPACK(
269 PASS_ARG_X_FLOAT(packed_argument),
270 & unpacked_argument,
271 EXPM1_CLASS_TO_ACTION_MAP,
272 PASS_RET_X_FLOAT(packed_result)
273 OPT_EXCEPTION_INFO);
274
275 if (0 > fp_class)
276 RETURN_X_FLOAT(packed_result);
277
278 constants = EXP_CONSTANT_TABLE_ADDRESS;
279 scale = UX_EXP_REDUCE( &unpacked_argument, &reduced_argument, constants);
280 if (scale == 0)
281 {
282 /*
283 ** abs(reduced_argument) < ln2/2. computing expm1(x) as
284 ** exp(x) - 1, could result in a serve loss of significance,
285 ** so use a direct polynomial evaluation instead. We use the
286 ** low EXP_COEF_ARRAY_DEGREE - 1 terms of the exp polynomial.
287 ** This has the side effect that the exponent field of the
288 ** result is 1 to small.
289 */
290 EVALUATE_RATIONAL(
291 &reduced_argument,
292 (FIXED_128 *) &constants[EXP_COEF_INDEX],
293 constants[EXP_DEGREE_INDEX] - 1,
294 NUMERATOR_FLAGS(POST_MULTIPLY),/* Post multiply by x */
295 &unpacked_result);
296 UX_INCR_EXPONENT(&unpacked_result, 1);
297 }
298 else
299 {
300 /*
301 ** Compute expm1(x) = exp(x) - 1. Since |scale| >= 1,
302 ** exp(x) <= 1/sqrt(2) and exp(x) >= sqrt(2)
303 */
304 EVALUATE_RATIONAL(
305 &reduced_argument,
306 (FIXED_128 *) &constants[EXP_COEF_INDEX],
307 constants[EXP_DEGREE_INDEX],
308 NUMERATOR_FLAGS(STANDARD),
309 &unpacked_result);
310 UX_INCR_EXPONENT(&unpacked_result, scale);
311
312 ADDSUB(
313 &unpacked_result,
314 UX_ONE,
315 SUB | NO_NORMALIZATION | MAGNITUDE_ONLY,
316 &unpacked_result
317 );
318 }
319
320 PACK(
321 &unpacked_result,
322 PASS_RET_X_FLOAT(packed_result),
323 NOT_USED,
324 EXPM1_OVERFLOW
325 OPT_EXCEPTION_INFO);
326
327 RETURN_X_FLOAT(packed_result);
328
329 }
330
331 /*
332 ** F_EXP10_NAME is the user level packed x-float exp10 routine
333 */
334
335 #undef F_ENTRY_NAME
336 #define F_ENTRY_NAME F_EXP10_NAME
337
X_X_PROTO(F_ENTRY_NAME,packed_result,packed_argument)338 X_X_PROTO(F_ENTRY_NAME, packed_result, packed_argument)
339 {
340 WORD fp_class;
341 UX_FLOAT unpacked_argument, unpacked_result;
342 EXCEPTION_INFO_DECL
343 DECLARE_X_FLOAT(packed_result)
344
345 INIT_EXCEPTION_INFO;
346 fp_class = UNPACK(
347 PASS_ARG_X_FLOAT(packed_argument),
348 & unpacked_argument,
349 EXP_CLASS_TO_ACTION_MAP,
350 PASS_RET_X_FLOAT(packed_result)
351 OPT_EXCEPTION_INFO );
352
353 if (0 > fp_class)
354 RETURN_X_FLOAT(packed_result);
355
356 UX_EXP_COMMON( &unpacked_argument, &unpacked_result,
357 EXP10_CONSTANT_TABLE_ADDRESS);
358
359 PACK(
360 &unpacked_result,
361 PASS_RET_X_FLOAT(packed_result),
362 EXP_UNDERFLOW,
363 EXP_OVERFLOW
364 OPT_EXCEPTION_INFO );
365
366 RETURN_X_FLOAT(packed_result);
367
368 }
369
370 /*
371 ** UX_HYPERBOLIC is the core processing for hyperbolic function of an unpacked
372 ** argument. Depending on the evaluation flags to UX_HYPERBOLIC, it computes
373 ** one of sinh, cosh, sinhcosh or tanh. In order to promote "efficiency" and
374 ** clarity, then evaluation flags are divided into three separate fields
375 ** containing (somewhat redundant) evaluation information. One field contains
376 ** the function to be evaluated (SINH, COSH, SINHCOSH or TANH); one field
377 ** contains the appropriate evaluation flags for EVALUATION_RATIONAL; and
378 ** one field containing the opcode to be used by the ADDSUB routine
379 */
380
381 #define __FLAGS(i,w,p) (((i) >> (p)) & MAKE_MASK(w,0))
382
383 #define EVAL_RATIONAL_POS 0
384 #define EVAL_RATIONAL_WIDTH (2*NUM_DEN_FIELD_WIDTH + 3)
385 #define EVAL_RATIONAL_FLAGS(i) __FLAGS(i,EVAL_RATIONAL_WIDTH,EVAL_RATIONAL_POS)
386
387 #define SINH_EVAL ( NUMERATOR_FLAGS( SQUARE_TERM | POST_MULTIPLY ) | SKIP)
388 #define COSH_EVAL ( SKIP | DENOMINATOR_FLAGS(SQUARE_TERM))
389 #define TANH_EVAL ( NUMERATOR_FLAGS( SQUARE_TERM | POST_MULTIPLY ) | \
390 DENOMINATOR_FLAGS(SQUARE_TERM) )
391 #define SINHCOSH_EVAL ( TANH_EVAL | NO_DIVIDE )
392
393 #define ADDSUB_POS (EVAL_RATIONAL_WIDTH + EVAL_RATIONAL_POS)
394 #define ADDSUB_WIDTH 2
395 #define ADDSUB_FLAGS(i) __FLAGS(i, ADDSUB_WIDTH, ADDSUB_POS)
396
397 #define FUNC_CODE_POS (ADDSUB_POS + ADDSUB_WIDTH)
398 #define SINH (1 << FUNC_CODE_POS)
399 #define COSH (2 << FUNC_CODE_POS)
400 #define SINHCOSH (4 << FUNC_CODE_POS)
401 #define TANH (8 << FUNC_CODE_POS)
402
403 #define EVAL_FLAGS(f,r,a) ( (f) | ((r) << EVAL_RATIONAL_POS) | \
404 ((a) << ADDSUB_POS))
405
406 #define UX_HYPERBOLIC __INTERNAL_NAME(ux_hyperbolic__)
407
408 void
UX_HYPERBOLIC(UX_FLOAT * unpacked_argument,WORD evaluation_flags,UX_FLOAT * unpacked_result)409 UX_HYPERBOLIC( UX_FLOAT * unpacked_argument, WORD evaluation_flags,
410 UX_FLOAT * unpacked_result)
411 {
412 UX_EXPONENT_TYPE scale;
413 UX_SIGN_TYPE sign;
414 UX_FLOAT reduced_argument, tmp[2];
415
416 /*
417 ** save sign of input and its absolute value before performing
418 ** argument reduction, x = I*ln2 + z, |z| < ln2/2. Note that
419 ** if this is a cosh(x) evaluation, we treat the sign as positive.
420 */
421
422 sign = G_UX_SIGN(unpacked_argument);
423 P_UX_SIGN(unpacked_argument, 0);
424 sign = ( evaluation_flags & COSH ) ? 0 : sign;
425 scale = UX_EXP_REDUCE( unpacked_argument, &reduced_argument,
426 EXP_CONSTANT_TABLE_ADDRESS);
427
428 /*
429 ** if scale == 0, then abs(x) < ln2/2 ==> sinh(x) or tanh(x) may have
430 ** a loss of significance if computed via the definition, so compute
431 ** by polynomial instead. Otherwise, we compute exp(z) and
432 ** exp(-z) as cosh(z) + sinh(z) and cosh(z) - sinh(z) respectively.
433 ** So, if scale == 0, used the passed in evaluation flags, otherwise
434 ** Force a SINHCOSH evaluation.
435 */
436
437 EVALUATE_RATIONAL(
438 &reduced_argument,
439 SINHCOSH_COEF_ARRAY,
440 SINHCOSH_COEF_ARRAY_DEGREE,
441 (scale == 0) ? EVAL_RATIONAL_FLAGS(evaluation_flags) :
442 SINHCOSH_EVAL,
443 unpacked_result );
444
445 if (scale)
446 {
447 /*
448 ** We want to compute sinh(x)/cosh(x) = (exp(x) -/+ exp(-x))/2.
449 ** Begin by computing exp(z) and exp(-z) and then scale them
450 ** to get exp(x)/2 and exp(-x)/2.
451 */
452 ADDSUB(
453 &unpacked_result[1], /* cosh(z) */
454 &unpacked_result[0], /* sinh(z) */
455 ADD_SUB | NO_NORMALIZATION,
456 &tmp[0] /* exp(z):exp(-z)*/
457 );
458
459 UX_INCR_EXPONENT(&tmp[0], (scale - 1));
460 UX_DECR_EXPONENT(&tmp[1], (scale + 1));
461
462 /*
463 ** Now add/sub exp(x)/2 and exp(-x)/2 to get sinh/cosh, if this
464 ** is a tanh evaluation, do the divide
465 */
466
467 ADDSUB(
468 &tmp[0], /* exp(x)/2 */
469 &tmp[1], /* exp(-x)/2 */
470 ADDSUB_FLAGS(evaluation_flags) | MAGNITUDE_ONLY | NO_NORMALIZATION,
471 &unpacked_result[0] /* sinh(x)/cosh(x) */
472 );
473
474 if (evaluation_flags & TANH)
475 DIVIDE(&unpacked_result[0], &unpacked_result[1], FULL_PRECISION,
476 &unpacked_result[0]);
477 }
478
479 P_UX_SIGN(unpacked_result, sign);
480 }
481
482 /*
483 ** C_UX_HYPERBOLIC is the common processing routine for the hyperbolic
484 ** routines: sinh, cosh, sinhcosh and tanh. It unpacks the input argument,
485 ** calls UX_HYPERBOLIC to computes sinh, cosh, sinhcosh or tanh, and packs the
486 ** results.
487 */
488
489 #define C_UX_HYPERBOLIC __INTERNAL_NAME(C_ux_hyperbolic__)
490
491 static void
C_UX_HYPERBOLIC(_X_FLOAT * packed_result,_X_FLOAT * packed_argument,U_WORD const * class_to_action_map,WORD evaluation_flags,WORD overflow_code OPT_EXCEPTION_INFO_DECLARATION)492 C_UX_HYPERBOLIC( _X_FLOAT * packed_result, _X_FLOAT * packed_argument,
493 U_WORD const * class_to_action_map, WORD evaluation_flags,
494 WORD overflow_code OPT_EXCEPTION_INFO_DECLARATION )
495 {
496 WORD fp_class;
497 UX_FLOAT unpacked_argument, unpacked_result[2];
498
499 fp_class = UNPACK(
500 packed_argument,
501 &unpacked_argument,
502 class_to_action_map,
503 &packed_result[0]
504 OPT_EXCEPTION_INFO_ARGUMENT );
505
506 if (0 > fp_class)
507 { /* If this is a SINHCOSH evaluation, write second result */
508 if (evaluation_flags & SINHCOSH)
509 {
510 (void) UNPACK(
511 packed_argument,
512 &unpacked_argument,
513 COSH_CLASS_TO_ACTION_MAP,
514 &packed_result[1]
515 OPT_EXCEPTION_INFO_ARGUMENT );
516 }
517 return;
518 }
519
520 UX_HYPERBOLIC(
521 &unpacked_argument,
522 evaluation_flags,
523 &unpacked_result[0]);
524
525 PACK(
526 &unpacked_result[0],
527 packed_result,
528 NOT_USED,
529 overflow_code
530 OPT_EXCEPTION_INFO_ARGUMENT );
531
532 if (evaluation_flags & SINHCOSH)
533 /* This was a sinhcosh evaluation */
534 PACK(
535 &unpacked_result[1],
536 &packed_result[1],
537 NOT_USED,
538 COSH_OVERFLOW
539 OPT_EXCEPTION_INFO_ARGUMENT );
540 }
541
542 /*
543 ** F_SINH_NAME, F_COSH_NAME, F_SINHCOSH_NAME and F_TANH_NAME are the packed
544 ** x-float sinh, cosh, sinhcosh and tanh routines. Each of these routines
545 ** simply invokes the common routine C_UX_HYPERBOLIC to unpack its arguments,
546 ** compute the result and pack it.
547 */
548
549 #undef F_ENTRY_NAME
550 #define F_ENTRY_NAME F_SINH_NAME
551
X_X_PROTO(F_ENTRY_NAME,packed_result,packed_argument)552 X_X_PROTO(F_ENTRY_NAME, packed_result, packed_argument)
553 {
554 EXCEPTION_INFO_DECL
555 DECLARE_X_FLOAT(packed_result)
556
557 INIT_EXCEPTION_INFO;
558 C_UX_HYPERBOLIC(
559 PASS_RET_X_FLOAT(packed_result),
560 PASS_ARG_X_FLOAT(packed_argument),
561 SINH_CLASS_TO_ACTION_MAP,
562 EVAL_FLAGS( SINH, SINH_EVAL, SUB ),
563 PACKED_ARG_IS_NEG(packed_argument) ? SINH_NEG_OVERFLOW : SINH_OVERFLOW
564 OPT_EXCEPTION_INFO );
565
566 RETURN_X_FLOAT(packed_result);
567
568 }
569
570 #undef F_ENTRY_NAME
571 #define F_ENTRY_NAME F_COSH_NAME
572
X_X_PROTO(F_ENTRY_NAME,packed_result,packed_argument)573 X_X_PROTO(F_ENTRY_NAME, packed_result, packed_argument)
574 {
575 EXCEPTION_INFO_DECL
576 DECLARE_X_FLOAT(packed_result)
577
578 INIT_EXCEPTION_INFO;
579 C_UX_HYPERBOLIC(
580 PASS_RET_X_FLOAT(packed_result),
581 PASS_ARG_X_FLOAT(packed_argument),
582 COSH_CLASS_TO_ACTION_MAP,
583 EVAL_FLAGS( COSH, COSH_EVAL, ADD),
584 COSH_OVERFLOW
585 OPT_EXCEPTION_INFO );
586
587 RETURN_X_FLOAT(packed_result);
588
589 }
590
591 #undef F_ENTRY_NAME
592 #define F_ENTRY_NAME F_SINHCOSH_NAME
593
RR_X_PROTO(F_ENTRY_NAME,packed_result0,packed_result1,packed_argument)594 RR_X_PROTO(F_ENTRY_NAME, packed_result0, packed_result1, packed_argument)
595 {
596 EXCEPTION_INFO_DECL
597 _X_FLOAT packed_result[2];
598
599 INIT_EXCEPTION_INFO;
600 C_UX_HYPERBOLIC(
601 packed_result, /*PASS_RET_X_FLOAT(packed_result)*/
602 PASS_ARG_X_FLOAT(packed_argument),
603 SINH_CLASS_TO_ACTION_MAP,
604 EVAL_FLAGS( SINHCOSH, SINHCOSH_EVAL, SUB_ADD),
605 PACKED_ARG_IS_NEG(packed_argument) ? SINH_NEG_OVERFLOW : SINH_OVERFLOW
606 OPT_EXCEPTION_INFO );
607
608 *packed_result0 = packed_result[0];
609 *packed_result1 = packed_result[1];
610
611 }
612
613 #undef F_ENTRY_NAME
614 #define F_ENTRY_NAME F_TANH_NAME
615
X_X_PROTO(F_ENTRY_NAME,packed_result,packed_argument)616 X_X_PROTO(F_ENTRY_NAME, packed_result, packed_argument)
617 {
618 EXCEPTION_INFO_DECL
619 DECLARE_X_FLOAT(packed_result)
620
621 C_UX_HYPERBOLIC(
622 PASS_RET_X_FLOAT(packed_result),
623 PASS_ARG_X_FLOAT(packed_argument),
624 TANH_CLASS_TO_ACTION_MAP,
625 EVAL_FLAGS( TANH, TANH_EVAL, SUB_ADD),
626 NOT_USED
627 OPT_EXCEPTION_INFO );
628 RETURN_X_FLOAT(packed_result);
629 }
630
631
632 #if defined(MAKE_INCLUDE)
633
634 @divert -append divertText
635
636 precision = ceil(UX_PRECISION/8) + 4;
637
638 # undef TABLE_NAME
639
640 START_TABLE;
641
642 TABLE_COMMENT("exp class-to-action-mapping");
643 PRINT_CLASS_TO_ACTION_TBL_DEF( "EXP_CLASS_TO_ACTION_MAP\t");
644 PRINT_64_TBL_ITEM( CLASS_TO_ACTION_DISP(5) +
645 CLASS_TO_ACTION( F_C_SIG_NAN, RETURN_QUIET_NAN, 0) +
646 CLASS_TO_ACTION( F_C_QUIET_NAN, RETURN_VALUE, 0) +
647 CLASS_TO_ACTION( F_C_POS_INF, RETURN_ERROR, 3) +
648 CLASS_TO_ACTION( F_C_NEG_INF, RETURN_ERROR, 2) +
649 CLASS_TO_ACTION( F_C_POS_DENORM, RETURN_VALUE, 1) +
650 CLASS_TO_ACTION( F_C_NEG_DENORM, RETURN_VALUE, 1) +
651 CLASS_TO_ACTION( F_C_POS_ZERO, RETURN_VALUE, 1) +
652 CLASS_TO_ACTION( F_C_NEG_ZERO , RETURN_VALUE, 1) );
653
654 TABLE_COMMENT("expm1 class-to-action-mapping");
655 PRINT_CLASS_TO_ACTION_TBL_DEF( "EXPM1_CLASS_TO_ACTION_MAP");
656 PRINT_64_TBL_ITEM( CLASS_TO_ACTION_DISP(4) +
657 CLASS_TO_ACTION( F_C_SIG_NAN, RETURN_QUIET_NAN, 0) +
658 CLASS_TO_ACTION( F_C_QUIET_NAN, RETURN_VALUE, 0) +
659 CLASS_TO_ACTION( F_C_POS_INF, RETURN_VALUE, 0) +
660 CLASS_TO_ACTION( F_C_NEG_INF, RETURN_NEGATIVE, 1) +
661 CLASS_TO_ACTION( F_C_POS_DENORM, RETURN_VALUE, 0) +
662 CLASS_TO_ACTION( F_C_NEG_DENORM, RETURN_VALUE, 0) +
663 CLASS_TO_ACTION( F_C_POS_ZERO, RETURN_VALUE, 0) +
664 CLASS_TO_ACTION( F_C_NEG_ZERO , RETURN_VALUE, 0) );
665
666 TABLE_COMMENT("sinh class-to-action-mapping");
667 PRINT_CLASS_TO_ACTION_TBL_DEF( "SINH_CLASS_TO_ACTION_MAP");
668 PRINT_64_TBL_ITEM( CLASS_TO_ACTION_DISP(3) +
669 CLASS_TO_ACTION( F_C_QUIET_NAN, RETURN_VALUE, 0) +
670 CLASS_TO_ACTION( F_C_POS_INF, RETURN_VALUE, 0) +
671 CLASS_TO_ACTION( F_C_NEG_INF, RETURN_VALUE, 0) +
672 CLASS_TO_ACTION( F_C_POS_DENORM, RETURN_VALUE, 0) +
673 CLASS_TO_ACTION( F_C_NEG_DENORM, RETURN_VALUE, 0) +
674 CLASS_TO_ACTION( F_C_POS_ZERO, RETURN_VALUE, 0) +
675 CLASS_TO_ACTION( F_C_NEG_ZERO , RETURN_VALUE, 0) );
676
677 TABLE_COMMENT("cosh class-to-action-mapping");
678 PRINT_CLASS_TO_ACTION_TBL_DEF( "COSH_CLASS_TO_ACTION_MAP");
679 PRINT_64_TBL_ITEM( CLASS_TO_ACTION_DISP(2) +
680 CLASS_TO_ACTION( F_C_SIG_NAN, RETURN_QUIET_NAN, 0) +
681 CLASS_TO_ACTION( F_C_QUIET_NAN, RETURN_VALUE, 0) +
682 CLASS_TO_ACTION( F_C_POS_INF, RETURN_VALUE, 0) +
683 CLASS_TO_ACTION( F_C_NEG_INF, RETURN_NEGATIVE, 0) +
684 CLASS_TO_ACTION( F_C_POS_DENORM, RETURN_VALUE, 1) +
685 CLASS_TO_ACTION( F_C_NEG_DENORM, RETURN_VALUE, 1) +
686 CLASS_TO_ACTION( F_C_POS_ZERO, RETURN_VALUE, 1) +
687 CLASS_TO_ACTION( F_C_NEG_ZERO , RETURN_VALUE, 1) );
688
689
690 TABLE_COMMENT("tanh class-to-action-mapping");
691 PRINT_CLASS_TO_ACTION_TBL_DEF( "TANH_CLASS_TO_ACTION_MAP");
692 PRINT_64_TBL_ITEM( CLASS_TO_ACTION_DISP(1) +
693 CLASS_TO_ACTION( F_C_SIG_NAN, RETURN_QUIET_NAN, 0) +
694 CLASS_TO_ACTION( F_C_QUIET_NAN, RETURN_VALUE, 0) +
695 CLASS_TO_ACTION( F_C_POS_INF, RETURN_VALUE, 1) +
696 CLASS_TO_ACTION( F_C_NEG_INF, RETURN_NEGATIVE, 1) +
697 CLASS_TO_ACTION( F_C_POS_DENORM, RETURN_VALUE, 0) +
698 CLASS_TO_ACTION( F_C_NEG_DENORM, RETURN_VALUE, 0) +
699 CLASS_TO_ACTION( F_C_POS_ZERO, RETURN_VALUE, 0) +
700 CLASS_TO_ACTION( F_C_NEG_ZERO , RETURN_VALUE, 0) );
701
702 TABLE_COMMENT("Data for the class to action mappings");
703 PRINT_U_TBL_ITEM( /* data 1 */ ONE );
704 PRINT_U_TBL_ITEM( /* data 2 */ EXP_OF_NEG_INF );
705 PRINT_U_TBL_ITEM( /* data 3 */ EXP_OF_INF );
706
707 /*
708 ** Create the "table" of exp constants. The table includes the constants
709 ** for the argument reduction, the degree of the polynomial and the
710 ** polynomial coefficients.
711 */
712
713 TABLE_COMMENT("Constant structure for exp based evaluations");
714 PRINT_UX_FRACTION_DIGIT_TBL_ADEF("EXP_CONSTANT_TABLE_ADDRESS");
715
716 save_precision = precision;
717 precision = ceil(2*UX_PRECISION/8);
718
719 ln2 = log(2);
720
721 precision = save_precision;
722
723 TABLE_COMMENT("High digits of 1/ln2, ln2 and binary exponent of ln2");
724 exp_cons_base_offset = MP_BIT_OFFSET;
725 ln2_hi = bround(ln2, BITS_PER_UX_FRACTION_DIGIT_TYPE);
726 tmp = bround(bldexp(1/ln2, BITS_PER_UX_FRACTION_DIGIT_TYPE - 2),
727 BITS_PER_UX_FRACTION_DIGIT_TYPE);
728 PRINT_UX_FRACTION_DIGIT_TBL_ITEM(tmp);
729 PRINT_UX_FRACTION_DIGIT_TBL_ITEM(
730 bldexp(ln2, BITS_PER_UX_FRACTION_DIGIT_TYPE) );
731 PRINT_UX_FRACTION_DIGIT_TBL_ITEM(0);
732
733 TABLE_COMMENT("ln2_lo = ln2 - ln2_hi in unpacked form");
734 PRINT_UX_TBL_ITEM( ln2 - ln2_hi);
735
736 /*
737 ** Compute polynomial coefficient for exp and expm1. Get coefficients
738 ** for expm1 and prepend a 1 to the front of the list
739 */
740
__expm1(x)741 function __expm1(x)
742 {
743 if (x == 0)
744 return 1.;
745 else
746 return expm1(x)/x;
747 }
748
749 save_precision = precision;
750 precision = ceil(UX_PRECISION/8) + 8;
751
752 max_arg = ln2/2;
753
754 remes(REMES_FIND_POLYNOMIAL + REMES_RELATIVE_WEIGHT, -max_arg,
755 max_arg, __expm1, UX_PRECISION, °ree, &ux_rational_coefs);
756
757 precision = save_precision;
758
759 for (i = degree + 1; i > 0; /* NULL */ )
760 ux_rational_coefs[i] = ux_rational_coefs[--i];
761 ux_rational_coefs[0] = 1;
762
763 #define __INDEX(z,b) ((z - b)/BITS_PER_UX_FRACTION_DIGIT_TYPE)
764 TABLE_COMMENT("Polynomial degree");
765 printf("#define EXP_DEGREE_INDEX\t\t%i\n",
766 __INDEX(MP_BIT_OFFSET, exp_cons_base_offset));
767 PRINT_UX_FRACTION_DIGIT_TBL_ITEM(degree+1);
768 TABLE_COMMENT("Fixed point coefficients for exp/expm1 evaluation");
769 printf("#define EXP_COEF_INDEX\t\t\t%i\n",
770 __INDEX(MP_BIT_OFFSET, exp_cons_base_offset));
771 print_ux_rational_coefs(degree + 1, 0, 0);
772
773 TABLE_COMMENT("1 in unpacked format");
774 PRINT_UX_TBL_ADEF_ITEM( "UX_ONE\t\t\t", 1);
775
776 /*
777 ** Create the "table" of exp10 constants. The layout is the same as for
778 ** the exp constants.
779 */
780
781
782 TABLE_COMMENT("Constant structure for exp10 based evaluations");
783 PRINT_UX_FRACTION_DIGIT_TBL_ADEF("EXP10_CONSTANT_TABLE_ADDRESS");
784
785 save_precision = precision;
786 precision = ceil(2*UX_PRECISION/8);
787
788 ln2_ov_ln10 = log(2)/log(10);
789
790 precision = save_precision;
791
792 TABLE_COMMENT(
793 "High digits of ln10/ln2, ln2/ln10 and binary exponent of ln2/ln10");
794 exp_cons_base_offset = MP_BIT_OFFSET;
795 ln2_ov_ln10_hi = bround(ln2_ov_ln10, BITS_PER_UX_FRACTION_DIGIT_TYPE);
796 tmp = bround(bldexp(1/ln2_ov_ln10, BITS_PER_UX_FRACTION_DIGIT_TYPE - 2),
797 BITS_PER_UX_FRACTION_DIGIT_TYPE);
798 PRINT_UX_FRACTION_DIGIT_TBL_ITEM(tmp);
799 PRINT_UX_FRACTION_DIGIT_TBL_ITEM(
800 bldexp(ln2_ov_ln10, BITS_PER_UX_FRACTION_DIGIT_TYPE + 1) );
801 PRINT_UX_FRACTION_DIGIT_TBL_ITEM(-1);
802
803 TABLE_COMMENT("ln2_ov_ln10_lo = ln2 - ln2_ov_ln10__hi in unpacked form");
804 PRINT_UX_TBL_ITEM( ln2_ov_ln10 - ln2_ov_ln10_hi);
805
806 /*
807 ** Compute polynomial coefficient for exp10.
808 */
809
__exp10(x)810 function __exp10(x)
811 {
812 return exp(x*log(10));
813 }
814
815 save_precision = precision;
816 precision = ceil(UX_PRECISION/8) + 8;
817
818 max_arg = ln2_ov_ln10/2;
819
820 remes(REMES_FIND_POLYNOMIAL + REMES_RELATIVE_WEIGHT, -max_arg,
821 max_arg, __exp10, UX_PRECISION, °ree, &ux_rational_coefs);
822
823 precision = save_precision;
824
825 TABLE_COMMENT("Polynomial degree");
826 PRINT_UX_FRACTION_DIGIT_TBL_ITEM(degree);
827 TABLE_COMMENT("Fixed point coefficients for exp10 evaluation");
828 print_ux_rational_coefs(degree, 0, 0);
829
830 /*
831 ** Now get sinh and cosh coefficients in the same array
832 */
833
__cosh(x)834 function __cosh(x) { return cosh(x); }
835
__sinh(x)836 function __sinh(x)
837 {
838 if (x == 0)
839 return 1.;
840 else
841 return sinh(x)/x;
842 }
843
844 save_precision = precision;
845 precision = ceil(UX_PRECISION/8) + 8;
846
847 max_arg = ln2/2;
848
849 remes(REMES_FIND_POLYNOMIAL + REMES_RELATIVE_WEIGHT + REMES_SQUARE_ARG,
850 0, max_arg, __sinh, UX_PRECISION, °ree, &ux_rational_coefs);
851
852 remes(REMES_FIND_POLYNOMIAL + REMES_RELATIVE_WEIGHT + REMES_SQUARE_ARG,
853 0, max_arg, __cosh, UX_PRECISION, &tmp_degree, &tmp_coefs);
854
855 for (i = 0; i <= tmp_degree; i++)
856 ux_rational_coefs[i + degree + 1] = tmp_coefs[i];
857
858 TABLE_COMMENT("Fixed point coefficients for sinh/cosh evaluation");
859 PRINT_FIXED_128_TBL_ADEF("SINHCOSH_COEF_ARRAY\t");
860 degree = print_ux_rational_coefs(degree, tmp_degree, 0);
861 PRINT_WORD_DEF("SINHCOSH_COEF_ARRAY_DEGREE", degree);
862
863
864 END_TABLE;
865
866 @end_divert
867 @eval my $tableText; \
868 my $outText = MphocEval( GetStream( "divertText" ) ); \
869 my $defineText = Egrep( "#define", $outText, \$tableText ); \
870 $outText = "$tableText\n\n$defineText"; \
871 my $headerText = GetHeaderText( STR(BUILD_FILE_NAME), \
872 "Definitions and constants exponential" . \
873 " and hyperbolic routines", __FILE__ ); \
874 print "$headerText\n\n$outText\n";
875
876
877 #endif
878