I don't know how to calculate the square root of a floating point number,Can you help me?

The floating point number is reptrsented as this:
FLOATING POINT REPRESENTATION
1byte Exponent and 3bytes Signed Mantissa.

Thank you very much.

Use Newton-Raphson to successively improve an initial guess. So if the number you want to square root is x and the initial guess is n then the improved guess is given by:

n' = (n + x / n) / 2

Repeat this a few times (say 10) or until n and n' vary by less than your required tolerance.

An easy way to create the initial guess is to use input number but shift the exponent to the right by one bit (e.g. divide the exponent by two).

_________________
Andrew Jacobs
6502 & PIC Stuff - http://www.obelisk.me.uk/
Cross-Platform 6502/65C02/65816 Macro Assembler - http://www.obelisk.me.uk/dev65/
Open Source Projects - https://github.com/andrew-jacobs

Thanks.

But I'm a beginner,can anybody give me a source code?

First you need the software to add, multiply and divide floating point numbers. All standard and wellknown but not to a "beginner". Subroutines for this in for example all BASIC implementations. If you are a beginner and just want to learn, start with something more elementary. If this software really is needed then then it should be written by a "non-beginner" with a working understanding of math!

Anyway, look in the "Code" section to learn about multiplication and division with 6502! Good luck!

I want some source code to complete a statistics program.Can anybody give me a source code? Or some detaile reference.
I have already writed a floating point add/sub/mul/div sub-procedure. They can process single precision floating point numbers.

I might point out that floating-point (FP) is seldom necessary outside the world of calculators. FP adds a lot of complexity and overhead, and slows the computer down a lot. I used to think it was necessary for most non-integer numbers and certainly for things like trig and log functions, and so I wrote my own FP routines for a product I developed in the late 1980's using a 65c02.

Later when I was introduced to the Forth programming language and read Leo Brodie's excellent books "Starting Forth" and later "Thinking Forth", I found how fixed-point and scaled-integer arithmetic could do a lot of things that I had thought required FP. Still, I was skeptical. However through many years of experience, I have seen how fixed-point and scaled-integer arithmetic can take care of most needs. (FP can indeed be done in Forth, and many Forths come standard with it. But Forth is partly about power through simplicity, which is why use of FP is discouraged when unnecessary.) As an example, if you have a number you want to multiply by PI (3.14...), you can multiply by 355 and divide this intermediate result by 113. The end result is like multiplying by PI, accurate to 8 digits.

Let's say you want the square root of 3.564837E-9 (which is approximately 5.971E-5). You can take the square root of 35648370 which is 1E16 times as much, get the result 5971, and understand that the decimal point is moved over 8 places. (Half of the 16 in "1E16" is 8. The "half" comes from X^(1/2), which is square root.) This means that the .00005971 is represented as 5971. You, the programmer, determine where the decimal point will be, using the available precision to achieve the desired range in fixed-point.

Scaled integer is similar to fixed-point, but the "decimal point" may be in the middle of a digit, so to speak. When we use scaled-integer arithmetic for trig functions where 16-bit cells are standard precision, we make the 360-degree circle to be 65,536. So what happens when you add an angle and there's overflow? It just wraps around to the same angle, perfectly! Then the count for one radian is 65,536 divided by 2PI, or 10,430 (28BE in hex), which can be used in trig calculations like to take the sine of an angle.

In reality, you would normally do this in binary, not decimal. Even though the 6502 has a direct decimal arithmetic mode, computers (including the 6502) work more efficiently in binary. You will typically convert to and from decimal for human I/O, but do the internals in binary.

There is definitely a place for FP; but insistence upon it for most applications indicates inexperience.

For the acual calculation of the square root, I've used the r=(x+r^2)/2r iterative method, which converges rather quickly. I have not tried BitWise's method (above). There is even a simple way to get a square root with no multiplying or dividing, but it is _extremely_ slow!!

Once you have the basic functions defined in assembly like * / + and - , it is normally beneficial to do the rest in a higher-level language. Since most of the processor's time will be spent in the multiply and divide building-block routines anyway, you won't gain much by trying to do the more complex functions in assembly. IOW, it will be pretty rare that you'll find a square-root, sin, cos, tan, log, FFT, etc. function defined in assembly.

Um, your formula:

r = (x + r^2) / 2r

is the same as BitWise's:

n' = (n + x / n) / 2

since n and r are the same thing.

Anyway, for source code, there is:

1. The DP18 (a BCD FP package) square root routine in Apple Assembly Line, linked to in the Source Code Repository (the 11-1984 and 10-1984 issues discuss the square root routine).

2. The KIMATH square root routine, in Simple Microcomputers and Trainers -> MOS Technology KIM-1

3. The 6502 FP routines by Wozniak & Rankin in the Source Code Repository -> Floating Point Math

4. The EhBASIC square root function (SQR) starting at LAB_SQR, linked to in the Homebuilt Projects on The Web -> Homebuilt 6502 Software Projects

#1 & #2 use the Newton's method formula above (KIMATH calls it Heron's method, but it's the same thing).

#3 does not have a square root function, but you can calculate the square root with:

Square root of X = EXP (LOG (X) / 2)

EXP(number) is e^number. LOG is the natural log, not log base 10. In fact, any pair of N^number and log base N functions will work, as long as N is the same in both.

#4 uses a integer square root routine on the mantissa and divides the exponent by two, elegantly handling odd exponents. In other words:

SQRT (mantissa * base ^ exponent) = SQRT (mantissa) * base ^ (exponent / 2)

#4 is the fastest of the three techniques (a lot of cycles can be squeezed out of the actual routine, and the temporary RAM usage can be reduced as well), then Newton's method, and EXP(LOG) is the slowest. You have to pay attention to round-off errors with EXP(LOG), but it's the least amount of additional code when EXP and LOG are already available.

Which is what EhBASIC does. Just remember to adjust the mantissa calculation if the exponent is odd.

Lee.

http://www.cs.berkeley.edu/~wkahan/ieee754status/reciprt.pdf gives an algorithm with no divison. The method for obtaining the initial approximation is particularly good: no table is required (although you'll get away with fewer iterations if you use one), and the floating point input can be treated as an integer (no need to separate the exponent).

The algorithm is for rsqrt, but it's trivial to get sqrt from that.

The initial approximation without a table is good to a few bits. Some applications won't need any refinement at all, which makes rsqrt an integer shift and subtract, and sqrt the same plus a floating point multiply.

A couple of things I should also mention.

First, for source code of EXP(LOG) method, see the disassembly of SYM 1 BASIC on Lee Davison's (EhBASIC) web site.

Second, note that EhBASIC doesn't round the square root of the mantissa. The computed 24-bit mantissa is always <= the exact square root value. For example, when the mantissa is $900003, the integer that EhBASIC takes the square root of is $900003000000, which is closer to $C00002 ^ 2 (= $900003000004) than $C00001 ^ 2 (= $900001800001), but EhBASIC returns $C00001. The square root can be properly rounded by setting FAC1_r (the FAC1 rounding byte) to $80 when the remainder > the root.