The case of the divide by zero

If you have any questions on programming, this is the place to ask them, whether you're a newbie or an experienced programmer. Discussion on programming in general is also welcome. We will help you with programming homework, but we will not do your work for you! Any porting requests must be made in Developmental Ideas.
Post Reply
User avatar
Newbie
Insane DCEmu
Insane DCEmu
Posts: 166
Joined: Sat Jul 27, 2013 1:16 pm
Has liked: 0
Been liked: 0

The case of the divide by zero

Post by Newbie » Sun Mar 26, 2017 7:06 pm

I could not obtain an error when I divide a float value by zero ....

I send a little C code that I use to make the test.

If you run it : you obtain an exception (quite normal).

If you replace in the source code the word "int" by "float" and run it : no error occurs.

It seems that for FPU unit, there is no CPU exception handling by default.

I try in assembly level to force SR.FD to 1 and just after execute a float instruction : no error occurs.

It is not what hardware manual says :
SR.jpg

Code: Select all


#ifndef KOS
#define KOS
#include <kos.h>
#endif

#pragma GCC push_options
#pragma GCC optimize ("O1")

int Get(uint32 code)
{
	switch(code)
	{
		case 1: 
		return 1.458;

		case 2: 
		return 2.878;

		case 3: 
		return 3.598;
		
		default :
		return 0.0;
	}
}

int TestDivideZero()
{
   volatile int x = 784.2541;
      
   volatile int y = Get(789);
   
   return x / y;  
}

int main(int argc, char* argv[])
{	
	volatile int result = TestDivideZero();
	
	return 0;
}
How strange and disturbing ...
User avatar
bogglez
Moderator
Moderator
Posts: 576
Joined: Sun Apr 20, 2014 9:45 am
Has liked: 0
Been liked: 0

Re: The case of the divide by zero

Post by bogglez » Sun Mar 26, 2017 7:27 pm

Just a wild guess, but that's probably GCC's "fault"
-mdiv=strategy
‘call-table’
Calls a library function that uses a lookup table for small divisors and the div1 instruction with case distinction for larger divisors. Division by zero calculates an unspecified result and does not trap. This is the default for SH4. Specifying this for targets that do not have dynamic shift instructions will default to call-div1.
https://gcc.gnu.org/onlinedocs/gcc-4.7. ... tions.html
Wiki & tutorials: http://dcemulation.org/?title=Development
Wiki feedback: viewtopic.php?f=29&t=103940
My libgl playground (not for production): https://bitbucket.org/bogglez/libgl15
My lxdream fork (with small fixes): https://bitbucket.org/bogglez/lxdream
User avatar
BlueCrab
The Crabby Overlord
The Crabby Overlord
Posts: 5383
Joined: Mon May 27, 2002 11:31 am
Location: Sailing the Skies of Arcadia
Has liked: 1 time
Been liked: 11 times
Contact:

Re: The case of the divide by zero

Post by BlueCrab » Sun Mar 26, 2017 7:58 pm

@Newbie: Are you sure that GCC is not just optimizing away the division entirely since the result is never actually being used? Have you looked at the assembly output and made sure that it is actually doing the division at all?
bogglez wrote:Just a wild guess, but that's probably GCC's "fault"
-mdiv=strategy
‘call-table’
Calls a library function that uses a lookup table for small divisors and the div1 instruction with case distinction for larger divisors. Division by zero calculates an unspecified result and does not trap. This is the default for SH4. Specifying this for targets that do not have dynamic shift instructions will default to call-div1.
https://gcc.gnu.org/onlinedocs/gcc-4.7. ... tions.html
That's only for integer division, I think. div0u/div0s and div1 only apply to integer division, after all. Floating point stuff should always be done by way of fdiv.
User avatar
Newbie
Insane DCEmu
Insane DCEmu
Posts: 166
Joined: Sat Jul 27, 2013 1:16 pm
Has liked: 0
Been liked: 0

Re: The case of the divide by zero

Post by Newbie » Mon Mar 27, 2017 11:12 am

Are you sure that GCC is not just optimizing away the division entirely since the result is never actually being used? Have you looked at the assembly output and made sure that it is actually doing the division at all?
Well, I thought about this possibility.

So to not be "corrupted" by how GCC could translate C code in asm, I made a second test in pure "volatile asm".

In this test, I get the SR value (with privileged instruction) , I make bit 15 to 1 to disable FPU, I refresh SR value (with privileged instruction), I execute a classic FPU instruction like FADD straight after.

When I run the test, it executes smoothly : but according to hardware manual it should trigger an exception ( a crash with dump of registers ...)

So it is very strange ...
User avatar
bogglez
Moderator
Moderator
Posts: 576
Joined: Sun Apr 20, 2014 9:45 am
Has liked: 0
Been liked: 0

Re: The case of the divide by zero

Post by bogglez » Mon Mar 27, 2017 1:58 pm

BlueCrab wrote:
bogglez wrote:Just a wild guess, but that's probably GCC's "fault"
-mdiv=strategy
‘call-table’
Calls a library function that uses a lookup table for small divisors and the div1 instruction with case distinction for larger divisors. Division by zero calculates an unspecified result and does not trap. This is the default for SH4. Specifying this for targets that do not have dynamic shift instructions will default to call-div1.
https://gcc.gnu.org/onlinedocs/gcc-4.7. ... tions.html
That's only for integer division, I think. div0u/div0s and div1 only apply to integer division, after all. Floating point stuff should always be done by way of fdiv.
Oops, you're right about that.

But wasn't there some libgcc function that's called to set up the fpu at the start of the program? Maybe it was turned off there.
Wiki & tutorials: http://dcemulation.org/?title=Development
Wiki feedback: viewtopic.php?f=29&t=103940
My libgl playground (not for production): https://bitbucket.org/bogglez/libgl15
My lxdream fork (with small fixes): https://bitbucket.org/bogglez/lxdream
nymus
DC Developer
DC Developer
Posts: 955
Joined: Tue Feb 11, 2003 4:12 pm
Location: In a Dream
Has liked: 0
Been liked: 1 time

Re: The case of the divide by zero

Post by nymus » Mon Mar 27, 2017 3:06 pm

I could be wrong, but I think KOS ignores FPU exceptions.

perhaps this default handler in irq.c explains it:

Code: Select all

/* Default FPU exception handler (can't seem to turn these off) */
static void irq_def_fpu(irq_t src, irq_context_t *context) {
    (void)src;
    context->pc += 2;
}
IEEE754 floating point generates exceptions for operations that might appear "normal." I remember being taught that instead of comparing to zero, one should compare to a delta because fp calculations are usually "inexact"... which is one of the exceptions specified in the standard.

IEEE754 also cannot offer exact representation for some numbers e.g. 33.3 will print 32.9999 so this might trigger an inexact exception. In graphics, where geometry is used a lot, I believe there are situations where infinity and zero should be generated but might not be or might be generated through truncation and rounding.

The FPU still writes sensible values to the target registers, e.g. divide by zero will write infinity which has a standard representation so you can compare to that and correctly conclude that whatever operation was performed was bad.

It would seem that exception handling is dependent on the user/programmer, e.g if you are comfortable with six digits of decimal precision, then you can ignore over/underflow caused by using 32-bit instead of 64-bit math.

just a naive example, but I also learned recently that using numbers less than 1 increases precision e.g if you want to multiply 550.750 * 150.250, you should multiply 0.550750 * 0.150250 then scale the result to reduce propagated errors.

Code: Select all

#include <stdio.h>
#include <stdlib.h>
#include <errno.h>

int main(int argc, char *argv[])
{
	errno = 0;
	float first = strtof(argv[1], NULL);
	if(0 != errno)
	{
		printf("\nBad input\n");
		exit(EXIT_FAILURE);
	}
	float second = strtof(argv[2], NULL);
	if(0 != errno)
	{
		printf("\nBad input\n");
		exit(EXIT_FAILURE);
	}

	printf("\n first (%f) - second (%f) = %f\n", first, second, first-second);

	printf("\n first (%f) / second (%f) = %f\n", first, second, first / second);
}
/* output:

$ ./fpu 360.333 180.333

 first (360.333008) - second (180.332993) = 180.000015

 first (360.333008) / second (180.332993) = 1.998154

*/
behold the mind
inspired by Dreamcast
nymus
DC Developer
DC Developer
Posts: 955
Joined: Tue Feb 11, 2003 4:12 pm
Location: In a Dream
Has liked: 0
Been liked: 1 time

Re: The case of the divide by zero

Post by nymus » Mon Mar 27, 2017 4:12 pm

I just thought I'd try some pseudocode to imagine how one might go about handling fpu exceptions... there could be errors.

Code: Select all

// fpscr cause and enable masks in one. will write cause (ignored?) and enable to fpscr when enabling
// will be compared to cause field when handling
// manual says cause is for current exception and
// flag for exception generated the last time it was cleared
// I'm just using cause for current exception

enum fpu_exception_source
{
	E_FPU_ERROR = 0x00020000,
	V_INVALID_OP = 0x00010800,
	Z_DIVIDE_BY_ZERO = 0x00008400,
	O_OVERFLOW = 0x00004200,
	U_UNDERFLOW = 0x00002100,
	I_INEXACT = 0x00001080
};

static void enable_fpu_exceptions(irq_context_t *context, int src)
{
	context->fpscr |= src;
}

static void irq_def_fpu(irq_t src, irq_context_t *context)
{
	int fpu_disabled = 0x800;
	int slot_fpu_disabled = 0x820;
	if((fpu_disabled == src) || (slot_fpu_disabled == src))
	{
		arch_panic("fpu is disabled");
	}

	if(0x120 != src)
	{
		arch_panic("incorrect handler");
	}
	int fpscr_cause = context->fpscr & 0x0003f000;
	switch(fpscr_cause)
	{
		case Z_DIVIDE_BY_ZERO:
			arch_panic("divide by zero");
			break;
		default:
			arch_panic("unhandled fpu error");
			break;
	}
	context->pc += 2;
}
behold the mind
inspired by Dreamcast
User avatar
bogglez
Moderator
Moderator
Posts: 576
Joined: Sun Apr 20, 2014 9:45 am
Has liked: 0
Been liked: 0

Re: The case of the divide by zero

Post by bogglez » Mon Mar 27, 2017 6:20 pm

nymus wrote: just a naive example, but I also learned recently that using numbers less than 1 increases precision e.g if you want to multiply 550.750 * 150.250, you should multiply 0.550750 * 0.150250 then scale the result to reduce propagated errors.
It's true that numbers close to 0 have more precision, but I think you're a bit confused about the rest.
This is a great read: https://docs.oracle.com/cd/E19957-01/80 ... dberg.html
Wiki & tutorials: http://dcemulation.org/?title=Development
Wiki feedback: viewtopic.php?f=29&t=103940
My libgl playground (not for production): https://bitbucket.org/bogglez/libgl15
My lxdream fork (with small fixes): https://bitbucket.org/bogglez/lxdream
nymus
DC Developer
DC Developer
Posts: 955
Joined: Tue Feb 11, 2003 4:12 pm
Location: In a Dream
Has liked: 0
Been liked: 1 time

Re: The case of the divide by zero

Post by nymus » Mon Mar 27, 2017 8:35 pm

Yes! That's the document I've seen refrerenced before. It doesn't make any more sense today than the couple of times I've looked at it before :D. However, being "confused about the rest" seems a bit broad... I certainly have an over-simplified view by comparison, but I wouldn't say I'm completely wrong?

Anyway, I noted the section about exception handling states that the default is to deliver a result and continue.
Exceptions, Flags and Trap Handlers

When an exceptional condition like division by zero or overflow occurs in IEEE arithmetic, the default is to deliver a result and continue. Typical of the default results are NaN for 0/0 and , and for 1/0 and overflow.
A couple of interesting links:
http://wiki.seas.harvard.edu/geos-chem/ ... ath_issues
http://www.cprogramming.com/tutorial/fl ... ation.html
http://davenewson.com/posts/2013/unity- ... cales.html

All the best.
behold the mind
inspired by Dreamcast
User avatar
bogglez
Moderator
Moderator
Posts: 576
Joined: Sun Apr 20, 2014 9:45 am
Has liked: 0
Been liked: 0

Re: The case of the divide by zero

Post by bogglez » Tue Mar 28, 2017 8:15 am

You mentioned 550.750 * 150.250 should be calculated as 0.550750 * 0.150250 instead, but that will give you exactly the same result.
The reason is that the numbers are stored in a normalized format, i.e. 0.mm * 2*exp. mm is the mantissa and its size depends on the data type you use (float vs double). The base is always 2 but you can change the exponent.

So assuming a mantissa with two bits you could have the following (I'm using binary numbers for the mantissa here since it's base 2)

Code: Select all

a = 0.10 * 2^0 = 0.10
b = 0.11 * 2^0 = 0.11
but
c = 0.10 * 2^1 = 1.00
d = 0.11 * 2^1 = 1.10
so
Distance between b and a: b-a = 0.11 - 0.10 = 0.01
Distance between d and c: d-c = 1.10 - 1.00 = 0.10 (the step is bigger, so we lost precision for bigger numbers here!)
Multiplication is actually not that much of a problem with floating point. The real problem is with addition and subtraction.

Code: Select all

#include <stdio.h>

int main() {
	float a = 100000000.0f;
	float b =         0.5f;
	float c = a+b;

	printf("  %f\n"
				 "+ %f\n"
				 "= %f\n",
				 a, b, c
	);
}

Code: Select all

  100000000.000000
+ 0.500000
= 100000000.000000
Wiki & tutorials: http://dcemulation.org/?title=Development
Wiki feedback: viewtopic.php?f=29&t=103940
My libgl playground (not for production): https://bitbucket.org/bogglez/libgl15
My lxdream fork (with small fixes): https://bitbucket.org/bogglez/lxdream
nymus
DC Developer
DC Developer
Posts: 955
Joined: Tue Feb 11, 2003 4:12 pm
Location: In a Dream
Has liked: 0
Been liked: 1 time

Re: The case of the divide by zero

Post by nymus » Tue Mar 28, 2017 10:10 am

Thanks! I wasn't sure which part in the post was inaccurate. That actually starts to make fp math clearer.

So, wrt the OP, there is little risk of arriving at zero by repeated division / multiplation (of values that fit) since the mantissa stays roughly the same and the exponent will hold the "largeness" so to speak?

For addition and subtraction, the problem is that the mantissa is actually being shifted around as happens with integer math?

Close enough?
behold the mind
inspired by Dreamcast
User avatar
bogglez
Moderator
Moderator
Posts: 576
Joined: Sun Apr 20, 2014 9:45 am
Has liked: 0
Been liked: 0

Re: The case of the divide by zero

Post by bogglez » Tue Mar 28, 2017 10:31 am

That will do. :)
Wiki & tutorials: http://dcemulation.org/?title=Development
Wiki feedback: viewtopic.php?f=29&t=103940
My libgl playground (not for production): https://bitbucket.org/bogglez/libgl15
My lxdream fork (with small fixes): https://bitbucket.org/bogglez/lxdream
User avatar
Newbie
Insane DCEmu
Insane DCEmu
Posts: 166
Joined: Sat Jul 27, 2013 1:16 pm
Has liked: 0
Been liked: 0

Re: The case of the divide by zero

Post by Newbie » Tue Mar 28, 2017 7:22 pm

The manual says clearly that there is a dedicated "divide by zero" exception for FPU.
As I do not achieve to trigger "divide by zero", I decide to move to the level up.

I said in a earlier post that in SH4 the SR register (Status register) contains a bit in position 15 called FD.
This FD bit is "FPU disable bit" : when FD = 1 there is no FPU instruction access. Besides, if we try to execute
an FPU instruction when FD bit is 1 in SR, the "general FPU disable" exception is triggered.
FD: FPU disable bit (cleared to 0 by a reset)
FD = 1: An FPU instruction causes a general FPU disable exception, and if the FPU instruction is in a delay
slot, a slot FPU disable exception is generated. (FPU instructions: H'F*** instructions, LDC(.L)/STS(.L)
instructions for FPUL/FPSCR)
In the manual the "general FPU disable" is described as follow :
General FPU Disable Exception
- Source: Decoding of an FPU instruction* not in a delay slot with SR.FD =1
- Transition address: VBR + H'0000 0100
- Transition operations:
The PC and SR contents for the instruction at which this exception occurred are saved in SPC and SSR.
Exception code H'800 is set in EXPEVT. The BL, MD, and RB bits are set to 1 in SR, and a branch is made to
PC = VBR + H'0100.
So I try to register a simple handler.

Code: Select all

irq_handler savedHandler = irq_get_handler(EXC_GENERAL_FPU);
irq_set_handler(EXC_GENERAL_FPU, MyHandler);
Set bit FD of SR to 1.

Code: Select all

__asm__ __volatile__("MOV #1, R10");      //R10 = 1
__asm__ __volatile__("SHLL8 R10");
__asm__ __volatile__("SHLL2 R10");
__asm__ __volatile__("SHLL2 R10");
__asm__ __volatile__("SHLL2 R10");
__asm__ __volatile__("STC SR, R11");      //SR -> R11
__asm__ __volatile__("MOV R11, R12");     //R11 -> R12
__asm__ __volatile__("OR R10, R11");      //R11 = R11 | R10
__asm__ __volatile__("LDC R11, SR");      //R11 -> SR
And execute any FPU instruction.

Code: Select all

__asm__ __volatile__("FADD FR0, FR1");
But the code produces no errors ....

I suspect none FPU exceptions is triggered (and not divide by zero too).

There is something wrong somewhere but i am stuck at this point ...
nymus
DC Developer
DC Developer
Posts: 955
Joined: Tue Feb 11, 2003 4:12 pm
Location: In a Dream
Has liked: 0
Been liked: 1 time

Re: The case of the divide by zero

Post by nymus » Wed Mar 29, 2017 3:48 am

I've tested disabling the fpu and it resets my Dreamcast right away, probably because of fp code in the kernel when interrupts are already blocked.

I think you need one more shift (shll) to set bit 15. shll8 shll2 shll2 shll2 moves 1 to position 14.

Code: Select all

// these handlers are not even called because cpu resets
static void exc_fpu(irq_t src, irq_context_t *context)
{
	if((EXC_GENERAL_FPU == src) || (EXC_SLOT_FPU == src))
	{
		asm __volatile__
		(
			"stc	sr, r0\n"
			"mov	#0x40, r1\n"
			"shll8	r1\n"
			"shll	r1\n"
			"not	r1, r1\n"
			"and	r1, r0\n"
			"ldc	r0, sr\n"
		);
		printf("\nFPU was disabled but is now enabled\n");
		context->pc += 2;
		return;
	}

	if(EXC_FPU != src)
	{
		printf("FPU exception handler called for other exception");
		arch_exit();
	}

	(void)src;
	context->pc += 2;
}

int main(int argc, char *arv[])
{
	printf("\nsaving old fpu handlers\n");
	irq_handler saved_general_fpu = irq_get_handler(EXC_GENERAL_FPU);
	irq_handler saved_slot_fpu = irq_get_handler(EXC_SLOT_FPU);

	printf("\nsetting custom fpu handler\n");
	irq_set_handler(EXC_GENERAL_FPU, exc_fpu);
	irq_set_handler(EXC_SLOT_FPU, exc_fpu);
	irq_set_handler(EXC_FPU, exc_fpu);

	printf("\ntesting fpu op\n");
	asm __volatile__
	(
		"fldi0	fr0\n"
	);
	printf("\ndisabling fpu\n");
	asm __volatile__
	(
		"stc	sr, r0\n"
		"mov	#0x40, r1\n"
		"shll8	r1\n"
		"shll	r1\n"
		"or		r1, r0\n"
		"ldc	r0, sr\n"
	);
/* not reached
	asm __volatile__
	(
		"fldi0	fr0\n"
	);
*/
}

behold the mind
inspired by Dreamcast
User avatar
Newbie
Insane DCEmu
Insane DCEmu
Posts: 166
Joined: Sat Jul 27, 2013 1:16 pm
Has liked: 0
Been liked: 0

Re: The case of the divide by zero

Post by Newbie » Wed Mar 29, 2017 7:57 am

I think you need one more shift (shll) to set bit 15. shll8 shll2 shll2 shll2 moves 1 to position 14
At first, I shifted the bits correctly but when I executed the code, I saw the machine make a reset.

I thought it was my fault and I had by error set a reserved bit (normally zero) in SR before making my LDC. So I removed a SHLL.

In fact, now I realize that my code was good and that the system reaction is wackily.

So it' s a "KOS Architectural dependence" ? We could not use "FPU exception handling" with KOS ?
Chilly Willy
DC Developer
DC Developer
Posts: 414
Joined: Thu Aug 20, 2009 11:00 am
Has liked: 0
Been liked: 2 times

Re: The case of the divide by zero

Post by Chilly Willy » Wed Mar 29, 2017 10:56 am

As said in earlier posts, you'll have to modify kos to do your own fp exception handler. The built-in one isn't designed for user exception code... at least, not for fp exceptions.
nymus
DC Developer
DC Developer
Posts: 955
Joined: Tue Feb 11, 2003 4:12 pm
Location: In a Dream
Has liked: 0
Been liked: 1 time

Re: The case of the divide by zero

Post by nymus » Wed Mar 29, 2017 12:40 pm

The following code seems to allow some fpu exception handling. I hope it helps.

Code: Select all

#include <kos.h>

typedef enum fpu_exception_source
{
	E_FPU_ERROR = 0x00020000,
	V_INVALID_OP = 0x00010800,
	Z_DIVIDE_BY_ZERO = 0x00008400,
	O_OVERFLOW = 0x00004200,
	U_UNDERFLOW = 0x00002100,
	I_INEXACT = 0x00001080,
	CAUSE_ENABLE_MASK = 0x0003ff80
} fpu_exc_src;

static void exc_fpu(irq_t src, irq_context_t *context)
{
	dbgio_printf("\nhandling fpu exception\n");

	if(EXC_FPU != src)
	{
		arch_panic("FPU exception handler called for other exception");

	}

	unsigned cause = context->fpscr & CAUSE_ENABLE_MASK;
	switch(cause)
	{
		case Z_DIVIDE_BY_ZERO:
			arch_panic("\ndivide by zero\n");
			break;
		default:
			arch_panic("\nunhandled fpu error\n");
			break;
	}
	
	context->pc += 2;
}

static void exc_fpu_enable(fpu_exc_src src)
{
	unsigned fpscr;

	asm __volatile__
	(
		"sts	fpscr, %0\n\t" : "=r"(fpscr)
	);

	fpscr |= src;

//	printf("\nwriting %x to fpscr\n", fpscr);

	asm __volatile__
	(
		"lds	%0, fpscr\n\t" : : "r"(fpscr)
	);
}

KOS_INIT_FLAGS(INIT_DEFAULT);

int main(int argc, char *arv[])
{
	printf("\nsetting custom fpu handler\n");
	irq_disable();
	irq_handler saved_exc_fpu = irq_get_handler(EXC_FPU);

	irq_set_handler(EXC_FPU, exc_fpu);
	irq_enable();

	printf("\ntesting fpu op\n");
	asm __volatile__
	(
		"fldi0	fr0\n\t"
	);
/*
	printf("\ndisabling fpu\n\t");
	asm __volatile__
	(
		"stc	sr, r0\n\t"
		"mov	#0x40, r1\n\t"
		"shll8	r1\n\t"
		"shll	r1\n\t"
		"or		r1, r0\n\t"
		"ldc	r0, sr\n\t"
	);

	asm __volatile__
	(
		"fldi0	fr0\n\t"
	);
*/
	printf("\nenabling divide by zero\n");
	exc_fpu_enable(Z_DIVIDE_BY_ZERO);

	float pi = 3.14159f;
	float zero;

	asm __volatile__
	(
		"fldi0	%0\n\t" : "=f"(zero)
	);
	printf("\n%f / %f = %f\n", pi, zero, pi/zero);

	printf("\ngoodbye\n");

	(void)saved_exc_fpu;

	return 0;
}
behold the mind
inspired by Dreamcast
User avatar
Newbie
Insane DCEmu
Insane DCEmu
Posts: 166
Joined: Sat Jul 27, 2013 1:16 pm
Has liked: 0
Been liked: 0

Re: The case of the divide by zero

Post by Newbie » Tue Apr 11, 2017 8:09 am

Thanks !
Post Reply