| 0,0 → 1,145 |
|---|
| #include<libc/asm.h> |
| |
| /* Copyright (C) 1999 DJ Delorie, see COPYING.DJ for details */ |
| .data |
| pinf: |
| .long 0xFF800000 |
| |
| NaN: |
| .long 0xFFC00000 |
| |
| temp: |
| .long 0, 0 |
| |
| onethird: |
| .long 1431655765 |
| |
| two54: |
| .long 0x5A800000 |
| |
| a0: |
| .float +1.87277957900533E+00 |
| |
| a1: |
| .float -1.87243905326548E+00 |
| |
| a2: |
| .float +1.60286399719912E+00 |
| |
| a3: |
| .float -7.46198924594210E-01 |
| |
| a4: |
| .float +1.42994392730009E-01 |
| |
| b0: |
| .float 14. |
| |
| b1: |
| .float -7. |
| |
| b2: |
| .float +2. |
| |
| one9th: |
| .tfloat +0.11111111111111111111 |
| |
| .text |
| MK_C_SYM(cbrt) |
| |
| movl 8(%esp), %eax |
| movl %eax, %ecx /* Save sign */ |
| |
| andl $0x7FFFFFFF, %eax /* fabs */ |
| movl %eax, 8(%esp) |
| |
| cmpl $0x7FF00000, %eax /* Control flows straight through for */ |
| jae abarg /* normal args: 0 < fabs(x) < +inf */ |
| testl $0x7FF00000, %eax |
| jz verysmall |
| |
| mull onethird |
| addl $0x2A9F7893, %edx |
| movl %edx, temp+4 /* First approximation good */ |
| /* to 5.5 bits */ |
| |
| have55: |
| fldl 4(%esp) |
| fld1 |
| fdivp /* recip */ |
| |
| fldl temp /* Load approximation */ |
| /* 4rd-order minimax to 24 bits */ |
| fld %st(0) /* x x recip */ |
| fmul %st(1) /* x^2 x recip */ |
| fmul %st(1) /* x^3 x recip */ |
| fmul %st(2) /* y x recip */ |
| fld %st(0) /* y y x recip */ |
| fmuls a4 /* P1' y x recip */ |
| fadds a3 /* P1 y x recip */ |
| fmul %st(1) /* P2' y x recip */ |
| fadds a2 /* P2 y x recip */ |
| fmul %st(1) /* P3' y x recip */ |
| fadds a1 /* P3 y x recip */ |
| fmulp /* P4' x recip */ |
| fadds a0 /* P4 x recip */ |
| fmulp /* x' recip */ |
| /* 2nd-order Taylor to 64 bits */ |
| fld %st(0) /* x x recip */ |
| fmul %st(1) /* x^2 x recip */ |
| fmul %st(1) /* x^3 x recip */ |
| fmul %st(2) /* y x recip */ |
| ffree %st(2) /* y x */ |
| fld %st(0) /* y y x */ |
| fmuls b2 |
| fadds b1 |
| fmulp /* ccc x */ |
| fadds b0 /* P(y) x */ |
| fmulp /* x'' */ |
| fldt one9th |
| fmulp |
| |
| cleanup: /* Restore sign */ |
| testl %ecx, %ecx |
| jns end |
| fchs |
| |
| end: |
| ret |
| |
| verysmall: /* Exponent is 0 */ |
| movl 8(%esp), %eax |
| testl %eax, %eax |
| jnz denormal |
| movl 4(%esp), %eax |
| testl %eax, %eax |
| jz special /* x = 0 */ |
| |
| denormal: |
| fldl 4(%esp) |
| fmuls two54 /* Multiply by 2^54 to normalize */ |
| fstpl temp |
| |
| movl temp+4, %eax |
| mull onethird |
| addl $0x297F7893, %edx /* Undo 2^54 multiplier */ |
| movl %edx, temp+4 /* First approximation to 5.5 bits */ |
| movl $0, temp |
| |
| jmp have55 |
| |
| abarg: /* x = inf, or NaN */ |
| testl $0x000FFFFF, %eax |
| jnz badarg |
| movl 4(%esp), %eax |
| testl %eax, %eax |
| jz special |
| |
| badarg: /* arg is negative or NaN */ |
| movl $1, C_SYM(errno) |
| flds NaN |
| ret |
| |
| special: |
| fldl 4(%esp) /* x = 0 or inf: just load x */ |
| jmp cleanup |