Home | History | Annotate | Line # | Download | only in fpe
fpu_exp.c revision 1.6 
      1 /*	$NetBSD: fpu_exp.c,v 1.6 2013/04/20 03:06:19 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.6 2013/04/20 03:06:19 isaki Exp $");
     36 
     37 #include "fpu_emulate.h"
     38 
     39 /*
     40  * fpu_exp.c: defines fpu_etox(), fpu_etoxm1(), fpu_tentox(), and fpu_twotox();
     41  */
     42 
     43 /*
     44  *                  x^2   x^3   x^4
     45  * exp(x) = 1 + x + --- + --- + --- + ...
     46  *                   2!    3!    4!
     47  */
     48 static struct fpn *
     49 fpu_etox_taylor(struct fpemu *fe)
     50 {
     51 	struct fpn res;
     52 	struct fpn x;
     53 	struct fpn s0;
     54 	struct fpn *s1;
     55 	struct fpn *r;
     56 	uint32_t k;
     57 
     58 	CPYFPN(&x, &fe->fe_f2);
     59 	CPYFPN(&s0, &fe->fe_f2);
     60 
     61 	/* res := 1 + x */
     62 	fpu_const(&fe->fe_f1, FPU_CONST_1);
     63 	r = fpu_add(fe);
     64 	CPYFPN(&res, r);
     65 
     66 	k = 2;
     67 	for (;; k++) {
     68 		/* s1 = s0 * x / k */
     69 		CPYFPN(&fe->fe_f1, &s0);
     70 		CPYFPN(&fe->fe_f2, &x);
     71 		r = fpu_mul(fe);
     72 
     73 		CPYFPN(&fe->fe_f1, r);
     74 		fpu_explode(fe, &fe->fe_f2, FTYPE_LNG, &k);
     75 		s1 = fpu_div(fe);
     76 
     77 		/* break if s1 is enough small */
     78 		if (ISZERO(s1))
     79 			break;
     80 		if (res.fp_exp - s1->fp_exp >= FP_NMANT)
     81 			break;
     82 
     83 		/* s0 := s1 for next loop */
     84 		CPYFPN(&s0, s1);
     85 
     86 		/* res += s1 */
     87 		CPYFPN(&fe->fe_f2, s1);
     88 		CPYFPN(&fe->fe_f1, &res);
     89 		r = fpu_add(fe);
     90 		CPYFPN(&res, r);
     91 	}
     92 
     93 	CPYFPN(&fe->fe_f2, &res);
     94 	return &fe->fe_f2;
     95 }
     96 
     97 /*
     98  * exp(x)
     99  */
    100 struct fpn *
    101 fpu_etox(struct fpemu *fe)
    102 {
    103 	struct fpn *fp;
    104 
    105 	if (ISNAN(&fe->fe_f2))
    106 		return &fe->fe_f2;
    107 	if (ISINF(&fe->fe_f2)) {
    108 		if (fe->fe_f2.fp_sign)
    109 			fpu_const(&fe->fe_f2, FPU_CONST_0);
    110 		return &fe->fe_f2;
    111 	}
    112 
    113 	if (fe->fe_f2.fp_sign == 0) {
    114 		/* exp(x) */
    115 		fp = fpu_etox_taylor(fe);
    116 	} else {
    117 		/* 1/exp(-x) */
    118 		fe->fe_f2.fp_sign = 0;
    119 		fp = fpu_etox_taylor(fe);
    120 
    121 		CPYFPN(&fe->fe_f2, fp);
    122 		fpu_const(&fe->fe_f1, FPU_CONST_1);
    123 		fp = fpu_div(fe);
    124 	}
    125 
    126 	return fp;
    127 }
    128 
    129 /*
    130  * exp(x) - 1
    131  */
    132 struct fpn *
    133 fpu_etoxm1(struct fpemu *fe)
    134 {
    135 	struct fpn *fp;
    136 
    137 	fp = fpu_etox(fe);
    138 
    139 	CPYFPN(&fe->fe_f1, fp);
    140 	/* build a 1.0 */
    141 	fp = fpu_const(&fe->fe_f2, FPU_CONST_1);
    142 	fe->fe_f2.fp_sign = !fe->fe_f2.fp_sign;
    143 	/* fp = f2 - 1.0 */
    144 	fp = fpu_add(fe);
    145 
    146 	return fp;
    147 }
    148 
    149 /*
    150  * 10^x = exp(x * ln10)
    151  */
    152 struct fpn *
    153 fpu_tentox(struct fpemu *fe)
    154 {
    155 	struct fpn *fp;
    156 
    157 	/* build a ln10 */
    158 	fp = fpu_const(&fe->fe_f1, FPU_CONST_LN_10);
    159 	/* fp = ln10 * f2 */
    160 	fp = fpu_mul(fe);
    161 
    162 	/* copy the result to the src opr */
    163 	CPYFPN(&fe->fe_f2, fp);
    164 
    165 	return fpu_etox(fe);
    166 }
    167 
    168 /*
    169  * 2^x = exp(x * ln2)
    170  */
    171 struct fpn *
    172 fpu_twotox(struct fpemu *fe)
    173 {
    174 	struct fpn *fp;
    175 
    176 	/* build a ln2 */
    177 	fp = fpu_const(&fe->fe_f1, FPU_CONST_LN_2);
    178 	/* fp = ln2 * f2 */
    179 	fp = fpu_mul(fe);
    180 
    181 	/* copy the result to the src opr */
    182 	CPYFPN(&fe->fe_f2, fp);
    183 
    184 	return fpu_etox(fe);
    185 }
    186