I mentioned in another post that I was converting a mainframe Fortran

IV program to F77 (and F90).

I found that running the converted program with the original check

data set that went with this, crashes on trying to take the square

root of a negative variable (to get the standard deviation as a

measure for selecting later optimum adjustments).

I was utterly lost to discover why all the many computational results

checked exactly to some 6 places, EXCEPT for one standard error and

its significance (both unused later) and the program crashing on one

square root operation before terminating correctly with a final table.

Obviously the mainframe did not crash and the PC did on the same data!

After some days (!!) I inserted before the SQRT line equivalent to

this

std=SQRT(var)

with the line

IF (var.LT.0.d0) var = -var

Problem solved. (And I hope SQRT(0.0) gives a zero result!)

I deduce that mainframe runtime replaces the subject variable of a

square root operation with an absolute value and continues, wheras the

PC runtime informs the user of the error of trying to take the root of

a negative value and stops.

Aaagh!

I have seen this many times. Note that six places isn't enough,

if the result is negative in the seventh place.

When the standard deviation is zero, or very close to zero

it is possible for rounding in the variance to give a

negative result. Probably better to set it to zero

instead of the sqrt() of the absolute value, but it will be

very small in any case.

See http://www.yqcomputer.com/ #Definition_and_calculation

The popular standard deviation algorithm uses an identity that is

true mathematically, but not always numerically. That is, the

mean of the squares minus the square of the mean should be positive

but with some types of rounding, it isn't.

-- glen

On 2008-11-05 19:33:03 -0400, Terence < XXXX@XXXXX.COM > said:

Did you first try printing out the offending value? In floating point

so see small magnitudes. I seriously doubt that any mainframe or PC

would silently take the absolute value of a square root. Much more

likely that the rounding effects have changed.

Calculating a standard deviation that is small is an absolutely classical

demonstration of the computational instability of the moment based formulas

that filled statistics text books from the days of hand calculators. Much

better to do a two pass or provisional means algorithm. Trivial test examples

that are exact fits tend to demonstrate these problems.

Did you first try printing out the offending value? In floating point

so see small magnitudes. I seriously doubt that any mainframe or PC

would silently take the absolute value of a square root. Much more

likely that the rounding effects have changed.

Calculating a standard deviation that is small is an absolutely classical

demonstration of the computational instability of the moment based formulas

that filled statistics text books from the days of hand calculators. Much

better to do a two pass or provisional means algorithm. Trivial test examples

that are exact fits tend to demonstrate these problems.

(snip)

x87 might just leave the original value on the stack.

The design is that you would do the checking in software before

using the x87 floating point instructions, but when generated inline

the checking is often not done. I have seen it for sin(1e100),

I wouldn't be surprised for sqrt(-1.23).

But yes, different rounding is more likely. Note, for example

that S/360 (and successors) truncate the quotient on floating

point divide, where many others round. That could go either

way in the standard deviation case, though.

-- glen

Well, yes, there was a 100 by100 matrix inversion, and a triple matrix

multiplication, and values wandered from about 0.4 up to 10**11, so it

might not be considered to be a stable calculation.

But the input is real-world data.

I was limited to REAL*8.

I accept the good point that 0.0 is a better default than the SQRT of

the absolute value, since the sample standard deviations SHOULD be

positive.

Anyway it works now, but I wonder just how many more pre-IBM360

algorithms have real-time problmes?? My best guide to "practical

methods" is the UK 1961 National Physica Laboratory White Paper on

Modern Computing Methods - often meaning hand-cranked calculators!

multiplication, and values wandered from about 0.4 up to 10**11, so it

might not be considered to be a stable calculation.

But the input is real-world data.

I was limited to REAL*8.

I accept the good point that 0.0 is a better default than the SQRT of

the absolute value, since the sample standard deviations SHOULD be

positive.

Anyway it works now, but I wonder just how many more pre-IBM360

algorithms have real-time problmes?? My best guide to "practical

methods" is the UK 1961 National Physica Laboratory White Paper on

Modern Computing Methods - often meaning hand-cranked calculators!

> std=SQRT(var) >> >> with the line >> IF (var.LT.0.d0) var = -var> >> > Problem solved. (And I hope SQRT(0.0) gives a zero result!)> >> > I deduce that mainframe runtime replaces the subject variable of > > square root operation with an absolute value and continues, wheras the> > PC runtime informs the user of the error of trying to take the root of> > a negative value and stops.> >> > Aaagh!

IIRC Forran I and II did this rather than crash. Likewise, early IBM-

PC Fortran was designed so that dividing by zero was not a fatal

error, instead the computation would march on.

-- E

A recent web search took me to a page from Cambridge about their EDSAC

1 and 2 computers. Amazing what they accomplished with vacuum tubes

and mercury delay lines. I bet that most people in the USA know little

about UK computing history (or jet aviation, etc.).

Interesting to see that programming was done in machine or assembly

code, with paper tape or mag tape.

- E

Oh Yes!

I worked for a time on guidance for a large UK missile, at a London

office.

We used the Mercury computer for trajectory calculations.

First thing, was to walk into the computer with a box of double

triodes valves, unplug and replace any of the thousands with no

filaments visibly lit (while in standby mode).

Then exit, go to run mode, feed the next target paper tape and cross

fingers for a 20 minute no-fail run, collect the output tape, feed it

into the teleprinter and go off with the printout.

The "screen" was just an oscilloscope for the engineers.

A few years later later I was installing IBM 1620 (addition and

multiply was table look-up, division was a subtraction routine) paper

tape machines for those who couldn't afford a 704 (and what went with

it), and had come to the limit of plug-board capacity relay

calculators (though I did a lot of plug-board computations on water

control strategy for the Thames Water Board, using 1930-1950 data. Now

we do it with the Out-of-Kilter algorithm.

Remember Remingtons?.

I'd like to ask the same as Gordon, what happens when you try to insert a print

statement and look at the offending value. Is it a negative number that is

almost zero within roundoff error, or is it some totally ridiculous negative

number ?

You might even find that inserting a write statement masks the problem, in which

case something more serious could going wrong.

On 2008-11-06 00:37:50 -0400, Terence < XXXX@XXXXX.COM > said:

Some folks would take the view that the negative square root is

a hint that there are problems to be solved. Not shoved under the rug

or papered over or forgotten or otherwise ignored.

Why not take a look at Ake Bjorck's SIAM book on Least Squares. It

is much much much more modern than NPL's book on Modern Computing

Methods published by HMSO in the 1960s. A lot of its stuff is based

on fixed point with double length accumulation, the fi_2 mode.

Some folks would take the view that the negative square root is

a hint that there are problems to be solved. Not shoved under the rug

or papered over or forgotten or otherwise ignored.

Why not take a look at Ake Bjorck's SIAM book on Least Squares. It

is much much much more modern than NPL's book on Modern Computing

Methods published by HMSO in the 1960s. A lot of its stuff is based

on fixed point with double length accumulation, the fi_2 mode.

> std=SQRT(var) >> >> with the line >> IF (var.LT.0.d0) var = -var> >> > Problem solved. (And I hope SQRT(0.0) gives a zero result!)> >> > I deduce that mainframe runtime replaces the subject variable of > > square root operation with an absolute value and continues, wheras the> > PC runtime informs the user of the error of trying to take the root of> > a negative value and stops.> >> > Aaagh!

As the others have said, most probably a rounding problem.

I came across the very same problem when a user tried to

determine the standard deviation of a set of equal numbers.

The "smart" one-pass algorithm I had used came up with a tiny

negative variance. And so the program crashed in the next

step.

As for this not happening in earlier computers: they used

a wide variety of implementations to do floating-point

computations. And so did early PC compilers (remember

the ones that did not have a mathematical coprocessor?).

I think I prefer the current, much more homogeneous

situation. But I can imagine your frustration.

Regards,

Arjen

Square roots of negative numbers can also happen for other reasons. The

usual formula for computing airspeed from instrumentation data involves

subtracting the static pressure from the total pressure and then taking

a square root. I've seen plenty of naive versions of such code that

regularly blow up when the aircraft is sitting on the ramp doing

preflight tests. The coders don't check for negative values because the

static pressure "can't" be greater than the total pressure (that and

because lots of coders are just plain careless about checking much of

anything). While the "true" static pressure can't be greater than the

total pressure, the measured data are another matter. Measured data are

imperfect. If the airplane is sitting on the ramp, it is quite easy for

the measured static pressure to be greater than the measured total

pressure, be it from sensor noise, imperfect calibration, or what.

I might add that the best answer in this case is zero; not the square

root of the absolute value. But one can't conclude that for all cases of

square roots of negatives. It depends on the particular application. For

computing airspeed, zero turns out to be the best answer.

--

Richard Maine | Good judgment comes from experience;

email: last name at domain . net | experience comes from bad judgment.

domain: summertriangle | -- Mark Twain

Well, one might comment on that from two perspectives.

1. The definition. Airspeed is a magnitude. It doesn't have direction.

If you are going backwards, that's still a positive airspeed (with an

angle of attack of 180 degrees - or you could call it an angle of

sidelsliip of 180 degrees, depending on how you normalize those things).

2. The practical. It can, in fact, quite easily happen on the ground

that the airplane is going backwards (wind relative). All it takes is a

little tailwind with the airplane sitting still. But the usual

arrangement of sensors won't pick that up. The total pressure sensor is

only good for modest angles of attack and sideslip. There exist sensors

that can handle more extreme cases, but that's not the normal

arrangement. (Well, maybe on helicopters, but I don't "do" them.) With

normal sensors, you pretty much won't pick up anything useful when the

plane is going backwards. It will show as total pressure being equal to

static (i.e. zero airspeed) plus or minus some noise... and the noise

will get you into that negative square root territory.

Anyway, back to the somewhat vaguely Fortran-related part, the square

root of the absolute value would not be the right answer. Zero would be

about as well as you could do with normal instrumentats. If you happened

to have instruments suitable to the case, then you wouldn't have a

negative value.

--

Richard Maine | Good judgment comes from experience;

email: last name at domain . net | experience comes from bad judgment.

domain: summertriangle | -- Mark Twain

1. Simulink square root of a negative

2. Error testing... division by zero and square root of negative number

3. Formula to make Negative Values Positive & Positive Values Negative?

5. I want to change the background color of cells, if value is a square root.

6. Compute rms(root mean squared) value of a matrix

7. Square root program doesn't produce the right value

8. Roots other than square roots in Emacs calc

9. square root function in excel (roots greater than 12)

10. "Matching" column A's values in column B

11. Finding the Value (of text and numbers) between N/A's

12. Negative Adjusted R-square

14. Large-Scale Non-Negative Linear Least Squares

15. Least squares optimizing for equal positive & negative error sums