| b_exp.c (1573) | b_exp.c (8870) |
|---|---|
| 1/* 2 * Copyright (c) 1985, 1993 3 * The Regents of the University of California. 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 7 * are met: 8 * 1. Redistributions of source code must retain the above copyright --- 24 unchanged lines hidden (view full) --- 33 34#ifndef lint 35static char sccsid[] = "@(#)exp.c 8.1 (Berkeley) 6/4/93"; 36#endif /* not lint */ 37 38/* EXP(X) 39 * RETURN THE EXPONENTIAL OF X 40 * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS) | 1/* 2 * Copyright (c) 1985, 1993 3 * The Regents of the University of California. 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 7 * are met: 8 * 1. Redistributions of source code must retain the above copyright --- 24 unchanged lines hidden (view full) --- 33 34#ifndef lint 35static char sccsid[] = "@(#)exp.c 8.1 (Berkeley) 6/4/93"; 36#endif /* not lint */ 37 38/* EXP(X) 39 * RETURN THE EXPONENTIAL OF X 40 * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS) |
| 41 * CODED IN C BY K.C. NG, 1/19/85; | 41 * CODED IN C BY K.C. NG, 1/19/85; |
| 42 * REVISED BY K.C. NG on 2/6/85, 2/15/85, 3/7/85, 3/24/85, 4/16/85, 6/14/86. 43 * 44 * Required system supported functions: | 42 * REVISED BY K.C. NG on 2/6/85, 2/15/85, 3/7/85, 3/24/85, 4/16/85, 6/14/86. 43 * 44 * Required system supported functions: |
| 45 * scalb(x,n) 46 * copysign(x,y) | 45 * scalb(x,n) 46 * copysign(x,y) |
| 47 * finite(x) 48 * 49 * Method: | 47 * finite(x) 48 * 49 * Method: |
| 50 * 1. Argument Reduction: given the input x, find r and integer k such | 50 * 1. Argument Reduction: given the input x, find r and integer k such |
| 51 * that | 51 * that |
| 52 * x = k*ln2 + r, |r| <= 0.5*ln2 . | 52 * x = k*ln2 + r, |r| <= 0.5*ln2 . |
| 53 * r will be represented as r := z+c for better accuracy. 54 * | 53 * r will be represented as r := z+c for better accuracy. 54 * |
| 55 * 2. Compute exp(r) by | 55 * 2. Compute exp(r) by |
| 56 * 57 * exp(r) = 1 + r + r*R1/(2-R1), 58 * where 59 * R1 = x - x^2*(p1+x^2*(p2+x^2*(p3+x^2*(p4+p5*x^2)))). 60 * 61 * 3. exp(x) = 2^k * exp(r) . 62 * 63 * Special cases: --- 74 unchanged lines hidden (view full) --- 138 /* return 2^k*[1+x+x*c/(2+c)] */ 139 z=x*x; 140 c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5)))); 141 return scalb(1.0+(hi-(lo-(x*c)/(2.0-c))),k); 142 143 } 144 /* end of x > lntiny */ 145 | 56 * 57 * exp(r) = 1 + r + r*R1/(2-R1), 58 * where 59 * R1 = x - x^2*(p1+x^2*(p2+x^2*(p3+x^2*(p4+p5*x^2)))). 60 * 61 * 3. exp(x) = 2^k * exp(r) . 62 * 63 * Special cases: --- 74 unchanged lines hidden (view full) --- 138 /* return 2^k*[1+x+x*c/(2+c)] */ 139 z=x*x; 140 c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5)))); 141 return scalb(1.0+(hi-(lo-(x*c)/(2.0-c))),k); 142 143 } 144 /* end of x > lntiny */ 145 |
| 146 else | 146 else |
| 147 /* exp(-big#) underflows to zero */ 148 if(finite(x)) return(scalb(1.0,-5000)); 149 150 /* exp(-INF) is zero */ 151 else return(0.0); 152 } 153 /* end of x < lnhuge */ 154 | 147 /* exp(-big#) underflows to zero */ 148 if(finite(x)) return(scalb(1.0,-5000)); 149 150 /* exp(-INF) is zero */ 151 else return(0.0); 152 } 153 /* end of x < lnhuge */ 154 |
| 155 else | 155 else |
| 156 /* exp(INF) is INF, exp(+big#) overflows to INF */ 157 return( finite(x) ? scalb(1.0,5000) : x); 158} 159 160/* returns exp(r = x + c) for |c| < |x| with no overlap. */ 161 162double __exp__D(x, c) 163double x, c; --- 19 unchanged lines hidden (view full) --- 183 z=x*x; 184 c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5)))); 185 c = (x*c)/(2.0-c); 186 187 return scalb(1.+(hi-(lo - c)), k); 188 } 189 /* end of x > lntiny */ 190 | 156 /* exp(INF) is INF, exp(+big#) overflows to INF */ 157 return( finite(x) ? scalb(1.0,5000) : x); 158} 159 160/* returns exp(r = x + c) for |c| < |x| with no overlap. */ 161 162double __exp__D(x, c) 163double x, c; --- 19 unchanged lines hidden (view full) --- 183 z=x*x; 184 c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5)))); 185 c = (x*c)/(2.0-c); 186 187 return scalb(1.+(hi-(lo - c)), k); 188 } 189 /* end of x > lntiny */ 190 |
| 191 else | 191 else |
| 192 /* exp(-big#) underflows to zero */ 193 if(finite(x)) return(scalb(1.0,-5000)); 194 195 /* exp(-INF) is zero */ 196 else return(0.0); 197 } 198 /* end of x < lnhuge */ 199 | 192 /* exp(-big#) underflows to zero */ 193 if(finite(x)) return(scalb(1.0,-5000)); 194 195 /* exp(-INF) is zero */ 196 else return(0.0); 197 } 198 /* end of x < lnhuge */ 199 |
| 200 else | 200 else |
| 201 /* exp(INF) is INF, exp(+big#) overflows to INF */ 202 return( finite(x) ? scalb(1.0,5000) : x); 203} | 201 /* exp(INF) is INF, exp(+big#) overflows to INF */ 202 return( finite(x) ? scalb(1.0,5000) : x); 203} |