Discussion:
what is the type of the result of a 128 bit integer square root?
(too old to reply)
jacob navia
2015-05-29 07:32:23 UTC
Permalink
Raw Message
The lcc-win compiler system features 128 bit integers. Now, the question
arises of the type of the square root of such numbers.

If lcc-win uses the highest precision available in hardware (what was
doing till last week) it uses then an 80 bit floatiung point number.

The square roots will lose precision!

A customer complained about that last week. I have developed then an
exact INTEGER square root that gives all 128 bit precision.

OK, but now the type of the result is a 128 bit integer, not a floating
point number.

Is this correct?

Since C99 and tgmath.h we can decide to what function the expression

sqrt(x)

expands to, according to the type of the argument. It is possible then
to expand the square root of an integer of 128 bits into a 128 integer
data type without any change to anything in C99.

lcc-win features a software 128 bit floating point. But even those
numbers would not have enough precision to keep the 128 bit precision of
a 128 bit integer. I would have to use the very extended qfloat format
that has 512 bits precision, to be able to match the 128 bit integer
precision. But that format is obviously very slow, by more than two
orders of magnitude!

I am doing now a sqrt(int128) in about 35 ns using several hardware
tricks. Using qfloat would mean multiplying that by around 1000.

BUT

If I do this for sqrt (keeping all precision) I should also modify some
other functions like log, for instance or maybe many others, I haven't
really investigated that part...

What should I do?

Note that since 128 bit integers are not in really widespread use, it
would be nice if the committee decided this quickly since there is now
no need to care about existing code...

Later it will be, well... too late

:-)
Francis Glassborow
2015-05-29 10:04:42 UTC
Permalink
Raw Message
Post by jacob navia
The lcc-win compiler system features 128 bit integers. Now, the question
arises of the type of the square root of such numbers.
If lcc-win uses the highest precision available in hardware (what was
doing till last week) it uses then an 80 bit floatiung point number.
The square roots will lose precision!
I think you need to address what the user really needs. 128-bit integers
are nice for some purposes (and unlimited precision are even nicer for
most of those). However most if not all these are purely interested in
integer arithmetic. In this context we have to ask what is meant by a
square root. In terms of integer only problems square roots have a very
limited value.

Note that you can use conversion from 128-bit integer to 80-bit floating
point to compute a first approximation of a square root and then a
simple integer algorithm (divide by candidate square root and take
median of divisor and quotient, repeat till stable) quickly converges to
the desired result)

So the real question is what extended floating point types are required.
Note that, in general, floating point computational problems are
satisfied by much more limited precision. The exception being in some
computational intense applications where extended guard digits are
needed to deal with the progressive degradation caused by rounding errors.

So what does your customer/client require? And is he happy to pay you to
provide something that has limited general utility.

Francis
m***@yahoo.co.uk
2015-05-29 11:22:13 UTC
Permalink
Raw Message
Post by jacob navia
The lcc-win compiler system features 128 bit integers. Now, the question
arises of the type of the square root of such numbers.
If lcc-win uses the highest precision available in hardware (what was
doing till last week) it uses then an 80 bit floating point number.
The square roots will lose precision!
Will it?
Post by jacob navia
A customer complained about that last week. I have developed then an
exact INTEGER square root that gives all 128 bit precision.
A 128-bit integer will have a square root that fits in a 64 bit integer (if the answer is integral). An 80-bit intel float has 64 bits of fraction, so will return *at least* as precise a result as your integer version. (The answer will be more precise when they start to have a fractional part).

I concur with Francis' advice to decide what it is you need, and then implement that. Note that the iteration he suggested:
estimate <- (estimate + square/estimate)/2
is VERY efficient. Once you are in the right ball park, it will double
the number of correct bits every iteration. That means that if you get an 80-bit result in hardware (with 64 bits of precision), a single run round
the loop in a slow software format will get you 128 bits of precision (if
your software format can manage that).
jacob navia
2015-05-29 16:06:38 UTC
Permalink
Raw Message
Post by m***@yahoo.co.uk
A 128-bit integer will have a square root that fits in a 64 bit integer (if the answer is integral).
An 80-bit intel float has 64 bits of fraction, so will return*at least*
as precise a result as your integer version.

(The answer will be more precise when they start to have a fractional part).

No.

Flloating point numbers are NOT evenly spaced in the real line. They are
more finely concentrated along the real axis in the neighbourhood of -1
to 1. This is basic floating point stuff...

For numbers around 10^25 or bigger the spacing makes inevitable rounding
errors.

Integers however, are spaced evenly along the real line...


Example:
sqrt(97656270007716649790232399025) yields
312500031982830, but the correct exact answer is
312500032012345.

The basic question I am asking is:

When the integer precision is bigger than the maximum floating point
precision, shouldn't the system give an integer answer?

Thanks
m***@yahoo.co.uk
2015-05-29 16:29:56 UTC
Permalink
Raw Message
Post by jacob navia
Post by m***@yahoo.co.uk
A 128-bit integer will have a square root that fits in a 64 bit integer (if the answer is integral).
An 80-bit intel float has 64 bits of fraction, so will return*at least*
as precise a result as your integer version.
(The answer will be more precise when they start to have a fractional part).
No.
I think you are mistaken.
Post by jacob navia
Floating point numbers are NOT evenly spaced in the real line. They are
more finely concentrated along the real axis in the neighbourhood of -1
to 1. This is basic floating point stuff...
Oh indeed. My point is that an 80 bit floating point number can represent
every 64 bit integer precisely. It can also represent fractional values for smaller numbers.
Post by jacob navia
For numbers around 10^25 or bigger the spacing makes inevitable rounding
errors.
Yes. (Although I make it about 2*10^19.)
Post by jacob navia
Integers however, are spaced evenly along the real line...
sqrt(97656270007716649790232399025) yields
312500031982830, but the correct exact answer is
312500032012345.
I suspect that the problem is that sqrt(97656270007716649790232399025)
is actually calling (in a C++ like syntax):
sqrt<long double>((long double)97656270007716649790232399025)

and you are losing the precision when the argument is cast to 80 bit
floating point. I was assuming that you claimed the function was behaving
as if prototyped

long double sqrt(long long arg)
Post by jacob navia
When the integer precision is bigger than the maximum floating point
precision, shouldn't the system give an integer answer?
That also means that sqrt(2LL) will return 1. My proposed signature
would return 1.414..., and could also return 312500032012345.0 for your
test case.


If you know that your integers are always large, or you know they are
perfect squares, it is easy enough to write an efficient function which
returns the value.
long long isqrt(long long arg) {
long long result = sqrt(arg)
long long delta;
while ((delta = (arg/result - result) / 2))
result += delta;
return result;
}
... should be close
Francis Glassborow
2015-05-29 16:33:27 UTC
Permalink
Raw Message
Post by jacob navia
Post by m***@yahoo.co.uk
A 128-bit integer will have a square root that fits in a 64 bit
integer (if the answer is integral).
An 80-bit intel float has 64 bits of fraction, so will return*at least*
as precise a result as your integer version.
(The answer will be more precise when they start to have a fractional part).
No.
Flloating point numbers are NOT evenly spaced in the real line. They are
more finely concentrated along the real axis in the neighbourhood of -1
to 1. This is basic floating point stuff...
For numbers around 10^25 or bigger the spacing makes inevitable rounding
errors.
Integers however, are spaced evenly along the real line...
sqrt(97656270007716649790232399025) yields
312500031982830, but the correct exact answer is
312500032012345.
Which reflect errors introduced by the mechanism for computing square
roots in floatinjg point (often via logarithms. But now that answer can
be refined using the algorithm you have been given.
Post by jacob navia
When the integer precision is bigger than the maximum floating point
precision, shouldn't the system give an integer answer?
Thanks
Come on Jacob, think about it.

Precision is roughly the number of value bits. The fact that floating
point is represented by a bit pattern that is notionally considered in
the range -1 to 1 + a multiplier is irrelevant. For example 65536 can be
competely and precisely represented by both a long and a long double

So what are you proposing is the square root of 8? 2, 3 or 2.82.....
The motive lies with the programmer. Until we know why she wants a
square root we do not know what we should do and this has nothing to do
with the available precision of built-in types.

Francis
jacob navia
2015-05-29 17:46:57 UTC
Permalink
Raw Message
Post by Francis Glassborow
So what are you proposing is the square root of 8? 2, 3 or 2.82.....
This depends on the current rounding mode.

lcc-win starts with round to nearest but this can be changed by the
user at any time.

The rounding mode is well defined in C99.
Keith Thompson
2015-05-29 18:19:01 UTC
Permalink
Raw Message
Post by jacob navia
Post by Francis Glassborow
So what are you proposing is the square root of 8? 2, 3 or 2.82.....
This depends on the current rounding mode.
lcc-win starts with round to nearest but this can be changed by the
user at any time.
The rounding mode is well defined in C99.
I think that *integer* square root is usually defined to yield
the integer part of the real result, truncated toward zero, so
sqrti(8) == 2, just as 8/3 yields 2, not 3. Rounding mode generally
applies to floating-point operations, not integer operations.

See, for example, https://en.wikipedia.org/wiki/Integer_square_root
--
Keith Thompson (The_Other_Keith) kst-***@mib.org <http://www.ghoti.net/~kst>
Working, but not speaking, for JetHead Development, Inc.
"We must do something. This is something. Therefore, we must do this."
-- Antony Jay and Jonathan Lynn, "Yes Minister"
jacob navia
2015-05-29 18:34:22 UTC
Permalink
Raw Message
Post by Keith Thompson
Post by jacob navia
Post by Francis Glassborow
So what are you proposing is the square root of 8? 2, 3 or 2.82.....
This depends on the current rounding mode.
lcc-win starts with round to nearest but this can be changed by the
user at any time.
The rounding mode is well defined in C99.
I think that *integer* square root is usually defined to yield
the integer part of the real result, truncated toward zero, so
sqrti(8) == 2, just as 8/3 yields 2, not 3. Rounding mode generally
applies to floating-point operations, not integer operations.
See, for example, https://en.wikipedia.org/wiki/Integer_square_root
Yes, you are right. That was a mistake

Kaz Kylheku
2015-05-29 12:59:55 UTC
Permalink
Raw Message
Post by jacob navia
The lcc-win compiler system features 128 bit integers. Now, the question
arises of the type of the square root of such numbers.
If lcc-win uses the highest precision available in hardware (what was
doing till last week) it uses then an 80 bit floatiung point number.
The square roots will lose precision!
A customer complained about that last week. I have developed then an
exact INTEGER square root that gives all 128 bit precision.
OK, but now the type of the result is a 128 bit integer, not a floating
point number.
Is this correct?
Correct with respect to what specification? If someone wants the square
root of 8 to produce 2.82... then an integer square root function isn't
correct for that requirement.
Post by jacob navia
Since C99 and tgmath.h we can decide to what function the expression
sqrt(x)
expands to, according to the type of the argument.
sqrt is expected to be floating-point, even if passed integer arguments. An
integer square root is something else, and should be given its own name, like
isqrt.

E.g. look at the ANSI Lisp language:

clisp -q
[1]> (isqrt 25)
5
[2]> (isqrt 8)
2
[3]> (sqrt 25)
5
[4]> (sqrt 8)
2.828427

Of course, here we also get integer type results when the square root is
exact, regardless of which we use: dynamic typing!

The point though is that a programmer might want either operation. Sometimes
the requirement might be 2, and sometimes 2.82.
m***@yahoo.co.uk
2015-05-29 14:16:30 UTC
Permalink
Raw Message
Post by Kaz Kylheku
clisp -q
[1]> (isqrt 25)
5
[2]> (isqrt 8)
2
[3]> (sqrt 25)
5
[4]> (sqrt 8)
2.828427
Of course, here we also get integer type results when the square root is
exact, regardless of which we use: dynamic typing!
The point though is that a programmer might want either operation.
Sometimes the requirement might be 2, and sometimes 2.82.
Interesting. Of course there are three possible definitions of isqrt.
Lisp is obviously using floor(sqrt(x)) {using C syntax}. One might
also want round(sqrt(x)) or ceil(sqrt(x))
Keith Thompson
2015-05-29 16:17:04 UTC
Permalink
Raw Message
Post by jacob navia
The lcc-win compiler system features 128 bit integers. Now, the
question arises of the type of the square root of such numbers.
If lcc-win uses the highest precision available in hardware (what was
doing till last week) it uses then an 80 bit floatiung point number.
The square roots will lose precision!
A customer complained about that last week. I have developed then an
exact INTEGER square root that gives all 128 bit precision.
OK, but now the type of the result is a 128 bit integer, not a
floating point number.
Is this correct?
The C standard library does not include an integer square root function.
If you want to provide such a function in a non-standard library, you
can define it any way you like; the standard has nothing to say about
it.

All the existing floating-point square root functions (sqrtf, sqrt,
sqrtl, and the sqrt macro in <tgmath.h>) return results of the same type
as the argument. Based on that, it probably makes sense for a 128-bit
integer square root function to return a 128-bit integer result. In
fact, I'd be surprised if it did anything else.
Post by jacob navia
Since C99 and tgmath.h we can decide to what function the expression
sqrt(x)
expands to, according to the type of the argument. It is possible then
to expand the square root of an integer of 128 bits into a 128 integer
data type without any change to anything in C99.
I wouldn't suggest expanding <tgmath.h>'s sqrt to operate on integers.
Instead, I'd suggest defining your own function(s) and/or macro(s) in an
implementation-specific header.

[...]
Post by jacob navia
Note that since 128 bit integers are not in really widespread use, it
would be nice if the committee decided this quickly since there is now
no need to care about existing code...
Later it will be, well... too late
The current standard allows for 128-bit integers, either by
making long long 128 bits or by providing extended integer types
and making [u]intmax_t 128 bits. I don't think that *requiring*
128-bit integer support (as C99 required 64-bit integer support)
is feasible; there are still systems where it would be difficult
and/or would result in very inefficient code. It would probably
make sense to require an implementation to define [u]int128_t *if
it can*, as it currently does for [u]int{8,16,32,64}_t.
--
Keith Thompson (The_Other_Keith) kst-***@mib.org <http://www.ghoti.net/~kst>
Working, but not speaking, for JetHead Development, Inc.
"We must do something. This is something. Therefore, we must do this."
-- Antony Jay and Jonathan Lynn, "Yes Minister"
Loading...