| .file "wm_sqrt.S" | |
| /*---------------------------------------------------------------------------+ | |
| | wm_sqrt.S | | |
| | | | |
| | Fixed point arithmetic square root evaluation. | | |
| | | | |
| | Copyright (C) 1992,1993,1995,1997 | | |
| | W. Metzenthen, 22 Parker St, Ormond, Vic 3163, | | |
| | Australia. E-mail billm@suburbia.net | | |
| | | | |
| | Call from C as: | | |
| | int wm_sqrt(FPU_REG *n, unsigned int control_word) | | |
| | | | |
| +---------------------------------------------------------------------------*/ | |
| /*---------------------------------------------------------------------------+ | |
| | wm_sqrt(FPU_REG *n, unsigned int control_word) | | |
| | returns the square root of n in n. | | |
| | | | |
| | Use Newton's method to compute the square root of a number, which must | | |
| | be in the range [1.0 .. 4.0), to 64 bits accuracy. | | |
| | Does not check the sign or tag of the argument. | | |
| | Sets the exponent, but not the sign or tag of the result. | | |
| | | | |
| | The guess is kept in %esi:%edi | | |
| +---------------------------------------------------------------------------*/ | |
| #include "exception.h" | |
| #include "fpu_emu.h" | |
| #ifndef NON_REENTRANT_FPU | |
| /* Local storage on the stack: */ | |
| #define FPU_accum_3 -4(%ebp) /* ms word */ | |
| #define FPU_accum_2 -8(%ebp) | |
| #define FPU_accum_1 -12(%ebp) | |
| #define FPU_accum_0 -16(%ebp) | |
| /* | |
| * The de-normalised argument: | |
| * sq_2 sq_1 sq_0 | |
| * b b b b b b b ... b b b b b b .... b b b b 0 0 0 ... 0 | |
| * ^ binary point here | |
| */ | |
| #define FPU_fsqrt_arg_2 -20(%ebp) /* ms word */ | |
| #define FPU_fsqrt_arg_1 -24(%ebp) | |
| #define FPU_fsqrt_arg_0 -28(%ebp) /* ls word, at most the ms bit is set */ | |
| #else | |
| /* Local storage in a static area: */ | |
| .data | |
| .align 4,0 | |
| FPU_accum_3: | |
| .long 0 /* ms word */ | |
| FPU_accum_2: | |
| .long 0 | |
| FPU_accum_1: | |
| .long 0 | |
| FPU_accum_0: | |
| .long 0 | |
| /* The de-normalised argument: | |
| sq_2 sq_1 sq_0 | |
| b b b b b b b ... b b b b b b .... b b b b 0 0 0 ... 0 | |
| ^ binary point here | |
| */ | |
| FPU_fsqrt_arg_2: | |
| .long 0 /* ms word */ | |
| FPU_fsqrt_arg_1: | |
| .long 0 | |
| FPU_fsqrt_arg_0: | |
| .long 0 /* ls word, at most the ms bit is set */ | |
| #endif /* NON_REENTRANT_FPU */ | |
| .text | |
| ENTRY(wm_sqrt) | |
| pushl %ebp | |
| movl %esp,%ebp | |
| #ifndef NON_REENTRANT_FPU | |
| subl $28,%esp | |
| #endif /* NON_REENTRANT_FPU */ | |
| pushl %esi | |
| pushl %edi | |
| pushl %ebx | |
| movl PARAM1,%esi | |
| movl SIGH(%esi),%eax | |
| movl SIGL(%esi),%ecx | |
| xorl %edx,%edx | |
| /* We use a rough linear estimate for the first guess.. */ | |
| cmpw EXP_BIAS,EXP(%esi) | |
| jnz sqrt_arg_ge_2 | |
| shrl $1,%eax /* arg is in the range [1.0 .. 2.0) */ | |
| rcrl $1,%ecx | |
| rcrl $1,%edx | |
| sqrt_arg_ge_2: | |
| /* From here on, n is never accessed directly again until it is | |
| replaced by the answer. */ | |
| movl %eax,FPU_fsqrt_arg_2 /* ms word of n */ | |
| movl %ecx,FPU_fsqrt_arg_1 | |
| movl %edx,FPU_fsqrt_arg_0 | |
| /* Make a linear first estimate */ | |
| shrl $1,%eax | |
| addl $0x40000000,%eax | |
| movl $0xaaaaaaaa,%ecx | |
| mull %ecx | |
| shll %edx /* max result was 7fff... */ | |
| testl $0x80000000,%edx /* but min was 3fff... */ | |
| jnz sqrt_prelim_no_adjust | |
| movl $0x80000000,%edx /* round up */ | |
| sqrt_prelim_no_adjust: | |
| movl %edx,%esi /* Our first guess */ | |
| /* We have now computed (approx) (2 + x) / 3, which forms the basis | |
| for a few iterations of Newton's method */ | |
| movl FPU_fsqrt_arg_2,%ecx /* ms word */ | |
| /* | |
| * From our initial estimate, three iterations are enough to get us | |
| * to 30 bits or so. This will then allow two iterations at better | |
| * precision to complete the process. | |
| */ | |
| /* Compute (g + n/g)/2 at each iteration (g is the guess). */ | |
| shrl %ecx /* Doing this first will prevent a divide */ | |
| /* overflow later. */ | |
| movl %ecx,%edx /* msw of the arg / 2 */ | |
| divl %esi /* current estimate */ | |
| shrl %esi /* divide by 2 */ | |
| addl %eax,%esi /* the new estimate */ | |
| movl %ecx,%edx | |
| divl %esi | |
| shrl %esi | |
| addl %eax,%esi | |
| movl %ecx,%edx | |
| divl %esi | |
| shrl %esi | |
| addl %eax,%esi | |
| /* | |
| * Now that an estimate accurate to about 30 bits has been obtained (in %esi), | |
| * we improve it to 60 bits or so. | |
| * | |
| * The strategy from now on is to compute new estimates from | |
| * guess := guess + (n - guess^2) / (2 * guess) | |
| */ | |
| /* First, find the square of the guess */ | |
| movl %esi,%eax | |
| mull %esi | |
| /* guess^2 now in %edx:%eax */ | |
| movl FPU_fsqrt_arg_1,%ecx | |
| subl %ecx,%eax | |
| movl FPU_fsqrt_arg_2,%ecx /* ms word of normalized n */ | |
| sbbl %ecx,%edx | |
| jnc sqrt_stage_2_positive | |
| /* Subtraction gives a negative result, | |
| negate the result before division. */ | |
| notl %edx | |
| notl %eax | |
| addl $1,%eax | |
| adcl $0,%edx | |
| divl %esi | |
| movl %eax,%ecx | |
| movl %edx,%eax | |
| divl %esi | |
| jmp sqrt_stage_2_finish | |
| sqrt_stage_2_positive: | |
| divl %esi | |
| movl %eax,%ecx | |
| movl %edx,%eax | |
| divl %esi | |
| notl %ecx | |
| notl %eax | |
| addl $1,%eax | |
| adcl $0,%ecx | |
| sqrt_stage_2_finish: | |
| sarl $1,%ecx /* divide by 2 */ | |
| rcrl $1,%eax | |
| /* Form the new estimate in %esi:%edi */ | |
| movl %eax,%edi | |
| addl %ecx,%esi | |
| jnz sqrt_stage_2_done /* result should be [1..2) */ | |
| #ifdef PARANOID | |
| /* It should be possible to get here only if the arg is ffff....ffff */ | |
| cmp $0xffffffff,FPU_fsqrt_arg_1 | |
| jnz sqrt_stage_2_error | |
| #endif /* PARANOID */ | |
| /* The best rounded result. */ | |
| xorl %eax,%eax | |
| decl %eax | |
| movl %eax,%edi | |
| movl %eax,%esi | |
| movl $0x7fffffff,%eax | |
| jmp sqrt_round_result | |
| #ifdef PARANOID | |
| sqrt_stage_2_error: | |
| pushl EX_INTERNAL|0x213 | |
| call EXCEPTION | |
| #endif /* PARANOID */ | |
| sqrt_stage_2_done: | |
| /* Now the square root has been computed to better than 60 bits. */ | |
| /* Find the square of the guess. */ | |
| movl %edi,%eax /* ls word of guess */ | |
| mull %edi | |
| movl %edx,FPU_accum_1 | |
| movl %esi,%eax | |
| mull %esi | |
| movl %edx,FPU_accum_3 | |
| movl %eax,FPU_accum_2 | |
| movl %edi,%eax | |
| mull %esi | |
| addl %eax,FPU_accum_1 | |
| adcl %edx,FPU_accum_2 | |
| adcl $0,FPU_accum_3 | |
| /* movl %esi,%eax */ | |
| /* mull %edi */ | |
| addl %eax,FPU_accum_1 | |
| adcl %edx,FPU_accum_2 | |
| adcl $0,FPU_accum_3 | |
| /* guess^2 now in FPU_accum_3:FPU_accum_2:FPU_accum_1 */ | |
| movl FPU_fsqrt_arg_0,%eax /* get normalized n */ | |
| subl %eax,FPU_accum_1 | |
| movl FPU_fsqrt_arg_1,%eax | |
| sbbl %eax,FPU_accum_2 | |
| movl FPU_fsqrt_arg_2,%eax /* ms word of normalized n */ | |
| sbbl %eax,FPU_accum_3 | |
| jnc sqrt_stage_3_positive | |
| /* Subtraction gives a negative result, | |
| negate the result before division */ | |
| notl FPU_accum_1 | |
| notl FPU_accum_2 | |
| notl FPU_accum_3 | |
| addl $1,FPU_accum_1 | |
| adcl $0,FPU_accum_2 | |
| #ifdef PARANOID | |
| adcl $0,FPU_accum_3 /* This must be zero */ | |
| jz sqrt_stage_3_no_error | |
| sqrt_stage_3_error: | |
| pushl EX_INTERNAL|0x207 | |
| call EXCEPTION | |
| sqrt_stage_3_no_error: | |
| #endif /* PARANOID */ | |
| movl FPU_accum_2,%edx | |
| movl FPU_accum_1,%eax | |
| divl %esi | |
| movl %eax,%ecx | |
| movl %edx,%eax | |
| divl %esi | |
| sarl $1,%ecx /* divide by 2 */ | |
| rcrl $1,%eax | |
| /* prepare to round the result */ | |
| addl %ecx,%edi | |
| adcl $0,%esi | |
| jmp sqrt_stage_3_finished | |
| sqrt_stage_3_positive: | |
| movl FPU_accum_2,%edx | |
| movl FPU_accum_1,%eax | |
| divl %esi | |
| movl %eax,%ecx | |
| movl %edx,%eax | |
| divl %esi | |
| sarl $1,%ecx /* divide by 2 */ | |
| rcrl $1,%eax | |
| /* prepare to round the result */ | |
| notl %eax /* Negate the correction term */ | |
| notl %ecx | |
| addl $1,%eax | |
| adcl $0,%ecx /* carry here ==> correction == 0 */ | |
| adcl $0xffffffff,%esi | |
| addl %ecx,%edi | |
| adcl $0,%esi | |
| sqrt_stage_3_finished: | |
| /* | |
| * The result in %esi:%edi:%esi should be good to about 90 bits here, | |
| * and the rounding information here does not have sufficient accuracy | |
| * in a few rare cases. | |
| */ | |
| cmpl $0xffffffe0,%eax | |
| ja sqrt_near_exact_x | |
| cmpl $0x00000020,%eax | |
| jb sqrt_near_exact | |
| cmpl $0x7fffffe0,%eax | |
| jb sqrt_round_result | |
| cmpl $0x80000020,%eax | |
| jb sqrt_get_more_precision | |
| sqrt_round_result: | |
| /* Set up for rounding operations */ | |
| movl %eax,%edx | |
| movl %esi,%eax | |
| movl %edi,%ebx | |
| movl PARAM1,%edi | |
| movw EXP_BIAS,EXP(%edi) /* Result is in [1.0 .. 2.0) */ | |
| jmp fpu_reg_round | |
| sqrt_near_exact_x: | |
| /* First, the estimate must be rounded up. */ | |
| addl $1,%edi | |
| adcl $0,%esi | |
| sqrt_near_exact: | |
| /* | |
| * This is an easy case because x^1/2 is monotonic. | |
| * We need just find the square of our estimate, compare it | |
| * with the argument, and deduce whether our estimate is | |
| * above, below, or exact. We use the fact that the estimate | |
| * is known to be accurate to about 90 bits. | |
| */ | |
| movl %edi,%eax /* ls word of guess */ | |
| mull %edi | |
| movl %edx,%ebx /* 2nd ls word of square */ | |
| movl %eax,%ecx /* ls word of square */ | |
| movl %edi,%eax | |
| mull %esi | |
| addl %eax,%ebx | |
| addl %eax,%ebx | |
| #ifdef PARANOID | |
| cmp $0xffffffb0,%ebx | |
| jb sqrt_near_exact_ok | |
| cmp $0x00000050,%ebx | |
| ja sqrt_near_exact_ok | |
| pushl EX_INTERNAL|0x214 | |
| call EXCEPTION | |
| sqrt_near_exact_ok: | |
| #endif /* PARANOID */ | |
| or %ebx,%ebx | |
| js sqrt_near_exact_small | |
| jnz sqrt_near_exact_large | |
| or %ebx,%edx | |
| jnz sqrt_near_exact_large | |
| /* Our estimate is exactly the right answer */ | |
| xorl %eax,%eax | |
| jmp sqrt_round_result | |
| sqrt_near_exact_small: | |
| /* Our estimate is too small */ | |
| movl $0x000000ff,%eax | |
| jmp sqrt_round_result | |
| sqrt_near_exact_large: | |
| /* Our estimate is too large, we need to decrement it */ | |
| subl $1,%edi | |
| sbbl $0,%esi | |
| movl $0xffffff00,%eax | |
| jmp sqrt_round_result | |
| sqrt_get_more_precision: | |
| /* This case is almost the same as the above, except we start | |
| with an extra bit of precision in the estimate. */ | |
| stc /* The extra bit. */ | |
| rcll $1,%edi /* Shift the estimate left one bit */ | |
| rcll $1,%esi | |
| movl %edi,%eax /* ls word of guess */ | |
| mull %edi | |
| movl %edx,%ebx /* 2nd ls word of square */ | |
| movl %eax,%ecx /* ls word of square */ | |
| movl %edi,%eax | |
| mull %esi | |
| addl %eax,%ebx | |
| addl %eax,%ebx | |
| /* Put our estimate back to its original value */ | |
| stc /* The ms bit. */ | |
| rcrl $1,%esi /* Shift the estimate left one bit */ | |
| rcrl $1,%edi | |
| #ifdef PARANOID | |
| cmp $0xffffff60,%ebx | |
| jb sqrt_more_prec_ok | |
| cmp $0x000000a0,%ebx | |
| ja sqrt_more_prec_ok | |
| pushl EX_INTERNAL|0x215 | |
| call EXCEPTION | |
| sqrt_more_prec_ok: | |
| #endif /* PARANOID */ | |
| or %ebx,%ebx | |
| js sqrt_more_prec_small | |
| jnz sqrt_more_prec_large | |
| or %ebx,%ecx | |
| jnz sqrt_more_prec_large | |
| /* Our estimate is exactly the right answer */ | |
| movl $0x80000000,%eax | |
| jmp sqrt_round_result | |
| sqrt_more_prec_small: | |
| /* Our estimate is too small */ | |
| movl $0x800000ff,%eax | |
| jmp sqrt_round_result | |
| sqrt_more_prec_large: | |
| /* Our estimate is too large */ | |
| movl $0x7fffff00,%eax | |
| jmp sqrt_round_result |