Author: Gerd Isenberg
Date: 01:10:57 09/17/03
Go up one level in this thread
On September 16, 2003 at 20:06:58, Dieter Buerssner wrote:
>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.
Of course you are right, i will write my own "long" arithmetic class with
scalable precision to fresh up the basics ;-)
Gerd
>
>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.