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: log.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/* double log (double A)
38   The log() function returns the natural logarithm of A.
39 */
40
41#define	FL_M1	0xbf800000	/* -1.0	*/
42#define	LN_2	0x3f317218	/* ln(2.0) = 0.69314718	*/
43
44#define	rC3	YH
45#define	rC2	YL		/* adiw,sbiw are used	*/
46#define	rC1	r17
47#define	rC0	r16
48#define	rCE	r15
49
50FUNCTION log
51
52.L_nf:	brts	.L_nan		; branch, if -Inf
53	XJMP	_U(__fp_mpack)	; pass as is: NaN --> NaN, +Inf --> +Inf
54.L_nan:	XJMP	_U(__fp_nan)
55.L_min:	set
56	XJMP	_U(__fp_inf)
57
58ENTRY log
59	XCALL	_U(__fp_splitA)
60	brcs	.L_nf		; !isfinite(A)
61	tst	rA3
62	breq	.L_min		; -Inf
63	brts	.L_nan		; log(negative)
64
65	push	rC3
66	push	rC2
67	push	rC1
68	push	rC0
69	push	rCE
70  ; normalize (if Subnormal)
71	mov	rC2, rA3	; place exponent to pair rC3.rC2
72	clr	rC3
73	tst	rA2
74	brmi	2f
751:	sbiw	rC2, 1
76	lsl	rA0
77	rol	rA1
78	rol	rA2
79	brpl	1b
802:
81  ; calculate log(1+x) for x:  -0x0.20p+0 <= x < 0.E0p+0
82	/* Different tables are used:
83	     .L_tlow  - for small values: -0x0.20p+0 <= x < + 0x0.18p+0
84	     .L_thigh - for others	*/
85	ldi	rB0,  lo8(FL_M1)
86	ldi	rB1,  hi8(FL_M1)
87	ldi	rB2, hlo8(FL_M1)
88	ldi	rB3, hhi8(FL_M1)
89	ldi	rA3, 0x3f
90	cpi	rA2, 0x98
91	brlo	3f
92	cpi	rA2, 0xE0
93	brlo	4f
94	adiw	rC2, 1
95	andi	rA2, ~0x80
963:	XCALL	_U(__addsf3)
97	ldi	ZL, lo8(.L_tlow)
98	ldi	ZH, hi8(.L_tlow)
99	rjmp	5f
1004:	XCALL	_U(__addsf3)
101	ldi	ZL, lo8(.L_thigh)
102	ldi	ZH, hi8(.L_thigh)
1035:	XCALL	_U(__fp_powser)
104
105	X_movw	rC0, rA0
106	X_movw	rA0, rC2	; rA1.rA0 = exponent (possible negative)
107	X_movw	rC2, rA2
108	mov	rCE, rAE
109  ; calculate log() of exponent
110	subi	rA0, 127	; 127 is exponent field of 1.0
111	sbc	rA1, r1
112	asr	rA1
113	rol	rA1		; C = rA1 sign bit
114	sbc	rA2, rA2
115	sbc	rA3, rA3
116	XCALL	_U(__floatsisf)
117	ldi	rB0,  lo8(LN_2)
118	ldi	rB1,  hi8(LN_2)
119	ldi	rB2, hlo8(LN_2)
120	ldi	rB3, hhi8(LN_2)
121	XCALL	_U(__mulsf3x)
122
123	mov	rBE, rCE
124	X_movw	rB0, rC0
125	X_movw	rB2, rC2
126
127	pop	rCE
128	pop	rC0
129	pop	rC1
130	pop	rC2
131	pop	rC3
132
133	XCALL	_U(__addsf3x)
134	XJMP	_U(__fp_round)
135ENDFUNC
136
137	PGM_SECTION
138
139.L_tlow:		; __fp_powser() table for small values
140	.byte	8
141	.byte        0x00,0x00,0x00,0xbe        ; -0.1250000000
142	.byte   0x92,0x24,0x49,0x12,0x3e        ;  0.1428571429
143	.byte   0xab,0xaa,0xaa,0x2a,0xbe        ; -0.1666666667
144	.byte   0xcd,0xcc,0xcc,0x4c,0x3e        ;  0.2000000000
145	.byte   0x00,0x00,0x00,0x80,0xbe        ; -0.2500000000
146	.byte   0xab,0xaa,0xaa,0xaa,0x3e        ;  0.3333333333
147	.byte   0x00,0x00,0x00,0x00,0xbf        ; -0.5000000000
148	.byte   0x00,0x00,0x00,0x80,0x3f        ;  1.0000000000
149	.byte   0x00,0x00,0x00,0x00,0x00        ;  0.0000000000
150
151.L_thigh:		; __fp_powser() table for other values
152	.byte	8
153	.byte        0x41,0x78,0xd3,0xbb        ; -0.0064535442
154	.byte   0x43,0x87,0xd1,0x13,0x3d        ;  0.0360884937
155	.byte   0x19,0x0e,0x3c,0xc3,0xbd        ; -0.0953293897
156	.byte   0x42,0x82,0xad,0x2b,0x3e        ;  0.1676540711
157	.byte   0x68,0xec,0x82,0x76,0xbe        ; -0.2407338084
158	.byte   0xd9,0x8f,0xe1,0xa9,0x3e        ;  0.3317990258
159	.byte   0x4c,0x80,0xef,0xff,0xbe        ; -0.4998741238
160	.byte   0x01,0xc4,0xff,0x7f,0x3f        ;  0.9999964239
161	.byte   0x00,0x00,0x00,0x00,0x00        ;  0.0000000000
162
163	.end
164
165#endif /* !defined(__AVR_TINY__) */
166