1/* Copyright (c) 2002  Michael Stumpf  <mistumpf@de.pepperl-fuchs.com>
2   Copyright (c) 2006  Dmitry Xmelkov
3   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 are met:
7
8   * Redistributions of source code must retain the above copyright
9     notice, this list of conditions and the following disclaimer.
10   * Redistributions in binary form must reproduce the above copyright
11     notice, this list of conditions and the following disclaimer in
12     the documentation and/or other materials provided with the
13     distribution.
14   * Neither the name of the copyright holders nor the names of
15     contributors may be used to endorse or promote products derived
16     from this software without specific prior written permission.
17
18   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19   AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20   IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21   ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
22   LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
23   CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
24   SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25   INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
26   CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
27   ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
28   POSSIBILITY OF SUCH DAMAGE. */
29
30/* $Id: exp.S 2473 2015-04-09 08:10:22Z pitchumani $ */
31
32#if !defined(__AVR_TINY__)
33
34#include "fp32def.h"
35#include "asmdef.h"
36
37/* float exp (float x);
38     The exp() function returns the value of e (the base of natural
39     logarithms) raised to the power of x.
40 */
41
42#define	X2BIG	 0x43000000	/* exp(-X2BIG) < 0x00000001	*/
43#define	FL_1_LN2 0x3fb8aa3b	/* 1 / ln(2)			*/
44
45/* second arg for modf(), ldexp() functions.	*/
46#define	exp_lo	r20
47#define	exp_hi	r21
48
49FUNCTION exp
50
51.L_nf:	brne	.L_nan
52.L_tb:	brts	.L_zr
53	XJMP	_U(__fp_inf)
54.L_zr:	XJMP	_U(__fp_zero)
55.L_nan:	XJMP	_U(__fp_nan)
56
57ENTRY exp
58  ; split and analize A
59	XCALL	_U(__fp_splitA)
60	brcs	.L_nf		; A is not a finite number
61	cpi	rA3, hhi8(X2BIG << 1)
62	brsh	.L_tb		; A is too big (in absolute value)
63  ; save sign and jump to positive
64	bld	r0, 7
65	push	r0
66	clt
67  ; A = A / ln(2)
68	ldi	rB0,  lo8(FL_1_LN2)
69	ldi	rB1,  hi8(FL_1_LN2)
70	ldi	rB2, hlo8(FL_1_LN2 | 0x800000)		; hidden '1'
71	ldi	rB3, hhi8(FL_1_LN2 << 1)		; exponent
72	XCALL	_U(__mulsf3_pse)
73  ; split A into fraction and integral parts
74	push	r0
75	push	r0
76	push	r0
77	in	exp_lo, SPL_IO_ADDR
78	in	exp_hi, SPH_IO_ADDR
79	push	r0
80	XCALL	_U(modf)
81  ; calculate 2**(-x) for 0 <= x < 1
82	ldi	ZL, lo8(.L_table)
83	ldi	ZH, hi8(.L_table)
84	XCALL	_U(__fp_powser)
85  ; get integral part
86	pop	exp_lo
87	pop	exp_hi
88	pop	ZL
89	pop	ZH
90  ; cast to integer
91	asr	ZL
92	rol	ZL
93	rol	ZH
94	breq	2f
95	subi	ZH, 0x7e
96	ori	ZL, 0x80
97	clr	exp_lo
981:	lsl	ZL
99	rol	exp_lo
100	dec	ZH
101	brne	1b
102  ; negate and scale
103	neg	exp_lo
104	sbc	exp_hi, exp_hi
1052:	XCALL	_U(ldexp)
106  ; inverse for positive arg.
107	pop	r0
108	sbrs	r0, 7
109	XJMP	_U(inverse)
110	ret
111ENDFUNC
112
113	PGM_SECTION
114.L_table:
115	.byte	7
116	.byte	     0x63,0x42,0x36,0xb7	; -0.0000108635
117	.byte	0x9b,0xd8,0xa7,0x1a,0x39	;  0.0001474911
118	.byte	0x68,0x56,0x18,0xae,0xba	; -0.0013282400
119	.byte	0xab,0x55,0x8c,0x1d,0x3c	;  0.0096159779
120	.byte	0xb7,0xcc,0x57,0x63,0xbd	; -0.0555036542
121	.byte	0x6d,0xed,0xfd,0x75,0x3e	;  0.2402264688
122	.byte	0xf6,0x17,0x72,0x31,0xbf	; -0.6931471802
123	.byte	0x00,0x00,0x00,0x80,0x3f	;  1.0000000000
124	.end
125
126#endif /* !defined(__AVR_TINY__) */
127