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: pow.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 pow (float x, float y); 38 The pow() function returns the value of x raised to the power of y. 39 40 Args combinations: 41 x y pow(x,y) 42 -------------------------------- 43 +0 NaN NaN exp(log(x)*y) 44 +0 +0,-0 +1 45 +0 +Inf +0 exp(log(x)*y) 46 +0 -Inf +Inf exp(log(x)*y) 47 +0 y > 0 +0 exp(log(x)*y) 48 +0 y < 0 +Inf exp(log(x)*y) 49 50 -0 NaN NaN exp(log(x)*y) 51 -0 +0,-0 +1 52 -0 1,3,5... -0 -exp(log(x)*y) 53 -0 y > 0 +0 exp(log(x)*y) 54 -0 -1,-3,-5... -Inf -exp(log(x)*y) 55 -0 y < 0 +Inf exp(log(x)*y) 56 57 +1 NaN +1 58 +1 +Inf,-Inf +1 59 +1 else +1 exp(log(x)*y) 60 61 -1 +0,-0 +1 62 -1 1,3,5... -1 63 -1 -1,-3,-5... -1 64 -1 2,4,6... +1 65 -1 -2,-4,-6... +1 66 -1 else NaN exp(log(x)*y) 67 68 +Inf NaN NaN exp(log(x)*y) 69 +Inf +0,-0 +1 70 +Inf y > 0 +Inf exp(log(x)*y) 71 +Inf y < 0 +0 exp(log(x)*y) 72 73 -Inf NaN NaN exp(log(x)*y) 74 -Inf +0,-0 +1 75 -Inf 1,3,5... -Inf 76 -Inf y > 0 +Inf 77 -Inf -1,-3,-5... -0 78 -Inf y < 0 +0 79 80 NaN +0,-0 +1 81 NaN else NaN exp(log(x)*y) 82 83 (0,1) NaN NaN exp(log(x)*y) 84 (0,1) +0,-0 +1 exp(log(x)*y) 85 (0,1) +Inf +0 exp(log(x)*y) 86 (0,1) -Inf +Inf exp(log(x)*y) 87 88 (-1,0) NaN NaN 89 (-1,0) +0,-0 +1 90 (-1,0) +Inf +0 91 (-1,0) -Inf +Inf 92 (-1,0) nonintegral NaN 93 94 x > 1 NaN NaN 95 x > 1 +0,-0 +1 96 x > 1 +Inf +Inf 97 x > 1 -Inf +0 98 99 x < -1 NaN NaN 100 x < -1 +0,-0 +1 101 x < -1 +Inf +Inf 102 x < -1 -Inf +0 103 x < -1 nonintegral NaN 104 */ 105 106#define FL_1 0x3f800000 /* +1.0 */ 107 108ENTRY pow 109 ; ZH := exponent of y 110 X_movw ZL, rB2 111 lsl ZL 112 rol ZH 113 ; y == 0 ? 114 adiw ZL, 0 115 cpc rB0, r1 116 cpc rB1, r1 117 breq .L_one 118 ; preliminary check 119 cp rA0, r1 120 cpc rA1, r1 121 brne 0f ; skip a bit of comparisons 122 ; x == 1.0 ? 123 cpi rA2, hlo8(FL_1) 124 ldi rAE, hhi8(FL_1) 125 cpc rA3, rAE 126 breq .L_ret 127 ; x == -0.0 ? 128 set ; flag: nonintegral y is a legal value 129 cpi rA3, 0x80 130 cpc rA2, r1 131 breq .L_int 132 ; x == -Inf ? 133 cpi rA2, 0x80 134 ldi rAE, 0xff 135 cpc rA3, rAE 136 breq .L_int 137 ; x >= 0 ? 1380: tst rA3 139 brpl .L_pow 140 ; isinf(y) ? 141 cpi ZH, 0xff 142 cpc ZL, r1 143 cpc rB1, r1 144 cpc rB0, r1 145 breq .L_big 146 ; isintegral(y) ? 147 clt ; nonintegral y is not a legal value 148.L_int: 149 /* Now we have: 150 y is nonzero value 151 ZL == (rB2 << 1) 152 ZH == exponenta, ZH <= 254 */ 153 sec ; hidden bit 154 ror ZL ; This is incorrect for subnormals, no sense: 155 ; result would NaN. 156 ; ffs(). Next two loops are finite due to above 'sec'. 157 X_movw XL, rB0 158 ; Byte search loop. 1591: tst XL 160 brne 2f 161 mov XL, XH 162 mov XH, ZL 163 subi ZH, -8 164 brcs 1b 165 rjmp .L_noint ; mantisa too big 166 ; Bit search loop. 1672: subi ZH, -1 168 brcc .L_noint ; mantisa too big 169 lsr XL 170 brcc 2b 171 ; Check exponent, is y an integral value? 172 /* Example: 1.0 == 0x3f800000: 173 exponent: ZH := 0x7f 174 byte search: ZH += 2*8 --> 0x8f 175 bit search: ZH += 8 --> 0x97 */ 176 cpi ZH, 0x97 177 brlo .L_noint 178 breq 3f ; y % 2 == 1 179 cpi ZH, 0x97 + 24 180 brsh .L_noint 181 andi rA3, 0x7f ; y is integral, y % 2 == 0 1823: push rA3 183 rcall .L_pow 184 pop r0 185 sbrc r0, 7 186 subi rA3, 0x80 187.L_ret: 188 ret 189 ; y is not an integral number 190.L_noint: 191 brts .L_pow 192.L_nan: 193 XJMP _U(__fp_nan) 194.L_one: 195 ldi rA0, lo8(FL_1) 196 ldi rA1, hi8(FL_1) 197 ldi rA2, hlo8(FL_1) 198 ldi rA3, hhi8(FL_1) 199 ret 200 ; replace Inf --> big finite (to exclude '0 * Inf' for legal x == -1) 201.L_big: 202 ldi rB2, 0x7f 203.L_pow: 204 andi rA3, 0x7f 205 push rB3 206 push rB2 207 push rB1 208 push rB0 209 XCALL _U(log) 210 pop rB0 211 pop rB1 212 pop rB2 213 pop rB3 214 XCALL _U(__mulsf3) 215 XJMP _U(exp) 216ENDFUNC 217 218#endif /* !defined(__AVR_TINY__) */ 219