Home | History | Annotate | Line # | Download | only in fpe
fpu_exp.c revision 1.5.2.1
      1  1.5.2.1    yamt /*	$NetBSD: fpu_exp.c,v 1.5.2.1 2014/05/22 11:39:55 yamt Exp $	*/
      2      1.1  briggs 
      3      1.1  briggs /*
      4      1.1  briggs  * Copyright (c) 1995  Ken Nakata
      5      1.1  briggs  *	All rights reserved.
      6      1.1  briggs  *
      7      1.1  briggs  * Redistribution and use in source and binary forms, with or without
      8      1.1  briggs  * modification, are permitted provided that the following conditions
      9      1.1  briggs  * are met:
     10      1.1  briggs  * 1. Redistributions of source code must retain the above copyright
     11      1.1  briggs  *    notice, this list of conditions and the following disclaimer.
     12      1.1  briggs  * 2. Redistributions in binary form must reproduce the above copyright
     13      1.1  briggs  *    notice, this list of conditions and the following disclaimer in the
     14      1.1  briggs  *    documentation and/or other materials provided with the distribution.
     15      1.1  briggs  * 3. Neither the name of the author nor the names of its contributors
     16      1.1  briggs  *    may be used to endorse or promote products derived from this software
     17      1.1  briggs  *    without specific prior written permission.
     18      1.1  briggs  *
     19      1.1  briggs  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
     20      1.1  briggs  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
     21      1.1  briggs  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
     22      1.1  briggs  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
     23      1.1  briggs  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
     24      1.1  briggs  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
     25      1.1  briggs  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
     26      1.1  briggs  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
     27      1.1  briggs  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
     28      1.1  briggs  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
     29      1.1  briggs  * SUCH DAMAGE.
     30      1.1  briggs  *
     31      1.1  briggs  *	@(#)fpu_exp.c	10/24/95
     32      1.1  briggs  */
     33      1.2   lukem 
     34      1.2   lukem #include <sys/cdefs.h>
     35  1.5.2.1    yamt __KERNEL_RCSID(0, "$NetBSD: fpu_exp.c,v 1.5.2.1 2014/05/22 11:39:55 yamt Exp $");
     36  1.5.2.1    yamt 
     37  1.5.2.1    yamt #include <machine/ieee.h>
     38      1.1  briggs 
     39      1.1  briggs #include "fpu_emulate.h"
     40      1.1  briggs 
     41  1.5.2.1    yamt /* The number of items to terminate the Taylor expansion */
     42  1.5.2.1    yamt #define MAX_ITEMS	(2000)
     43  1.5.2.1    yamt 
     44      1.1  briggs /*
     45      1.1  briggs  * fpu_exp.c: defines fpu_etox(), fpu_etoxm1(), fpu_tentox(), and fpu_twotox();
     46      1.1  briggs  */
     47      1.1  briggs 
     48  1.5.2.1    yamt /*
     49  1.5.2.1    yamt  *                  x^2   x^3   x^4
     50  1.5.2.1    yamt  * exp(x) = 1 + x + --- + --- + --- + ...
     51  1.5.2.1    yamt  *                   2!    3!    4!
     52  1.5.2.1    yamt  */
     53  1.5.2.1    yamt static struct fpn *
     54  1.5.2.1    yamt fpu_etox_taylor(struct fpemu *fe)
     55  1.5.2.1    yamt {
     56  1.5.2.1    yamt 	struct fpn res;
     57  1.5.2.1    yamt 	struct fpn x;
     58  1.5.2.1    yamt 	struct fpn s0;
     59  1.5.2.1    yamt 	struct fpn *s1;
     60  1.5.2.1    yamt 	struct fpn *r;
     61  1.5.2.1    yamt 	uint32_t k;
     62  1.5.2.1    yamt 
     63  1.5.2.1    yamt 	CPYFPN(&x, &fe->fe_f2);
     64  1.5.2.1    yamt 	CPYFPN(&s0, &fe->fe_f2);
     65  1.5.2.1    yamt 
     66  1.5.2.1    yamt 	/* res := 1 + x */
     67  1.5.2.1    yamt 	fpu_const(&fe->fe_f1, FPU_CONST_1);
     68  1.5.2.1    yamt 	r = fpu_add(fe);
     69  1.5.2.1    yamt 	CPYFPN(&res, r);
     70  1.5.2.1    yamt 
     71  1.5.2.1    yamt 	k = 2;
     72  1.5.2.1    yamt 	for (; k < MAX_ITEMS; k++) {
     73  1.5.2.1    yamt 		/* s1 = s0 * x / k */
     74  1.5.2.1    yamt 		CPYFPN(&fe->fe_f1, &s0);
     75  1.5.2.1    yamt 		CPYFPN(&fe->fe_f2, &x);
     76  1.5.2.1    yamt 		r = fpu_mul(fe);
     77  1.5.2.1    yamt 
     78  1.5.2.1    yamt 		CPYFPN(&fe->fe_f1, r);
     79  1.5.2.1    yamt 		fpu_explode(fe, &fe->fe_f2, FTYPE_LNG, &k);
     80  1.5.2.1    yamt 		s1 = fpu_div(fe);
     81  1.5.2.1    yamt 
     82  1.5.2.1    yamt 		/* break if s1 is enough small */
     83  1.5.2.1    yamt 		if (ISZERO(s1))
     84  1.5.2.1    yamt 			break;
     85  1.5.2.1    yamt 		if (res.fp_exp - s1->fp_exp >= EXT_FRACBITS)
     86  1.5.2.1    yamt 			break;
     87  1.5.2.1    yamt 
     88  1.5.2.1    yamt 		/* s0 := s1 for next loop */
     89  1.5.2.1    yamt 		CPYFPN(&s0, s1);
     90  1.5.2.1    yamt 
     91  1.5.2.1    yamt 		/* res += s1 */
     92  1.5.2.1    yamt 		CPYFPN(&fe->fe_f2, s1);
     93  1.5.2.1    yamt 		CPYFPN(&fe->fe_f1, &res);
     94  1.5.2.1    yamt 		r = fpu_add(fe);
     95  1.5.2.1    yamt 		CPYFPN(&res, r);
     96  1.5.2.1    yamt 	}
     97  1.5.2.1    yamt 
     98  1.5.2.1    yamt 	CPYFPN(&fe->fe_f2, &res);
     99  1.5.2.1    yamt 	return &fe->fe_f2;
    100  1.5.2.1    yamt }
    101  1.5.2.1    yamt 
    102  1.5.2.1    yamt /*
    103  1.5.2.1    yamt  * exp(x)
    104  1.5.2.1    yamt  */
    105      1.1  briggs struct fpn *
    106      1.4     dsl fpu_etox(struct fpemu *fe)
    107      1.1  briggs {
    108  1.5.2.1    yamt 	struct fpn *fp;
    109  1.5.2.1    yamt 
    110  1.5.2.1    yamt 	if (ISNAN(&fe->fe_f2))
    111  1.5.2.1    yamt 		return &fe->fe_f2;
    112  1.5.2.1    yamt 	if (ISINF(&fe->fe_f2)) {
    113  1.5.2.1    yamt 		if (fe->fe_f2.fp_sign)
    114  1.5.2.1    yamt 			fpu_const(&fe->fe_f2, FPU_CONST_0);
    115  1.5.2.1    yamt 		return &fe->fe_f2;
    116  1.5.2.1    yamt 	}
    117  1.5.2.1    yamt 
    118  1.5.2.1    yamt 	if (fe->fe_f2.fp_sign == 0) {
    119  1.5.2.1    yamt 		/* exp(x) */
    120  1.5.2.1    yamt 		fp = fpu_etox_taylor(fe);
    121  1.5.2.1    yamt 	} else {
    122  1.5.2.1    yamt 		/* 1/exp(-x) */
    123  1.5.2.1    yamt 		fe->fe_f2.fp_sign = 0;
    124  1.5.2.1    yamt 		fp = fpu_etox_taylor(fe);
    125  1.5.2.1    yamt 
    126  1.5.2.1    yamt 		CPYFPN(&fe->fe_f2, fp);
    127  1.5.2.1    yamt 		fpu_const(&fe->fe_f1, FPU_CONST_1);
    128  1.5.2.1    yamt 		fp = fpu_div(fe);
    129  1.5.2.1    yamt 	}
    130  1.5.2.1    yamt 
    131  1.5.2.1    yamt 	return fp;
    132      1.1  briggs }
    133      1.1  briggs 
    134  1.5.2.1    yamt /*
    135  1.5.2.1    yamt  * exp(x) - 1
    136  1.5.2.1    yamt  */
    137      1.1  briggs struct fpn *
    138      1.4     dsl fpu_etoxm1(struct fpemu *fe)
    139      1.1  briggs {
    140  1.5.2.1    yamt 	struct fpn *fp;
    141  1.5.2.1    yamt 
    142  1.5.2.1    yamt 	fp = fpu_etox(fe);
    143  1.5.2.1    yamt 
    144  1.5.2.1    yamt 	CPYFPN(&fe->fe_f1, fp);
    145  1.5.2.1    yamt 	/* build a 1.0 */
    146  1.5.2.1    yamt 	fp = fpu_const(&fe->fe_f2, FPU_CONST_1);
    147  1.5.2.1    yamt 	fe->fe_f2.fp_sign = !fe->fe_f2.fp_sign;
    148  1.5.2.1    yamt 	/* fp = f2 - 1.0 */
    149  1.5.2.1    yamt 	fp = fpu_add(fe);
    150  1.5.2.1    yamt 
    151  1.5.2.1    yamt 	return fp;
    152      1.1  briggs }
    153      1.1  briggs 
    154  1.5.2.1    yamt /*
    155  1.5.2.1    yamt  * 10^x = exp(x * ln10)
    156  1.5.2.1    yamt  */
    157      1.1  briggs struct fpn *
    158      1.4     dsl fpu_tentox(struct fpemu *fe)
    159      1.1  briggs {
    160  1.5.2.1    yamt 	struct fpn *fp;
    161  1.5.2.1    yamt 
    162  1.5.2.1    yamt 	/* build a ln10 */
    163  1.5.2.1    yamt 	fp = fpu_const(&fe->fe_f1, FPU_CONST_LN_10);
    164  1.5.2.1    yamt 	/* fp = ln10 * f2 */
    165  1.5.2.1    yamt 	fp = fpu_mul(fe);
    166  1.5.2.1    yamt 
    167  1.5.2.1    yamt 	/* copy the result to the src opr */
    168  1.5.2.1    yamt 	CPYFPN(&fe->fe_f2, fp);
    169  1.5.2.1    yamt 
    170  1.5.2.1    yamt 	return fpu_etox(fe);
    171      1.1  briggs }
    172      1.1  briggs 
    173  1.5.2.1    yamt /*
    174  1.5.2.1    yamt  * 2^x = exp(x * ln2)
    175  1.5.2.1    yamt  */
    176      1.1  briggs struct fpn *
    177      1.4     dsl fpu_twotox(struct fpemu *fe)
    178      1.1  briggs {
    179  1.5.2.1    yamt 	struct fpn *fp;
    180  1.5.2.1    yamt 
    181  1.5.2.1    yamt 	/* build a ln2 */
    182  1.5.2.1    yamt 	fp = fpu_const(&fe->fe_f1, FPU_CONST_LN_2);
    183  1.5.2.1    yamt 	/* fp = ln2 * f2 */
    184  1.5.2.1    yamt 	fp = fpu_mul(fe);
    185  1.5.2.1    yamt 
    186  1.5.2.1    yamt 	/* copy the result to the src opr */
    187  1.5.2.1    yamt 	CPYFPN(&fe->fe_f2, fp);
    188  1.5.2.1    yamt 
    189  1.5.2.1    yamt 	return fpu_etox(fe);
    190      1.1  briggs }
    191