.mmregs
             .include "tab_ld8a.inc"
             .global  Inv_sqrt
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Compute 1/sqrt(L_x)
;Algorithm:                                                               |
; |                                                                           |
; |   The function 1/sqrt(L_x) is approximated by a table and linear          |
; |   interpolation.                                                          |
; |                                                                           |
; |   1- Normalization of L_x.                                                |
 ;|   2- If (30-exponent) is even then shift right once.                      |
 ;|   3- exponent = (30-exponent)/2  +1                                       |
 ;|   4- i = bit25-b31 of L_x,    16 <= i <= 63  ->because of normalization.  |
 ;|   5- a = bit10-b24                                                        |
 ;|   6- i -=16                                                               |
 ;|   7- L_y = tabsqr[i]<<16 - (tabsqr[i] - tabsqr[i+1]) * a * 2                 |
; |   8- L_y >>= exponent 
;Word32 Inv_sqrt(   /* (o) Q30 : output value   (range: 0<=val<1)           */
;  Word32 L_x       /* (i) Q0  : input value    (range: 0<=val<=7fffffff)   */
;)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
Inv_sqrt:
    frame -3    ;*sp(0)--exp,*sp(1)--i,*sp(2)--a
                 
    bc   Inv_sqrt_start,agt
    ld   #3fffh,16,a
    or   #0ffffh,a
    b    Inv_sqrt_end
Inv_sqrt_start:
    exp  a
    norm a
    ld   a,b
    st   t,*sp(0)
    ld   #30,16,a
    sub  *sp(0),16,a
    sth  a,*sp(0)
    and  #1000h,4,a
    xc   1,aeq
    ld   b,-1,b
    ld   *sp(0),-1,a
    add  #1,a             ;exp = add(exp, 1)
    stl  a,*sp(0)
    
    ld   b,-9,b
    sth  b,*sp(1)
    
    ld   b,-1,b
    and  #7fffh,b
    stl  b,*sp(2)
    
    ld   *sp(1),a
    sub  #16,a
    stlm a,ar0
    stm  #tabsqr,ar2
    ld   *ar2+0,a
    
    ld   *ar2,16,b
    sub  *ar2(1),16,b,a
    ld   a,-16,a
    stlm  a,t
    mas   *sp(2),b
    
    ld   *sp(0),a
    neg  a
    stl  a,*sp(0)
    ld   *sp(0),asm
    ld   b,asm,a
    
        
Inv_sqrt_end:
    frame  3
    ret
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    .end