~ubuntu-branches/ubuntu/maverick/avr-libc/maverick

« back to all changes in this revision

Viewing changes to libm/fplib/sqrt.S

  • Committer: Bazaar Package Importer
  • Author(s): Hakan Ardo
  • Date: 2009-10-31 11:52:10 UTC
  • mfrom: (1.2.2 upstream)
  • mto: This revision was merged to the branch mainline in revision 10.
  • Revision ID: james.westby@ubuntu.com-20091031115210-crjd42sn6ezrj52c
* New upstream relese (closes: #544030)
* Added lintian overrides (closes: #553265)

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
1
/* Copyright (c) 2002  Michael Stumpf  <mistumpf@de.pepperl-fuchs.com>
2
2
   Copyright (c) 2006  Dmitry Xmelkov
 
3
   Copyright (c) 2008  Ruud v Gessel
 
4
   
3
5
   All rights reserved.
4
6
 
5
7
   Redistribution and use in source and binary forms, with or without
27
29
   ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
28
30
   POSSIBILITY OF SUCH DAMAGE. */
29
31
 
30
 
/* $Id: sqrt.S,v 1.8 2007/01/14 15:13:10 dmix Exp $ */
31
 
 
32
32
#include "fp32def.h"
33
33
#include "asmdef.h"
34
34
 
56
56
  ; normalize, if A is subnormal
57
57
        sbrs    rA2, 7
58
58
        rcall   _U(__fp_norm2)
59
 
  ; calculate result exponent
60
 
        lsr     rB3
61
 
        ror     rA3
62
 
  ; expand A mantissa to rAE.rA2.rA1.rA0
63
 
        ldi     rAE, 0
64
 
        brcc    1f              ; after 'ror rA3'
65
 
        lsl     rA0
66
 
        rol     rA1
67
 
        rol     rA2
68
 
        rol     rAE
69
 
1:
70
 
 
71
 
#define mvb0    ZL
72
 
#define mvb1    ZH
73
 
#define mvb2    rBE
74
 
#define rem0    r16
75
 
#define rem1    r17
76
 
#define rem2    r0
77
 
#define rem3    rB3
78
 
  ; save temporary regs.
79
 
        push    rem1
80
 
        push    rem0
81
 
  ; init variables
82
 
        /* arg's mantissa:      rAE.rA2.rA1.rA0
83
 
           result:              rB2.rB2.rB0
84
 
           moving bit:          mvb2.mvb1.mvb0
85
 
           remain value:        rem3.rem2.rem1.rem0     */
86
 
        clr     r0
87
 
        X_movw  rB0, r0
88
 
        X_movw  rB2, r0
89
 
        X_movw  rem0, r0
90
 
        X_movw  mvb0, r0
91
 
        ldi     mvb2, 0x80
92
 
.Loop:
93
 
  ; rem += mvb
94
 
        add     rem0, mvb0
95
 
        adc     rem1, mvb1
96
 
        adc     rem2, mvb2
97
 
        adc     rem3, r1
98
 
  ; A -= rem
99
 
        sub     rA0, rem0
100
 
        sbc     rA1, rem1
101
 
        sbc     rA2, rem2
102
 
        sbc     rAE, rem3
103
 
        brsh    1f
104
 
  ; restore A
105
 
        add     rA0, rem0
106
 
        adc     rA1, rem1
107
 
        adc     rA2, rem2
108
 
        adc     rAE, rem3
109
 
  ; restore rem
110
 
        sub     rem0, mvb0
111
 
        sbc     rem1, mvb1
112
 
        sbc     rem2, mvb2
113
 
        sbc     rem3, r1
114
 
        rjmp    2f
115
 
  ; B += mvb
116
 
1:      add     rB0, mvb0
117
 
        adc     rB1, mvb1
118
 
        adc     rB2, mvb2
119
 
  ; rem += mvb
120
 
        add     rem0, mvb0
121
 
        adc     rem1, mvb1
122
 
        adc     rem2, mvb2
123
 
        adc     rem3, r1
124
 
  ; A <<= 1
125
 
2:      lsl     rA0
126
 
        rol     rA1
127
 
        rol     rA2
128
 
        rol     rAE
129
 
  ; while (mvb >>= 1)
130
 
        lsr     mvb2
131
 
        ror     mvb1
132
 
        ror     mvb0
133
 
        brcc    .Loop
134
 
  ; round
135
 
        cp      rem0, rA0
136
 
        cpc     rem1, rA1
137
 
        cpc     rem2, rA2
138
 
        cpc     rem3, rAE       ; C is set if rem < A
139
 
        adc     rB0, r1
140
 
        adc     rB1, r1
141
 
        adc     rB2, r1
142
 
  ; pop temporary regs.
143
 
        pop     rem0
144
 
        pop     rem1
145
 
  ; merge result and return
 
59
 
 
60
#define msk0    r0
 
61
#define msk1    r1
 
62
#define msk2    rBE
 
63
 
 
64
        clr     msk0            ; msk1=R1 already 0
 
65
        ldi     msk2, 0x60      ; Initial rotation mask   =
 
66
                                ;       01100000.00000000.00000000
 
67
        ldi     rB2, 0xa0
 
68
        X_movw  rB0, msk0       ; Initial developing root =
 
69
                                ;       10100000.00000000.00000000
 
70
 
 
71
/* TODO: Now the Avr-libs does not have an infrastructure to build and
 
72
   *test automaticaly* with both OPTIMIZE_SPEED definitions. So the
 
73
   one variant is enabled today only.
 
74
 */
 
75
#if     1 /* defined(OPTIMIZE_SPEED) && OPTIMIZE_SPEED  */
 
76
 
 
77
  ;** Optimized for speed (9 code words larger than size optimized,
 
78
  ; 67 less cycles in average)
 
79
 
 
80
        subi    rA2, 0x80
 
81
        lsr     rB3
 
82
        ror     rA3             ; Divide exponent by 2, C==>exponent was odd
 
83
        brcc    1f              ; Jump for even exponent in argument
 
84
        subi    rA2, lo8(-0x40) ; Initial remainder for odd exponent.
 
85
 
 
86
        ; Loop for upper 23 bits
 
87
.Loop:  lsl     rA0
 
88
        rol     rA1
 
89
        rol     rA2             ; Shift left remainder argument
 
90
        brcs    2f              ; C --> Bit is always 1 (rA * 2 gave C)
 
91
1:      cp      rB0, rA0
 
92
        cpc     rB1, rA1
 
93
        cpc     rB2, rA2        ; Does test value fit?
 
94
        brcc    3f              ; NC --> nope, bit is 0
 
95
2:      sub     rA0, rB0
 
96
        sbc     rA1, rB1
 
97
        sbc     rA2, rB2        ; Prepare remainder argument for next bits
 
98
        or      rB0, msk0
 
99
        or      rB1, msk1
 
100
        or      rB2, msk2       ; Set developing bit to 1
 
101
3:      lsr     msk2
 
102
        ror     msk1
 
103
        ror     msk0            ; Shift right mask, C --> end loop
 
104
        eor     rB0, msk0
 
105
        eor     rB1, msk1
 
106
        eor     rB2, msk2       ; Shift right test bit in developing root
 
107
        brcc    .Loop           ; Develop 23 bits of the sqrt
 
108
 
 
109
        ; Loop for bit 0 and rounding
 
110
.Loop1: lsl     rA0
 
111
        rol     rA1
 
112
        rol     rA2             ; Shift left remainder argument 
 
113
        brcs    4f              ; C--> Last bits always 1
 
114
        cp      rB0, rA0
 
115
        cpc     rB1, rA1
 
116
        cpc     rB2, rA2        ; Test for last bit 1
 
117
        brcc    5f              ; Nope, stays the same
 
118
4:      sbc     rA0, rB0        ; MUST BE SBC !!
 
119
        sbc     rA1, rB1
 
120
        sbc     rA2, rB2        ; Prepare remainder argument for next bit
 
121
        add     rB0, msk0
 
122
        adc     rB1, msk1
 
123
        adc     rB2, msk1       ; Add 1 to result
 
124
5:      com     msk2            ; ZF if second time
 
125
        brne    .Loop1          ; 1 for last bit, 1 for rounding
 
126
 
 
127
#else   /* vs. to OPTIMIZE_SPEED        */
 
128
  ;** Optimized for size (9 code words smaller than speed optimized,
 
129
  ; 67 more cycles in average)
 
130
 
 
131
#define tv      rAE
 
132
 
 
133
        clr     tv              ; Test value for end of loop
 
134
        subi    rA2, 0x40       ; Initial remainder for odd exponent
 
135
        lsr     rB3
 
136
        ror     rA3             ; Divide exponent by 2, C==>exponent was odd
 
137
        brcs    3f              ; Jump for odd exponent in argument
 
138
        subi    rA2, 0x40       ; Initial remainder for even exponent, C=0
 
139
 
 
140
  ; Loop for all 24 bits
 
141
.Loop:  brcc    2f              ; NC --> nope, bit is 0
 
142
        cp      msk2, tv        ; Only needed to get the proper rounding
 
143
                                ;   for ffffff
 
144
        sbc     rA0, rB0
 
145
        sbc     rA1, rB1
 
146
        sbc     rA2, rB2        ; Prepare remainder argument for next bits
 
147
        or      rB0, msk0
 
148
        or      rB1, msk1
 
149
        or      rB2, msk2       ; Set developing bit to 1
 
150
2:      lsr     msk2
 
151
        ror     msk1
 
152
        ror     msk0            ; Shift right mask, C --> end loop
 
153
        rol     tv              ; Bit 1 set if end of loop
 
154
        eor     rB0, msk0
 
155
        eor     rB1, msk1
 
156
        eor     rB2, msk2       ; Shift right test bit in developing root
 
157
3:      lsl     rA0
 
158
        rol     rA1
 
159
        rol     rA2             ; Shift left remainder argument (C used
 
160
                                ;   at .Loop)
 
161
        brcs    4f
 
162
        cp      rB0, rA0
 
163
        cpc     rB1, rA1
 
164
        cpc     rB2, rA2
 
165
 
 
166
4:      sbrs    tv, 1
 
167
        rjmp    .Loop
 
168
 
 
169
  ; Rounding
 
170
        adc     rB0, msk1
 
171
        adc     rB1, msk1
 
172
        adc     rB2, msk1       ; Rounded mantissa ready (msk1=0)
 
173
 
 
174
#endif  /* !OPTIMIZE_SPEED */
 
175
 
146
176
        X_movw  rA0, rB0
147
 
        mov     rA2, rB2
148
 
        subi    rA3, lo8(-127)          ; exponent bias
 
177
        mov     rA2, rB2        ; Copy to rA
 
178
 
 
179
        subi    rA3, lo8(-127)  ; exponent bias
149
180
        lsl     rA2
150
181
        lsr     rA3
151
182
        ror     rA2