Re: Calculating SQRT in ADA

Re: Calculating SQRT in ADA

[Chris Morgan]

"cmcrae" <cmcrae@orca.st.usm.edu> writes:

>   I am trying to calculate the SQRT in ADA.  Is there a function that I
> can access?

Look for sqrt in the Ada Language Reference Manual. There, that wasn't
too hard was it?
-- 
Chris Morgan <mihalis at ix.netcom.com                http://mihalis.net
     "I don't have time to read a 1200 page book. 
      I am afraid to even let one in my house." 
                                     - Philip Greenspun 

Calculating SQRT in ADA

[cmcrae]

  I am trying to calculate the SQRT in ADA.  Is there a function that I
can access?
--
Posted via Talkway - http://www.talkway.com
Surf Usenet at home, on the road, and by email -- always at Talkway.

Re: Calculating SQRT in ADA

[Marin David Condic]

Chris Morgan wrote:
> 
> "cmcrae" <cmcrae@orca.st.usm.edu> writes:
> 
> >   I am trying to calculate the SQRT in ADA.  Is there a function that I
> > can access?
> 
> Look for sqrt in the Ada Language Reference Manual. There, that wasn't
> too hard was it?

I believe the original poster is using Ada83 - although he has not been
explicit about this. In which case, there may be vendor supplied
routines, but not standard ones to be looked up in the RM.

If "cmcrae" were to provide the language version (83/95), the compiler
and OS, we might be able to provide some help. Without that info, the
best I can suggest is to a) check with the Public Ada Library and maybe
http://www.adahome.com/ for sources of preexisting code or b) look up an
algorithm and implement one from scratch. Its not that big a job. (I
used to have a Sqrt which I think used Newton-Raphson - I just don't
remember where I put it! ;-)

MDC

-- 
Marin David Condic
Real Time & Embedded Systems, Propulsion Systems Analysis
United Technologies, Pratt & Whitney, Large Military Engines
M/S 731-95, P.O.B. 109600, West Palm Beach, FL, 33410-9600
***To reply, remove "bogon" from the domain name.***

Re: Calculating SQRT in ADA

[John Herro]

Here's a generic square root routine in Ada 83, taken from my AdaTutor program:

generic
  type Dummy is digits <>;
function Sqrt(X :in Dummy) return Dummy;
function Sqrt(X :in Dummy) return Dummy is
  Guess : Dummy := X;
begin
  if X < 0.0 then
    raise Constraint_Error;
  end if;
  while X /= 0.0 and then
       abs(Guess*Guess/X - 1.0) > 3.0*Dummy'Epsilon loop
    Guess := (X/Guess + Guess) * 0.5;
  end loop;
  return Guess;
end Sqrt;

I tested this routine with Vax Ada 83, which has types Float, Long_Float, and
Long_Long_Float.  That last type gives at least 33 decimal digits of precision.
 I instantiated this function, as well as Text_IO.Float_IO, for all three
types.  When I tested with the three types, all displayed digits of the answers
were correct.  I hereby place this long and complex function :-) into the
public domain.  Enjoy.
- John Herro
You can download a shareware AdaTutor program at
http://members.aol.com/AdaTutor

Re: Calculating SQRT in ADA

[Niklas Holsti]

bob wrote:
 
   [snip]

> In Ada83 or Ada95, if the target is real, possibly the following:
> x : real;
> [snip]
> x = x**0.5;

Nope, for the standard "**" the right operand must be integer
(Integer'Base, actually).

> All of the Ada85 comilers I have used supply a "math" package
> containing sqrt also.

That agrees with my experience and would be the first place to
look.

Niklas Holsti
Working at but not speaking for Space Systems Finland Ltd.

Re: Calculating SQRT in ADA

[Howard W. LUDWIG]

robert_dewar@my-dejanews.com wrote:

> I disagree, it is generally confusing to talk about
> obsolete languages unless necessary. If people ask about
> Ada, it is reasonable to assume they mean Ada 95. If people
> really want to ask about Ada 83 (or Ada 80, or whatever),
> then they need to qualify the questions appropriately.

I would like the situation to be as you describe.  If someone asks specific
questions about tagged records, modular integers, etc. using appropriate
terminology, then it makes good sense to assume Ada 95, as Ada 83 did not
deal with such.  However, we seem to get many newbies who have just been
hired to maintain Ada 83 code, and many times it is apparent that either
the person is using Ada 83 (because the answer to the question is obvious
to the casual observer in Ada 95) or does not know where to find
documentation.  (Of course, we do encounter the scenario of people not
wanting to bother to look up documentation, as most of us know in Robert's
pulling his secret documentation out of his drawer to answer GNAT
questions.)

I suppose my main point is really that we should be more careful in how we
address questions of potential newbies (who may not even know there is an
Ada 95 in contrast to Ada 83).  Some of us have this tendency to jump on
such people, belittling their honest ignorance, and showing off our
superior knowledge.  Such responses make the true newbie a bit wary of Ada
based on an apparent elitist culture in the Ada community.  We don't need
that kind of negative marketing.

HWL

Re: Calculating SQRT in ADA

[Hans Marqvardsen]

John Herro wrote:
> 
> Here's a generic square root routine in Ada 83, taken from my AdaTutor program:
> 
> generic
>   type Dummy is digits <>;
> function Sqrt(X :in Dummy) return Dummy;
> function Sqrt(X :in Dummy) return Dummy is
>   Guess : Dummy := X;

    -- insert, just in case the compiler fail to recognize 
    -- a constant
    Tolerance : constant Dummy := 3.0*Dummy'epsilon*X;

> begin
>   if X < 0.0 then
>     raise Constraint_Error;
>   end if;
>   while X /= 0.0 and then
>        abs(Guess*Guess/X - 1.0) > 3.0*Dummy'Epsilon loop

    -- replace the above with simpler and faster:
    while (Guess*Guess - X) > tolerance loop


>     Guess := (X/Guess + Guess) * 0.5;
>   end loop;
>   return Guess;
> end Sqrt;
> 
> I tested this routine with Vax Ada 83, which has types Float, Long_Float, and
> Long_Long_Float.  That last type gives at least 33 decimal digits of precision.
>  I instantiated this function, as well as Text_IO.Float_IO, for all three
> types.  When I tested with the three types, all displayed digits of the answers
> were correct.  I hereby place this long and complex function :-) into the
> public domain.  Enjoy.
> - John Herro
> You can download a shareware AdaTutor program at
> http://members.aol.com/AdaTutor

Bravo, amazing that it is this simple to find a Square root
efficiently!

A few questions out of general interest:

Where doest the multiplier 3.0 come from ?

Cant you make the routine even simpler and faster, 
multiplying with X, after having established that X is 
not-negative? 
(At the same time avoiding a division 
and the need to test for X=0.0?)

Do you really need to use abs ? 
(Isn't always Guess >=  'true sqrt' anyway?)

Sincerely,
Hans Marqvardsen

Re: Calculating SQRT in ADA

[Hans Marqvardsen]

John Herro wrote:
> 
> Robert Dewar wrote:
> >>   while ...
> >>        abs(Guess*Guess/X - 1.0) > 3.0*Dummy'Epsilon loop
> > Why the epsilon test, this iteration will terminate to
> > exact equality in any reasonable floating-point model ...
> 
> I don't know if the Open Ada 83 compiler has a "reasonable" floating-point
> model, but I do know that if I test for exact equality, instantiate with type
> Float, and compute Sqrt(2.0) with that compiler, the program goes into an
> infinite loop.  I think that in general, you can't guarantee that Guess*Guess
> will ever exactly equal X.
> 

Clearly, in any floating-point model the square root of X can only
be represented  exactly  for very few values of X.

However any "reasonable" model may be expected to contain two 
neighbour numbers G and H such that

   G*G <= X <= H*H

Then one may expect convergence to either
G, or H, or to the alternating sequence G,H,G,H,...

Hence you need the epsilon test after all, 
but probably the constant 3.0 can safely be reduced.

Sincerely, 
Hans Marqvardsen

(In my previous posting, please disregard the idea of excluding abs)

Re: Calculating SQRT in ADA

[bob]

Several possibilities.

If you are using gnat Ada95, the following is a possibility:

[snip]
with text_io;
with Ada.Numerics.Generic_Elementary_Functions;
use text_io;

package body vectors is

  package math is new Ada.Numerics.Generic_Elementary_Functions(real);
  use math;
  package fio is new float_io(real);
  use fio;
[snip]
  -- Evaluate the magnitude of a vector 
  function mag(v : vector) return real is
  begin
    return sqrt(v(1)**2 + v(2)**2 + v(3)**2);
  end mag;
[snip]
end vectors;

In Ada83 or Ada95, if the target is real, possibly the following:
x : real;
[snip]
x = x**0.5;
[snip]

All of the Ada85 comilers I have used supply a "math" package containing
sqrt also.

cheers....bob

cmcrae <cmcrae@orca.st.usm.edu> wrote in article
<auWJ2.6969$2G6.3578@c01read02.service.talkway.com>...
>   I am trying to calculate the SQRT in ADA.  Is there a function that I
> can access?
> --
> Posted via Talkway - http://www.talkway.com
> Surf Usenet at home, on the road, and by email -- always at Talkway.
> 
> 

Re: Ada 83 - Sometimes still chosen

[Richard D Riehle]

In article <36F913E0.75F51763@lmco.com>,
	"Howard W. LUDWIG" <howard.w.ludwig@lmco.com> wrote:

>We need to be a little more careful in addressing questions of new people.
>The answer Chris gave _assumes_ Ada 95.  Many people, for whatever reason,
>still use Ada 83.  

 Howard's point needs to be emphasized.  We still get inquiries for
 Ada 83 training.  Not everyone is making a transition to Ada 95. A
 few compiler publishers are still selling a remarkable number of
 licenses for Ada 83 compilers.  Why is this?  

 One reason is that Ada 95 is not ported to all the environments used
 in real applications.  One platform that comes to mind is the HP
 real time computer for which only an Alsys Ada 83 compiler is available.
 Another is the MIL-STD 1750A, still widely used in space applications
 and military aircraft.  There are, of course, many others.  For example,
 though it may be difficult to believe, there are still designers who 
 specify "older" 8086 and 80286 processors for certain applications.

 In a few situations, it is a matter of "authorization" to use the new
 Ada standard.  This is a reflection of the wonderful world of bureaucracy
 where a contract specifies MIL-STD 1815A and the program manager refuses
 to allow any deviation from the contract.  

 There are some sensible reasons for sticking to Ada 83.  Here in Palo
 Alto there is a company developing long-lived space applications. Ada
 83 has been used for satellites that will survive through the first
 quarter of the next century.  The code has been written in Ada 83 and
 Ada 83 will continue to be used for maintaining that code for the life
 span of the satellite.  Long-lived software is one of Ada's reasons for
 existence.  Robert Dewar often stresses that Ada is intended to me more
 readable than writeable.  In the year 2030, when some programmer is 
 trying to code a program for uplink to a satellite to compensate for
 a component destroyed by electrostatic discharge, that programmer will
 need to be able to read the code written in 1997 by some other programmer
 who is either dead or retired.  It will still be Ada 83 code.

 Richard Riehle
 richard@adaworks.com
 http://www.adaworks.com

Re: Ada 83 - Sometimes still chosen

[Richard D Riehle]

In article <7dbvd6$56q$1@nnrp1.dejanews.com>,
	robert_dewar@my-dejanews.com wrote:

in response to my observation about the continuation of Ada 83,

>And almost certainly it will still be Ada 95 code, seeing
>as Ada 95 is almost exactly upwards compatible with Ada 83

 On some projects some programming managers feel it is appropriate 
 to freeze even the version of the compiler used to develop some 
 long-lived software product.  That is, once the satellite is 
 deployed, all updates will be made using exactly the same language
 and version of the compiler used for the original development.  

 Richard Riehle
 richard@adaworks.com
 http://www.adaworks.com
 

Re: Ada 83 - Sometimes still chosen

[Richard D Riehle]

In article <7dbvl2$5bl$1@nnrp1.dejanews.com>,
	robert_dewar@my-dejanews.com wrote:

>Sure, we know that some projects get stuck using obsolete
>hardware, software, operating systems, compilers etc, and
>that is understandable (and not necessarily an
>inappropriate decision).

 I am uneasy with the selection of the word "obsolete."  A
 technology is obsolete if it is not the best choice for the
 solution to a particular problem.  Sometimes a newer CPU
 architecture might offer some improvement in performance,
 but the price/performance ratio might make the obsolete 
 processor a better option.

 I return to the example of the 1750A.  This processor, by all
 standards for earth-bound computing is obsolete.  For a variety
 of technical reasons, it is still the optimal choice for many 
 space applications.  Sometimes the economics of an 8086 can make 
 it the correct choice.  In many circumstances we find some of 
 our clients still choosing "obsolete eight-bit processors because
 they are excellent at small,single threaded tasks -- and one 
 can populate a circuit board with them for incredible throughput.

 Ada 83 is the only Ada compiler available for many of these "obsolete"
 processors.  But they are not obsolete.  It is not likely anyone will
 bother to write a production quality Ada 95 compiler for an 8086, but
 some designers will continue to choose that processor because it offers
 better economics for the application being designed. 

>But in comp.lang.ada, I think we do not want to be
>emphasizing this kind of stuck-in-the-dark-ages view.

 I know you like to tease us with hyperbole from time to time
 Robert.  My own antiquity is such that I do long for those
 simpler times when I could repair my own car, map the entire
 memory of my program, and see the actual printing on a printed
 circuit board.  Sigh.  I suppose those were truly dark ages.

>Ada 95 is a modern, recently designed language, and it
>is important to always have the message that Ada is alive
>and well (not just that ancient projects forced to use it
>a long time ago are still alive and well).

 I agree with your observations about Ada 95.  I am also delighted
 with the progress we have seen you and your colleagues at ACT
 make with it.  I wonder if those of us, living in the dark ages,
 will see a time when Ada 95 will be available for some of those
 platforms you consider obsolete.  I suspect not.  It is a matter
 of economics.  Consequently, when an older architecture is chosen,
 for perfectly good technical reasons, it would be appropriate to
 use a compiler that already exists.  Often, this will be Ada 83.

>Sure, if someone from one of these projects needs to know
>something, we can probably provide the answer here to
>questions about Ada 83, provided the question specifies
>this (for that matter, this is probably the only place you
>could still get Jovial questions answered :-)

 I realize you were not writing off those who have a genuine need
 for Ada 83 information.  Also, there are few people who subscribe
 to this channel who are more able to answer such questions.  I 
 simply got a little over-exercised from your choice of the word
 "obsolete."  Perhaps I took it to personally.  :-)

 Richard Riehle
 richard@adaworks.com
 http://www.adaworks.com 

Re: Ada 83 - Sometimes still chosen

[Larry Kilgallen]

In article <7dc9ld$2u6@dfw-ixnews9.ix.netcom.com>, Richard D Riehle <laoXhai@ix.netcom.com> writes:
> In article <7dbvl2$5bl$1@nnrp1.dejanews.com>,
> 	robert_dewar@my-dejanews.com wrote:
> 
>>Sure, we know that some projects get stuck using obsolete
>>hardware, software, operating systems, compilers etc, and
>>that is understandable (and not necessarily an
>>inappropriate decision).
> 
>  I am uneasy with the selection of the word "obsolete."  A
>  technology is obsolete if it is not the best choice for the
>  solution to a particular problem.  Sometimes a newer CPU
>  architecture might offer some improvement in performance,
>  but the price/performance ratio might make the obsolete 
>  processor a better option.

Price/performance is irrelevant if the price has already been
paid and the performance is adequate. The quality of the new
version of my bagel-counting program is going to be better
than the previous version mainly due to increased sophistication
of the bagel-counting algorithm, not due to compiler attributes.

There is a _great_ deal of value to sticking with the same tools
if repeatable results are important.

Larry Kilgallen

Re: Ada 83 - Sometimes still chosen

[Marin David Condic]

Richard D Riehle wrote:
> 
> In article <7dbvd6$56q$1@nnrp1.dejanews.com>,
>         robert_dewar@my-dejanews.com wrote:
> 
> in response to my observation about the continuation of Ada 83,
> 
> >And almost certainly it will still be Ada 95 code, seeing
> >as Ada 95 is almost exactly upwards compatible with Ada 83
> 
>  On some projects some programming managers feel it is appropriate
>  to freeze even the version of the compiler used to develop some
>  long-lived software product.  That is, once the satellite is
>  deployed, all updates will be made using exactly the same language
>  and version of the compiler used for the original development.
> 
I'd say this is true of most mission critical applications which require
any significant amount of verification. And the freezing point will come
far sooner than system deployment. The usual question goes like this:
"If you change the compiler version, host platform, etc. and feed it the
source code you gave me last week, can you guarantee that I'll get
exactly the same bits out this time as I did then?" The answer is
obviously "No." (It gets worse that that when you consider that it is
possible to have the same compiler version, OS version, etc. and run on
two different days with different loads on the system and because of
dynamic memory situations, the optimizer might give you different object
code. This actually happened to us once.)

Once you go to the effort of expensive verification testing, you need to
be able to reproduce the image precisely without disturbing anything you
did not intend to modify.

-- 
Marin David Condic
Real Time & Embedded Systems, Propulsion Systems Analysis
United Technologies, Pratt & Whitney, Large Military Engines
M/S 731-95, P.O.B. 109600, West Palm Beach, FL, 33410-9600
***To reply, remove "bogon" from the domain name.***

Re: Calculating SQRT in ADA

[David C. Hoos, Sr.]

Hans Marqvardsen wrote in message <36F8D9B7.1F8D@ddre.dk>...
>John Herro wrote:
>>
>> Robert Dewar wrote:
>> >>   while ...
>> >>        abs(Guess*Guess/X - 1.0) > 3.0*Dummy'Epsilon loop
>> > Why the epsilon test, this iteration will terminate to
>> > exact equality in any reasonable floating-point model ...
>>
>> I don't know if the Open Ada 83 compiler has a "reasonable" floating-point
>> model, but I do know that if I test for exact equality, instantiate with type
>> Float, and compute Sqrt(2.0) with that compiler, the program goes into an
>> infinite loop.  I think that in general, you can't guarantee that Guess*Guess
>> will ever exactly equal X.
>>
>
>Clearly, in any floating-point model the square root of X can only
>be represented  exactly  for very few values of X.
>
>However any "reasonable" model may be expected to contain two
>neighbour numbers G and H such that
>
>   G*G <= X <= H*H
>
>Then one may expect convergence to either
>G, or H, or to the alternating sequence G,H,G,H,...
>
>Hence you need the epsilon test after all,
>but probably the constant 3.0 can safely be reduced.

I take issue with all of the previous answers on at least one point in each.

First, the test should not be for Guess * Guess / X = 1.0, because this is far too expensive of CPU cycles and unnecessary.  The
correct test, as Hans points out is for equality of successive estimates, or alternating.  This permits the solution to be obtained
with one division, one multiplication, one addition, and one or two floating point comparisons per iteration.

So, no epsilon test _per se_ is required.

As Robert correctly points out the choice of the initial estimate is important.  One solution to this is to examine the exponent of
the value of which the root is desired (argument), and start with an estimate of the root having an exponent half that of the
argument.

For example, of the machine radix is 2, then if the argument has an odd exponent, one can denormalize the number by increasing the
exponent by one (so it can now be evenly divisible by 2), and halving the mantissa.
This puts the mantissa of the denormalized argument on the range 0.25 .. 1.0, so its square root will be the correct value for the
mantissa of the square root, and the exponent of the square root is simply half of the (possibly incremented) exponent of the
argument.

This approach, however, is machine-dependent, but we actually found this approach faster than any other, on a machine whose
OS-provided math library (which was used by all language compilers) provided abysmal performance.

The square root of the mantissa can be closely estimated using a Pade' approximation, such that generally no more than two
Newton-Raphson iterations are required for convergence.

I became very concerned about this sort of thing about thirty years ago when CPU time cost me 40 cents per second on a GE 635.

Here is an implementation which is machine-independent, but suffers from the initial estimate problem.  It will, however,
unconditionally converge.


   generic
      type Real is digits <>;
   function Sqrt (X : real) return Real;

   function Sqrt (X : real) return real is
      Previous_Estimate : Real := X;
      New_Estimate : Real;
      Increasing : Boolean;
      Reversed : Boolean := False;
      Alternating : Boolean := False;
   begin
      if X < 0.0 then
         raise Constraint_Error;
      end if;
      if X = 0.0 then
         return 0.0;
      end if;
      New_Estimate :=
        0.5 * (X / Previous_Estimate + Previous_Estimate);
      Increasing := New_Estimate > Previous_Estimate;
      while  not alternating and then
        New_Estimate /= Previous_estimate loop
         if not Reversed then
            Reversed :=
              Increasing /= (New_Estimate > Previous_Estimate);
         else
            Alternating :=
              Increasing = (New_Estimate > Previous_Estimate);
         end if;
         Previous_Estimate := New_Estimate;
         New_Estimate :=
           0.5 * (X / Previous_Estimate + Previous_Estimate);
      end loop;
      return New_Estimate;
   end Sqrt;

Re: Calculating SQRT in ADA

[Larry Kilgallen]

In article <36F9AA62.24118202@lmco.com>, "Howard W. LUDWIG" <howard.w.ludwig@lmco.com> writes:

> I suppose my main point is really that we should be more careful in how we
> address questions of potential newbies (who may not even know there is an
> Ada 95 in contrast to Ada 83).  Some of us have this tendency to jump on
> such people, belittling their honest ignorance, and showing off our
> superior knowledge.  Such responses make the true newbie a bit wary of Ada
> based on an apparent elitist culture in the Ada community.  We don't need
> that kind of negative marketing.

Not to excuse any of such behaviour, but comp.lang.ada is much less
harsh than many newsgroups in such regards.

Larry Kilgallen

Re: Calculating SQRT in ADA

[robert_dewar]

In article <19990324201959.00800.00000708@ngol04.aol.com>,
  johnherro@aol.com (John Herro) wrote:
> Robert Dewar wrote:
> >>   while ...
> >>        abs(Guess*Guess/X - 1.0) > 3.0*Dummy'Epsilon
loop
> > Why the epsilon test, this iteration will terminate to
> > exact equality in any reasonable floating-point model
...
>
> I don't know if the Open Ada 83 compiler has a
"reasonable" floating-point
> model, but I do know that if I test for exact equality,
instantiate with type
> Float, and compute Sqrt(2.0) with that compiler, the
program goes into an
> infinite loop.  I think that in general, you can't
guarantee that Guess*Guess
> will ever exactly equal X.


Of course not!

The proper way for testing for convergence is to see if
the new guess is the same as the last guess!

-----------== Posted via Deja News, The Discussion Network ==----------
http://www.dejanews.com/       Search, Read, Discuss, or Start Your Own    

Re: Calculating SQRT in ADA

[robert_dewar]

In article <19990324201959.00800.00000708@ngol04.aol.com>,
  johnherro@aol.com (John Herro) wrote:
> I admit that in cases where execution time must be
> minimized, the extra programming effort is called for,
> but in general my first rule is to minimize
> program complexity instead.

The best way to minimize effort in computing square roots
is to use a library routine written by someone who knows
how to write things efficiently.

The example here is

(a) subject to doing too many operations per iteration
(b) takes more time that necessary per iteration, due to
      the very inefficient termination test
(c) is not fully last bit accurate

Sqrt is an easy routine to write, and any library routine
will undoubtedly do the right things with respect to the
above three points!

-----------== Posted via Deja News, The Discussion Network ==----------
http://www.dejanews.com/       Search, Read, Discuss, or Start Your Own    

Re: Calculating SQRT in ADA

[David C. Hoos, Sr.]

Mark A Biggar wrote in message <36FAA3DF.42C31CF@lmco.com>...
>Hans Marqvardsen wrote:
>>
>> John Herro wrote:
>
>No,no,no,no!
>
>This algorithm has quadratic convergence.  This means that each time
through the
>loop doubles the number of correct bits in the answer.  So with IEEE double
>precession float you need loop at most 6 times, and even with extended
float
>7 times, so don't test against a torerance at all, just loop 6 or 7 times
>and be done with it.
>
No,no,no,no! (to coin a phrase)

The number of correct bits doubles with every iteration, only if you are
counting place holding bits as well as mantissa bits.  So,
if you take the square root of 2.0e100, for example, 173 iterations are
required if the argument is used as the initial estimate.

similar number is required for 0.5e-100, but only six iterations are
required for 2.0e00 or 0.5e00.

This is why as Robert stated first, and I reiterated in my earlier post, the
coice of the initial estimate is very important, if
you're concerned about execution time.

If the machine representation is known, an improved initial estimate can be
made using the halving the exponent method I outlined
in my earlier post.

Re: Calculating SQRT in ADA

[robert_dewar]

In article <36F913E0.75F51763@lmco.com>,
  howard.w.ludwig@lmco.com wrote:
>
> --------------F93A017E5BB78033596BF5CF
> Content-Type: text/plain; charset=us-ascii
> Content-Transfer-Encoding: 7bit
>
> We need to be a little more careful in addressing
> questions of new people.
> The answer Chris gave _assumes_ Ada 95.  Many people, for
> whatever reason,
> still use Ada 83.  It is clear from another post by the
> original poster that
> Ada 83 is needed.  A more appropriate, complete answer is
> something like the
> following:

I disagree, it is generally confusing to talk about
obsolete languages unless necessary. If people ask about
Ada, it is reasonable to assume they mean Ada 95. If people
really want to ask about Ada 83 (or Ada 80, or whatever),
then they need to qualify the questions appropriately.

-----------== Posted via Deja News, The Discussion Network ==----------
http://www.dejanews.com/       Search, Read, Discuss, or Start Your Own    

Re: Calculating SQRT in ADA

[robert_dewar]

In article <19990324125518.01416.00000432@ngol08.aol.com>,
  johnherro@aol.com (John Herro) wrote:
> Here's a generic square root routine in Ada 83, taken
from my AdaTutor program:
>
> generic
>   type Dummy is digits <>;
> function Sqrt(X :in Dummy) return Dummy;
> function Sqrt(X :in Dummy) return Dummy is
>   Guess : Dummy := X;
> begin
>   if X < 0.0 then
>     raise Constraint_Error;
>   end if;
>   while X /= 0.0 and then
>        abs(Guess*Guess/X - 1.0) > 3.0*Dummy'Epsilon loop
>     Guess := (X/Guess + Guess) * 0.5;
>   end loop;
>   return Guess;
> end Sqrt;

Why the epsilon test, this iteration will terminate to
exact equality in any reasonable floating-point model,
and you cannot guarantee the 3*epsilon convergence if
the floating-point model is silly anyway.

-----------== Posted via Deja News, The Discussion Network ==----------
http://www.dejanews.com/       Search, Read, Discuss, or Start Your Own    

Re: Calculating SQRT in ADA

[robert_dewar]

In article <19990324125518.01416.00000432@ngol08.aol.com>,
  johnherro@aol.com (John Herro) wrote:
> Here's a generic square root routine in Ada 83, taken
from my AdaTutor program:
>
> generic
>   type Dummy is digits <>;
> function Sqrt(X :in Dummy) return Dummy;
> function Sqrt(X :in Dummy) return Dummy is
>   Guess : Dummy := X;

Notice also that no one would actually program a square
root this way (with such a poor initial guess). You need
to do a table lookup to get the initial guess reasonably
close, otherwise you waste several iterations getting near
the target value.

-----------== Posted via Deja News, The Discussion Network ==----------
http://www.dejanews.com/       Search, Read, Discuss, or Start Your Own    

Re: Ada 83 - Sometimes still chosen

[robert_dewar]

In article <7dbcj3$e0l@dfw-ixnews7.ix.netcom.com>,
  Richard D Riehle <laoXhai@ix.netcom.com> wrote:
  It will still be Ada 83 code.


And almost certainly it will still be Ada 95 code, seeing
as Ada 95 is almost exactly upwards compatible with Ada 83

-----------== Posted via Deja News, The Discussion Network ==----------
http://www.dejanews.com/       Search, Read, Discuss, or Start Your Own    

Re: Ada 83 - Sometimes still chosen

[robert_dewar]

In article <7dbcj3$e0l@dfw-ixnews7.ix.netcom.com>,
  Richard D Riehle <laoXhai@ix.netcom.com> wrote:
> In article <36F913E0.75F51763@lmco.com>,
> 	"Howard W. LUDWIG" <howard.w.ludwig@lmco.com>
wrote:
>
> >We need to be a little more careful in addressing
> >questions of new people.
> >The answer Chris gave _assumes_ Ada 95.  Many people,
> >for whatever reason,
> >still use Ada 83.
>
>  Howard's point needs to be emphasized.

Here is why I think it is a bad idea by default to always
talk about Ada 83.

Sure, we know that some projects get stuck using obsolete
hardware, software, operating systems, compilers etc, and
that is understandable (and not necessarily an
inappropriate decision).

But in comp.lang.ada, I think we do not want to be
emphasizing this kind of stuck-in-the-dark-ages view.
Ada 95 is a modern, recently designed language, and it
is important to always have the message that Ada is alive
and well (not just that ancient projects forced to use it
a long time ago are still alive and well).

Sure, if someone from one of these projects needs to know
something, we can probably provide the answer here to
questions about Ada 83, provided the question specifies
this (for that matter, this is probably the only place you
could still get Jovial questions answered :-)

-----------== Posted via Deja News, The Discussion Network ==----------
http://www.dejanews.com/       Search, Read, Discuss, or Start Your Own    

Re: Calculating SQRT in ADA

[John Herro]

Robert Dewar wrote:
>>   while ...
>>        abs(Guess*Guess/X - 1.0) > 3.0*Dummy'Epsilon loop
> Why the epsilon test, this iteration will terminate to
> exact equality in any reasonable floating-point model ...

I don't know if the Open Ada 83 compiler has a "reasonable" floating-point
model, but I do know that if I test for exact equality, instantiate with type
Float, and compute Sqrt(2.0) with that compiler, the program goes into an
infinite loop.  I think that in general, you can't guarantee that Guess*Guess
will ever exactly equal X.

Robert also wrote:
>>   Guess : Dummy := X;
> Notice also that no one would actually program a square
> root this way (with such a poor initial guess). You need
> to do a table lookup to get the initial guess reasonably
> close, otherwise you waste several iterations getting near
> the target value.

I agree that this poor initial guess causes the program to require several
extra iterations, and that a table lookup for the initial guess would shorten
the run time.  However, my philosophy is, wherever possible, to minimize
program complexity rather than execution time, to save programming and
debugging effort.  Computers are supposed to save people time.
Convergence is very rapid.  If we want to compute Sqrt(9.0) and the initial
guess is 9.0, successive guesses, to nine significant digits, are  5.00000000,
3.40000000, 3.02352941, 3.00009155, 3.00000000.
I admit that in cases where execution time must be minimized, the extra
programming effort is called for, but in general my first rule is to minimize
program complexity instead.
- John Herro
You can download a shareware AdaTutor program at
http://members.aol.com/AdaTutor

Re: Calculating SQRT in Ada

[John Herro]

Robert Dewar said:
> The proper way for testing for convergence is to see if
> the new guess is the same as the last guess!

Hans Marqvardsen said:
>>   while X /= 0.0 and then
>>        abs(Guess*Guess/X - 1.0) > 3.0*Dummy'Epsilon loop
> Where does the multiplier 3.0 come from?

There are some cases where the loop will alternate between two guesses, so
testing that the new guess is the same as the last guess won't guarantee
convergence with every implementation.  However, it might be possible to
terminate when the guess is the same as EITHER the last guess OR the
second-to-last guess.  I can't prove to myself that the loop would always
terminate in this case, but it seems likely that it will.

I've found that with a multiplier of 1.0, the loop doesn't always terminate. 
With a multiplier of 10.0, not all displayed digits of the answer are correct
when I simply use Put in an instantiation of Float_IO, without specifying Aft,
Exp, etc.  So I chose a multiplier that's approximately the geometric mean of
1.0 and 10.0.  With the multiplier of 3.0, the loop seems always to terminate,
and all displayed digits are correct.  However, Robert Dewar is right when he
says

> The example here ... is not fully last bit accurate
> ... and any library routine will undoubtedly do the
> right things ...

It's interesting to note that I have yet to see an Ada 83 compiler (except
Ada/Ed), that DIDN'T come with a math package that included a square root
routine.
- John Herro
http://members.aol.com/AdaTutor

Re: Calculating SQRT in ADA

[Tom Moran]

I found "Software Manual for the Elementary Functions", by Cody and
Waite, ISBN 0-13-822064-6 quite helpful, though perhaps there's
something newer out.  As I recall, a fixed point 386 asm version with
an initial estimate only needed something like a hard coded three or
four iterations and was quite fast.

Re: Ada 83 - Sometimes still chosen

[robert_dewar]

In article <7dc9ld$2u6@dfw-ixnews9.ix.netcom.com>,
  Richard D Riehle <laoXhai@ix.netcom.com> wrote:
> In article <7dbvl2$5bl$1@nnrp1.dejanews.com>,
> 	robert_dewar@my-dejanews.com wrote:
>  I am uneasy with the selection of the word "obsolete."
>  A technology is obsolete if it is not the best choice
>  for the solution to a particular problem.

I disagree, I would say that sometimes the use of obsolete
equipment is the right solution. I just bought a 1987 Ford
F250 truck, definitely an obsolete piece of equipment, but
just fine for my purposes.

>  I return to the example of the 1750A.  This processor,
>  by all standards for earth-bound computing is obsolete.
>  For a variety of technical reasons, it is still the
>  optimal choice for many space applications.

Actually we see very few projects selecting the 1750A for
new projects. There are other alternatives these days, such
as the ERC.

>  Sometimes the economics of an 8086 can make it the
> correct choice.

Very rarely these days. A 68K is almost always the better
choice, and is far more often used, which is why there are
Ada 95 compilers available for the 68K.

These days, it is quite inexpensive to generate a new
Ada 95 compiler, and if there are systems for which Ada 95
is not used, it is an indicator that very few new projects
are choosing the hardware in question (yes, I realize that
there is a chicken and egg problem, but it is minimal given
the low cost of porting an Ada 95 compiler). If no Ada 95
compiler is available for a given system, it almost
certainly means that the customer demand for such a
compiler is minimal or non-existant.

My fundamental point here is that the general impression
of the community is that Ada 83 was a failure. Yes, that
is hyperobole of course, but on the other hand, at this
stage, Ada 83 is pretty creaky.

Just yesterday, I talked to someone doing research into
object oriented component design. He was trying to fit
into C++, and having various troubles. He made a list of
deficiencies of C++, all fixed in Ada 95. He then was
playing with a modified Java to accomodate his ideas, all
of which could have been used directly in Ada 95.

When I asked him about Ada, he replied "Oh, does Ada have
object oriented features, I didn't know that".

Now suppose he were to ask on this list if Ada has object
oriented features. I would far rather he get an unambiguous
yes, rather than a:

  "well it depends which version of Ada you are using"

type response, which could easily be interpreted to mean
"not really".

The number of people asking about Ada 83 on this list is
small. The number of people asking about Ada 83 on this
list who do not know that they are asking about Ada 83 and
thus do not say is even smaller.

So of the people asking questions about Ada (without
specifying which version), the overwhelming majority are
talking about Ada 95, and it is a definite disservice to
tell people that Ada 83 is still in wide use and that they
must worry about the Ada 83 answer as well as the Ada 95
answer. This would be like asking a question on the C++
group, and having responders worry that you might really be
asking about C!

-----------== Posted via Deja News, The Discussion Network ==----------
http://www.dejanews.com/       Search, Read, Discuss, or Start Your Own    

Re: Ada 83 - Sometimes still chosen

[robert_dewar]

In article <7dc7tu$2uk@dfw-ixnews7.ix.netcom.com>,
  Richard D Riehle <laoXhai@ix.netcom.com> wrote:
>  On some projects some programming managers feel it is
>  appropriate to freeze even the version of the compiler
>  used to develop some long-lived software product.  That
>  is, once the satellite is deployed, all updates will be
>  made using exactly the same language and version of the
>  compiler used for the original development.

Richard, of *course* I know this (I am after all working as
an Ada compiler vendor these days! :-)

Yes, of course some projects are still using Ada 83, and
there are even a few new projects using Ada 83 (though this
is very rapidly dwindling).

Absolutely no one that I know of considers this
unreasonable, for all the reasons you state.

But that has nothing whatsoever to do with the point at
hand, which, please remember, is the following.

If someone asks a question about Ada on comp.lang.ada, do
we assume that they are asking about Ada 95, a language
that has now been around for four years, or do we by
default worry them about the possibility of different
answers for the now obsolete Ada 83 (obsolete *is* the
right word when a new ISO standard supplants an old one).

I think it is very important NOT to confuse answers by
referring to shortcomings in Ada 83. If people are working
on a project using Ada 83, they should know what they are
doing, and ask the question appropriately. Sure, I can see
that students might not know, but they are almost *sure*
to be using Ada 95.

I see no reason to confuse the discussions and answers on
comp.lang.ada for the sake of the small number of people
working on large projects under the conditions you describe
who also do not know enough to specify Ada 83 when they ask
a question!

Robert Dewar



-----------== Posted via Deja News, The Discussion Network ==----------
http://www.dejanews.com/       Search, Read, Discuss, or Start Your Own    

Re: Calculating SQRT in Ada

[David C. Hoos, Sr.]

John Herro wrote in message
<19990325232830.01123.00000079@ngol02.aol.com>...
>Robert Dewar said:
>> The proper way for testing for convergence is to see if
>> the new guess is the same as the last guess!
>
>Hans Marqvardsen said:
>>>   while X /= 0.0 and then
>>>        abs(Guess*Guess/X - 1.0) > 3.0*Dummy'Epsilon loop
>> Where does the multiplier 3.0 come from?
>
>There are some cases where the loop will alternate between two guesses, so
>testing that the new guess is the same as the last guess won't guarantee
>convergence with every implementation.  However, it might be possible to
>terminate when the guess is the same as EITHER the last guess OR the
>second-to-last guess.  I can't prove to myself that the loop would always
>terminate in this case, but it seems likely that it will.
>
John -- Did you look at my convergence test? It always works.

Re: Calculating SQRT in Ada

[John Herro]

David C. Hoos, Sr. wrote
> while  not alternating and then
>  New_Estimate /= Previous_estimate loop

> John -- Did you look at my convergence test?
> It always works.

Your routine, although a bit more complicated than my original post, is easy to
understand, and it does indeed look like it will always converge.
- John Herro
http://members.aol.com/AdaTutor

Re: Calculating SQRT in ADA

[Ole-Hjalmar Kristensen]

robert_dewar@my-dejanews.com writes:

> In article <19990324201959.00800.00000708@ngol04.aol.com>,
>   johnherro@aol.com (John Herro) wrote:
> > I admit that in cases where execution time must be
> > minimized, the extra programming effort is called for,
> > but in general my first rule is to minimize
> > program complexity instead.
> 
> The best way to minimize effort in computing square roots
> is to use a library routine written by someone who knows
> how to write things efficiently.
> 
> The example here is
> 
> (a) subject to doing too many operations per iteration
> (b) takes more time that necessary per iteration, due to
>       the very inefficient termination test
> (c) is not fully last bit accurate
> 
> Sqrt is an easy routine to write, and any library routine
> will undoubtedly do the right things with respect to the
> above three points!
> 
> -----------== Posted via Deja News, The Discussion Network ==----------
> http://www.dejanews.com/       Search, Read, Discuss, or Start Your Own    


Your faith in library routines is touching :-)

-- 
E pluribus Unix

Re: Ada 83 - Sometimes still chosen

[Tarjei Tj�stheim Jensen]

robert_dewar@my-dejanews.com wrote:
>So of the people asking questions about Ada (without
>specifying which version), the overwhelming majority are
>talking about Ada 95, and it is a definite disservice to
>tell people that Ada 83 is still in wide use and that they
>must worry about the Ada 83 answer as well as the Ada 95
>answer. This would be like asking a question on the C++
>group, and having responders worry that you might really be
>asking about C!


I think it would be more apropriate to use ANSI C versus K&R C as an example.

Anyway. I wholehartedly agree that Ada 95 is Ada. Ada 83 is Ada 83 and is as
obsolete as K&R C is obsolete in the C world.

Greetings,

Re: Calculating SQRT in ADA

[als0045]

[-- Attachment #1: Type: text/plain, Size: 536 bytes --]

Nice to see so much interest in a numerics thread!
In practice it is good to avoid flt. point "/"  ... it's usually
10 to 50 times slower than "*".  One trick is to calculate
Inverse_Sqrt followed by one call to "/".  THere are
no "/"'s in the Newton-Raphson calculation for Inverse_Sqrt.

Inverse_Sqrt is extremely useful indepepently of Sqrt.
The following Inverse_Sqrt is usually faster than 1.0 / Sqrt on the
machines I tested, (but 1.0 / Inverse_Sqrt is not necessarily as fast or
as accurate as the machine's Sqrt).

--  Jonathan

[-- Attachment #2: sqr2 --]
[-- Type: text/plain, Size: 9004 bytes --]

--
--  Uses Newton-Raphson to get Y = A**(-1/N):
--    Inverse_Nth_Root =  Y_k+1 = Y_k + (1 - A * Y_k**N) * Y_k / N.
--  where  N is 2.
--
generic 

  type Real is digits <>;
  --  Anything up to Real'Digits = 24 digits is OK.  This is not checked.

package invsqrt is
     
  function Inverse_Sqrt (Arg : Real) return Real;
  pragma Inline(Inverse_Sqrt);

end invsqrt;
package body invsqrt is

  function Inverse4_Sqrt_At (X : Real) return Real;
--pragma Inline (Inverse4_Sqrt_At);

  ---------------------
  -- Inverse_Sqrt_At --
  ---------------------
  --
  -- Fit quartic spline to 2 pts. on Inverse_Sqrt curve.
  -- Should be good to 10+ bits on [1,4]. IF  No_of_Iterations := 3, THEN
  -- GOOD ENOUGH FOR  Real'Digits = 24.
  --   
  function Inverse4_Sqrt_At (X : Real) return Real is
     X3, X2 : Real;
   --X4, X3, X2 : Real;
     Y : Real;
     a1 : constant Real :=  2.05658516425871E+00;
     b1 : constant Real := -1.89528168315569E+00;
     c1 : constant Real :=  1.16136225073854E+00;
     d1 : constant Real := -3.70446006083701E-01;
     e1 : constant Real :=  4.74502518096651E-02;

     a2 : constant Real :=  1.45422531573500E+00;
     b2 : constant Real := -6.70083265209039E-01;
     c2 : constant Real :=  2.05301780727833E-01;
     d2 : constant Real := -3.27431103706594E-02;
     e2 : constant Real :=  2.09702467647666E-03;

  begin 

    X2 := X*X; 
    X3 := X2*X; 
    --X4 := X2*X2; 

    IF X > 2.0 THEN  

    -- Y   := a2 + X*(b2 + X*(c2 + X*(d2 + X*e2)));
    -- Y   := a2 + b2*X + c2*X2 + d2*X3 + e2*X4;  -- fastest on dec alpha
       Y   := a2 + b2*X + c2*X2 + (d2 + e2*X)*X3 ;

    ELSE

    -- Y   := a1 + X*(b1 + X*(c1 + X*(d1 + X*e1)));
    -- Y   := a1 + b1*X + c1*X2 + d1*X3 + e1*X4;
       Y   := a1 + b1*X + c1*X2 + (d1 + e1*X)*X3;

    END IF;

    return Y;
 
  end Inverse4_Sqrt_At;

  ---------------------
  -- Inverse_Sqrt_At --
  ---------------------
  --
  -- Fit cubic spline to 2 pts.
  -- Should be good to 7+ bits on [1,4].  IF  No_of_Iterations := 3, THEN
  -- ALMOST GOOD ENOUGH FOR  Real'Digits = 18, SO USE IF Real'Digits = 15.
  --   
  function Inverse3_Sqrt_At (X : Real) return Real is
     X3, X2 : Real;
     Y : Real;
     a1 : constant Real :=  1.84286198063662E+00;  
     b1 : constant Real := -1.29816848576085E+00;
     c1 : constant Real :=  5.43101858353914E-01;
     d1 : constant Real := -8.92682549640998E-02;

     a2 : constant Real :=  1.30310020329903E+00;
     b2 : constant Real := -4.58971869702085E-01;
     c2 : constant Real :=  9.60077517292674E-02;
     d2 : constant Real := -7.89027355372562E-03;

  begin 

    X2 := X*X; 
    X3 := X2*X; 

    IF X > 2.0 THEN  

     Y   := a2 + b2*X + c2*X2 + d2*X3;

    ELSE

     Y   := a1 + b1*X + c1*X2 + d1*X3;

    END IF;

    return Y;
 
  end Inverse3_Sqrt_At;

  pragma Inline (Inverse3_Sqrt_At);

  ------------------
  -- Inverse_Sqrt --
  ------------------
  --
  --  Uses Newton-Raphson to get Y = A**(-1/N):
  --    Inverse_Nth_Root =  Y_k+1 = Y_k + (1 - A * Y_k**N) * Y_k / N.
  --
  --     N : constant Integer := 2;
  --
  --IF Real'Digits < 25 THEN  
  --   No_of_Iterations := 3; -- Should be good to >80 bits.
  --ELSE
  --   No_of_Iterations := 4; --  Should be good to >160 bits.
  --END IF;
  --   
  function Inverse_Sqrt (Arg : Real) return Real is
   
     Scaled_Arg_Over_N  : Real;
     Scaled_Arg   : Real := Arg;
     Scale_Factor : Real := 1.0;
     Y_0          : Real;
     N            : constant Integer := 2;
     One_Over_N   : constant Real := 0.5;
     One_Plus_One_Over_N  : constant Real := 1.0 + One_Over_N;

     No_of_Iterations  : constant Integer := 3;

  begin 
   
   IF Arg  <= 0.0 THEN
      raise Constraint_error;
   END IF;
   
   IF Arg > 1.0 THEN

   IF Scaled_Arg > 256.0 THEN
      Scaled_Arg   := Scaled_Arg * 0.0625**2; 
      Scale_Factor := Scale_Factor * 0.0625;
      while Scaled_Arg > 256.0 loop
         Scaled_Arg   := Scaled_Arg * 0.0625**2; 
         Scale_Factor := Scale_Factor * 0.0625;
      end loop;
   END IF;

   IF Scaled_Arg > 16.0 THEN
      Scaled_Arg   := Scaled_Arg * 0.0625; 
      Scale_Factor := Scale_Factor * 0.25;
   END IF;

   IF Scaled_Arg > 4.0 THEN
      Scaled_Arg   := Scaled_Arg * 0.25; 
      Scale_Factor := Scale_Factor * 0.5;
   END IF;

   ELSE

   IF Scaled_Arg < 0.0625**2 THEN
      Scaled_Arg   := Scaled_Arg * 256.0; 
      Scale_Factor := Scale_Factor * 16.0;
      while Scaled_Arg < 0.0625**2 loop
         Scaled_Arg   := Scaled_Arg * 256.0; 
         Scale_Factor := Scale_Factor * 16.0;
      end loop;
   END IF;

   IF Scaled_Arg < 0.0625 THEN
      Scaled_Arg   := Scaled_Arg * 16.0; 
      Scale_Factor := Scale_Factor * 4.0;
   END IF;

   --IF Scaled_Arg < 0.25 THEN
   --   Scaled_Arg   := Scaled_Arg * 4.0; 
   --   Scale_Factor := Scale_Factor * 2.0;
   --END IF;

   IF Scaled_Arg < 0.25 THEN
      Scaled_Arg   := Scaled_Arg * 16.0; 
      Scale_Factor := Scale_Factor * 4.0;
   ELSE
      Scaled_Arg   := Scaled_Arg * 4.0; 
      Scale_Factor := Scale_Factor * 2.0;
   END IF;

   END IF;
  
     Scaled_Arg_Over_N := One_Over_N * Scaled_Arg; 

     Y_0 := Inverse4_Sqrt_At (Scaled_Arg); 
 
     for j in 1 .. No_of_Iterations-1 loop
        Y_0 := Y_0 * (One_Plus_One_Over_N - Scaled_Arg_Over_N * Y_0*Y_0); 
     end loop;
     
     -- last iteration (better accuracy, esp. at large N):
     Y_0 := Y_0 + Y_0 * (One_Over_N - Scaled_Arg_Over_N * Y_0*Y_0);

     return Y_0 * Scale_Factor;
      
  end Inverse_Sqrt;

end invsqrt;
with Text_IO;  use Text_IO;
with invsqrt; 
procedure sqr1tst1 is

 type Real is digits 15;
 
 Max : Real := 0.0;
 X, Y : Real := 0.0;

 Package FLIO is new Float_IO (Real);
 use FLIO;

 Package invsqr is new invsqrt (Real);
 use invsqr;

begin

  X :=  1.00000000;

  Max := -1.0E20;
  for I in 1..10000000 loop
  
    Y := 1.0 - Inverse_SQRT(X*X) * X;
    IF Abs Y > Max THEN    
       Max := Abs Y;
    END IF;
    X := X + 2.0*1024.0**2;
  end loop;
  
  new_line; put("Max error is:  "); put(Max); New_line;
  
  X :=  0.000000001;

  Max := -1.0E20;
  for I in 1..10000000 loop
  
    X := X + 1000.1;
    Y := 1.0 - Inverse_SQRT(X*X) * X;
    IF Abs Y > Max THEN    
       Max := Abs Y;
    END IF;
  end loop;
  
  new_line; put("Max error is:  "); put(Max); New_line;
  
  X :=  0.000000001;

  Max := -1.0E20;
  for I in 1..10000000 loop
  
    X := X + 1.1;
    Y := 1.0 - Inverse_SQRT(X*X) * X;
    IF Abs Y > Max THEN    
       Max := Abs Y;
    END IF;
  end loop;
  
  new_line; put("Max error is:  "); put(Max); New_line;
  
  X :=  0.000000001;

  Max := -1.0E20;
  for I in 1..10000000 loop
  
    X := X + 0.001;
    Y := 1.0 - Inverse_SQRT(X*X) * X;
    IF Abs Y > Max THEN    
       Max := Abs Y;
    END IF;
  end loop;
  
  new_line; put("Max error is:  "); put(Max); New_line;
  
  X :=  0.000000001;

  Max := -1.0E20;
  for I in 1..10000000 loop
  
    X := X + 0.00001;
    Y := 1.0 - Inverse_SQRT(X*X) * X;
    IF Abs Y > Max THEN    
       Max := Abs Y;
    END IF;
  end loop;
  
  new_line; put("Max error is:  "); put(Max); New_line;
  
  X :=  0.1;

  Max := -1.0E20;
  for I in 1..10000000 loop
  
    X := X + 0.0000001;
    Y := 1.0 - Inverse_SQRT(X*X) * X;
    IF Abs Y > Max THEN    
       Max := Abs Y;
    END IF;
  end loop;
  
  new_line; put("Max error is:  "); put(Max); New_line;

  X :=  0.0000000001111;

  Max := -1.0E20;
  for I in 1..10000000 loop
  
    X := X + 0.01111;
    Y := 1.0 - Inverse_SQRT(X**2) * X;
    IF Abs Y > Max THEN    
       Max := Abs Y;
    END IF;
  end loop;
  
  new_line; put("Max error is:  "); put(Max); New_line;
  
  X :=  0.0000000001111;

  Max := -1.0E20;
  for I in 1..10000000 loop
  
    X := X + 0.0001111;
    Y := 1.0 - Inverse_SQRT(X**2) * X;
    IF Abs Y > Max THEN    
       Max := Abs Y;
    END IF;
  end loop;
  
  new_line; put("Max error is:  "); put(Max); New_line;
  
  X :=  0.01;

  Max := -1.0E20;
  for I in 1..10000000 loop
  
    X := X + 0.000001111;
    Y := 1.0 - Inverse_SQRT(X**2) * X;
    IF Abs Y > Max THEN    
       Max := Abs Y;
    END IF;
  end loop;
  
  new_line; put("Max error is:  "); put(Max); New_line;
  
  X :=  0.0;

  Max := -1.0E20;
  for I in 1..10000000 loop
  
    X := X + 0.00000000000000111;
    Y := 1.0 - Inverse_SQRT(X**2) * X;
    IF Abs Y > Max THEN    
       Max := Abs Y;
    END IF;
  end loop;
  
  new_line; put("Max error is:  "); put(Max); New_line;
  
end;
with Ada.Numerics.Generic_Elementary_Functions;
with Text_IO;  use Text_IO;
with invsqrt; 
procedure sqr1tst2 is

 type Integer32 is range -2**31..2**31-1;
 
 type Real is digits 15;
 
 X, Y : Real := 0.0;

 package Mathz is new Ada.Numerics.Generic_Elementary_Functions(Real);
 use Mathz;

 Package FLIO is new Float_IO (Real);
 use FLIO;

 Package invsqr is new invsqrt (Real);
 use invsqr;

begin

  X :=  0.1;
  Y :=  0.0;

  for I in 1..100000000 loop
  
    X := X + 1.0e-7;
    Y := Y + Inverse_SQRT(X);
   -- Y := Y + 1.0 / SQRT(X);

  end loop;
  
  new_line;  put(Y); New_line;
  
end;

Re: Calculating SQRT in ADA

[Howard W. LUDWIG]

David C. Hoos, Sr. wrote:

> Mark A Biggar wrote in message <36FAA3DF.42C31CF@lmco.com>...
> >Hans Marqvardsen wrote:
> >>
> >> John Herro wrote:
> >
> >No,no,no,no!
> >
> >This algorithm has quadratic convergence.  This means that each time
> through the
> >loop doubles the number of correct bits in the answer.  So with IEEE double
> >precession float you need loop at most 6 times, and even with extended
> float
> >7 times, so don't test against a torerance at all, just loop 6 or 7 times
> >and be done with it.
> >
> No,no,no,no! (to coin a phrase)
>
> The number of correct bits doubles with every iteration, only if you are
> counting place holding bits as well as mantissa bits.  So,
> if you take the square root of 2.0e100, for example, 173 iterations are
> required if the argument is used as the initial estimate.
>
> similar number is required for 0.5e-100, but only six iterations are
> required for 2.0e00 or 0.5e00.

No, no, no, no! (to plagiarize a coined phrase)

It is a common oversimplified assumption that Newton[-Raphson] has quadratic
convergence.  There are many well-known cases where the method does not converge
at all.  Even when the method does converge, the convergence can be linear.
Quadratic convergence applies when the quadratic term is negligible compared to
the linear term in the Taylor series expansion of the function in question.
When the initial guess is far off, "negligible" does not describe the situation.

To find the actual  a = sqrt(c), where c is known, involves calculating
    x(n+1) = 0.5 * (x(n) + c/x(n)).
Let x(n) = a + d(n), where d(n) measures the deviation of the estimate x(n) from
the actual value a.  Then
    d(n+1) = x(n+1) - a
           = 0.5 * (a + d(n) + a**2/(a + d(n))) - a
           = 0.5 * d(n)**2/(a + d(n)).
When d(n) is small compared to a, the denominator behaves more and more like a,
so the deviation is indeed proportional to d(n)**2.  However, when d(n) is large
compared to a, the denominator for some iterations looks more like d(n),
canceling out one of the d(n) factors in the numerator, so the deviation is cut
each time by a factor of 2, thus linear convergence, until the iterations get to
the place of having the deviation be of the order of a, at which time the
convergence transitions to quadratic.

In the case of the example given [c = 2E100, a =1.4142...E50], starting with an
initial guess of c requires 166 iterations of halving the deviation each time to
get to the point of the deviation being smaller than the actual root a.  Then
quadratic convergence kicks in, apparently requiring only another 7 iterations
to achieve the ultimate desired result. If we simply counted place-holder bits,
roughly 167, plus the desired 53 bits of precision, then we have about 220 bits
to deal with.  _Assuming_ the first bit to be fine and doubling the number of
good bits each iteration, only 8 iterations would be needed (2^8 = 256), nowhere
close to the 173 iterations quoted above.  This shows that you cannot assume
that the first bit is good.  If you want to do a "back of the envelope"
calculation based on quadratic convergence, the 2E100 case using 2E100 as the
initial guess would say your first guess has less than 1E-50 bit good, in order
to double the number of good bits each iteration and to obtain 53 good bits
after 173 iterations.  A better way of looking at it is that the place-holder
bits are lopped off at the rate of one per iteration for the first 166
iterations.

HWL

Re: Calculating SQRT in ADA

[Marin David Condic]

cmcrae wrote:
> 
>   I am trying to calculate the SQRT in ADA.  Is there a function that I
> can access?
> --
> Posted via Talkway - http://www.talkway.com
> Surf Usenet at home, on the road, and by email -- always at Talkway.

I certainly hope that the degeneration of this thread into a detailed
discussion of how to compute SQRT has not discouraged you from using Ada
:-) I hope you ultimately found a satisfactory answer. (Either the
numerics appendix in A.5 for Ada95 or compiler specific - post the brand
and platform & maybe we can help.) 

I sometimes wonder if questions like this aren't just tossed out to see
how many people will run off on a tangent? :-) "Hey everybody! I'm
looking for the best sorting algorithm available in Ada. Can somebody
tell me what is the best algorithm?" :-))

MDC
-- 
Marin David Condic
Real Time & Embedded Systems, Propulsion Systems Analysis
United Technologies, Pratt & Whitney, Large Military Engines
M/S 731-95, P.O.B. 109600, West Palm Beach, FL, 33410-9600
***To reply, remove "bogon" from the domain name.***

Re: Calculating SQRT in ADA

[als0045]

als0045 wrote:
> 
> One trick is to calculate
> Inverse_Sqrt followed by one call to "/". 
No, dummy, just multiply your 1 / Sqrt(x) by x.

Come to think of it, One_Over_Sqrt(x) is a much better name for it.

--     Jonathan Parker 
--     Department of Applied Mathematics and Theoretical Physics
--     The Queen's University of Belfast

Re: Ada 83 - Sometimes still chosen

[Tom Moran]

>sometimes the use of obsolete
>equipment is the right solution
As I think somebody pointed out, a better word in this case is
"obsolescent" meaning "becoming obsolete".

Re: Ada 83 - Sometimes still chosen

[Tom Moran]

>If they fail to specify the 
> version of Ada...
If an elementary school student asks you "is the moon round" the
appropriate answer is "yes".  But if a spacecraft designer asks you
the same question, the answer is something like "approximately, but
with such and such deviations".  To answer a question it really is
necessary to make some assumptions about where the questioner is
coming from.  Assuming they are simply too stupid to have looked in
the Ada 95 LRM may not be an appropriate assumption, though teachers
of beginning programming may be biased toward that assumption.

Re: Calculating SQRT in ADA

[David C. Hoos, Sr.]

Marin David Condic wrote in message <36FBCEE2.5DCA764F@pwfl.com>...
>cmcrae wrote:
>>
>>   I am trying to calculate the SQRT in ADA.  Is there a function that I
>> can access?
>> --
>> Posted via Talkway - http://www.talkway.com
>> Surf Usenet at home, on the road, and by email -- always at Talkway.
>
>I certainly hope that the degeneration of this thread into a detailed
>discussion of how to compute SQRT has not discouraged you from using Ada
>:-) I hope you ultimately found a satisfactory answer. (Either the
>numerics appendix in A.5 for Ada95 or compiler specific - post the brand
>and platform & maybe we can help.)
>
>I sometimes wonder if questions like this aren't just tossed out to see
>how many people will run off on a tangent? :-) "Hey everybody! I'm
>looking for the best sorting algorithm available in Ada. Can somebody
>tell me what is the best algorithm?" :-))
>
At the risk of continuing on a tangent ... although the point for Sqrt _per se_ is moot, discussions of quadratic convergence, etc.,
is pertinent to many problems for which the standard libraries don't provide solutions.

To hopefully wrap this up, Ada95 provides a mechanism for cheaply generating the first estimate for square root, which had to be
done in a machine-specific way in Ada83.  The new Ada95 attributes 'Exponent and 'Compose provide the means, as demonstrated in the
example code which follows.

This implementation takes the best of all of the "roll your own" solutions discussed in this thread, and provides a solution in less
than six iterations.

It should be pointed out, however, that the gnat implementation which imports the C library's "sqrt" function is about ten times
faster than this one, on Solaris 7, x86

with Ada.Text_Io;
procedure Test_Sqrt is
   generic
      type Real is digits <>;
   function Sqrt (X : real) return Real;
   function Sqrt (X : real) return real is
      Previous_Previous_Estimate : Real :=
        Real'Compose (0.5, (Real'Exponent (X) + 1) / 2);
      Previous_Estimate          : Real;
      New_Estimate               : Real;
   begin
      if X < 0.0 then
         raise Constraint_Error;
      end if;

      -- Special case Sqrt (0.0) to preserve possible minus sign per
      -- IEEE?? This seeems strange -- why would we ever return a
      -- negative result for a real square root?
      if X = 0.0 then
         return X;
      end if;

      Previous_Estimate := 0.5 *
        (X / Previous_Previous_Estimate + Previous_Previous_Estimate);
      New_Estimate := 0.5 *
        (X / Previous_Estimate + Previous_Estimate);

      while  New_Estimate /= Previous_Estimate and then
        New_Estimate /= Previous_Previous_Estimate loop
         Previous_Previous_Estimate := Previous_Estimate;
         Previous_Estimate := New_Estimate;
         New_Estimate :=
           0.5 * (X / Previous_Estimate + Previous_Estimate);
      end loop;

      return New_Estimate;

   end Sqrt;

   function Long_Float_Sqrt is new Sqrt (Long_Float);

   Root : Long_Float;
   begin
      Ada.Text_Io.Put_Line
        (Long_Float'Image (Long_Float_Sqrt (2.0)));
      Ada.Text_Io.Put_Line
        (Long_Float'Image (Long_Float_Sqrt (0.5)));
      Ada.Text_Io.Put_Line
        (Long_Float'Image (Long_Float_Sqrt (2.0e100)));
      Ada.Text_Io.Put_Line
        (Long_Float'Image (Long_Float_Sqrt (0.5e-100)));
      Ada.Text_Io.Put_Line
        (Long_Float'Image (Long_Float_Sqrt (4.0)));
      Ada.Text_Io.Put_Line
        (Long_Float'Image (Long_Float_Sqrt (0.25)));
   end Test_Sqrt;

Re: Ada 83 - Sometimes still chosen

[Larry Kilgallen]

In article <36fbe743.5372929@news.pacbell.net>, tmoran@bix.com (Tom Moran) writes:
>>If they fail to specify the 
>> version of Ada...
> If an elementary school student asks you "is the moon round" the
> appropriate answer is "yes".  But if a spacecraft designer asks you
> the same question, the answer is something like "approximately, but
> with such and such deviations".

Whereas through internet newsgroups you can hear that the moon is closer
to a sphere than to a circle :-)

Larry Kilgallen

Re: Calculating SQRT in ADA

[bob]

Niklas Holsti <nholsti@icon.fi> wrote in article
<36F94393.5D1B317D@icon.fi>...
> bob wrote:
>  
>    [snip]
> 
> > In Ada83 or Ada95, if the target is real, possibly the following:
> > x : real;
> > [snip]
> > x = x**0.5;
> 
> Nope, for the standard "**" the right operand must be integer
> (Integer'Base, actually).

Agree for the standard. However many supply packages for the real**real.
One that comes to mind immediately is
"Ada.Numerics.Generic_Elementary_Functions in gnat.

with Ada.Numerics.Generic_Elementary_Functions;
with text_io;
use  text_io;

procedure sqrt is
  subtype real is long_float;
  package math is new Ada.Numerics.Generic_Elementary_Functions(real);
  use math;
  package fio is new text_io.float_io(real);
  use fio;

  x : real := 2.3;
begin
  put(" The sqrt of ");
  put(x);
  put(" is ");
  put(x**0.5);
  new_line;
end sqrt;

ada sqrt.adb
gcc -c -g sqrt.adb
gnatbind -x sqrt.ali
gnatlink -g sqrt.ali
linux:~/atest> sqrt
 The sqrt of  2.30000000000000E+00 is  1.51657508881031E+00

> 
> > All of the Ada85 comilers I have used supply a "math" package
> > containing sqrt also.
> 
> That agrees with my experience and would be the first place to
> look.
> 
> Niklas Holsti
> Working at but not speaking for Space Systems Finland Ltd.
> 

cheers

Re: Ada 83 - Sometimes still chosen

[Richard D Riehle]

In article <7dfqr1$g5j$1@nnrp1.dejanews.com>,
	robert_dewar@my-dejanews.com wrote:

>I disagree, I would say that sometimes the use of obsolete
>equipment is the right solution. I just bought a 1987 Ford
>F250 truck, definitely an obsolete piece of equipment, but
>just fine for my purposes.

 I suddenly have this image of Robert in a tattered baseball
 cap and faded bluejeans, a little spot of grease on the tip of his
 nose from tinkering under the hood of his pickup truck.  Those
 who know Robert might also be amused by this image.  Next, I
 expect him to reveal that he has abandoned the pipe organ for
 an avocation as a country and western singer with a steel string
 guitar.  :-)

 Seriously, though, we simply have different ideas of the meaning
 of obsolete.  I am typing this on one of my office computers, a
 Pentium 233.  It is of course, obsolete by one meaning.  It is not
 obsolete for this purpose.  The 1987 Ford is not obsolete if it
 does the job intended.

>Actually we see very few projects selecting the 1750A for
>new projects. There are other alternatives these days, such
>as the ERC.

 I just recently had yet another inquiry from a [prospective]
 client who has selected the 1750A for a space application. It
 turns out that the alternatives to 1750A are not universally
 trusted by developers of communications satellites.  BTW, is
 there an Ada 95 for ERC?
 
>These days, it is quite inexpensive to generate a new
>Ada 95 compiler, and if there are systems for which Ada 95
>is not used, 

 Some developers of space applications are conservative.
 If a technology is proven over time, has been successful
 on other projects, and personnel are experienced with
 it, the risk of staying with that technology might be deemed 
 less severe.  This is the well-known "If it aint broke don't
 fix it" viewpoint.  A surprising number of program managers
 adopt this peculiar stategy when designing for safety-critical
 software.  

>it is an indicator that very few new projects
>are choosing the hardware in question (yes, I realize that
>there is a chicken and egg problem, but it is minimal given
>the low cost of porting an Ada 95 compiler). If no Ada 95
>compiler is available for a given system, it almost
>certainly means that the customer demand for such a
>compiler is minimal or non-existant.

 I interpret this to mean that very few projects are inquiring
 of ACT about the availability of a GNAT solution.  I wonder if
 a compiler publisher such as DDCI is having the same experience?
 The last I heard, they were still selling a lot of Ada 83 compilers
 to their international customer base.  

>My fundamental point here is that the general impression
>of the community is that Ada 83 was a failure. Yes, that
>is hyperobole of course, but on the other hand, at this
>stage, Ada 83 is pretty creaky.

 Creaky?  I hear the passenger's door of Robert's Ford truck opening
 on a winter morning.  I hear the engine cranking its first rotation
 of the morning.  

 The success of Ada 83 on real projects would not lead me to use
 the word "creaky" to describe it.  I am, as are you, an advocate
 of transitioning to Ada 95 whenever feasible.  I am delighted
 with the improvements in the standard and agree that Ada is the
 current Ada 95.  On the other hand, many of those improvements are
 superfluous in the design of working software.  For the space
 applications mentioned earlier, one of the few really useful 
 improvements is unsigned numerics.  Many communication satellite
 designs avoid tasking, have no use for polymorphism, and are generally 
 conservative in the use of many other Ada features.    

>
>Just yesterday, I talked to someone doing research into
>object oriented component design. He was trying to fit
>into C++, and having various troubles. He made a list of
>deficiencies of C++, all fixed in Ada 95. He then was
>playing with a modified Java to accomodate his ideas, all
>of which could have been used directly in Ada 95.

 The list of deficiencies must have been long and mournful.
 I have decided that, since there is not perfect programming
 language, such comparisons are useful only for reinforcing
 a previously held viewpoint. 

>When I asked him about Ada, he replied "Oh, does Ada have
>object oriented features, I didn't know that".

 We need to do more to get the word out on the new Ada standard.
 Notice that few computer conferences have any participants from
 the Ada community.  An exception to this is the annual TOOLS USA
 conference which has been accepting Ada tutorials and speakers for
 several years now - along with other languages.  

>So of the people asking questions about Ada (without
>specifying which version), the overwhelming majority are
>talking about Ada 95, and it is a definite disservice to
>tell people that Ada 83 is still in wide use and that they
>must worry about the Ada 83 answer as well as the Ada 95
>answer. 

 Ada 83 is still in wide use.  It may be in wider use than Ada 95.
 Those who are using it also have questions.  Most of those who use
 Ada 83 do not frequent this list.  They are in the trenches grinding
 out code.  Many are new to Ada.  When they do find this forum and 
 need to get the answer to a question, they could be unaware of the
 existence of a new standard (as was you C++ colleague).  All they
 want is to get a question answered. If they fail to specify the 
 version of Ada, it is because they are still tilling the fields 
 originally sown with MIL-STD 1815, the only Ada they know.

 Sometimes, because so many of us are active in the current version
 of the language, it is hard to remember that so many are still 
 struggling with the issues of the previous standard.  

 It will not cause any permanent psychological damage for us to be
 charitable toward those gentle jungle folk seeking help from
 an unrehabilitated crocodile.

This would be like asking a question on the C++
>group, and having responders worry that you might really be
>asking about C!

 HmmMMMMmmmmm. Does C++ have templates? Depends of what version
 of C++ you are using.  Does C++ have something like an Ada 
 package?  It depends on which version of C++ you are using and
 how you understand "namespace."   

 I agree with most of your observations, Robert.  Ada 95 is Ada.
 Not everyone programming in Ada 83 knows this.

 Richard Riehle
 richard@adaworks.com
 http://www.adaworks.com

Re: Ada 83 - Sometimes still chosen

[robert_dewar]

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #1: Type: text/plain, Size: 393 bytes --]

In article <7dg5f0$fia1@ftp.kvaerner.com>,
  "Tarjei Tj�stheim Jensen" <tarjei.jensen@kvaerner.no>
wrote:
> I think it would be more apropriate to use ANSI C versus
> K&R C as an example.

No, that's too narrow a gap to be comparable ...


-----------== Posted via Deja News, The Discussion Network ==----------
http://www.dejanews.com/       Search, Read, Discuss, or Start Your Own    

Re: Calculating SQRT in ADA

[robert_dewar]

In article <umqzp50qstq.fsf@maestro.clustra.com>,
  Ole-Hjalmar Kristensen <ohk@maestro.clustra.com> wrote:
> Your faith in library routines is touching :-)

I have *far* more faith in library routines than in math
routines cooked up by programmers who do not know how to
compute elementary functions!

I must say that in a *lot* of experience, I have never seen
a library sqrt function that was bad. As I said, this is a
particularly easy one (that does not mean it cannot be
messed up by people whose knowledge of numerical algorithms
is minimal :-)

-----------== Posted via Deja News, The Discussion Network ==----------
http://www.dejanews.com/       Search, Read, Discuss, or Start Your Own    

Re: Ada 83 - Sometimes still chosen

[Matthew Heaney]

Richard D Riehle <laoXhai@ix.netcom.com> writes:

>  We need to do more to get the word out on the new Ada standard.

This is what I've been trying to do with the ACM patterns list.

<mailto:patterns@acm.org>
<http://www.acm.org/archives/patterns.html>

I have converted every C++ example in the GoF book to Ada95.  Hopefully
this will lower the entry barrier for those Ada programmers just
learning pattern technology, and allow programmers in other languages to
see how it's done in Ada.

You can subscribe to the patterns list by sending the message (body) 


subscribe patterns <your full name>


to the ACM mailing-list server.

<mailto:listserv@acm.org>


Matt

Re: Ada 83 - Sometimes still chosen

[Tarjei Tj�stheim Jensen]

robert_dewar@my-dejanews.com wrote in message
<7dhje8$26s$1@nnrp1.dejanews.com>...
>In article <7dg5f0$fia1@ftp.kvaerner.com>,
>  "Tarjei Tj stheim Jensen" <tarjei.jensen@kvaerner.no>
>wrote:
>> I think it would be more apropriate to use ANSI C versus
>> K&R C as an example.
>
>No, that's too narrow a gap to be comparable ...


That may very well be so. But the C++ / C example is flawed as both languages
officially lead their own lives (possibly aware of each other). C++ does not
obsolete C.  The example would be apropriate if Ada 83 and Ada 95 were on
different standards tracks.

Greetings,

Re: Ada 83 - Sometimes still chosen

[Marin David Condic]

Tom Moran wrote:
> 
> >If they fail to specify the
> > version of Ada...
> If an elementary school student asks you "is the moon round" the
> appropriate answer is "yes".  But if a spacecraft designer asks you

The critical point is that it absolutely never pays to be rude to those
who have not been rude to you. Some poor soul posts an innocent enough
question to this newsgroup and gets flogged with "RTFM!" or "Why didn't
you tell us it was Ada83, you moron!" and he's really going to want to
come back here and as another question, right? ;-)

There are more than enough people out there who are ignorant of The
Great Ada83/Ada95 Debate. Ignorance is not a crime. Ignorance is not
stupidity. Ignorance is not shameful. Ignorance is just ignorance and we
should *politely* try to enlighten those who suffer from the condition.

MDC
-- 
Marin David Condic
Real Time & Embedded Systems, Propulsion Systems Analysis
United Technologies, Pratt & Whitney, Large Military Engines
M/S 731-95, P.O.B. 109600, West Palm Beach, FL, 33410-9600
***To reply, remove "bogon" from the domain name.***

Re: Ada 83 - Sometimes still chosen

[Tarjei Tj�stheim Jensen]

Marin David Condic wrote:
 >There are more than enough people out there who are ignorant of The
>Great Ada83/Ada95 Debate. Ignorance is not a crime. Ignorance is not
>stupidity. Ignorance is not shameful. Ignorance is just ignorance and we
>should *politely* try to enlighten those who suffer from the condition.


I second, third and fourth that.



Greegings,

Re: Calculating SQRT in ADA

[Robert I. Eachus]

robert_dewar@my-dejanews.com wrote:

> I have *far* more faith in library routines than in math
> routines cooked up by programmers who do not know how to
> compute elementary functions!
> 
> I must say that in a *lot* of experience, I have never seen
> a library sqrt function that was bad. As I said, this is a
> particularly easy one (that does not mean it cannot be
> messed up by people whose knowledge of numerical algorithms
> is minimal :-)
 
   Robert, your memory is going as you are getting older. ;-)  In the
sixties and seventies, there were some extremely atrocious
implementations in language run-time libraries.  The worst horror was
the URAND for the IBM 360 which caused dozens of peer-reviewed papers to
be withdrawn.  But I also remember many trigonometric libraries which
used polynomical approximations that were good to only six decimal
places--even for double precision.  Notice that even one N-R iteration
would have fixed that problem.  (The justification! was that the same
library function was used for both single and double precision, but most
single precision users wanted speed.  Of course, in many cases I was
able to write faster library functions that were much more accurate for
the double case, and gave LSB results for the single precision case.

   Of course, as Robert Dewar knows, the IEEE floating-point arithmetic
project did a lot to fix the situation, and to get trig, log, and
exponential functions officially in the hardware.

   Nowadays it is very hard to improve on the "hardware" IEEE
instructions, even when they involve some software emulation.  Back to
the original topic, the best way to compute square roots is not to use
NR, but to use the same algorithm you learned in school (with
calculating the exponent correctly for floating-point values).  On most
RISC hardware it is slightly faster than one floating-point divide, so
Newton-Rhapson can't compete.  (Of course, the hardware instruction, if
it is really done in hardware, is faster still.)

Re: Calculating SQRT in ADA

[robert_dewar]

In article <36FFF83A.BE789C93@mitre.org>,
  "Robert I. Eachus" <eachus@mitre.org> wrote:
>    Robert, your memory is going as you are getting older.
> ;-)  In the
> sixties and seventies, there were some extremely
> atrocious
> implementations in language run-time libraries.  The
> worst horror was
> the URAND for the IBM 360 which caused dozens of
> peer-reviewed papers to
> be withdrawn.  But I also remember many trigonometric
> libraries which
> used polynomical approximations that were good to only
> six decimal
> places--even for double precision.

It would be nice to see chapter and verse cited for this
(perhaps this horrible behavioor comes from the Honeywell
compiler that you always like to cite even though it never
saw the light of day :-)

In fact the IBM trigonometric routines were excellent on
both the 7094 and the 360. Indeed it is from Hirondo Kuki
(I hope I remember the spelling right, this is from a long
time ago on the 7094) that I learned about how math library
routines were written.

I will say again, I *far* better trust library routines
written by someone who knows what they are doing (which
ever since 7094 days tends to be the case, given the high
standards set by HK), than routines written by people who
just don't know math stuff.

And please note, this discussion is about elementary
functions, your reference to URAND, though it reflects
your own particular interest in random number generation,
is not germane to the thread on square root!

-----------== Posted via Deja News, The Discussion Network ==----------
http://www.dejanews.com/       Search, Read, Discuss, or Start Your Own    

Re: Calculating SQRT in ADA

[robert_dewar]

In article <36FFF83A.BE789C93@mitre.org>,
  "Robert I. Eachus" <eachus@mitre.org> wrote:
>    Nowadays it is very hard to improve on the "hardware"
> IEEE
> instructions,

I am not sure what this refers to, what "hardware" IEEE
instructions are you referring to. Certainly IEEE does not
include elementary functions except for sqrt, and this is
of course NOT hardware on most machines.

> even when they involve some software
> emulation.

"some software emulation" = done by a software library
routine for most architectures for sqrt, and there are no
other operations defined in IEEE.

  Back to
> the original topic, the best way to compute square roots
> is not to use
> NR, but to use the same algorithm you learned in school
> (with calculating the exponent correctly for
> floating-point values).  On most
> RISC hardware it is slightly faster than one
> floating-point divide, so
> Newton-Rhapson can't compete.  (Of course, the hardware
> instruction, if
> it is really done in hardware, is faster still.)

Sure, a sqrt in hardware can be as fast as a divide, since
a very similar algorithm can be used. But I challenge your
initial statement here. Please cough up code on a specific
machine to justify the statement that you can do a sqrt in
floating-point divide time.

-----------== Posted via Deja News, The Discussion Network ==----------
http://www.dejanews.com/       Search, Read, Discuss, or Start Your Own    

Re: Calculating SQRT in ADA

[Robert I. Eachus]

robert_dewar@my-dejanews.com wrote:

> I am not sure what this refers to, what "hardware" IEEE
> instructions are you referring to. Certainly IEEE does not
> include elementary functions except for sqrt, and this is
> of course NOT hardware on most machines.

   News to me.  A lot of the current processor architectures emulate
some of the trigonmetric and trancendental functions with microcode or
hardware traps to library routines, but there are still available
members of these families that implement the instructions in hardware. 
For example, in the 68000 family, the earlier processor families had all
the instructions in hardware in the 68881 and 68882 coprocessor chips. 
The 68040 implemented many floating point instructions in hardware and
emulated others, in the 68060, almost all of the instructions other than
the basic floating point operations are done in emulation libraries.

> Sure, a sqrt in hardware can be as fast as a divide, since
> a very similar algorithm can be used. But I challenge your
> initial statement here. Please cough up code on a specific
> machine to justify the statement that you can do a sqrt in
> floating-point divide time.

  Well the first manual I grabbed off the shelf surprised me slightly:
on the 68881, the floating point sqare root took two cycles more than a
divide--out of about 130. (The exact number of cycles depends on
register modes.)  Of course the cost of loading the second operand for
divide takes longer than two clocks--four to 40 depending on source and
memory speed.

  The FSQRT has been part of the SPARC architecture since version 7, it
is implemented in hardware on almost all chipsets.   I don't have timing
tables handy, but I have tested several SPARC processors where FSQRT is
faster than FDIV.  As above, the speed advantage comes from only having
one operand more than anything.  Loading FP registers, especially from
memory, costs.  (Of course, YMMV, but I was more concerned with cases
where I was calculating for a large set of points, so the high speed
caches didn't much effect the data loading.)

  For the integer case, a div.l takes about 90 clocks on a 68020, while
the corresponding square root algorithm takes 16 iterations through a
loop:

     L0:   MOVE.L (operand),D5;
           TRAPMI                       ; Error if operand is negative
           MOVEQ #15,D4;
           MOVEQ #0,D6;
           MOVEQ #1,D7;
      L1:  ROR #1,D6                    ; First rotation has no effect.
           ROR #2,D7                    ; First rotation results
in                                                                                                                                                                                                                                                                                                                                
; #8000000
           CMP.L D6,D5;
           BLT L2;
           ADD D7,D6;
           SUB D6,D5;
      L2:  DBF D4,L1;

     I'm not sure I have this correct from memory, but it is close.  The
version I used unwound the loop, used a 64-bit operand, and did a BFFFO
to skip leading zeros.  I needed to do SQRT(X*X+Y*Y) fast, again for
lots of points.

Re: Calculating SQRT in ADA

[Robert I. Eachus]

robert_dewar@my-dejanews.com wrote:

> It would be nice to see chapter and verse cited for this
> (perhaps this horrible behavioor comes from the Honeywell
> compiler that you always like to cite even though it never
> saw the light of day :-)

   The Ada/SIL compiler was both used internally at Honeywell and for
training military personnel.  But you are right that it was never sold
commercially.
 
> In fact the IBM trigonometric routines were excellent on
> both the 7094 and the 360. Indeed it is from Hirondo Kuki
> (I hope I remember the spelling right, this is from a long
> time ago on the 7094) that I learned about how math library
> routines were written.

   Yes, the math libraries on the 7094 were excellent.  The original
libraries on the 360 series were not, but IBM improved them rapidly. 
However, for the Apollo program, NASA continued to use 4 7094s that were
supposed to have been replaced by 360's:  IBM could get the math right,
or the necessary speed but not both.  (To be fair, it was both the
shorter word lenght and the trucaction semantics on the 360 that made it
unsuitable for this application.  IBM did finally come up with the
360/91 that could make the grade.)

   However, the first IBM PCs had an abominable math library for BASIC. 
(Can you say Microsoft kiddies?  I knew you could.)  On the other hand,
on the 68000 family, Microsoft BASIC used the excellent MOtorola
libraries.

> And please note, this discussion is about elementary
> functions, your reference to URAND, though it reflects
> your own particular interest in random number generation,
> is not germane to the thread on square root!

   Remember the Savage benchmark?  It is not surprising to find no
cumulative error with it nowadays.  But even ten years ago, it wasn't
that unusual to be much more concerned about the error introduced than
about the speed performance.

Re: Calculating SQRT in ADA

[robert_dewar]

In article <3705555E.5572A782@mitre.org>,
  "Robert I. Eachus" <eachus@mitre.org> wrote:
>   For the integer case, a div.l takes about 90 clocks on
> a 68020

We are NOT talking about the integer case!
We are NOT talking about obsolete CISC machines!

We are talking about modern RISC processors with fast
floating-point units. Look at the start of this thread,
it is clear that the focus is floating-point.

Now, backup your original statement that NR is not the
fastest approach in this environment.

-----------== Posted via Deja News, The Discussion Network ==----------
http://www.dejanews.com/       Search, Read, Discuss, or Start Your Own