VBCC - What's going on here?

[PeterK]

My C programming knowledge has almost faded away in the last 12 years and was never good, but you could try something to display the hexadecimal values of the floating point results like this (correct it if it's wrong):

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

int main(int argc, char const *argv[]) {
    double a = M_PI / 2;
    double x = cos(a);
    double y = sin(a);
    ulong *a1 = (ulong *)&a;
    ulong *a2 = (ulong *)&a+4;
    ulong *x1 = (ulong *)&x;
    ulong *x2 = (ulong *)&x+4;
    ulong *y1 = (ulong *)&y;
    ulong *y2 = (ulong *)&y+4;
    printf("%f: %f, %f\n", a, x, y);
    printf("a = %08lx %08lx, x = %08lx %08lx, y = %08lx %08lx\n", a1, a2, x1, x2, y1 ,y2);
    return 0;
}

[asymetrix]

Are you sure its not a precision error or rounding error. I read somewhere a while back to not use float and stick to double for more consistent / accurate maths.
What is the accuracy of cos/sin functions and does the storing of the variables crop the accuracy ?

We dont know the implementation used to calculate Cos/Sin - during the internal calculation, temp variable storage/math could have errors.

Scientists /mathematicians NEED accuracy and reliability - without it - its just not useful tool.
AKA standards.

Casting is NEVER good. Your just asking for trouble.
Compiling on different architecture systems would lead to very different results as all machines store/process the data differently.
very unsafe practice.

[alkis]

It could be either
a) Not original library, altered. (Cause you know...internet)
b) emulation side-effect bug

Probably someone with 2.04 has to run this http://eab.abime.net/showpost.php?p=...1&postcount=18 on real hardware. If it works, we need the md5sum of mathieeedoubbas.library. If it doesn't work then original system had the bug too, which I think is unlikely.

@PeterK
Something like that?

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

void print(char *name, double *d, int size) {
  char buf[80];
  char *c = (char *)d;
  char *s = &buf[0];

  for (int i = 0; i < size; ++i, s+=2, ++c)
    sprintf(s, "%02X", *c);
  *s = '\0';

  printf("%s %f %s\n", name, *d, buf);
}

int main(int argc, char const *argv[]) {
    double a = M_PI / 2;
    double x = cos(a);
    double y = sin(a);

    printf("%f: %f, %f\n", a, x, y);
    printf("%f\n", cos(M_PI / 2));
    printf("%f\n", sin(M_PI / 2));
    print("a", &a, sizeof(double));
    print("x", &x, sizeof(double));
    print("y", &y, sizeof(double));
    return 0;
}

[PeterK]

@alkis

No, I don't know what you are trying now, but all I would like to do is to print the hexadecimal values of a, x and y as they are stored in memory like

a = 1.570796 = $3FF921FB $54442D18
x = 0.000000 = $00000000 $00000000
y = 1.000000 = $3FF00000 $00000000

... but I don't want to get in touch with a C compiler now. I hate C !
Can you compile and test my code for me? Thx

@asymetrix
casting is only used here to keep the compiler quiet. I'm just trying to assign pointers to the memory locations of the double values of a, x and y. I'm not trying to convert the contents of a, x or y in any way.

[Photon]

I have ruled out emulator settings, and I think we can rule out the emulator also. With the only difference of putting the exe on a WB2.1 disk instead of a 2.04 disk, under WB2.1 it shows the correct value. Tested with and without SetPatch.

But this only means that you need WB 2.1 to run the generated code and get the expected result. So that we run the same executable doesn't mean anything, cos we swapped out the OS.

If we suspect libs, then it only takes one line of code to see if sin(M_PI/2) gives the expected result on WB 2.04. This could be done in any language that uses sin(). It could be done several times with different precision.

[PeterK]

The bad code is not in sin() directly, because the mathieeedoubtrans.library is the same in 2.04, 2.1 and 3.1, isn't it? The bad code comes from a function in mathieeedoubas.library. And that's what we still want to find out. It could be Add(), Sub(), Mul() or Div(), but I won't suspect these basic functions to be buggy. My guess is the Fix() function. Or maybe a bad Cmp(), who knows?

[alkis]

Ok, to make things even worse...

I compiled (beebo's gcc 6) the 'show me the hex bytes' thing. The values (aside from endianess) are the same on amiga and linux.

The thing is, executable from gcc displays correct values with "bad" and with "good" library.

I remember seeing Beebo mentioning something like you need to open math-single library even if you are to use only the double or you get incorrect values...let me try to find that reference.

Found it: https://github.com/bebbo/gcc/issues/...ment-366931667

[PeterK]

Found at least 2 bugs in the v37 Double Cmp() function. That could be used by Sin(), by Printf() or other functions, but most likely the buggy Cmp() is used by Printf() and that gives a wrong decision.

The thing is, executable from gcc displays correct values with "bad" and with "good" library.

GCC could have a different Printf() code which won't use Cmp() or Fix() from the mathlibs. It's easy to write a Cmp function with integer instructions only.

Update:
Found the same 2 Cmp() bugs also in v38.2, which is still used by OS 3.1! (this is fixed in v45 and perhaps OS 3.1.4)

A fix for v38.2 is attached to post #55 below.

So, we have to search for the culprit again ... later ...

[hmn]

Originally Posted by PeterK

My C programming knowledge has almost faded away in the last 12 years and was never good, but you could try something to display the hexadecimal values of the floating point results like this (correct it if it's wrong):

    double a = M_PI / 2;
    [...]
    ulong *a1 = (ulong *)&a;
    ulong *a2 = (ulong *)&a+4;

That is wrong Pointer addition does not work in bytes(*) - so &a+4 points to four sizeof(double) after the address of a.

EDIT: (*) unless you have a byte pointer of course...

[PeterK]

Originally Posted by hmn

That is wrong Pointer addition does not work in bytes(*) - so &a+4 points to four sizeof(double) after the address of a.
...

So what now ?
ulong *a2 = ((ulong *)&a)+1;

[deimos]

Maybe the problem is only showing up in printf - playing with the format string can cause the problem to disappear.

The following code should print three lines of 0 and 1, each line formatted slightly differently, but still printing the same values. the first and third lines are correct, but the second one isn't:

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

int main(int argc, char const *argv[]) {
    double y = sin(M_PI / 2);
    double x = cos(M_PI / 2);
    printf("%0.0f %0.0f\n", x, y);
    printf("%f %f\n", x, y);
    printf("%7.3f %7.3f\n", x, y);
    return 0;
}

gives:

0 1
0.000000 0.000000
0.00 1.000

[deimos]

And this shows that it only happens when I have more then one %f in my printf:

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

int main(int argc, char const *argv[]) {
    double x = cos(M_PI / 2);
    double y = sin(M_PI / 2);

    printf("zero zero then zero zero\n");
    printf("%f %f\n", x, x);
    printf("%f ", x); printf("%f\n", x);
    printf("%0.3f %0.3f\n", x, x);
    printf("\n");

    printf("zero one then zero one\n");
    printf("%f %f\n", x, y);
    printf("%f ", x); printf("%f\n", y);
    printf("%0.3f %0.3f\n", x, y);
    printf("\n");

    printf("one one then one one\n");
    printf("%f %f\n", y, y);
    printf("%f ", y); printf("%f\n", y);
    printf("%0.3f %0.3f\n", y, y);
    return 0;
}

Only the first line in the middle group is incorrect.

[deimos]

And the problem doesn't happen with variables that have been set to 1.0, only those that are very very close to 1.0 like the result of sin(M_PI / 2).

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

int main(int argc, char const *argv[]) {
    double x = 0.0;
    double y = 1.0;
    
    printf("zero one zero\n");
    printf("%f %f %f\n", x, y, x);
    printf("\n");

    printf("one zero one\n");
    printf("%f %f %f\n", y, x, y);
    printf("\n");
    
    x = cos(M_PI / 2);
    y = sin(M_PI / 2);
    
    printf("zero one zero\n");
    printf("%f %f %f\n", x, y, x);
    printf("\n");

    printf("one zero one\n");
    printf("%f %f %f\n", y, x, y);
    printf("\n");
    
    return 0;
}

[deimos]

Correction to above - it's not 1.0 that's the problem, it's 0.0.

The result of cos(M_PI / 2) is poisonous to printf.

The following code gives

-1.000000 0.000000 1.000000 2.000000 3.000000
-1.000000 0.000000 0.000000 0.000000 0.000000

Once that poisonous zero gets eaten by printf, everything it does is broken.

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

int main(int argc, char const *argv[]) {
    double neg_one = sin(-M_PI / 2);
    double zero = cos(M_PI / 2);
    double really_zero = 0.0;
    double one = sin(M_PI / 2);
    double two = 2.0;
    double three = 3.0;
    
    printf("%f %f %f %f %f\n", neg_one, really_zero, one, two, three);  // good
    printf("%f %f %f %f %f\n", neg_one, zero, one, two, three); // bad

    return 0;
}

[PeterK]

Btw, I found a mathieeedoubbas.library v39.1 on my harddisk, which is exactly a fix for the Cmp() bug in v38.2 that I mentioned in my previous post. Since I could not find the source for this patch from 1999, I've recreated it now. It has nothing to do with the printf() problem, but it could be useful for OS 3.1.

[alkis]

@deimos
I just saw your makefile. You don't link with math library. No -lmieee

Try this
vc bug.c -o bug -lmieee

[deimos]

Because I'm targeting 1.3, so I have -lm13 instead. -lmieee fails - No such file or directory.

Originally Posted by alkis

@deimos
I just saw your makefile. You don't link with math library. No -lmieee

Try this
vc bug.c -o bug -lmieee

[deimos]

Changing everything to target 2+ and adding -lmieee instead of -lm13 doesn't fix the problem.

Originally Posted by deimos

Because I'm targeting 1.3, so I have -lm13 instead. -lmieee fails - No such file or directory.

[deimos]

Perhaps I'll just do this:

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

double f(double a) {
    if (a > -1e-9 && a < 1e-9)
        return 0.0;
    return a;
}

int main(int argc, char const *argv[]) {
    double neg_one = sin(-M_PI / 2);
    double zero = cos(M_PI / 2);
    double really_zero = 0.0;
    double one = sin(M_PI / 2);
    double two = 2.0;
    double three = 3.0;
    
    printf("%f %f %f %f %f\n", neg_one, really_zero, one, two, three);  // good
    printf("%f %f %f %f %f\n", neg_one, zero, one, two, three); // bad
    printf("%f %f %f %f %f\n", f(neg_one), f(zero), f(one), f(two), f(three)); // fudged

    return 0;
}

[deimos]

Printing enough decimal places also makes the problem go away.

printf("%0.13f\n%0.13f\n%0.13f\n%0.13f\n%0.13f\n\n", neg_one, zero, one, two, three); // good
printf("%0.12f\n%0.12f\n%0.12f\n%0.12f\n%0.12f\n\n", neg_one, zero, one, two, three); // bad