FPN math bug

From https://cboard.cprogramming.com/c-programming/179473-i-cant-spot-my-mistake.html
Welp the new more thorough tests revealed a bug as I was expecting in FPN multiplication, division remains to be seen as the tests now stop on the 1st error, I'm not getting any ideas presently so will just post the code & output, upload the code and pray either God gives me an idea for fixing it or someone else gets given the idea and they post it.
Code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
int_t alup__mul( alup_t _NUM, alup_t _VAL, void *_cpy, void *_tmp )
{
	if ( alup_floating( _NUM ) || alup_floating( _VAL ) )
	{
		int_t ret;
		alup_t _DST, _SRC, _DMAN, _SMAN;
		ssize_t exp, dexp, sexp, dbias, sbias, dbits, sbits, bits, size;
		bool_t dneg, sneg;
		
		/* Ensure dealing with just floats, impossible for both to take _tmp
		 * given the above if statment */
		
		if ( alup_floating( _NUM ) )
		{
			_DST = _NUM;
		}
		else
		{
			bits = _NUM.upto - _NUM.from;
			size = BITS2SIZE( bits );
			alup_init_floating( _DST, _tmp, size );
			alup_mov( _DST, _NUM );
		}
		
		if ( alup_floating( _VAL ) )
		{
			_SRC = _VAL;
		}
		else
		{
			bits = _VAL.upto - _VAL.from;
			size = BITS2SIZE( bits );
			alup_init_floating( _SRC, _tmp, size );
			alup_mov( _SRC, _VAL );
		}
		
		dneg = alup_below0( _DST );
		sneg = alup_below0( _SRC );
		alub_set( _DST.data, _DST.upto - 1, dneg != sneg );
		
		dexp = alup_get_exponent( _DST );
		sexp = alup_get_exponent( _SRC );
		
		if ( !dexp || !sexp )
		{
			_NUM.upto--;
			alup_set( _NUM, 0 );
			return 0;
		}
		
		dbias = alup_get_exponent_bias( _DST );
		dexp -= dbias;
		
		sbias = alup_get_exponent_bias( _DST );
		sexp -= sbias;
		
		exp = dexp + sexp;
				
		alup_init_mantissa( _DST, _DMAN );
		alup_init_mantissa( _SRC, _SMAN );
		
		alub_set( _DMAN.data, _DMAN.upto, 1 );
		alub_set( _SMAN.data, _SMAN.upto, 1 );
		_DMAN.upto++;
		_SMAN.upto++;
		
		dbits = _DMAN.upto - _DMAN.from;
		sbits = _SMAN.upto - _SMAN.from;
		bits = LOWEST( dbits, sbits );
		
		_DMAN.from = _DMAN.upto - EITHER( dexp >= 0 && dexp < dbits, dexp, dbits );
		_SMAN.from = _SMAN.upto - EITHER( sexp >= 0 && sexp < sbits, sexp, sbits );
		
		ret = alup__mul_int2int( _DMAN, _SMAN, _cpy );
		
		/* We mangled this so restore it now */
		alup_set_exponent( _SRC, sexp + sbias );
		
		/* Multiplication prep may have changed this value, ensure is what it
		 * started as before using it any more */
		_DMAN.from = _DST.from;
		
		if ( exp > dbias )
		{	
			/* Set infinity */
			_DMAN.upto--;
			alup_set_exponent( _DST, (dbias << 1) | 1 );
			alup_set( _DMAN, 0 );
		}
		else
		{
			/* Normalise */
			alub_t final = alup_final_one( _DMAN );
			size_t pos = (_DMAN.upto - final.bit) - 1;
			
			if ( ret == EOVERFLOW )
			{
				exp++;
				alup__shr( _DMAN, pos );
				alub_set( _DMAN.data, _DMAN.upto - 2, 1 );
			}
			
			alup_set_exponent( _DST, exp + dbias );
		}
		
		if ( alup_floating( _NUM ) )
			return 0;
			
		return alup_mov( _NUM, _DST );
	}
	
	return alup__mul_int2int( _NUM, _VAL, _cpy );
}
...
const int per_func = LDBL_MANT_DIG;
bool_t stop_checks = false;

START_TEST( test_alup__math_floating )
{
	if ( !stop_checks )
	{
		ck_assert( alu_upto(alu) > 0 );
		size_t func = _i / per_func;
		double tmp_value1, tmp_value2
			, _value1 = _i
			, _value2 = _i % per_func
			, expect = _value1
			, result = _value1;
		alup_t _RESULT, _VALUE2;
		
		alup_init_floating( _RESULT, &result, sizeof(double) );
		alup_init_floating( _VALUE2, &_value2, sizeof(double) );
		
		expect = flt_math[func]( _value1, _value2 );
		alup__math[func]( _RESULT, _VALUE2, &tmp_value1, &tmp_value2 );
		
		if ( memcmp( &result, &expect, sizeof(double) ) != 0 )
		{
			alup_t _EXP, _MAN, _EXPECT, _VALUE1;
			
			alup_init_floating( _EXPECT, &expect, sizeof(double) );
			alup_init_floating( _VALUE1, &_value1, sizeof(double) );
			alup_init_exponent( _VALUE1, _EXP );
			alup_init_mantissa( _VALUE1, _MAN );
			
			alu_printf
			(
				"%f %c %f = %f, got %f, exp_dig = %zu, man_dig = %zu"
				, _value1
				, op_math[func]()
				, _value2
				, expect
				, result
				, _EXP.upto - _EXP.from
				, _MAN.upto - _MAN.from
			);
			
			alup_print( _VALUE1, 0, 1 );
			alup_print( _VALUE2, 0, 1 );
			alup_print( _EXPECT, 0, 1 );
			alup_print( _RESULT, 0, 1 );
			stop_checks = true;
		}
		
		ck_assert( memcmp( &result, &expect, sizeof(double) ) == 0 );
	}
}
END_TEST

Output:
1
2
3
4
5
6
7
8
9
10
make check.run
...
Running suite(s): ALU
../tests/check_alu.c:1151: test_alup__math_floating_fn() 131.000000 * 3.000000 = 393.000000, got 262.000000, exp_dig = 11, man_dig = 52
../tests/check_alu.c:1163: test_alup__math_floating_fn() _VALUE1 = 0 +0007 10000000110 0000011000000000000000000000000000000000000000000000
../tests/check_alu.c:1164: test_alup__math_floating_fn() _VALUE2 = 0 +0001 10000000000 1000000000000000000000000000000000000000000000000000
../tests/check_alu.c:1165: test_alup__math_floating_fn() _EXPECT = 0 +0008 10000000111 1000100100000000000000000000000000000000000000000000
../tests/check_alu.c:1166: test_alup__math_floating_fn() _RESULT = 0 +0008 10000000111 0000011000000000000000000000000000000000000000000000
99%: Checks: 1815, Failures: 1, Errors: 0
...
Never mind, fixed the code like this:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
int_t alup__mul( alup_t _NUM, alup_t _VAL, void *_cpy, void *_tmp )
{
	if ( alup_floating( _NUM ) || alup_floating( _VAL ) )
	{
		int_t ret;
		alup_t _DST, _SRC, _DMAN, _SMAN;
		ssize_t exp, dexp, sexp, dbias, sbias, dbits, sbits, bits, size;
		bool_t dneg, sneg;
		ssize_t ddiff = 0, sdiff = 0;
		size_t man_upto;
		alub_t final, prior;
		
		/* Ensure dealing with just floats, impossible for both to take _tmp
		 * given the above if statment */
		
		if ( alup_floating( _NUM ) )
		{
			_DST = _NUM;
		}
		else
		{
			bits = _NUM.upto - _NUM.from;
			size = BITS2SIZE( bits );
			alup_init_floating( _DST, _tmp, size );
			alup_mov( _DST, _NUM );
		}
		
		if ( alup_floating( _VAL ) )
		{
			_SRC = _VAL;
		}
		else
		{
			bits = _VAL.upto - _VAL.from;
			size = BITS2SIZE( bits );
			alup_init_floating( _SRC, _tmp, size );
			alup_mov( _SRC, _VAL );
		}
		
		dneg = alup_below0( _DST );
		sneg = alup_below0( _SRC );
		
		dexp = alup_get_exponent( _DST );
		sexp = alup_get_exponent( _SRC );
		
		if ( !dexp || !sexp )
		{
			_NUM.upto--;
			alup_set( _NUM, 0 );
			alub_set( _DST.data, _DST.upto - 1, dneg != sneg );
			return 0;
		}
		
		dbias = alup_get_exponent_bias( _DST );
		dexp -= dbias;
		
		sbias = alup_get_exponent_bias( _DST );
		sexp -= sbias;
		
		exp = dexp + sexp;
				
		alup_init_mantissa( _DST, _DMAN );
		alup_init_mantissa( _SRC, _SMAN );
		
		dbits = _DMAN.upto - _DMAN.from;
		sbits = _SMAN.upto - _SMAN.from;
		bits = LOWEST( dbits, sbits );
		
		ddiff = EITHER( dexp >= 0 && dexp < dbits, dexp, dbits );
		sdiff = EITHER( sexp >= 0 && sexp < sbits, sexp, sbits );
		prior = alub( _DST.data, _DST.upto - 1 );
		*(prior.ptr) &= ~(prior.mask);
		alup_set_exponent( _DST, 0 );
		
		alub_set( _DST.data, _DMAN.upto, 1 );
		alub_set( _SRC.data, _SMAN.upto, 1 );
		
		man_upto = _DMAN.upto;
		
		_SMAN.from = _SMAN.upto - sdiff;
		_DMAN.upto = _DST.upto;
		_SMAN.upto++;
		
		alup__shr( _DMAN, dbits - ddiff );
		prior = alup_final_one( _DMAN );
		
		ret = alup__mul_int2int( _DMAN, _SMAN, _cpy );
		
		/* We mangled this so restore it now */
		alup_set_exponent( _SRC, sexp + sbias );
		
		if ( exp >= dbias )
		{	
			/* Set infinity */
			alup_init_mantissa( _DST, _DMAN );
			
			alup_set( _DMAN, 0 );
			alup_set_exponent( _DST, (dbias << 1) | 1 );
		}
		else
		{	
			final = alup_final_one( _DMAN );
			_DMAN.upto = man_upto;
			
			if ( final.bit >= prior.bit )
			{
				/* Normalise */
				size_t mov = _DMAN.upto - final.bit;
				size_t add = final.bit - prior.bit;
				
				exp = dexp + add;
				alup__shl( _DMAN, mov );
			}
			else if ( ret == EOVERFLOW )
			{
				exp++;
				alup__shr( _DMAN, 1 );
				alub_set( _DMAN.data, _DMAN.upto - 1, 1 );
			}
			
			alup_set_exponent( _DST, exp + dbias );
		}
		
		alub_set( _DST.data, _DST.upto - 1, dneg != sneg );
		
		if ( alup_floating( _NUM ) )
			return 0;
			
		return alup_mov( _NUM, _DST );
	}
	
	return alup__mul_int2int( _NUM, _VAL, _cpy );
}
Topic archived. No new replies allowed.