Multi precision floating point

[jacob navia]

I was specifically referring to jacob's qfloat precision which he claims
is 352-bits. I presume this is mantissa bits, including the hidden bit,
but I'm not sure. I doubt there are many compilers which provide long
double of this size.

Phil

This is mantissa bits. Gives around 107 decimal digits...

 

[CBFalconer]

(e-mail address removed) wrote:
.... snip ...


I was specifically referring to jacob's qfloat precision which he
claims is 352-bits. I presume this is mantissa bits, including the
hidden bit, but I'm not sure. I doubt there are many compilers
which provide long double of this size.

Those are off topic here, where we deal only with standard C. Try
comp.programming and/or comp.compilers.lcc.

 

[Philip Potter]

Those are off topic here, where we deal only with standard C. Try
comp.programming and/or comp.compilers.lcc.

Topicality drifts, and deep-thread topicality violations are not as
damaging as offtopic new threads.

Besides, if you read the entire thread you will learn that I was
advocating the OP use Standard C.

 

[christian.bau]

Does there exist some C compiler supporting quadruple precision?
The wikipedia page says that "long double" can be the 80-bits of
the x86 or sometimes the 128 bits of quadruple precision. What are
those "remarkably few" compilers? Also, how could we interface
fortran and C for quadruple precision?

Buy a used Macintosh with a G4 or G5 processor. The gcc compiler that
ships with it supports 128 bit long double. Each long double operation
is implemented using about ten ordinary floating point operations on
the average.

 

[Flash Gordon]

santosh wrote, On 29/11/07 11:05:

(e-mail address removed) wrote:


jacob was part way through implementing his lcc-win32 for Linux (x86)
when licensing restrictions forced him to freeze the effort. But he may
be able to still help you out if you can establish a commercial license
with him.

He can't do that until he has sorted out the licensing problems. Note
that this is NOT a complaint about Jacob since he did not suggest this.

 

[rosewater]

santosh wrote, On 29/11/07 11:05:

He can't do that until he has sorted out the licensing problems. Note
that this is NOT a complaint about Jacob since he did not suggest this.

The root cause of Navia's licensing problems is his own selfishness.
If he released his compiler under a free license like the GPL, then
his Linux version would be releasable.

Of course, there's still the question of why anyone would want to use
a compiler maintained by someone with no interest in C, but only his C+
+-like extensions, and who exposes his abysmal lack of knowledge about
the language every day in this group.

 

[jacob navia]

santosh wrote, On 29/11/07 11:05:

He can't do that until he has sorted out the licensing problems. Note
that this is NOT a complaint about Jacob since he did not suggest this.

Yeah I have been working again trying to use the system assembler with

1) writing the assembler output into /tmp/something66677.s
2) calling system("as /tmp/something66677.s");

But now there are some problems with the output as it seems,
and the compiler doesn't compile itself any more.

Is it because the assembler is different? It looks like.
Or is it because since I was assuming the lcc assembler behavior
in the compiler and that doesn't work any more?

What a waste of time...

 

[Kenny McCormack]

....

The root cause of Navia's licensing problems is his own selfishness.
If he released his compiler under a free license like the GPL, then
his Linux version would be releasable.

Of course, there's still the question of why anyone would want to use
a compiler maintained by someone with no interest in C, but only his C+
+-like extensions, and who exposes his abysmal lack of knowledge about
the language every day in this group.

I'd like to introduce a new metric - called the Heathfield-Navia delta -
that is defined as the number of responses between when Navia's name is
mentioned, in a CLC article or response, and when Navia gets slammed by
Heathfield (either directly or via a sock puppet).

 

[user923005]

Yes that is somewhat more what I had in mind.
Actually, the intel fortran compiler version 10 has support for
quadruple precision of IEEE754r.

Does there exist some C compiler supporting quadruple precision?

The QD package found here:
http://crd.lbl.gov/~dhbailey/mpdist/
has 32 and 64 digit floating point. It's C++ and not C, though.

The wikipedia page says that "long double" can be the 80-bits of
the x86 or sometimes the 128 bits of quadruple precision. What are
those "remarkably few" compilers? Also, how could we interface
fortran and C for quadruple precision?

The QD package contains an interface for Fortran. Generally one would
expect a DLL under windows and a shared library under UNIX variants to
provide an interface to a special purpose library.

That's the extend of my problem. I have an existing C library
depending on lapack/blas, which I need to extend to computing
in higher precision, and quadruple might not be enough.

If quadruple precision is not enough, then your matrices must be
awfully close to singular. I smell trouble.

Otherwise, thanks to all and to Jacob for lcc-win, which I did
not knew about. But this work is supposed to be on linux x86
computers.

The lcc-win numeric package is just the Cephes library from Moshier.
You can use it with any other compiler as well. See the qlib link
from this site:
http://www.moshier.net/#Cephes

 

[user923005]

The root cause of Navia's licensing problems is his own selfishness.
If he released his compiler under a free license like the GPL, then
his Linux version would be releasable.

I guess that this is not possible, since his compiler was derived from
LCC, which does not have a license compatible with GPL. At any rate,
whatever licence he chooses should be his own choice -- within
whatever limits are possible.

Of course, there's still the question of why anyone would want to use
a compiler maintained by someone with no interest in C, but only his C+
+-like extensions, and who exposes his abysmal lack of knowledge about
the language every day in this group.

I think is qfloat extension is quite interesing. I would like it if
every C compiler had a Cephes compatible interface for all functions
and for all of the float sizes in the Cephes lib, including qfloat.
As far as Jacob's C knowledge goes, he is clearly above average. He's
a bit combative at times (which may be reactionary) and I think that
sometimes he deliberately pretends not to understand something, but
everyone who posts here has their own little quirks (myself definitely
included).

IMO-YMMV.

 

[Keith Thompson]

On Nov 29, 1:23 pm, (e-mail address removed) wrote: [...]

The root cause of Navia's licensing problems [...]

I guess that this is not possible, since [...]

Please don't feed the troll. (e-mail address removed) has done
nothing but flame jacob navia; whatever you might think of jacob, this
rosewater person's trolling is not worth a response.

 

[mathieu.dutour]

Does there exist some C compiler supporting quadruple precision?

Have you investigated the excellent Intel C/C++ compilers? They should
provide 128 bit floating point support.

<http://www.intel.com/cd/software/products/asmo-na/eng/compilers/28413...>

Yes, after looking at the forum there is a compiler option
for doing just that, like in fortran.
We will experiment with that.

Then you need something like MPFR.

<http://en.wikipedia.org/wiki/MPFR>
<http://www.mpfr.org/>

The big problem of MPFR is that it is slow even if it faster than
the competition. Otherwise, yes it is very flexible and well designed.

My question is now off topic, but how hard is it do write explicitly
an octuple precision floating point arithmetic (+, x, /) in C.
For integers, this has been done upon the "long long" type.

 

[jacob navia]

The lcc-win numeric package is just the Cephes library from Moshier.
You can use it with any other compiler as well. See the qlib link
from this site:
http://www.moshier.net/#Cephes

This is not exactly true.
1) I rewrote all the basic operations in assembler. Speed is 50%
increased compared to the original versions.

2) I rewrote all the C99 library, adapting all functions from the
original library, increasing precision when possible

3) Using operator overloading, I integrated this library
like a basic type into the language.

 

[Nick Keighley]

(e-mail address removed) wrote:

If you can, think about using tested existing packages instead of
writing yet another implementation. Something like GMP is a good
candidate

note: GMP uses the full GPL (or did last time I looked) this may
limit where you can use it.

<duck>

<snip>

--
Nick Keighley

Quantum Boggum Sort:
Q1. use a source of quantum noise (eg. radioactive decay) to
randomly permutate an array.
Q2. if the array is not ordered, destroy the universe (*)
Q3. if you reached this step your universe has sorted the array
in O(n) time.
(*) [100] this is left as an exercise

 

[Keith Thompson]

This is not exactly true.
1) I rewrote all the basic operations in assembler. Speed is 50%
increased compared to the original versions.

2) I rewrote all the C99 library, adapting all functions from the
original library, increasing precision when possible

3) Using operator overloading, I integrated this library
like a basic type into the language.

Where the phrase "the language" refers to something other than C.

 

[jacob navia]

Where the phrase "the language" refers to something other than C.

No. All this extensions are compatible with the C standard as
I have explained thousand times in this newsgroup.

Note that the C standard does NOT forbid extensions

 

[Flash Gordon]

jacob navia wrote, On 30/11/07 21:32:

No. All this extensions are compatible with the C standard as
I have explained thousand times in this newsgroup.

Note that the C standard does NOT forbid extensions

Well, your compiler does not conform to the C standard *without* the
-ansic flag (for example declaring atof in math.h), so does it support
operator overloading *with* the -ansic flag? If not then it compiles C
when the -ansic flag is specified and something other than C that
supports operator overloading with it.

 

[user923005]

No. All this extensions are compatible with the C standard as
I have explained thousand times in this newsgroup.

Note that the C standard does NOT forbid extensions

By the same logic, this is C code, since GCC compiles it:

PROGRAM STEFACC
C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994
C To accompany the text:
C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd
Ed, 1992
C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
C This free software is complements of the author.
C
C Algorithm 2.7 (Steffensen's Acceleration).
C Section 2.5, Aitken's Process & Steffensen's & Muller's Methods,
Page 96
C
PARAMETER(Delta=5E-6,Epsilon=5E-6,Max=99)
INTEGER Cond,K
REAL Dp,P0,P3
CHARACTER ANS*1
EXTERNAL F,F1
10 CALL INPUT(P0)
CALL STEFFEN(F,F1,P0,Delta,Epsilon,Max,P3,Dp,Cond,K)
CALL RESULTS(P0,P3,Dp,Cond,K)
WRITE(9,*)' '
WRITE(9,*)'WANT TO TRY ANOTHER STARTING VALUE ? <Y/N> '
READ(9,'(A)') ANS
IF (ANS.EQ.'Y' .OR. ANS.EQ.'y') GOTO 10
STOP
END

REAL FUNCTION F(X)
REAL X
F=X*X*X-3*X+2
RETURN
END

REAL FUNCTION F1(X)
REAL X
F1=3*X*X-3
RETURN
END

SUBROUTINE PRINTFUN
WRITE(9,*)'F(X) = X*X*X-3*X+2'
RETURN
END

SUBROUTINE STEFFEN(F,F1,P0,Delta,Epsilon,Max,P3,Dp,Cond,K)
PARAMETER(Small=1E-20)
INTEGER Cond,K,Max
REAL D1,D2,Delta,DF0,DF1,Dp,Epsilon,P0,P1,P2,P3,Pold,RelErr,Y3
EXTERNAL F,F1
Pold=P0
K=0
Cond=0
P3=P0
P2=P0+1
P1=P0+2
100 continue
IF ((K.LT.Max).AND.(Cond.EQ.0)) then
P0=P3
DF0=F1(P0)
IF (DF0.NE.0) THEN
P1=P0 - F(P0)/DF0
ELSE
Cond=1
Dp=P3-P2
P3=P0
ENDIF
DF1=F1(P1)
IF (DF1.EQ.0) THEN
Cond=1
Dp=P1-P0
P3=P1
ELSE
P2=P1 - F(P1)/DF1
D1=(P1-P0)*(P1-P0)
D2=P2-2*P1+P0
IF (D2.EQ.0) THEN
Cond=1
Dp=P2-P1
P3=P2
ELSE
P3=P0-D1/D2
Dp=P3-P2
ENDIF
Y3=F(P3)
RelErr=ABS(Dp)/(ABS(P3)+Small)
IF (RelErr.LT.Delta) then
Cond=2
endif
IF (ABS(Y3).LT.Epsilon) then
Cond=3
endif
IF ((Cond.EQ.2).AND.(ABS(Y3).LT.Epsilon)) then
Cond=4
endif
ENDIF
K=K+1
WRITE(9,1000) K,P3,P3
goto 100
ENDIF
P0=Pold
PAUSE
RETURN
1000 FORMAT(I2,4X,F15.7,4X,F15.7)
END

SUBROUTINE XSTEFFEN(F,F1,P0,Delta,Epsilon,Max,P3,Dp,Cond,K)
C This subroutine uses simulated WHILE loop(s).
PARAMETER(Small=1E-20)
INTEGER Cond,K,Max
REAL D1,D2,Delta,DF0,DF1,Dp,Epsilon,P0,P1,P2,P3,Pold,RelErr,Y3
EXTERNAL F,F1
Pold=P0
K=0
Cond=0
P3=P0
P2=P0+1
P1=P0+2
10 IF ((K.LT.Max).AND.(Cond.EQ.0)) THEN
P0=P3
DF0=F1(P0)
IF (DF0.NE.0) THEN
P1=P0 - F(P0)/DF0
ELSE
Cond=1
Dp=P3-P2
P3=P0
ENDIF
DF1=F1(P1)
IF (DF1.EQ.0) THEN
Cond=1
Dp=P1-P0
P3=P1
ELSE
P2=P1 - F(P1)/DF1
D1=(P1-P0)*(P1-P0)
D2=P2-2*P1+P0
IF (D2.EQ.0) THEN
Cond=1
Dp=P2-P1
P3=P2
ELSE
P3=P0-D1/D2
Dp=P3-P2
ENDIF
Y3=F(P3)
RelErr=ABS(Dp)/(ABS(P3)+Small)
IF (RelErr.LT.Delta) then
Cond=2
endif
IF (ABS(Y3).LT.Epsilon) then
Cond=3
endif
IF ((RelErr.LT.Delta).AND.(ABS(Y3).LT.Epsilon)) then
Cond=4
endif
ENDIF
K=K+1
WRITE(9,1000) K,P3,P3
GOTO 10
ENDIF
P0=Pold
PAUSE
RETURN
1000 FORMAT(I2,4X,F15.7,4X,F15.7)
END

SUBROUTINE INPUT(P0)
INTEGER I
REAL P0
DO 10 I=1,18
WRITE(9,*)' '
10 CONTINUE
WRITE(9,*)' '
WRITE(9,*)'STEFFENSEN`S ACCELERATION OF THE NEWTON-RAPHSON'
WRITE(9,*)' '
WRITE(9,*)'METHOD IS USED TO FIND A ZERO OF THE FUNCTION'
WRITE(9,*)' '
CALL PRINTFUN
WRITE(9,*)' '
WRITE(9,*)'ONE INITIAL APPROXIMATION P0 IS NEEDED.'
WRITE(9,*)' '
WRITE(9,*)'ENTER P0 = '
READ(9,*) P0
WRITE(9,*)' '
RETURN
END

SUBROUTINE RESULTS(P0,P3,Dp,Cond,K)
INTEGER Cond,I,K
REAL P0,P3,Dp
DO 10 I=1,15
WRITE(9,*)' '
10 CONTINUE
WRITE(9,*)'STEFFENSEN`S ACCELERATION OF THE NEWTON-RAPHSON'
WRITE(9,*)' '
WRITE(9,*)'METHOD WAS USED TO FIND A ZERO OF THE FUNCTION'
WRITE(9,*)' '
CALL PRINTFUN
WRITE(9,*)' '
WRITE(9,*)'STARTING WITH THE APPROXIMATION P0 =',P0
WRITE(9,*)' '
WRITE(9,*)'AFTER ',K,' ITERATIONS AN APPROXIMATE VALUE OF THE
ZERO
+ IS'
WRITE(9,*)' '
WRITE(9,*)' P =',P3
WRITE(9,*)' '
WRITE(9,*)' DP =',ABS(Dp),' IS AN ESTIMATE OF THE ACCURACY.'
WRITE(9,*)' '
WRITE(9,*)' F(',P3,' ) =',F(P3)
WRITE(9,*)' '
IF (F(P3).EQ.0) THEN
WRITE(9,*)'THE COMPUTED FUNCTION VALUE IS EXACTLY ZERO!'
WRITE(9,*)' '
ENDIF
IF (Cond.EQ.0) THEN
WRITE(9,*)'CONVERGENCE IS DOUBTFUL BECAUSE'
WRITE(9,*)' '
WRITE(9,*)'THE MAXIMUM NUMBER OF ITERATIONS WAS EXCEEDED.'
ELSEIF (Cond.EQ.1) THEN
WRITE(9,*)'CONVERGENCE IS DOUBTFUL. DIVISION BY ZERO WAS
ENCOUNT
+ERED.'
ELSEIF (Cond.EQ.2) THEN
WRITE(9,*)'THE APPROXIMATION P IS WITHIN THE DESIRED
TOLERANCE.'
ELSEIF (Cond.EQ.3) THEN
WRITE(9,*)'THE FUNCTION VALUE F(P) IS WITHIN THE DESIRED
TOLERAN
+CE.'
ELSEIF (Cond.EQ.4) THEN
WRITE(9,*)'THE APPROXIMATION P AND THE FUNCTION VALUE '
WRITE(9,*)'F(P) ARE BOTH WITHIN THE DESIRED TOLERANCES.'
ENDIF
RETURN
END

 

[Ian Collins]

jacob navia wrote, On 30/11/07 21:32:

Well, your compiler does not conform to the C standard *without* the
-ansic flag (for example declaring atof in math.h), so does it support
operator overloading *with* the -ansic flag? If not then it compiles C
when the -ansic flag is specified and something other than C that
supports operator overloading with it.

Just like gcc?

 

[CBFalconer]

Where the phrase "the language" refers to something other than C.

Since the result is a library to go with the non-C capabilities of
lcc-win32, that doesn't matter, and is of negligible interest on
c.l.c. In fact, probably off-topic.