/[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.6 by joerg_wunsch, Thu Jun 30 22:31:07 2005 UTC revision 1.6.2.1 by joerg_wunsch, Sat Nov 19 21:30:20 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 38  Line 39 
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"
# Line 50  Line 71 
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
   
   

Legend:
Removed from v.1.6  
changed lines
  Added in v.1.6.2.1

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