Return to wm_sqrt.S CVS log

Up to [MW Coherent from dump] / coherent / b / kernel / emulator

Annotation of coherent/b/kernel/emulator/wm_sqrt.S, revision 1.1

1.1     ! root        1:        .file   "wm_sqrt.S"
        !             2: /*---------------------------------------------------------------------------+
        !             3:  |  wm_sqrt.S                                                                |
        !             4:  |                                                                           |
        !             5:  | Fixed point arithmetic square root evaluation.                            |
        !             6:  |                                                                           |
        !             7:  | Copyright (C) 1992    W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
        !             8:  |                       Australia.  E-mail [email protected]    |
        !             9:  |                                                                           |
        !            10:  | Call from C as:                                                           |
        !            11:  |   void wm_sqrt(FPU_REG *n)                                                |
        !            12:  |                                                                           |
        !            13:  +---------------------------------------------------------------------------*/
        !            14: 
        !            15: /*---------------------------------------------------------------------------+
        !            16:  |  wm_sqrt(FPU_REG *n)                                                      |
        !            17:  |    returns the square root of n in n.                                     |
        !            18:  |                                                                           |
        !            19:  |  Use Newton's method to compute the square root of a number, which must   |
        !            20:  |  be in the range  [1.0 .. 4.0),  to 64 bits accuracy.                     |
        !            21:  |  Does not check the sign or tag of the argument.                          |
        !            22:  |  Sets the exponent, but not the sign or tag of the result.                |
        !            23:  |                                                                           |
        !            24:  |  The guess is kept in %esi:%edi                                           |
        !            25:  +---------------------------------------------------------------------------*/
        !            26: 
        !            27: #include "exception.h"
        !            28: #include "fpu_asm.h"
        !            29: 
        !            30: 
        !            31: .data
        !            32: /*
        !            33:        Local storage:
        !            34:  */
        !            35:        .align 4,0
        !            36: accum_2:
        !            37:        .long   0               // ms word
        !            38: accum_1:
        !            39:        .long   0
        !            40: accum_0:
        !            41:        .long   0
        !            42: 
        !            43: sq_2:
        !            44:        .long   0               // ms word
        !            45: sq_1:
        !            46:        .long   0
        !            47: sq_0:
        !            48:        .long   0
        !            49: 
        !            50: .text
        !            51:        .align 2,144
        !            52: 
        !            53: .globl wm_sqrt
        !            54: 
        !            55: wm_sqrt:
        !            56:        pushl   %ebp
        !            57:        movl    %esp,%ebp
        !            58:        pushl   %esi
        !            59:        pushl   %edi
        !            60: //     pushl   %ebx
        !            61: 
        !            62:        movl    PARAM1,%esi
        !            63: 
        !            64:        movl    SIGH(%esi),%eax
        !            65:        movl    SIGL(%esi),%ecx
        !            66:        xorl    %edx,%edx
        !            67: 
        !            68: // We use a rough linear estimate for the first guess..
        !            69: 
        !            70:        cmpl    EXP_BIAS,EXP(%esi)
        !            71:        jnz     L10
        !            72: 
        !            73:        shrl    $1,%eax                 /* arg is in the range  [1.0 .. 2.0) */
        !            74:        rcrl    $1,%ecx
        !            75:        rcrl    $1,%edx
        !            76: 
        !            77: L10:
        !            78: // From here on, n is never accessed directly again until it is
        !            79: // replaced by the answer.
        !            80: 
        !            81:        movl    %eax,sq_2               // ms word of n
        !            82:        movl    %ecx,sq_1
        !            83:        movl    %edx,sq_0
        !            84: 
        !            85:        shrl    $1,%eax
        !            86:        addl    $0x40000000,%eax
        !            87:        movl    $0xaaaaaaaa,%ecx
        !            88:        mull    %ecx
        !            89:        shll    $1,%edx                 /* max result was 7fff... */
        !            90:        testl   $0x80000000,%edx        /* but min was 3fff... */
        !            91:        jnz     no_adjust
        !            92: 
        !            93:        movl    $0x80000000,%edx        /* round up */
        !            94: 
        !            95: no_adjust:
        !            96:        movl    %edx,%esi       // Our first guess
        !            97: 
        !            98: /* We have now computed (approx)   (2 + x) / 3, which forms the basis
        !            99:    for a few iterations of Newton's method */
        !           100: 
        !           101:        movl    sq_2,%ecx       // ms word
        !           102: 
        !           103: // From our initial estimate, three iterations are enough to get us
        !           104: // to 30 bits or so. This will then allow two iterations at better
        !           105: // precision to complete the process.
        !           106: 
        !           107: // Compute  (g + n/g)/2  at each iteration (g is the guess).
        !           108:        shrl    $1,%ecx         // Doing this first will prevent a divide
        !           109:                                // overflow later.
        !           110: 
        !           111:        movl    %ecx,%edx
        !           112:        divl    %esi
        !           113:        shrl    $1,%esi
        !           114:        addl    %eax,%esi
        !           115: 
        !           116:        movl    %ecx,%edx
        !           117:        divl    %esi
        !           118:        shrl    $1,%esi
        !           119:        addl    %eax,%esi
        !           120: 
        !           121:        movl    %ecx,%edx
        !           122:        divl    %esi
        !           123:        shrl    $1,%esi
        !           124:        addl    %eax,%esi
        !           125: 
        !           126: // Now that an estimate accurate to about 30 bits has been obtained,
        !           127: // we improve it to 60 bits or so.
        !           128: 
        !           129: // The strategy from now on is to compute new estimates from
        !           130: //      guess := guess + (n - guess^2) / (2 * guess)
        !           131: 
        !           132: // First, find the square of the guess
        !           133:        movl    %esi,%eax
        !           134:        mull    %esi
        !           135: // guess^2 now in %edx:%eax
        !           136: 
        !           137:        movl    sq_1,%ecx
        !           138:        subl    %ecx,%eax
        !           139:        movl    sq_2,%ecx       // ms word of normalized n
        !           140:        sbbl    %ecx,%edx
        !           141:        jnc     l40
        !           142: 
        !           143: // subtraction gives a negative result
        !           144: // negate the reult before division
        !           145:        notl    %edx
        !           146:        notl    %eax
        !           147:        addl    $1,%eax
        !           148:        adcl    $0,%edx
        !           149: 
        !           150:        divl    %esi
        !           151:        movl    %eax,%ecx
        !           152: 
        !           153:        movl    %edx,%eax
        !           154:        divl    %esi
        !           155:        jmp     l45
        !           156: 
        !           157: l40:
        !           158:        divl    %esi
        !           159:        movl    %eax,%ecx
        !           160: 
        !           161:        movl    %edx,%eax
        !           162:        divl    %esi
        !           163: 
        !           164:        notl    %ecx
        !           165:        notl    %eax
        !           166:        addl    $1,%eax
        !           167:        adcl    $0,%ecx
        !           168: 
        !           169: l45:
        !           170:        sarl    $1,%ecx         // divide by 2
        !           171:        rcrl    $1,%eax
        !           172: 
        !           173:        movl    %eax,%edi
        !           174:        addl    %ecx,%esi
        !           175: 
        !           176: // Now the square root has been computed to better than 60 bits
        !           177: 
        !           178: // Find the square of the guess
        !           179:        movl    %edi,%eax               // ls word of guess
        !           180:        mull    %edi
        !           181:        movl    %edx,accum_0
        !           182: 
        !           183:        movl    %esi,%eax
        !           184:        mull    %esi
        !           185:        movl    %edx,accum_2
        !           186:        movl    %eax,accum_1
        !           187: 
        !           188:        movl    %edi,%eax
        !           189:        mull    %esi
        !           190:        addl    %eax,accum_0
        !           191:        adcl    %edx,accum_1
        !           192:        adcl    $0,accum_2
        !           193: 
        !           194:        movl    %esi,%eax
        !           195:        mull    %edi
        !           196:        addl    %eax,accum_0
        !           197:        adcl    %edx,accum_1
        !           198:        adcl    $0,accum_2
        !           199: 
        !           200: // guess^2 now in accum_2:accum_1:accum_0
        !           201: 
        !           202:        movl    sq_1,%eax               // get normalized n
        !           203:        subl    %eax,accum_1
        !           204:        movl    sq_2,%eax               // ms word of normalized n
        !           205:        sbbl    %eax,accum_2
        !           206:        jnc     l60
        !           207: 
        !           208: // subtraction gives a negative result
        !           209: // negate the reult before division
        !           210:        notl    accum_0
        !           211:        notl    accum_1
        !           212:        notl    accum_2
        !           213:        addl    $1,accum_0
        !           214:        adcl    $0,accum_1
        !           215: 
        !           216: #ifdef PARANOID
        !           217:        adcl    $0,accum_2      // This must be zero
        !           218:        jz      l51
        !           219: 
        !           220:        pushl   EX_INTERNAL|0x207
        !           221:        call    EXCEPTION
        !           222: 
        !           223: l51:
        !           224: #endif PARANOID
        !           225: 
        !           226:        movl    accum_1,%edx
        !           227:        movl    accum_0,%eax
        !           228:        divl    %esi
        !           229:        movl    %eax,%ecx
        !           230: 
        !           231:        movl    %edx,%eax
        !           232:        divl    %esi
        !           233: 
        !           234:        sarl    $1,%ecx         // divide by 2
        !           235:        rcrl    $1,%eax
        !           236: 
        !           237:        // round the result
        !           238:        addl    $0x80000000,%eax
        !           239:        adcl    $0,%ecx
        !           240: 
        !           241:        addl    %ecx,%edi
        !           242:        adcl    $0,%esi
        !           243: 
        !           244:        jmp     l65
        !           245: 
        !           246: l60:
        !           247:        movl    accum_1,%edx
        !           248:        movl    accum_0,%eax
        !           249:        divl    %esi
        !           250:        movl    %eax,%ecx
        !           251: 
        !           252:        movl    %edx,%eax
        !           253:        divl    %esi
        !           254: 
        !           255:        sarl    $1,%ecx         // divide by 2
        !           256:        rcrl    $1,%eax
        !           257: 
        !           258:        // round the result
        !           259:        addl    $0x80000000,%eax
        !           260:        adcl    $0,%ecx
        !           261: 
        !           262:        subl    %ecx,%edi
        !           263:        sbbl    $0,%esi
        !           264: 
        !           265: l65:
        !           266:        movl    PARAM1,%ecx
        !           267: 
        !           268:        movl    %edi,SIGL(%ecx)
        !           269:        movl    %esi,SIGH(%ecx)
        !           270: 
        !           271:        movl    EXP_BIAS,EXP(%ecx)      /* Result is in  [1.0 .. 2.0) */
        !           272: 
        !           273: //     popl    %ebx
        !           274:        popl    %edi
        !           275:        popl    %esi
        !           276:        leave
        !           277:        ret

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.