Author: Dieter Buerssner
Date: 17:06:58 09/16/03
Go up one level in this thread
On September 16, 2003 at 19:13:14, Gerd Isenberg wrote:
>On September 16, 2003 at 18:38:12, Dieter Buerssner wrote:
>
>>On September 16, 2003 at 17:48:32, Gerd Isenberg wrote:
>>
>>Hi Gerd,
>>only a comment to the later part.
>>
>>>void findGoodReciprocalConstants()
>>>{
>>> double errorMargin = 0.0000001; // !? 0 is also fine here with double
>>
>>Yes. At least until the result of the devision will be rather small, 0.0 will
>>work equally well.
>>
>>> double c = 0x100000000;
>>> c = c * c - 1.0; // 2**64-1
>>
>>Note, that this number is no exact double precision floating point number. You
>>will end up with 2**64. It is an exact extended precision floating point number
>>on x86 (80 bits with 64 bit mantissa. GNU-C uses this as type long double).
>>
>
>Hi Dieter,
>
>i see the problems - quick and dirty on the fly to find some constants with MSC.
>What about VB, or better FORTRAN.
I never used VB, and no FORTRAN for a long time. I have used Fortran, that could
use 128 bit floating point (I think on Sun). That would be enough. I am not
sure, wheter 80 bit floating point is enough (perhaps the last bit can be wrong,
then, in some cases). You could try GNU C on the PC. Last I checked, long double
was the 80 bit type. Also old versions of MSVC (which probably was not called
MSVC yet) used 80 bit long double. Newer MSVC uses the same type for double and
long double.
>>> for (int i = 1; i < 0x10000; i += 2)
>>> {
>>> double fdiv = (c + i)/i;
>>
>>The smallest result, you can get here is around 2**48. This leaves 5 more bits
>>of precision in type double. So all the results will have a smallest fraction of
>>2**-5. So you will never get numbers like integer-part.001. You can get
>>integer-part.00000 or integer-part.03125. At least when I did not do an off by
>>one error. Even if you make it 2**-6 or 2**-7, 0.0 in the error margin is the
>>same, as 0.0000001.
>>
>>> unsigned __int64 idiv = fdiv;
>>> double rdiv = (double)idiv;
>>> if ( fdiv >= rdiv - errorMargin && fdiv <= rdiv + errorMargin)
>>
>>Actually, with a decent compiler and good floating point implementation, you
>>could use
>>
>> if (rdiv == fdiv)
>>
>>here. This is one advantage of well defined IEEE floating point arithmetics.
>>
>>> printf("0x%08x%08x for mod(%5d)\n", (int)(idiv>>32), (int)idiv, i);
>>> }
>>>}
>>>
>>>For division/modulos by odd 64-bit constant less than 4096 a correction
>>>is not necessary.
>>
>>OTOH, I don't understand this conclusion. I think, you have missed, that double
>>precision floating point math is not exact enough here.
>>
>
>hmm... no inductive prove - only a feeling, excuse my math skills.
I think, it has nothing to do with math skills. Only with the details of
floating point math.
>If you try this "dirty" one with MSC:
>
>void findGoodReciprocalConstants()
>{
> double c = 0x100000000;
> c = c * c - 1.0;
> int count = 1;
> for (int i = 2; i < 0x10000; i += 1)
> {
> double fdiv = (c + i)/i;
> unsigned __int64 idiv = fdiv;
> double rdiv = (double)idiv;
> if ( fdiv == rdiv )
> printf("0x%08x%08x for mod(%5d) %d\n",
> (int)(idiv>>32), (int)idiv, i, ++count);
> }
>}
>
>and look to the trailing digits of the reciprocal constants:
>
>0x8000000000000000 for mod( 2) 2
>0x5555555555555400 for mod( 3) 3
>0x4000000000000000 for mod( 4) 4
>0x3333333333333400 for mod( 5) 5
>0x2aaaaaaaaaaaaa00 for mod( 6) 6
>0x2492492492492400 for mod( 7) 7
>0x2000000000000000 for mod( 8) 8
>0x1c71c71c71c71c00 for mod( 9) 9
>0x1999999999999a00 for mod( 10) 10
In my other post, you can see, that last 8 bits are wrong for 3, 5, 7, 9. They
are calso wrong for the even numbers.
Cheers,
Dieter
This page took 0 seconds to execute
Last modified: Thu, 15 Apr 21 08:11:13 -0700
Current Computer Chess Club Forums at Talkchess. This site by Sean Mintz.