mersenneforum.org  

Go Back   mersenneforum.org > Extra Stuff > Programming

Reply
 
Thread Tools
Old 2005-07-27, 14:35   #1
R.D. Silverman
 
R.D. Silverman's Avatar
 
Nov 2003

11101001001002 Posts
Thumbs up Bizarre!!!! Alert!!!!

Here is a snippet of code: (called in a loop that increments e)

if (e < lower_midpt) f_lower = round_double1(lm1 * (double)e + lc1 + .5);
else f_lower = round_double1(lm2 * (double)e + lc2 + .5);

if (e < lower_midpt) printf("LESS: \n");
else printf("MORE: \n");

if (e < lower_midpt) printf("LESS: %lf, %d\n",lower_midpt, e);
else printf("MORE: %lf %d\n", lower_midpt, e);

=========================================================
Here is some output

LESS:
LESS: -36.319540, -42
old_flower = 52
LESS:
LESS: -36.319540, -41
old_flower = 44
LESS:
LESS: -36.319540, -40
old_flower = 35
LESS:
LESS: -36.319540, -39
old_flower = 27
LESS:
LESS: -36.319540, -38
old_flower = 18
LESS:
LESS: -36.319540, -37
old_flower = 10
MORE:
MORE: -36.319540 -36
old_flower = 4
LESS:
MORE: -36.319540 -35
old_flower = 4
LESS:
MORE: -36.319540 -34
old_flower = 4
LESS:
MORE: -36.319540 -33
old_flower = 4
LESS:
MORE: -36.319540 -32
old_flower = 3
LESS:
MORE: -36.319540 -31
old_flower = 3
LESS:
MORE: -36.319540 -30
old_flower = 3
LESS:
MORE: -36.319540 -29
old_flower = 3
LESS:
MORE: -36.319540 -28


Would anyone care to explain this? As e increments from -36 to -35
suddenly the first comparison against midpt (-36.3) says -35 is less
than -36.3. But the second comparison gives a different answer.

I am not making this up!!!

Bob
R.D. Silverman is offline   Reply With Quote
Old 2005-07-27, 14:45   #2
R.D. Silverman
 
R.D. Silverman's Avatar
 
Nov 2003

22·5·373 Posts
Thumbs up

Quote:
Originally Posted by R.D. Silverman

Would anyone care to explain this? As e increments from -36 to -35
suddenly the first comparison against midpt (-36.3) says -35 is less
than -36.3. But the second comparison gives a different answer.



Bob
BTW, I know the answer.
R.D. Silverman is offline   Reply With Quote
Old 2005-07-27, 15:45   #3
akruppa
 
akruppa's Avatar
 
"Nancy"
Aug 2002
Alexandria

1001101000112 Posts
Default

What does the round_double1() function do? Is that the asm function to convert a double to an int we discussed in the other thread?

My initial guess is that the compiler keeps lower_midpt on top of the cpu stack and, since it doesn't know that round_double1() works on that stack, rearranges the opcodes so that the first if() ends up comparing e to the argument of the lower_midpt() function.

I'd need to take a look at the asm output for this code snippet.

Alex
akruppa is offline   Reply With Quote
Old 2005-07-27, 16:08   #4
R.D. Silverman
 
R.D. Silverman's Avatar
 
Nov 2003

746010 Posts
Default

Quote:
Originally Posted by akruppa
What does the round_double1() function do? Is that the asm function to convert a double to an int we discussed in the other thread?

My initial guess is that the compiler keeps lower_midpt on top of the cpu stack and, since it doesn't know that round_double1() works on that stack, rearranges the opcodes so that the first if() ends up comparing e to the argument of the lower_midpt() function.

I'd need to take a look at the asm output for this code snippet.

Alex
Yes. round_double1() is the function discussed earlier.

I am not sure if your explanation is right, or if somehow fistp resets
the FPU control word somehow..

Adding finit to round_double1() fixes the problem. I wish there were
a way to avoid it, since finit is expensive.

I am going to try popping the stack by fistp rather than fist

Bob
R.D. Silverman is offline   Reply With Quote
Old 2005-07-27, 16:18   #5
akruppa
 
akruppa's Avatar
 
"Nancy"
Aug 2002
Alexandria

246710 Posts
Default

Oh. right. I hadn't noticed. If this is the code you're using

__inline int iceil(double a)
{
int d;
const double h = 0.5000001f;
_asm
{
fld a
fadd h
fist d
};
return d;
}

it's surprising it works at all. If you load a value onto the stack with FLD, you need to get rid of it again, so use FISTP. Otherwise you'll just keep pushing things onto the stack and any value that was on it will never see the light of day again. The FIST would be good if you had a way to tell the compiler that the value is staying there, which is easy with the gcc asm() block - I don't know how to do it with VC, tough. Also, make the _asm_ block volatile so the compiler keeps it in a single block - reordering might lead to trouble if the compiler does not know about what is going on on the FPU stack.

I REALLY recommend using the lrint() function instead.

Alex

Last fiddled with by akruppa on 2005-07-27 at 16:19
akruppa is offline   Reply With Quote
Old 2005-07-27, 16:21   #6
R.D. Silverman
 
R.D. Silverman's Avatar
 
Nov 2003

22·5·373 Posts
Thumbs up

Quote:
Originally Posted by akruppa
Oh. right. I hadn't noticed. If this is the code you're using

__inline int iceil(double a)
{
int d;
const double h = 0.5000001f;
_asm
{
fld a
fadd h
fist d
};
return d;
}

it's surprising it works at all. If you load a value onto the stack with FLD, you need to get rid of it again, so use FISTP. Otherwise you'll just keep pushing things onto the stack and any value that was on it will never see the light of day again. The FIST would be good if you had a way to tell the compiler that the value is staying there, which is easy with the gcc asm() block - I don't know how to do it with VC, tough. Also, make the _asm_ block volatile so the compiler keeps it in a single block - reordering might lead to trouble if the compiler does not know about what is going on on the FPU stack.

I REALLY recommend using the lrint() function instead.

Alex

Popping the is indeed correct. As for lrint(), Microsoft VC++ doesn't
seem to know about it. At least, it is not indexed under 'help'' at all.
R.D. Silverman is offline   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
2 holes in bizarre theorem about composite Mersenne Numbers wildrabbitt Math 120 2016-09-29 21:52
Alert Sound Cricage Information & Answers 3 2016-01-21 16:49
Bizarre problem with gnfs-lasieve4I15e fivemack Factoring 1 2007-11-27 17:12
Spam Alert (all over again) S485122 Forum Feedback 8 2006-11-10 15:51
SPAMBOT ALERT! Andi47 Lounge 25 2006-11-03 04:33

All times are UTC. The time now is 17:47.

Wed May 12 17:47:45 UTC 2021 up 34 days, 12:28, 1 user, load averages: 1.77, 2.19, 2.79

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.