home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / bsd_srcs / lib / libc / tahoe / stdlib / atof.s < prev    next >
Encoding:
Text File  |  1991-04-12  |  11.4 KB  |  364 lines

  1. /*
  2.  * Copyright (c) 1988 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms, with or without
  6.  * modification, are permitted provided that the following conditions
  7.  * are met:
  8.  * 1. Redistributions of source code must retain the above copyright
  9.  *    notice, this list of conditions and the following disclaimer.
  10.  * 2. Redistributions in binary form must reproduce the above copyright
  11.  *    notice, this list of conditions and the following disclaimer in the
  12.  *    documentation and/or other materials provided with the distribution.
  13.  * 3. All advertising materials mentioning features or use of this software
  14.  *    must display the following acknowledgement:
  15.  *    This product includes software developed by the University of
  16.  *    California, Berkeley and its contributors.
  17.  * 4. Neither the name of the University nor the names of its contributors
  18.  *    may be used to endorse or promote products derived from this software
  19.  *    without specific prior written permission.
  20.  *
  21.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  22.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  25.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  26.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  27.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  28.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  29.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  30.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  31.  * SUCH DAMAGE.
  32.  */
  33.  
  34. #if defined(LIBC_SCCS) && !defined(lint)
  35.     .asciz "@(#)atof.s    5.3 (Berkeley) 6/1/90"
  36. #endif /* LIBC_SCCS and not lint */
  37.  
  38. #include "DEFS.h"
  39.  
  40. /*
  41.  *    atof: convert ascii to floating
  42.  *
  43.  *    C usage:
  44.  *
  45.  *        double atof (s)
  46.  *        char *s;
  47.  *
  48.  *    Register usage:
  49.  *
  50.  *        r0-1:    value being developed
  51.  *        r2:    first section: pointer to the next character
  52.  *            second section: binary exponent
  53.  *        r3:    flags
  54.  *        r4:    first section: the current character
  55.  *            second section: scratch
  56.  *        r5:    the decimal exponent
  57.  *        r6-7:    scratch
  58.  */
  59.     .set    msign,0        # mantissa has negative sign
  60.     .set    esign,1        # exponent has negative sign
  61.     .set    decpt,2        # decimal point encountered
  62.  
  63. ENTRY(atof, R6|R7)
  64. /*
  65.  *    Initialization
  66.  */
  67.     clrl    r3        # All flags start out false
  68.     movl    4(fp),r2    # Address the first character
  69.     clrl    r5        # Clear starting exponent
  70. /*
  71.  *    Skip leading white space
  72.  */
  73. sk0:    movzbl    (r2),r4        # Fetch the next (first) character
  74.     incl    r2
  75.     cmpb    $' ,r4        # Is it blank?
  76.     beql    sk0        #   ...yes
  77.     cmpb    r4,$8        # 8 is lowest of white-space group
  78.     blss    sk1        # Jump if char too low to be white space
  79.     cmpb    r4,$13        # 13 is highest of white-space group
  80.     bleq    sk0        # Jump if character is white space
  81. sk1:
  82. /*
  83.  *    Check for a sign
  84.  */
  85.     cmpb    $'+,r4        # Positive sign?
  86.     beql    cs1        #   ... yes
  87.     cmpb    $'-,r4        # Negative sign?
  88.     bneq    cs2        #   ... no
  89.     orb2    $1<msign,r3    # Indicate a negative mantissa
  90. cs1:    movzbl    (r2),r4        # Skip the character
  91.     incl    r2
  92. cs2:
  93. /*
  94.  *    Accumulate digits, keeping track of the exponent
  95.  */
  96.     clrl    r1
  97.     clrl    r0        # Clear the accumulator
  98. ad0:    cmpb    r4,$'0        # Do we have a digit?
  99.     blss    ad4        #   ... no, too small
  100.     cmpb    r4,$'9
  101.     bgtr    ad4        #   ... no, too large
  102. /*
  103.  *    We got a digit.  Accumulate it
  104.  */
  105.     cmpl    r0,$214748364    # Would this digit cause overflow?
  106.     bgeq    ad1        #   ... yes
  107. /*
  108.  *    Multiply (r0,r1) by 10.  This is done by developing
  109.  *    (r0,r1)*2 in (r6,r7), shifting (r0,r1) left three bits,
  110.  *    and adding the two quadwords.
  111.  */
  112.     shlq    $1,r0,r6    # (r6,r7)=(r0,r1)*2
  113.     shlq    $3,r0,r0    # (r0,r1)=(r0,r1)*8
  114.     addl2    r7,r1        # Add low halves
  115.     adwc    r6,r0        # Add high halves
  116. /*
  117.  *    Add in the digit
  118.  */
  119.     subl2    $'0,r4        # Get the digit value
  120.     addl2    r4,r1        # Add it into the accumulator
  121.     adwc    $0,r0        # Possible carry into high half
  122.     brb    ad2        # Join common code
  123. /*
  124.  *    Here when the digit won't fit in the accumulator
  125.  */
  126. ad1:    incl    r5        # Ignore the digit, bump exponent
  127. /*
  128.  *    If we have seen a decimal point, decrease the exponent by 1
  129.  */
  130. ad2:    bbc    $decpt,r3,ad3    # Jump if decimal point not seen
  131.     decl    r5        # Decrease exponent
  132. ad3:
  133. /*
  134.  *    Fetch the next character, back for more
  135.  */
  136.     movzbl    (r2),r4        # Fetch
  137.     incl    r2
  138.     brb    ad0        # Try again
  139. /*
  140.  *    Not a digit.  Could it be a decimal point?
  141.  */
  142. ad4:    cmpb    r4,$'.        # If it's not a decimal point, either it's
  143.     bneq    ad5        #   the end of the number or the start of
  144.                 #   the exponent.
  145.     bbs    $decpt,r3,ad5
  146.     orb2    $1<decpt,r3    # If it IS a decimal point, we record that
  147.     brb    ad3        #   we've seen one, and keep collecting
  148.                 #   digits if it is the first one.
  149.             
  150. /*
  151.  *    Check for an exponent
  152.  */
  153. ad5:    clrl    r6        # Initialize the exponent accumulator
  154.  
  155.     cmpb    r4,$'e        # We allow both lower case e
  156.     beql    ex1        #   ... and ...
  157.     cmpb    r4,$'E        #   upper-case E
  158.     bneq    ex7
  159. /*
  160.  *    Does the exponent have a sign?
  161.  */
  162. ex1:    movzbl    (r2),r4        # Get next character
  163.     incl    r2
  164.     cmpb    r4,$'+        # Positive sign?
  165.     beql    ex2        #   ... yes ...
  166.     cmpb    r4,$'-        # Negative sign?
  167.     bneq    ex3        #   ... no ...
  168.     orb2    $1<esign,r3    # Indicate exponent is negative
  169. ex2:    movzbl    (r2),r4        # Grab the next character
  170.     incl    r2
  171. /*
  172.  *    Accumulate exponent digits in r6
  173.  */
  174. ex3:    cmpb    r4,$'0        # A digit is within the range
  175.     blss    ex4        # '0' through
  176.     cmpb    r4,$'9        # '9',
  177.     bgtr    ex4        # inclusive.
  178.     cmpl    r6,$214748364    # Exponent outrageously large already?
  179.     bgeq    ex2        #   ... yes
  180.     moval    (r6)[r6],r6    # r6 *= 5
  181.     movaw    -'0(r4)[r6],r6    # r6 = r6 * 2 + r4 - '0'
  182.     brb    ex2        # Go 'round again
  183. ex4:
  184. /*
  185.  *    Now get the final exponent and force it within a reasonable
  186.  *    range so our scaling loops don't take forever for values
  187.  *    that will ultimately cause overflow or underflow anyway.
  188.  *    A tight check on over/underflow will be done by ldexp.
  189.  */
  190.     bbc    $esign,r3,ex5    # Jump if exponent not negative
  191.     mnegl    r6,r6        # If sign, negate exponent
  192. ex5:    addl2    r6,r5        # Add given exponent to calculated exponent
  193.     cmpl    r5,$-100    # Absurdly small?
  194.     bgtr    ex6        #   ... no
  195.     movl    $-100,r5    #   ... yes, force within limit
  196. ex6:    cmpl    r5,$100        # Absurdly large?
  197.     blss    ex7        #   ... no
  198.     movl    $100,r5        #   ... yes, force within bounds
  199. ex7:
  200. /*
  201.  *    Our number has now been reduced to a mantissa and an exponent.
  202.  *    The mantissa is a 63-bit positive binary integer in r0,r1,
  203.  *    and the exponent is a signed power of 10 in r5.  The msign
  204.  *    bit in r3 will be on if the mantissa should ultimately be
  205.  *    considered negative.
  206.  *
  207.  *    We now have to convert it to a standard format floating point
  208.  *    number.  This will be done by accumulating a binary exponent
  209.  *    in r2, as we progressively get r5 closer to zero.
  210.  *
  211.  *    Don't bother scaling if the mantissa is zero
  212.  */
  213.     tstl    r1
  214.     bneq    1f
  215.     tstl    r0        # Mantissa zero?
  216.     jeql    exit        #   ... yes
  217.  
  218. 1:    clrl    r2        # Initialize binary exponent
  219.     tstl    r5        # Which way to scale?
  220.     bleq    sd0        # Scale down if decimal exponent <= 0
  221. /*
  222.  *    Scale up by "multiplying" r0,r1 by 10 as many times as necessary,
  223.  *    as follows:
  224.  *
  225.  *    Step 1: Shift r0,r1 right as necessary to ensure that no
  226.  *    overflow can occur when multiplying.
  227.  */
  228. su0:    cmpl    r0,$429496729    # Compare high word to (2**31)/5
  229.     blss    su1        # Jump out if guaranteed safe
  230.     shrq    $1,r0,r0    # Else shift right one bit
  231.     incl    r2        #    bump exponent to compensate
  232.     brb    su0        #    and go back to test again.
  233. /*
  234.  *    Step 2: Multiply r0,r1 by 5, by appropriate shifting and
  235.  *    double-precision addition
  236.  */
  237. su1:    shlq    $2,r0,r6    # (r6,r7) := (r0,r1) * 4
  238.     addl2    r7,r1        # Add low-order halves
  239.     adwc    r6,r0        #   and high-order halves
  240. /*
  241.  *    Step 3: Increment the binary exponent to take care of the final
  242.  *    factor of 2, and go back if we still need to scale more.
  243.  */
  244.     incl    r2        # Increment the exponent
  245.     decl    r5        # ...sobgtr r5,su0
  246.     bgtr    su0        #    and back for more (maybe)
  247.  
  248.     brb    cm0        # Merge to build final value
  249.  
  250. /*
  251.  *    Scale down.  We must "divide" r0,r1 by 10 as many times
  252.  *    as needed, as follows:
  253.  *
  254.  *    Step 0: Right now, the condition codes reflect the state
  255.  *    of r5.  If it's zero, we are done.
  256.  */
  257. sd0:    beql    cm0        # If finished, build final number
  258. /*
  259.  *    Step 1: Shift r0,r1 left until the high-order bit (not counting
  260.  *    the sign bit) is nonzero, so that the division will preserve
  261.  *    as much precision as possible.
  262.  */
  263.     tstl    r0        # Is the entire high-order half zero?
  264.     bneq    sd2        #   ...no, go shift one bit at a time
  265.     shlq    $30,r0,r0    #   ...yes, shift left 30,
  266.     subl2    $30,r2        #   decrement the exponent to compensate,
  267.                 #   and now it's known to be safe to shift
  268.                 #   at least once more.
  269. sd1:    shlq    $1,r0,r0    # Shift (r0,r1) left one, and
  270.     decl    r2        #   decrement the exponent to compensate
  271. sd2:    bbc    $30,r0,sd1    # If the high-order bit is off, go shift
  272. /*
  273.  *    Step 2: Divide the high-order part of (r0,r1) by 5,
  274.  *    giving a quotient in r1 and a remainder in r7.
  275.  */
  276. sd3:    movl    r0,r7        # Copy the high-order part
  277.     clrl    r6        # Zero-extend to 64 bits
  278.     ediv    $5,r6,r0,r6    # Divide (cannot overflow)
  279. /*
  280.  *    Step 3: Divide the low-order part of (r0,r1) by 5,
  281.  *    using the remainder from step 2 for rounding.
  282.  *    Note that the result of this computation is unsigned,
  283.  *    so we have to allow for the fact that an ordinary division
  284.  *    by 5 could overflow.  We make allowance by dividing by 10,
  285.  *    multiplying the quotient by 2, and using the remainder
  286.  *    to adjust the modified quotient.
  287.  */
  288.     addl3    $2,r1,r7    # Dividend is low part of (r0,r1) plus
  289.     adwc    $0,r6        #  2 for rounding plus
  290.                 #  (2**32) * previous remainder
  291.     ediv    $10,r6,r1,r7    # r1 := quotient, r7 := remainder.
  292.     addl2    r1,r1        # Make r1 result of dividing by 5
  293.     cmpl    r7,$5        # If remainder is 5 or greater,
  294.     blss    sd4        #   increment the adjustted quotient.
  295.     incl    r1
  296. /*
  297.  *    Step 4: Increment the decimal exponent, decrement the binary
  298.  *    exponent (to make the division by 5 into a division by 10),
  299.  *    and back for another iteration.
  300.  */
  301. sd4:    decl    r2        # Binary exponent
  302.     aoblss    $0,r5,sd2
  303. /*
  304.  *    We now have the following:
  305.  *
  306.  *    r0:    high-order half of a 64-bit integer
  307.  *    r1:    load-order half of the same 64-bit integer
  308.  *    r2:    a binary exponent
  309.  *
  310.  *    Our final result is the integer represented by (r0,r1)
  311.  *    multiplied by 2 to the power contained in r2.
  312.  *    We will transform (r0,r1) into a floating-point value,
  313.  *    set the sign appropriately, and let ldexp do the
  314.  *    rest of the work.
  315.  *
  316.  *    Step 1: if the high-order bit (excluding the sign) of
  317.  *    the high-order half (r0) is 1, then we have 63 bits of
  318.  *    fraction, too many to convert easily.  However, we also
  319.  *    know we won't need them all, so we will just throw the
  320.  *    low-order bit away (and adjust the exponent appropriately).
  321.  */
  322. cm0:    bbc    $30,r0,cm1    # jump if no adjustment needed
  323.     shrq    $1,r0,r0    # lose the low-order bit
  324.     incl    r2        # increase the exponent to compensate
  325. /*
  326.  *    Step 2: split the 62-bit number in (r0,r1) into two
  327.  *    31-bit positive quantities
  328.  */
  329. cm1:    shlq    $1,r0,r0    # put the high-order bits in r0
  330.                 #   and a 0 in the bottom of r1
  331.     shrl    $1,r1,r1    # right-justify the bits in r1
  332.                 #   moving 0 into the sign bit.
  333. /*
  334.  *    Step 3: convert both halves to floating point
  335.  */
  336.     cvld    r1
  337.     std    r6        # low-order part in r6-r7
  338.     cvld    r0
  339.     std    r0        # high-order part in r0-r1
  340. /*
  341.  *    Step 4: multiply the high order part by 2**31 and combine them
  342.  */
  343.     ldd    two31
  344.     muld    r0        # multiply
  345.     addd    r6        # combine
  346. /*
  347.  *    Step 5: if appropriate, negate the floating value
  348.  */
  349.     bbc    $msign,r3,cm2    # Jump if mantissa not signed
  350.     negd            # If negative, make it so
  351. /*
  352.  *    Step 6: call ldexp to complete the job
  353.  */
  354. cm2:    pushl    r2        # Put exponent in parameter list
  355.     pushd            #    and also mantissa
  356.     calls    $3,_ldexp    # go combine them
  357.  
  358. exit:
  359.     ret
  360.  
  361.     .align    2
  362. two31:    .long    0x50000000    # (=2147483648) 2 ** 31 in floating-point
  363.     .long    0        # so atof doesn't have to convert it
  364.