1*b7275c88Smartynas /* $OpenBSD: elog.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 /* xlog.c */
20*b7275c88Smartynas /* natural logarithm */
21*b7275c88Smartynas /* by Stephen L. Moshier. */
22*b7275c88Smartynas
23*b7275c88Smartynas #include "mconf.h"
24*b7275c88Smartynas #include "ehead.h"
25*b7275c88Smartynas
26*b7275c88Smartynas
27*b7275c88Smartynas
elog(x,y)28*b7275c88Smartynas void elog( x, y )
29*b7275c88Smartynas unsigned short *x, *y;
30*b7275c88Smartynas {
31*b7275c88Smartynas unsigned short xx[NE], z[NE], a[NE], b[NE], t[NE], qj[NE];
32*b7275c88Smartynas long ex;
33*b7275c88Smartynas int fex;
34*b7275c88Smartynas
35*b7275c88Smartynas
36*b7275c88Smartynas if( x[NE-1] & (unsigned short )0x8000 )
37*b7275c88Smartynas {
38*b7275c88Smartynas eclear(y);
39*b7275c88Smartynas mtherr( "elog", DOMAIN );
40*b7275c88Smartynas return;
41*b7275c88Smartynas }
42*b7275c88Smartynas if( ecmp( x, ezero ) == 0 )
43*b7275c88Smartynas {
44*b7275c88Smartynas einfin( y );
45*b7275c88Smartynas eneg(y);
46*b7275c88Smartynas mtherr( "elog", SING );
47*b7275c88Smartynas return;
48*b7275c88Smartynas }
49*b7275c88Smartynas if( ecmp( x, eone ) == 0 )
50*b7275c88Smartynas {
51*b7275c88Smartynas eclear( y );
52*b7275c88Smartynas return;
53*b7275c88Smartynas }
54*b7275c88Smartynas
55*b7275c88Smartynas /* range reduction: log x = log( 2**ex * m ) = ex * log2 + log m */
56*b7275c88Smartynas efrexp( x, &fex, xx );
57*b7275c88Smartynas /*
58*b7275c88Smartynas emov(x, xx );
59*b7275c88Smartynas ex = xx[NX-1] & 0x7fff;
60*b7275c88Smartynas ex -= 0x3ffe;
61*b7275c88Smartynas xx[NX-1] = 0x3ffe;
62*b7275c88Smartynas */
63*b7275c88Smartynas
64*b7275c88Smartynas /* Adjust range to 1/sqrt(2), sqrt(2) */
65*b7275c88Smartynas esqrt2[NE-1] -= 1;
66*b7275c88Smartynas if( ecmp( xx, esqrt2 ) < 0 )
67*b7275c88Smartynas {
68*b7275c88Smartynas fex -= 1;
69*b7275c88Smartynas emul( xx, etwo, xx );
70*b7275c88Smartynas }
71*b7275c88Smartynas esqrt2[NE-1] += 1;
72*b7275c88Smartynas
73*b7275c88Smartynas esub( eone, xx, a );
74*b7275c88Smartynas if( a[NE-1] == 0 )
75*b7275c88Smartynas {
76*b7275c88Smartynas eclear( y );
77*b7275c88Smartynas goto logdon;
78*b7275c88Smartynas }
79*b7275c88Smartynas eadd( eone, xx, b );
80*b7275c88Smartynas ediv( b, a, y ); /* store (x-1)/(x+1) in y */
81*b7275c88Smartynas
82*b7275c88Smartynas emul( y, y, z );
83*b7275c88Smartynas
84*b7275c88Smartynas emov( eone, a );
85*b7275c88Smartynas emov( eone, b );
86*b7275c88Smartynas emov( eone, qj );
87*b7275c88Smartynas do
88*b7275c88Smartynas {
89*b7275c88Smartynas eadd( etwo, qj, qj ); /* 2 * i + 1 */
90*b7275c88Smartynas emul( z, a, a );
91*b7275c88Smartynas ediv( qj, a, t );
92*b7275c88Smartynas eadd( t, b, b );
93*b7275c88Smartynas }
94*b7275c88Smartynas while( ((b[NE-1] & 0x7fff) - (t[NE-1] & 0x7fff)) < NBITS );
95*b7275c88Smartynas
96*b7275c88Smartynas
97*b7275c88Smartynas emul( b, y, y );
98*b7275c88Smartynas emul( y, etwo, y );
99*b7275c88Smartynas
100*b7275c88Smartynas logdon:
101*b7275c88Smartynas
102*b7275c88Smartynas /* now add log of 2**ex */
103*b7275c88Smartynas if( fex != 0 )
104*b7275c88Smartynas {
105*b7275c88Smartynas ex = fex;
106*b7275c88Smartynas ltoe( &ex, b );
107*b7275c88Smartynas emul( elog2, b, b );
108*b7275c88Smartynas eadd( b, y, y );
109*b7275c88Smartynas }
110*b7275c88Smartynas }
111