1 |
/* -*- Mode: Asm -*- */ |
/* -*- Mode: Asm -*- */ |
2 |
|
|
3 |
/* Copyright (c) 2002 Michael Stumpf <mistumpf@de.pepperl-fuchs.com> |
/* Copyright (c) 2002 Michael Stumpf <mistumpf@de.pepperl-fuchs.com> |
4 |
|
Copyright (c) 2005 Dmitry Xmelkov |
5 |
All rights reserved. |
All rights reserved. |
6 |
|
|
7 |
|
|
39 |
pow.S is part of FPlib V 0.3.0 ported to avr-as |
pow.S is part of FPlib V 0.3.0 ported to avr-as |
40 |
for details see readme.fplib |
for details see readme.fplib |
41 |
|
|
42 |
*---------------------------------------------------------------------------------------- |
*--------------------------------------------------------------------- |
43 |
* |
* |
44 |
* |
* A = pow(x,y) |
45 |
* A = pow(x,y) = exp( log(x)*y ) | x > 0 |
|
46 |
|
double pow (double x, double y) |
47 |
|
{ |
48 |
|
return |
49 |
|
!isfinite(x) || !isfinite(y) ? NaN |
50 |
|
: (x > 0) ? exp(log(x)*y) |
51 |
|
: (x == 0 && y > 0) ? 0 |
52 |
|
: (x == 0 && y <= 0) ? NaN |
53 |
|
: (is_int(y) && (y % 2)) ? sign(x)*exp(log(fabs(x))*y) |
54 |
|
: (is_int(y) && !(y % 2)) ? exp(log(fabs(x))*y) |
55 |
|
: NaN ; |
56 |
|
} |
57 |
|
|
58 |
|
Notes: |
59 |
|
* C99 says: |
60 |
|
"A domain error occurs if the result cannot be represented when x is |
61 |
|
zero and y is less than or equal to zero. A range error may occur." |
62 |
|
In this implementation pow(0,0)=NaN (shortest code). |
63 |
|
* To do not enlarge program size, the case y==0 is not analized |
64 |
|
separatly. |
65 |
|
* __fp_tstA, __fp_tstB will global in future. |
66 |
*/ |
*/ |
67 |
|
|
68 |
#include "gasava.inc" |
#include "gasava.inc" |
71 |
|
|
72 |
TEXT_SEG(fplib, pow) |
TEXT_SEG(fplib, pow) |
73 |
FUNCTION(pow) |
FUNCTION(pow) |
|
|
|
74 |
GLOBAL(pow) |
GLOBAL(pow) |
75 |
TST rB3 |
/* Check for NaN right here, before calling exp() as the latter |
76 |
BREQ 10f ; x ^ 0 = 1.0 |
is known to be buggy for NaN arguments. */ |
77 |
TST rA3 |
rcall __fp_tstB |
78 |
BRMI 1f |
brcs .L_nan ; isfinite(y) |
79 |
BREQ 10f |
rcall __fp_tstA |
80 |
|
brcs .L_ret ; !isfinite(x) |
81 |
PUSH rB3 |
brts .L_negx ; x < 0 |
82 |
PUSH rB2 |
brne .L_pow |
83 |
PUSH rB1 |
/* x == 0 */ |
84 |
PUSH rB0 |
rcall __fp_tstB |
85 |
RCALL _U(log) |
breq .L_nan ; y == 0 |
86 |
POP rB0 |
brtc .L_ret ; y > 0 --> return 0 |
87 |
POP rB1 |
.L_nan: |
88 |
POP rB2 |
rjmp __fp_nan |
89 |
POP rB3 |
|
90 |
RCALL _U(__mulsf3) |
.L_pow: |
91 |
RJMP _U(exp) |
push rB3 |
92 |
|
push rB2 |
93 |
1: |
push rB1 |
94 |
/* x < 0 -> no rational power allowed |
push rB0 |
95 |
* check if B is an integer number |
rcall log |
96 |
* if yes then pow(x,y) = |
pop rB0 |
97 |
*/ |
pop rB1 |
98 |
PUSH rA3 |
pop rB2 |
99 |
PUSH rA2 |
pop rB3 |
100 |
PUSH rA1 |
rcall __mulsf3 |
101 |
PUSH rA0 |
rjmp exp |
102 |
PUSH rB3 |
|
103 |
PUSH rB2 |
.L_negx: |
104 |
PUSH rB1 |
rcall __fp_tstB |
105 |
PUSH rB0 |
brcs .L_pow ; --> NaN, as x < 0 |
106 |
MOV rA3,rB3 |
breq .L_absx |
107 |
MOV rA2,rB2 |
/* Now we have: |
108 |
MOV rA1,rB1 |
y is normal, nonzero value |
109 |
MOV rA0,rB0 |
ZL == (rB2 << 1) |
110 |
RCALL _U(__fixsfsi) |
ZH == exponenta, ZH <= 254 */ |
111 |
|
sec ; hidden bit |
112 |
RCALL _U(__subsf3) |
ror ZL ; This is incorrect for subnormals, no sense: |
113 |
CP rA3,rT1c |
; result would NaN. |
114 |
CPC rA2,rT1c |
/* ffs(). Next two loops are finite due to above 'sec'. */ |
115 |
CPC rA1,rT1c |
X_movw XL, rB0 |
116 |
CPC rA0,rT1c |
/* Byte search loop. */ |
117 |
POP rB0 |
1: tst XL |
118 |
POP rB1 |
brne 2f |
119 |
POP rB2 |
mov XL, XH |
120 |
POP rB3 |
mov XH, ZL |
121 |
POP rA0 |
subi ZH, -8 |
122 |
POP rA1 |
brcs 1b |
123 |
POP rA2 |
rjmp .L_pow ; mantisa too big |
124 |
POP rA3 |
/* Bit search loop. */ |
125 |
BREQ 2f |
2: subi ZH, -1 |
126 |
ANDI rA3,0x7F |
brcc .L_pow ; mantisa too big |
127 |
SUBI rB3,0x80 |
lsr XL |
128 |
RCALL 1b |
brcc 2b |
129 |
SUBI rB3,0x80 |
/* Check exponent, is y an integral value? |
130 |
|
Example: 1.0 == 0x3f800000: |
131 |
2: |
__fp_tstB: ZH := 0x7f |
132 |
RJMP _U(__fp_nanEDOM) |
byte search: ZH += 2*8 --> 0x8f |
133 |
|
bit search: ZH += 8 --> 0x97 */ |
134 |
10: |
cpi ZH, 0x97 |
135 |
LDI rA3,0x3F |
brlo .L_nan |
136 |
LDI rA2,0x80 |
breq 3f ; y % 2 == 1 |
137 |
CLR rA1 |
.L_absx: |
138 |
CLR rA0 |
andi rA3, 0x7f ; y is integral, y % 2 == 0 |
139 |
RET |
3: push rA3 |
140 |
|
andi rA3, 0x7f |
141 |
|
rcall .L_pow |
142 |
|
pop ZL |
143 |
|
andi ZL, 0x80 |
144 |
|
eor rA3, ZL |
145 |
|
.L_ret: |
146 |
|
ret |
147 |
|
|
148 |
|
/* Test of 32-bits float value in rA0..rA3 registers. |
149 |
|
Return: |
150 |
|
C = 1 Z = 0 T = * , if NaN |
151 |
|
C = 0 Z = 1 T = 0 , if 0.0 (or -0.0) |
152 |
|
C = 0 Z = 0 T = 1 , if < 0.0 |
153 |
|
C = 0 Z = 0 T = 0 , if > 0.0 |
154 |
|
ZL = (rA2 << 1) |
155 |
|
ZH = exponent (without subnormals correction) |
156 |
|
Value and other registers are not changed. |
157 |
|
*/ |
158 |
|
__fp_tstA: |
159 |
|
clt |
160 |
|
X_movw ZL, rA2 |
161 |
|
lsl ZL |
162 |
|
rol ZH |
163 |
|
adiw ZL, 0 |
164 |
|
cpc rA0, __zero_reg__ |
165 |
|
cpc rA1, __zero_reg__ |
166 |
|
breq 1f ; C == 0, Z == 1 |
167 |
|
inc ZH |
168 |
|
sbci ZH, 1 ; C == 0/1, Z == 0 |
169 |
|
bst rA3, 7 ; sign bit |
170 |
|
1: ret |
171 |
|
|
172 |
|
/* Test of 32-bits float value in rB0..rB3 registers. |
173 |
|
Return: |
174 |
|
C = 1 Z = 0 T = * , if NaN |
175 |
|
C = 0 Z = 1 T = 0 , if 0.0 (or -0.0) |
176 |
|
C = 0 Z = 0 T = 1 , if < 0.0 |
177 |
|
C = 0 Z = 0 T = 0 , if > 0.0 |
178 |
|
ZL = (rB2 << 1) |
179 |
|
ZH = exponent (without subnormals correction) |
180 |
|
Value and other registers are not changed. |
181 |
|
*/ |
182 |
|
__fp_tstB: |
183 |
|
clt |
184 |
|
X_movw ZL, rB2 |
185 |
|
lsl ZL |
186 |
|
rol ZH |
187 |
|
adiw ZL, 0 |
188 |
|
cpc rB0, __zero_reg__ |
189 |
|
cpc rB1, __zero_reg__ |
190 |
|
breq 1f ; C == 0, Z == 1 |
191 |
|
inc ZH |
192 |
|
sbci ZH, 1 ; C == 0/1, Z == 0 |
193 |
|
bst rB3, 7 ; sign bit |
194 |
|
1: ret |
195 |
|
|
196 |
ENDFUNC |
ENDFUNC |
|
|
|
|
|
|