1*b7275c88Smartynas /* $OpenBSD: etanh.c,v 1.1 2011/07/02 18:11:01 martynas Exp $ */
2*b7275c88Smartynas
3*b7275c88Smartynas /*
4*b7275c88Smartynas * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5*b7275c88Smartynas *
6*b7275c88Smartynas * Permission to use, copy, modify, and distribute this software for any
7*b7275c88Smartynas * purpose with or without fee is hereby granted, provided that the above
8*b7275c88Smartynas * copyright notice and this permission notice appear in all copies.
9*b7275c88Smartynas *
10*b7275c88Smartynas * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11*b7275c88Smartynas * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12*b7275c88Smartynas * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13*b7275c88Smartynas * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14*b7275c88Smartynas * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15*b7275c88Smartynas * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16*b7275c88Smartynas * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17*b7275c88Smartynas */
18*b7275c88Smartynas
19*b7275c88Smartynas /* xtanh.c */
20*b7275c88Smartynas /* hyperbolic tangent check routine */
21*b7275c88Smartynas /* this subroutine is used by the exponential function routine */
22*b7275c88Smartynas /* by Stephen L. Moshier. */
23*b7275c88Smartynas
24*b7275c88Smartynas
25*b7275c88Smartynas
26*b7275c88Smartynas #include "ehead.h"
27*b7275c88Smartynas
28*b7275c88Smartynas
etanh(x,y)29*b7275c88Smartynas void etanh( x, y )
30*b7275c88Smartynas unsigned short *x, *y;
31*b7275c88Smartynas {
32*b7275c88Smartynas unsigned short e[NE], r[NE], j[NE], xx[NE], m2[NE];
33*b7275c88Smartynas short i, n;
34*b7275c88Smartynas long lj;
35*b7275c88Smartynas
36*b7275c88Smartynas emov( x, r );
37*b7275c88Smartynas r[NE-1] &= (unsigned short )0x7fff;
38*b7275c88Smartynas if( ecmp(r, eone) >= 0 )
39*b7275c88Smartynas {
40*b7275c88Smartynas /* tanh(x) = (exp(x) - exp(-x)) / (exp(x) + exp(-x))
41*b7275c88Smartynas * Note eexp() calls xtanh, but with an argument less than (1 + log 2)/2.
42*b7275c88Smartynas */
43*b7275c88Smartynas eexp( r, e );
44*b7275c88Smartynas ediv( e, eone, r );
45*b7275c88Smartynas esub( r, e, xx );
46*b7275c88Smartynas eadd( r, e, j );
47*b7275c88Smartynas ediv( j, xx, y );
48*b7275c88Smartynas return;
49*b7275c88Smartynas }
50*b7275c88Smartynas
51*b7275c88Smartynas emov( etwo, m2 );
52*b7275c88Smartynas eneg( m2 );
53*b7275c88Smartynas
54*b7275c88Smartynas n = NBITS/8; /* Number of terms to do in the continued fraction */
55*b7275c88Smartynas lj = 2 * n + 1;
56*b7275c88Smartynas ltoe( &lj, j );
57*b7275c88Smartynas
58*b7275c88Smartynas emov( j, e );
59*b7275c88Smartynas emul( x, x, xx );
60*b7275c88Smartynas
61*b7275c88Smartynas /* continued fraction */
62*b7275c88Smartynas for( i=0; i<n; i++)
63*b7275c88Smartynas {
64*b7275c88Smartynas ediv( e, xx, r );
65*b7275c88Smartynas eadd( m2, j, j );
66*b7275c88Smartynas eadd( r, j, e );
67*b7275c88Smartynas }
68*b7275c88Smartynas
69*b7275c88Smartynas ediv( e, x, y );
70*b7275c88Smartynas }
71