Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

  1. #include<libc/asm.h>
  2.  
  3. /* Copyright (C) 1999 DJ Delorie, see COPYING.DJ for details */
  4.         .data
  5. pinf:
  6.         .long   0xFF800000
  7.  
  8. NaN:
  9.         .long   0xFFC00000
  10.  
  11. temp:
  12.         .long   0, 0
  13.  
  14. onethird:
  15.         .long   1431655765
  16.  
  17. two54:
  18.         .long   0x5A800000
  19.  
  20. a0:
  21.         .float  +1.87277957900533E+00
  22.  
  23. a1:
  24.         .float  -1.87243905326548E+00
  25.  
  26. a2:
  27.         .float  +1.60286399719912E+00
  28.  
  29. a3:
  30.         .float  -7.46198924594210E-01
  31.  
  32. a4:
  33.         .float  +1.42994392730009E-01
  34.  
  35. b0:
  36.         .float  14.
  37.  
  38. b1:
  39.         .float  -7.
  40.  
  41. b2:
  42.         .float  +2.
  43.  
  44. one9th:
  45.         .tfloat +0.11111111111111111111
  46.  
  47.         .text
  48. MK_C_SYM(cbrt)
  49.  
  50.         movl    8(%esp), %eax
  51.         movl    %eax, %ecx              /* Save sign */
  52.  
  53.         andl    $0x7FFFFFFF, %eax       /* fabs */
  54.         movl    %eax, 8(%esp)
  55.  
  56.         cmpl    $0x7FF00000, %eax       /* Control flows straight through for */
  57.         jae     abarg                   /* normal args: 0 < fabs(x) < +inf */
  58.         testl   $0x7FF00000, %eax
  59.         jz      verysmall
  60.  
  61.         mull    onethird
  62.         addl    $0x2A9F7893, %edx
  63.         movl    %edx, temp+4            /* First approximation good */
  64.                                         /* to 5.5 bits */
  65.  
  66. have55:
  67.         fldl    4(%esp)
  68.         fld1
  69.         fdivp                           /* recip */
  70.  
  71.         fldl    temp                    /* Load approximation */
  72.                                         /* 4rd-order minimax to 24 bits */
  73.         fld     %st(0)                  /* x    x   recip       */
  74.         fmul    %st(1)                  /* x^2  x   recip       */
  75.         fmul    %st(1)                  /* x^3  x   recip       */
  76.         fmul    %st(2)                  /* y    x   recip       */
  77.         fld     %st(0)                  /* y    y   x     recip */
  78.         fmuls   a4                      /* P1'  y   x     recip */
  79.         fadds   a3                      /* P1   y   x     recip */
  80.         fmul    %st(1)                  /* P2'  y   x     recip */
  81.         fadds   a2                      /* P2   y   x     recip */
  82.         fmul    %st(1)                  /* P3'  y   x     recip */
  83.         fadds   a1                      /* P3   y   x     recip */
  84.         fmulp                           /* P4'  x   recip       */
  85.         fadds   a0                      /* P4   x   recip       */
  86.         fmulp                           /* x'   recip           */
  87.                                         /* 2nd-order Taylor to 64 bits */
  88.         fld     %st(0)                  /* x    x   recip */
  89.         fmul    %st(1)                  /* x^2  x   recip */
  90.         fmul    %st(1)                  /* x^3  x   recip */
  91.         fmul    %st(2)                  /* y    x   recip */
  92.         ffree   %st(2)                  /* y    x         */
  93.         fld     %st(0)                  /* y    y   x     */
  94.         fmuls   b2
  95.         fadds   b1
  96.         fmulp                           /* ccc  x */
  97.         fadds   b0                      /* P(y) x */
  98.         fmulp                           /* x''    */
  99.         fldt    one9th
  100.         fmulp
  101.  
  102. cleanup:                                /* Restore sign */
  103.         testl   %ecx, %ecx
  104.         jns     end
  105.         fchs
  106.  
  107. end:
  108.         ret
  109.  
  110. verysmall:                              /* Exponent is 0 */
  111.         movl    8(%esp), %eax
  112.         testl   %eax, %eax
  113.         jnz     denormal
  114.         movl    4(%esp), %eax
  115.         testl   %eax, %eax
  116.         jz      special                 /* x = 0 */
  117.        
  118. denormal:
  119.         fldl    4(%esp)
  120.         fmuls   two54                   /* Multiply by 2^54 to normalize */
  121.         fstpl   temp
  122.  
  123.         movl    temp+4, %eax
  124.         mull    onethird
  125.         addl    $0x297F7893, %edx       /* Undo 2^54 multiplier */
  126.         movl    %edx, temp+4            /* First approximation to 5.5 bits */
  127.         movl    $0, temp
  128.  
  129.         jmp     have55
  130.  
  131. abarg:                                  /* x = inf, or NaN */
  132.         testl   $0x000FFFFF, %eax
  133.         jnz     badarg
  134.         movl    4(%esp), %eax
  135.         testl   %eax, %eax
  136.         jz      special
  137.  
  138. badarg:                                 /* arg is negative or NaN */
  139.         movl    $1, C_SYM(errno)
  140.         flds    NaN
  141.         ret
  142.  
  143. special:
  144.         fldl    4(%esp)                 /* x = 0 or inf: just load x */
  145.         jmp     cleanup
  146.