Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.c > #138834 > unrolled thread

storing and retrieving a double may show a difference?

Back to article view | Back to comp.lang.c

Contents

  storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-20 23:25 -0800
    Re: storing and retrieving a double may show a difference? fir <profesor.fir@gmail.com> - 2018-12-21 00:56 -0800
      Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-21 07:02 -0800
        Re: storing and retrieving a double may show a difference? Bart <bc@freeuk.com> - 2018-12-21 15:33 +0000
          Re: storing and retrieving a double may show a difference? "Rick C. Hodgin" <rick.c.hodgin@gmail.com> - 2018-12-21 11:16 -0500
            Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-21 09:05 -0800
              Re: storing and retrieving a double may show a difference? fir <profesor.fir@gmail.com> - 2018-12-21 09:43 -0800
              Re: storing and retrieving a double may show a difference? "Rick C. Hodgin" <rick.c.hodgin@gmail.com> - 2018-12-21 12:57 -0500
              Re: storing and retrieving a double may show a difference? Ben Bacarisse <ben.usenet@bsb.me.uk> - 2018-12-21 20:06 +0000
        Re: storing and retrieving a double may show a difference? Ben Bacarisse <ben.usenet@bsb.me.uk> - 2018-12-21 17:10 +0000
          Re: storing and retrieving a double may show a difference? fir <profesor.fir@gmail.com> - 2018-12-21 09:17 -0800
            Re: storing and retrieving a double may show a difference? fir <profesor.fir@gmail.com> - 2018-12-21 09:20 -0800
            Re: storing and retrieving a double may show a difference? fir <profesor.fir@gmail.com> - 2018-12-21 09:24 -0800
          Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-21 21:46 -0800
          Re: storing and retrieving a double may show a difference? Manfred <noname@add.invalid> - 2018-12-22 21:41 +0100
            Re: storing and retrieving a double may show a difference? Ben Bacarisse <ben.usenet@bsb.me.uk> - 2018-12-22 21:08 +0000
              Re: storing and retrieving a double may show a difference? Manfred <noname@add.invalid> - 2018-12-23 20:34 +0100
                Re: storing and retrieving a double may show a difference? Tim Rentsch <txr@alumni.caltech.edu> - 2018-12-28 09:14 -0800
            Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-22 19:16 -0800
        Re: storing and retrieving a double may show a difference? fir <profesor.fir@gmail.com> - 2018-12-21 09:12 -0800
        Re: storing and retrieving a double may show a difference? James Kuyper <jameskuyper@alumni.caltech.edu> - 2018-12-21 09:23 -0800
          Re: storing and retrieving a double may show a difference? Ben Bacarisse <ben.usenet@bsb.me.uk> - 2018-12-21 20:03 +0000
            Re: storing and retrieving a double may show a difference? James Kuyper <jameskuyper@alumni.caltech.edu> - 2018-12-22 00:10 -0500
              Re: storing and retrieving a double may show a difference? Ben Bacarisse <ben.usenet@bsb.me.uk> - 2018-12-22 11:06 +0000
    Re: storing and retrieving a double may show a difference? mark.bluemel@gmail.com - 2018-12-21 01:00 -0800
      What type of problem is it? (Was: storing and retrieving a double may show a difference?) gazelle@shell.xmission.com (Kenny McCormack) - 2018-12-21 09:10 +0000
        Re: What type of problem is it? (Was: storing and retrieving a double may show a difference?) mark.bluemel@gmail.com - 2018-12-21 01:33 -0800
    Re: storing and retrieving a double may show a difference? Bart <bc@freeuk.com> - 2018-12-21 10:16 +0000
    Re: storing and retrieving a double may show a difference? Ben Bacarisse <ben.usenet@bsb.me.uk> - 2018-12-21 11:25 +0000
    Re: storing and retrieving a double may show a difference? James Kuyper <jameskuyper@alumni.caltech.edu> - 2018-12-21 08:06 -0500
      Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-21 07:13 -0800
        Re: storing and retrieving a double may show a difference? James Kuyper <jameskuyper@alumni.caltech.edu> - 2018-12-21 09:10 -0800
      Re: storing and retrieving a double may show a difference? Tim Rentsch <txr@alumni.caltech.edu> - 2018-12-28 08:16 -0800
    Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-21 08:09 -0800
    Re: storing and retrieving a double may show a difference? "Chris M. Thomasson" <invalid_chris_thomasson@invalid.invalid> - 2018-12-22 17:23 -0800
      Re: storing and retrieving a double may show a difference? "Chris M. Thomasson" <invalid_chris_thomasson@invalid.invalid> - 2018-12-22 17:58 -0800
    Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-22 19:30 -0800
      Re: storing and retrieving a double may show a difference? Richard Damon <Richard@Damon-Family.org> - 2018-12-22 22:48 -0500
        Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-22 20:10 -0800
          Re: storing and retrieving a double may show a difference? "Chris M. Thomasson" <invalid_chris_thomasson@invalid.invalid> - 2018-12-22 20:41 -0800
          Re: storing and retrieving a double may show a difference? Barry Schwarz <schwarzb@dqel.com> - 2018-12-22 22:02 -0800
            Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-23 02:43 -0800
              Re: storing and retrieving a double may show a difference? Ben Bacarisse <ben.usenet@bsb.me.uk> - 2018-12-23 12:38 +0000
                Re: storing and retrieving a double may show a difference? chanrasjid@gmail.com - 2018-12-23 05:12 -0800
                  Re: storing and retrieving a double may show a difference? fir <profesor.fir@gmail.com> - 2018-12-23 05:32 -0800
                    Re: storing and retrieving a double may show a difference? Ben Bacarisse <ben.usenet@bsb.me.uk> - 2018-12-23 14:55 +0000
              Re: storing and retrieving a double may show a difference? Richard Damon <Richard@Damon-Family.org> - 2018-12-23 07:43 -0500
      Re: storing and retrieving a double may show a difference? Manfred <noname@add.invalid> - 2018-12-23 20:42 +0100

Page 2 of 3 — ← Prev page 1 [2] 3  Next page →

#138882

On Friday, December 21, 2018 at 10:03:11 AM UTC-5, chanr...@gmail.com wrote:
...
> int width_A, height_A;// compute earlier
> 
> 
> // gaussian elimination for MLR, num_war,etc...
> 
> i 1, x 158161154320.539 != diagonal 158161154320.547
> test: mlr-v1.c:561: gauss: Assertion `0' failed.
> 
> I think the first answer abou fpu 80 bits and ram 64 bits may be the reason.
> 
> 
> double gauss() {
> 	int i, j, k;
> 	double x, y, diagonal[256];
> // augmented matrix A ;  variables x eq
> // initialize A for system of equations for gaussian elimination/ MLR.
> // must loop r nested within loop c else last row (r == num_var+1) not initialized.

A is presumably a pointer. Could you show us the declaration for A and the expression that is used to allocate memory for the array?

> 	for (i = 0; i < width_A; ++i) {
> 		// c  == eq
> 		for (j = i; j < height_A; ++j) {
> //r == variables
> 			x = 0;
> 			for (k = 0; k < num_obsv; ++k) {
> 				x += *(obsv + j * num_obsv + k) * *(obsv + i * num_obsv + k);
> 			}
> 			*(A + j * width_A + i) = x;
> 			if (i == j){
> 				diagonal[i] = x;
> 			}
> 			if (j > i && j < width_A) {
> 				// symmetry about diagonal; must; && needed, else assign overwrite earlier rows.
> 				*(A + i * width_A + j) = *(A + j * width_A + i);
> 			}
> 		}
> 	}
> //print_A("A before gauss");
> 
> //reduce UT upper triangle to 0
> 	for (i = 0; i < width_A; ++i) {
> 		// div col i by leading non-zero element == diag element a-ii
> 		x = *(A + i * width_A + i);
> 		//diagonal
> 		if (i > 0 && x != diagonal[i]){
> 			printf("i %d, x %.15g != diagonal %.15g\n",i, x ,diagonal[i]);
> 			assert(0);
> 		}
> 		for (j = i; j < height_A; ++j) {
> 			*(A + j * width_A + i) /= x;
> 		}
> 		assert(*(A + i * width_A + i) == 1.0);
> 		// now diag elem = 1;
> 		for (k = i + 1; k < width_A; ++k) {
> 			x = *(A + i * width_A + k); // mult
> 			for (j = i; j < height_A; ++j) {
> 				*(A + j * width_A + k) -= x * *(A + j * width_A + i);
> 			}
> 			assert(*(A + i * width_A + k) < 1.0e-300);
> 		}
> 	}
> 
> //print_A("A after gauss");
> // get soln B
> 	for (i = num_var; i >= 0; --i) {
> 		// i == col
> 		B[i] = *(A + (num_var + 1) * width_A + i);
> 		for (j = i + 1; j <= num_var; ++j) {
> 			B[i] -= B[j] * *(A + j * width_A + i);
> 		}
> 	}
> 
> 	/* SSE = ∑(y - y-predicted)^2
> 	 * SSY = ∑ (y - y-mean)^2
> 	 *
> 	 * R2 = 1 - SSE / SSY
> 	 */
> 
> // get predicted y obsv
> 	x = 0;
> 	for (i = 0; i < num_obsv; ++i) {
> 		y = B[0];
> 		for (j = 0; j < num_var; ++j) {
> 			y += B[j + 1] * *(obsv + (j + 1) * num_obsv + i);
> 		}
> 		x += pow(obsv_y[i] - y, 2); // SSE
> 	}
> 
> 	x = x / SSY; // 1 - R2
> 	R2 = 1 - x;
> 	//assert(R2 >= 0 && R2 <= 1.0);
> 	if (R2 < 0 || R2 > 1.0) {
> 		printf("R2 %g\n", R2);
> 		return HUGE_VAL;
> 	}
> 
> 	return x;
> }

I didn't notice any errors in this code, but I'm much less
reliable when I say "I see no errors" than I am when I say "I see
an error".

The key thing seems to be the values of A[width_A+1] and
diagonal[1]. If they start out equal (which should be tested),
then one or the other of them has to be changing during the
course of this program - which, as I understand the logic of your
code, shouldn't be happening. But it appears to be happening, so
you need to find out where. 

If you have access to decent debugger, I'd recommend waiting
until diagonal[1] has been set, then check the value of both
A[width_A+1] and diagonal[1]. If they match, then set a trace on
both objects to determine when one or the other changes.

[toc] | [prev] | [next] | [standalone]

#138888

James Kuyper <jameskuyper@alumni.caltech.edu> writes:
<snip>
> The key thing seems to be the values of A[width_A+1] and
> diagonal[1]. If they start out equal (which should be tested),

They start out equal by inspection of the code!

> then one or the other of them has to be changing during the
> course of this program - which, as I understand the logic of your
> code, shouldn't be happening. But it appears to be happening, so
> you need to find out where.

The logic of the program does exactly that -- it changes A[1][1] (and
other diagonal values) after the diagonal has been stored.  Now that
enough code has been posted, it's clearly just plain logic error.

-- 
Ben.

[toc] | [prev] | [next] | [standalone]

#138898

On 12/21/18 15:03, Ben Bacarisse wrote:
> James Kuyper <jameskuyper@alumni.caltech.edu> writes:
> <snip>
>> The key thing seems to be the values of A[width_A+1] and
>> diagonal[1]. If they start out equal (which should be tested),
> 
> They start out equal by inspection of the code!

Yes, well one of the things I've learned is that when your code
misbehaves, and you're not sure why, you shouldn't trust things that
seem guaranteed to be true by "inspection of the code". Since the code
is malfunctioning, at least one of those things isn't actually true.

>> then one or the other of them has to be changing during the
>> course of this program - which, as I understand the logic of your
>> code, shouldn't be happening. But it appears to be happening, so
>> you need to find out where.
> 
> The logic of the program does exactly that -- it changes A[1][1] (and
> other diagonal values) after the diagonal has been stored.  Now that
> enough code has been posted, it's clearly just plain logic error.

As I said, I'm far less reliable when I say "I don't see an error".

[toc] | [prev] | [next] | [standalone]

#138902

James Kuyper <jameskuyper@alumni.caltech.edu> writes:

> On 12/21/18 15:03, Ben Bacarisse wrote:
>> James Kuyper <jameskuyper@alumni.caltech.edu> writes:
>> <snip>
>>> The key thing seems to be the values of A[width_A+1] and
>>> diagonal[1]. If they start out equal (which should be tested),
>> 
>> They start out equal by inspection of the code!
>
> Yes, well one of the things I've learned is that when your code
> misbehaves, and you're not sure why, you shouldn't trust things that
> seem guaranteed to be true by "inspection of the code". Since the code
> is malfunctioning, at least one of those things isn't actually true.

I disagree.  I find that inspecting the code is often the simplest way
to find errors.  But even when it does not, it can save time by ruling
out (or at least putting to the bottom of your list) numerous
hypotheses.  Here, the diagonal element of the array and the element of
the 1D diagonal array are set from the same variable, x, on consecutive
lines.  Adding assert(a == b) after a = x; b = x; seems a waste of time
to me, even if a and b are more complex.  I might do that when I start
doubting everything, but it will be way down my list.

Now, I have come across compiler bugs that mean you can't believe what
you see in the code.  That leads you down into a rabbit hole where you
don't know what to believe, but fortunately, that is very rare.

<snip>
-- 
Ben.

[toc] | [prev] | [next] | [standalone]

#138838

On Friday, 21 December 2018 07:25:55 UTC, chanr...@gmail.com  wrote:
> All values are doubles.
> 		
>         x = *(A + c * width_A + c);
> 	diagonal_A[c] = x;
>         
> When diagonal is stored and tested for equality immediately, it tests equal:
> 
> c 1, x 103997435021.915, diagonal 103997435021.915
> diff 0
> 
> But when the diagonal was stored much earlier, an after some blocks of codes, x and diagonal will fail the equality test. 
> for (c = 0; c < width_A; ++c) {
> 		x = *(A + c * width_A + c);
> 		//diagonal_A[c] = x;
> 		if (x != diagonal_A[c]) {
> 			printf("c %d, x %.15g != diagonal %.15g\n", c, x, diagonal_A[c]);
> 			printf("diff %.15g\n", fabs(x - diagonal_A[c]));
> 			assert(0);
> 		}
> 
> c 1, x 103997435021.915 != diagonal 103997435022.047
> diff 0.131515502929688
> 
> Why a slight difference?

It would be more helpful to post a complete (minimal) program which
demonstrates the problem, rather than two out of context snippets.

If the problem is as simple to reproduce as you suggest, I'd expect
that you could create such a program quickly and simply. Often the
exercise of doing so will help you debug the problem for yourself.

[toc] | [prev] | [next] | [standalone]

#138839 — What type of problem is it? (Was: storing and retrieving a double may show a difference?)

In article <ef3b5e4e-316d-465e-afbb-d0db02389b2f@googlegroups.com>,
 <mark.bluemel@gmail.com> wrote:
...
>> Why a slight difference?
>
>It would be more helpful to post a complete (minimal) program which
>demonstrates the problem, rather than two out of context snippets.
>
>If the problem is as simple to reproduce as you suggest, I'd expect
>that you could create such a program quickly and simply. Often the
>exercise of doing so will help you debug the problem for yourself.

I think that posts of this sort are not of type:

    1) Please help me debug my program.

but are, rather, of the form:

    2) Is there a systemic reason why this happens?

In this case, it appears that the answer to question 2) is "No".
And, thus, that the snarky net followup is: There's a bug in your code.

Note that for many questions involving floating point math in C, that are
of the form: Why don't these two numbers compare equal?
the answer is:

    3) See:
	http://what/ever/the/path/is/to/thatArticle/AboutFloatingPointMath.html

which explains why floating point math is not exact on binary computers.

And that answer (answer 3)) is not in any way dependent on OP posting code.

-- 
"Only a genius could lose a billion dollars running a casino."
"You know what they say: the house always loses."
"When life gives you lemons, don't pay taxes."
"Grab 'em by the p***y!"

[toc] | [prev] | [next] | [standalone]

#138840 — Re: What type of problem is it? (Was: storing and retrieving a double may show a difference?)

On Friday, 21 December 2018 09:10:33 UTC, Kenny McCormack  wrote:
> In article <ef3b5e4e-316d-465e-afbb-d0db02389b2f@googlegroups.com>,
>  <mark.bluemel@gmail.com> wrote:
> ...
> >> Why a slight difference?
> >
> >It would be more helpful to post a complete (minimal) program which
> >demonstrates the problem, rather than two out of context snippets.
> >
> >If the problem is as simple to reproduce as you suggest, I'd expect
> >that you could create such a program quickly and simply. Often the
> >exercise of doing so will help you debug the problem for yourself.
> 
> I think that posts of this sort are not of type:
> 
>     1) Please help me debug my program.
> 
> but are, rather, of the form:
> 
>     2) Is there a systemic reason why this happens?
> 
> In this case, it appears that the answer to question 2) is "No".
> And, thus, that the snarky net followup is: There's a bug in your code.
> 
> Note that for many questions involving floating point math in C, that are
> of the form: Why don't these two numbers compare equal?
> the answer is:
> 
>     3) See:
> 	http://what/ever/the/path/is/to/thatArticle/AboutFloatingPointMath.html
> 
> which explains why floating point math is not exact on binary computers.
> 
> And that answer (answer 3)) is not in any way dependent on OP posting code.

I wondered about posting answer 3, but if the situation is as the OP seems to 
indicate, then answer 3 can't apply.
Their question seems to indicate that the values are calculated twice, in the
same way, from the same data. Given that the OP seems to believe the data is
changing when left alone (!) I think we need more information 

[toc] | [prev] | [next] | [standalone]

#138842

On 21/12/2018 07:25, chanrasjid@gmail.com wrote:
> All values are doubles.
> 		
>          x = *(A + c * width_A + c);
> 	diagonal_A[c] = x;
>          
> When diagonal is stored and tested for equality immediately, it tests equal:
> 
> c 1, x 103997435021.915, diagonal 103997435021.915
> diff 0
> 
> But when the diagonal was stored much earlier, an after some blocks of codes, x and diagonal will fail the equality test.
> for (c = 0; c < width_A; ++c) {
> 		x = *(A + c * width_A + c);
> 		//diagonal_A[c] = x;
> 		if (x != diagonal_A[c]) {
> 			printf("c %d, x %.15g != diagonal %.15g\n", c, x, diagonal_A[c]);
> 			printf("diff %.15g\n", fabs(x - diagonal_A[c]));
> 			assert(0);
> 		}
> 
> c 1, x 103997435021.915 != diagonal 103997435022.047
> diff 0.131515502929688
> 
> Why a slight difference?
> 

How was diagonal_A[1] set up earlier?

Are all float variables actually 'double' rather than 'float'?

The discrepancy in your values (about 1 in 1e12) is too small to be due 
to using exclusively float, and too large to be due to slight 
differences from order of calculation, which usually affect the last few 
bits. The difference here is around a dozen bits.


-- 
bart

[toc] | [prev] | [next] | [standalone]

#138845

chanrasjid@gmail.com writes:

> All values are doubles.

Not, I hope, c and width_A.

>         x = *(A + c * width_A + c);

It's more normal to write this as x = A[c * with_A + c], and in C99 you
can even write A[c][c] provided A is declared in the right way.  (You
might want to ask about that once your current issue is resolved.  C99
has a number of things that help when writing numerical programs.)

As to your issue...  Are you sure width_A is correct and that A (and
diagonal_A) points to the right sized block of data?  Those would be the
first places I'd look, particularly as the problem first shows up when
c==1.  At the very least you need to post the declarations of these and
the code that "sets them up".

Can you, in fact, post the whole program?

<snip>
> But when the diagonal was stored much earlier, an after some blocks of
> codes, x and diagonal will fail the equality test.
>
> for (c = 0; c < width_A; ++c) {
> 		x = *(A + c * width_A + c);
> 		//diagonal_A[c] = x;
> 		if (x != diagonal_A[c]) {
> 			printf("c %d, x %.15g != diagonal %.15g\n", c, x, diagonal_A[c]);
> 			printf("diff %.15g\n", fabs(x - diagonal_A[c]));
> 			assert(0);
> 		}
>
> c 1, x 103997435021.915 != diagonal 103997435022.047
> diff 0.131515502929688
>
> Why a slight difference?

Rather than assert(0) when something is wrong, you can write

  assert(A[c * width_A + c] == diagonal_A[c]);

at the top of the loop body and leave it there so as to catch any
recurrence of this issue, even after this problem is fixed.

-- 
Ben.

[toc] | [prev] | [next] | [standalone]

#138852

On 12/21/18 02:25, chanrasjid@gmail.com wrote:
> All values are doubles.
> 		
>         x = *(A + c * width_A + c);
> 	diagonal_A[c] = x;
>         
> When diagonal is stored and tested for equality immediately, it tests equal:
> 
> c 1, x 103997435021.915, diagonal 103997435021.915
> diff 0
> 
> But when the diagonal was stored much earlier, an after some blocks of codes, x and diagonal will fail the equality test. 
> for (c = 0; c < width_A; ++c) {
> 		x = *(A + c * width_A + c);
> 		//diagonal_A[c] = x;
> 		if (x != diagonal_A[c]) {
> 			printf("c %d, x %.15g != diagonal %.15g\n", c, x, diagonal_A[c]);
> 			printf("diff %.15g\n", fabs(x - diagonal_A[c]));
> 			assert(0);
> 		}
> 
> c 1, x 103997435021.915 != diagonal 103997435022.047
> diff 0.131515502929688
> 
> Why a slight difference?

I can make a guess, but seeing a complete program that demonstrates this
problem would make it a lot easier to investigate this issue. However,
some basic information would help:

What type is diagonal_A? I'm going to assume that diagonal_A is an array
of double. The macros, format strings, and math library functions used
below will all require appropriate changes if it's either float or long
double.

Please add the following lines to your code, and tell us what is
displayed by the printf() lines below:

#include <float.h>
#include <math.h>
#include <stdio.h>

/* insert the following code into some appropriate
 * function's body: */
{
    printf("__STDC_VERSION__:%ld\n",
          __STDC_VERSION__);
    printf("DBL_EPSILON:%g\n",
        DBL_EPSILON);

/* Insert the following code after the initial
 * calculation of diagonal_A[1]: */
#if __STDC_VERSION__ > 199901L
    #if __STDC_IEC_559__ == 1
        printf("Conforms to Annex F.\n");
    #else
        printf("Might not conform to Annex F.\n");
    #endif

    printf("FLT_EVAL_METHOD:%d\n",
        FLT_EVAL_METHOD);
    printf("diagonal_A[%d]=%.*g nextafter=%.*g\n",
         c, DECIMAL_DIG,
         diagonal_A[c], DECIMAL_DIG,
         nextafter(diagonal_A[c], DBL_MAX));
    printf("sizeof(double):%zu "
        "sizeof(double_t):%zu\n",
        sizeof(double), sizeof(double_t));
#endif
}

[toc] | [prev] | [next] | [standalone]

#138861

On Friday, December 21, 2018 at 9:07:00 PM UTC+8, James Kuyper wrote:
> On 12/21/18 02:25, chanrasjid@gmail.com wrote:
> > All values are doubles.
> > 		
> >         x = *(A + c * width_A + c);
> > 	diagonal_A[c] = x;
> >         
> > When diagonal is stored and tested for equality immediately, it tests equal:
> > 
> > c 1, x 103997435021.915, diagonal 103997435021.915
> > diff 0
> > 
> > But when the diagonal was stored much earlier, an after some blocks of codes, x and diagonal will fail the equality test. 
> > for (c = 0; c < width_A; ++c) {
> > 		x = *(A + c * width_A + c);
> > 		//diagonal_A[c] = x;
> > 		if (x != diagonal_A[c]) {
> > 			printf("c %d, x %.15g != diagonal %.15g\n", c, x, diagonal_A[c]);
> > 			printf("diff %.15g\n", fabs(x - diagonal_A[c]));
> > 			assert(0);
> > 		}
> > 
> > c 1, x 103997435021.915 != diagonal 103997435022.047
> > diff 0.131515502929688
> > 
> > Why a slight difference?
> 
> I can make a guess, but seeing a complete program that demonstrates this
> problem would make it a lot easier to investigate this issue. However,
> some basic information would help:
> 
> What type is diagonal_A? I'm going to assume that diagonal_A is an array
> of double. The macros, format strings, and math library functions used
> below will all require appropriate changes if it's either float or long
> double.
> 
> Please add the following lines to your code, and tell us what is
> displayed by the printf() lines below:
> 
> #include <float.h>
> #include <math.h>
> #include <stdio.h>
> 
> /* insert the following code into some appropriate
>  * function's body: */
> {
>     printf("__STDC_VERSION__:%ld\n",
>           __STDC_VERSION__);
>     printf("DBL_EPSILON:%g\n",
>         DBL_EPSILON);
> 
> /* Insert the following code after the initial
>  * calculation of diagonal_A[1]: */
> #if __STDC_VERSION__ > 199901L
>     #if __STDC_IEC_559__ == 1
>         printf("Conforms to Annex F.\n");
>     #else
>         printf("Might not conform to Annex F.\n");
>     #endif
> 
>     printf("FLT_EVAL_METHOD:%d\n",
>         FLT_EVAL_METHOD);
>     printf("diagonal_A[%d]=%.*g nextafter=%.*g\n",
>          c, DECIMAL_DIG,
>          diagonal_A[c], DECIMAL_DIG,
>          nextafter(diagonal_A[c], DBL_MAX));
>     printf("sizeof(double):%zu "
>         "sizeof(double_t):%zu\n",
>         sizeof(double), sizeof(double_t));
> #endif
> }

__STDC_VERSION__:201710
DBL_EPSILON:2.22045e-16
Conforms to Annex F.
FLT_EVAL_METHOD:0
diagonal[23]=6.91576149128030239731e-310 nextafter=6.91576149128035180387e-310
sizeof(double):8 sizeof(double_t):8
i 1, x 158161154320.539 != diagonal 158161154320.547
Conforms to Annex F.
FLT_EVAL_METHOD:0
diagonal_A[1]=158161154320.547302246 nextafter=158161154320.547332764
sizeof(double):8 sizeof(double_t):8
test: mlr-v1.c:603: full_gauss: Assertion `0' failed.

[toc] | [prev] | [next] | [standalone]

#138875

On Friday, December 21, 2018 at 10:13:46 AM UTC-5, chanr...@gmail.com wrote:
> On Friday, December 21, 2018 at 9:07:00 PM UTC+8, James Kuyper wrote:
> > On 12/21/18 02:25, chanrasjid@gmail.com wrote:
> > > All values are doubles.
> > > 		
> > >         x = *(A + c * width_A + c);
> > > 	diagonal_A[c] = x;
> > >         
> > > When diagonal is stored and tested for equality immediately, it tests equal:
> > > 
> > > c 1, x 103997435021.915, diagonal 103997435021.915
> > > diff 0
> > > 
> > > But when the diagonal was stored much earlier, an after some blocks of codes, x and diagonal will fail the equality test. 
> > > for (c = 0; c < width_A; ++c) {
> > > 		x = *(A + c * width_A + c);
> > > 		//diagonal_A[c] = x;
> > > 		if (x != diagonal_A[c]) {
> > > 			printf("c %d, x %.15g != diagonal %.15g\n", c, x, diagonal_A[c]);
> > > 			printf("diff %.15g\n", fabs(x - diagonal_A[c]));
> > > 			assert(0);
> > > 		}
> > > 
> > > c 1, x 103997435021.915 != diagonal 103997435022.047
> > > diff 0.131515502929688
> > > 
> > > Why a slight difference?
> > 
> > I can make a guess, but seeing a complete program that demonstrates this
> > problem would make it a lot easier to investigate this issue. However,
> > some basic information would help:

I'm glad I didn't make a guess - your results are inconsistent
with the guess I would have made.

> > What type is diagonal_A? I'm going to assume that diagonal_A is an array
> > of double. The macros, format strings, and math library functions used
> > below will all require appropriate changes if it's either float or long
> > double.
> > 
> > Please add the following lines to your code, and tell us what is
> > displayed by the printf() lines below:
> > 
> > #include <float.h>
> > #include <math.h>
> > #include <stdio.h>
> > 
> > /* insert the following code into some appropriate
> >  * function's body: */
> > {
> >     printf("__STDC_VERSION__:%ld\n",
> >           __STDC_VERSION__);
> >     printf("DBL_EPSILON:%g\n",
> >         DBL_EPSILON);
> > 
> > /* Insert the following code after the initial
> >  * calculation of diagonal_A[1]: */
> > #if __STDC_VERSION__ > 199901L
> >     #if __STDC_IEC_559__ == 1
> >         printf("Conforms to Annex F.\n");
> >     #else
> >         printf("Might not conform to Annex F.\n");
> >     #endif
> > 
> >     printf("FLT_EVAL_METHOD:%d\n",
> >         FLT_EVAL_METHOD);
> >     printf("diagonal_A[%d]=%.*g nextafter=%.*g\n",
> >          c, DECIMAL_DIG,
> >          diagonal_A[c], DECIMAL_DIG,
> >          nextafter(diagonal_A[c], DBL_MAX));
> >     printf("sizeof(double):%zu "
> >         "sizeof(double_t):%zu\n",
> >         sizeof(double), sizeof(double_t));
> > #endif
> > }
> 
> __STDC_VERSION__:201710
> DBL_EPSILON:2.22045e-16

==2^-52

> Conforms to Annex F.
> FLT_EVAL_METHOD:0
> diagonal[23]=6.91576149128030239731e-310 nextafter=6.91576149128035180387e-310
> sizeof(double):8 sizeof(double_t):8
> i 1, x 158161154320.539 != diagonal 158161154320.547
> Conforms to Annex F.
> FLT_EVAL_METHOD:0

That rules out the 80bit/64bit argument. If intermediate results
were being kept in 80 bit registers, but being stored in 64 bit
objects, then FLT_EVAL_METHOD should be 2, not 0.
 
> diagonal_A[1]=158161154320.547302246 nextafter=158161154320.547332764

The difference you're seeing are about 0.018, If it was an
80bit/64 bit issue, the differences should have been much
smaller. They should have been comparable to the difference
between consecutive representable 64-bit values, which in this
case is .000030518, or 2^-15.

> sizeof(double):8 sizeof(double_t):8

If they were using an 80bit intermediate type, then
sizeof(double_t) should have been 10, not 8.

> test: mlr-v1.c:603: full_gauss: Assertion `0' failed.

[toc] | [prev] | [next] | [standalone]

#139171

James Kuyper <jameskuyper@alumni.caltech.edu> writes:

> [question about unexpected program behavior]
>
> Please add the following lines to your code, and tell us what is
> displayed by the printf() lines below:
>
> [...]
>
>     printf("__STDC_VERSION__:%ld\n",
>           __STDC_VERSION__);

Use of __STDC_VERSION__ should be protected by

    #if  defined __STDC_VERSION__
       ...
    #endif

if the code might be compiled in a C89/C90 mode (since it need not
be defined under the C90 standard, and indeed is not defined by
either gcc or clang in such cases).


> /* Insert the following code after the initial
>  * calculation of diagonal_A[1]: */
> #if __STDC_VERSION__ > 199901L

Presumably this test was meant to be 

    #if __STDC_VERSION__ > 199900L

or

    #if __STDC_VERSION__ >= 199901L

so the delimited block will be compiled in C99 and later,
rather than being compiled only for post-C99.

[toc] | [prev] | [next] | [standalone]

#138870

My processor is core2 quad; Arch Linux.

[toc] | [prev] | [next] | [standalone]

#138929

On 12/20/2018 11:25 PM, chanrasjid@gmail.com wrote:
> All values are doubles.
> 		
>          x = *(A + c * width_A + c);
> 	diagonal_A[c] = x;
>          
> When diagonal is stored and tested for equality immediately, it tests equal:
> 
> c 1, x 103997435021.915, diagonal 103997435021.915
> diff 0
> 
> But when the diagonal was stored much earlier, an after some blocks of codes, x and diagonal will fail the equality test.
> for (c = 0; c < width_A; ++c) {
> 		x = *(A + c * width_A + c);
> 		//diagonal_A[c] = x;
> 		if (x != diagonal_A[c]) {
> 			printf("c %d, x %.15g != diagonal %.15g\n", c, x, diagonal_A[c]);
> 			printf("diff %.15g\n", fabs(x - diagonal_A[c]));
> 			assert(0);
> 		}
> 
> c 1, x 103997435021.915 != diagonal 103997435022.047
> diff 0.131515502929688
> 
> Why a slight difference?
> 

Should the follow hyper contrived program ever print "Ah Shucks!"?
______________________________
#include <stdio.h>
#include <math.h>


typedef double ct_real;


struct ct_real_buf
{
     unsigned char bytes[sizeof(ct_real)];
};


union ct_holder
{
     ct_real target;
     struct ct_real_buf buf;
};


void ct_hex_printf(
     FILE* fout,
     unsigned char* buf,
     size_t buf_sz
) {
     for (size_t i = 0; i < buf_sz; i++)
     {
         fprintf(fout, "%02x", buf[i]);
     }
}


int main(void)
{
     union ct_holder a;
     a.target = sqrt(2.);

     union ct_holder b = { .0 };

     b.buf = a.buf; // store a double

     printf("a = %f, hex = ", a.target);
     ct_hex_printf(stdout, a.buf.bytes, sizeof(a.buf.bytes));

     printf("\nb = %f, hex = ", b.target);
     ct_hex_printf(stdout, a.buf.bytes, sizeof(a.buf.bytes));

     if (a.target != b.target) // load, compare
     {
         printf("\n\nAh Shucks!\n");
     }

     else
     {
         printf("\n\nEQUAL!\n");
     }

     return 0;
}
______________________________


Fwiw, I get:
______________________________
a = 1.414214, hex = cd3b7f669ea0f63f
b = 1.414214, hex = cd3b7f669ea0f63f

EQUAL!
______________________________


;^)

[toc] | [prev] | [next] | [standalone]

#138930

On 12/22/2018 5:23 PM, Chris M. Thomasson wrote:
> On 12/20/2018 11:25 PM, chanrasjid@gmail.com wrote:
>> All values are doubles.
[...]
> 
> Should the follow hyper contrived program ever print "Ah Shucks!"?
> ______________________________
[...]
> int main(void)
> {
>      union ct_holder a;
>      a.target = sqrt(2.);
> 
>      union ct_holder b = { .0 };
> 
>      b.buf = a.buf; // store a double
> 
>      printf("a = %f, hex = ", a.target);
>      ct_hex_printf(stdout, a.buf.bytes, sizeof(a.buf.bytes));
> 



>      printf("\nb = %f, hex = ", b.target);
>      ct_hex_printf(stdout, a.buf.bytes, sizeof(a.buf.bytes));
[...]
                              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

I should probably use b, instead of a here wrt the hex output:

printf("\nb = %f, hex = ", b.target);
ct_hex_printf(stdout, b.buf.bytes, sizeof(b.buf.bytes));

Sorry about that.
_________________________
[...]
int main(void)
{
     union ct_holder a;
     a.target = sqrt(2.);

     union ct_holder b = { .0 };

     b.buf = a.buf; // store a double

     printf("a = %f, hex = ", a.target);
     ct_hex_printf(stdout, a.buf.bytes, sizeof(a.buf.bytes));

     printf("\nb = %f, hex = ", b.target);
     ct_hex_printf(stdout, b.buf.bytes, sizeof(b.buf.bytes)); // <---

     if (a.target != b.target) // load, compare
     {
         printf("\n\nAh Shucks!\n");
     }

     else
     {
         printf("\n\nEQUAL!\n");
     }

     return 0;
}
_________________________

[toc] | [prev] | [next] | [standalone]

#138933

Just a simple C question now.

In C, is it a proper programming thing to do a comparison of values for x, a double as: If (x == 0.0) or if (x == 1.0).

Ben just mentioned ==, !=, etc are "just" comparisons, so it seems OK to compare floating values. I think (x/x) will evaluate to/compare == 1.0 unless x == 0.

[toc] | [prev] | [next] | [standalone]

#138934

On 12/22/18 10:30 PM, chanrasjid@gmail.com wrote:
> Just a simple C question now.
> 
> In C, is it a proper programming thing to do a comparison of values for x, a double as: If (x == 0.0) or if (x == 1.0).
> 
> Ben just mentioned ==, !=, etc are "just" comparisons, so it seems OK to compare floating values. I think (x/x) will evaluate to/compare == 1.0 unless x == 0.
> 

Equality comparisons with floating point have the issue that floating
point numbers and calculations are really approximations, so the results
may not be exactly what you expect.

For instance, x + y - x - y might not be exactly 0. and
1. / 3 * 3 may not be exactly 1.0

Test for floating point equality does work, it will check if the two
number are equal, it is just that if some of the numbers are the results
of a calculation, the number you get may not be exactly what you think
you should have gotten.

The one spot where comparisons can get really unexpected results is when
you get NAN values involved, as if x = NAN, then x == x is false.

[toc] | [prev] | [next] | [standalone]

#138935

On Sunday, December 23, 2018 at 11:49:00 AM UTC+8, Richard Damon wrote:
> On 12/22/18 10:30 PM, chanrasjid@gmail.com wrote:
> > Just a simple C question now.
> > 
> > In C, is it a proper programming thing to do a comparison of values for x, a double as: If (x == 0.0) or if (x == 1.0).
> > 
> > Ben just mentioned ==, !=, etc are "just" comparisons, so it seems OK to compare floating values. I think (x/x) will evaluate to/compare == 1.0 unless x == 0.
> > 
> 
> Equality comparisons with floating point have the issue that floating
> point numbers and calculations are really approximations, so the results
> may not be exactly what you expect.
> 
> For instance, x + y - x - y might not be exactly 0. and
> 1. / 3 * 3 may not be exactly 1.0
> 
> Test for floating point equality does work, it will check if the two
> number are equal, it is just that if some of the numbers are the results
> of a calculation, the number you get may not be exactly what you think
> you should have gotten.
> 
> The one spot where comparisons can get really unexpected results is when
> you get NAN values involved, as if x = NAN, then x == x is false.

I think I understand what you mean by "testing for equality works" even for floating values. If x = 1.23 and y = 1.23, then any sane compiler should (must) give x == y as true; my values x and y are exact values. But there is no exact value for 1.0 / 3.0. 

[toc] | [prev] | [next] | [standalone]

#138937

On 12/22/2018 8:10 PM, chanrasjid@gmail.com wrote:
> On Sunday, December 23, 2018 at 11:49:00 AM UTC+8, Richard Damon wrote:
>> On 12/22/18 10:30 PM, chanrasjid@gmail.com wrote:
>>> Just a simple C question now.
>>>
>>> In C, is it a proper programming thing to do a comparison of values for x, a double as: If (x == 0.0) or if (x == 1.0).
>>>
>>> Ben just mentioned ==, !=, etc are "just" comparisons, so it seems OK to compare floating values. I think (x/x) will evaluate to/compare == 1.0 unless x == 0.
>>>
>>
>> Equality comparisons with floating point have the issue that floating
>> point numbers and calculations are really approximations, so the results
>> may not be exactly what you expect.
>>
>> For instance, x + y - x - y might not be exactly 0. and
>> 1. / 3 * 3 may not be exactly 1.0
>>
>> Test for floating point equality does work, it will check if the two
>> number are equal, it is just that if some of the numbers are the results
>> of a calculation, the number you get may not be exactly what you think
>> you should have gotten.
>>
>> The one spot where comparisons can get really unexpected results is when
>> you get NAN values involved, as if x = NAN, then x == x is false.
> 
> I think I understand what you mean by "testing for equality works" even for floating values. If x = 1.23 and y = 1.23, then any sane compiler should (must) give x == y as true; my values x and y are exact values. But there is no exact value for 1.0 / 3.0.
> 

Testing withing the same compiler, program and hardware should go ahead 
and work:

double a = 1. / 3.;
double b = a;

double c = fabs(a - b);

c should be damn small, if at all!

Also, a and b should have the same bit pattern when stored as bytes, so 
they should compare as equal without an epsilon. Yikes!

[toc] | [prev] | [next] | [standalone]

Page 2 of 3 — ← Prev page 1 [2] 3  Next page →

Back to top | Article view | comp.lang.c

csiph-web