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