/[avr-libc]/avr-libc/libm/fplib/pow.S
ViewVC logotype

Diff of /avr-libc/libm/fplib/pow.S

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.1.4.2.1  
changed lines
  Added in v.1.4.2.2

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26