flush to zero denormalized floating point numbers

[veselic]

Dear all,

in numerical applications with freepascal I wish the program to set to zero
each  denormalized small number. I hope this is possible,
but do no know how. In answering please consider that I am not expert
programmer.

Thanks in advance, Kresimir.

[Blaazen]

What is "denormalized small number" ?

[rvk]

What is "denormalized small number" ?

Just what I was thinking. So I googled it.

It's the gap between zero and underflow (Where underflow is: "The term arithmetic underflow (or "floating point underflow", or just "underflow") is a condition in a computer program where the result of a calculation is a smaller number than the computer can actually store in memory.").

So my question would be: how would you calculate something that can't be stored in memory? Won't this already truncate to zero?

@veselic: Could you give a programming example of what you want, exactly.

[engkin]

It simply can not be represented according to the normalized format. For instance, a normalized single precision number can be written as 1.xxxxE(-127..+127). Try to fit 1/1E127 in that format.

If I did not make a mistake, this procedure tells, beside basic parts, if a single precision number is normalized or not:

Code: [Select]

procedure ShowNumber(s: single);overload;
type
  TBits23 = 0..(1 shl 23) - 1;
  TBits8 = 0..(1 shl 8) - 1;

  TSingleRec=bitpacked record
    Fraction: TBits23; {1. implied}
    Exponent: TBits8;
    Sign: boolean;
  end;
var
  r: TSingleRec absolute s;
begin
  writeLn('Number: ', s);
  writeLn('Negative: ', r.Sign);
  writeLn('Exponent: ',r.Exponent,', Biased: ',r.Exponent-127);
  writeLn('Fraction: ',r.Fraction);
  writeLn('Normalized: ',(r.Exponent<>0));  //<------
  writeLn('');
end;

Edit:
The attached image is from "IA-32 Intel Architecture Developer’s Manual - Volume 1 - Basic Architecture".

[rvk]

Try to fit 1/1E127 in that format.

And assiging 1/1e127 to a single will result in the single having value 0. So there is no need to "flush it to zero" like veselic asked for.

(Or is there a chance of overflow arithmetic errors?)

[Nitorami]

I guess what Veselic wants is what is called "Disabling denormal floats at the code level" in wikipedia, http://en.wikipedia.org/wiki/Denormal_number

It requires to set FPU flags but I do not know how. In the compiler source, I found code to initialise the FPU

Code: [Select]

softfloat_exception_mask := float_flag_underflow or float_flag_inexact or float_flag_denormal;

We can probably change this mask but I have no clue how and what then happens.

[engkin]

Nope, these are for sw emulation. The equivalent ones are in Math unit. SetExceptionMask & exDenormalized for instance.

[engkin]

For SSE instructions the following MXCSR bits:
"Flush to zero" is bit 15.
"Denormals are zeros" is bit 6.

[Nitorami]

For the 8087, there is a control word and a status word.
The status word tells you when a denormalised operation has occured.
With the control word, you can tell the FPU not to do denormalised calculations. Instead, the FPU will then raise an exception which you need to handle.

In the attached example, the FPU is set NOT to denormalise (flag FPU_Denormal not set), therefore an exception will occur when repeatedly dividing a number by 2.

Note: I am not an expert on this and just cobbled it together using the source code and a FPU manual, so I do not know what side effects this may have. If you are not an "expert programmer", as you say, you should probably better leave this alone. You will now have to handle the exception to prevent your program crashing etc.

Is there a specific reason you need this or have you just read that denormalised numbers may be processed slower ?

Code: [Select]

const
  FPU_Invalid = 1;
  FPU_Denormal = 2;
  FPU_DivisionByZero = 4;
  FPU_Overflow = 8;
  FPU_Underflow = $10;
  FPU_StackUnderflow = $20;
  FPU_StackOverflow = $40;
  FPU_ExceptionMask = $ff;

//  fpucw : word = $1300 or FPU_StackUnderflow or FPU_Underflow or FPU_Denormal;
  fpucw : word = $1300 or FPU_StackUnderflow or FPU_Underflow;

Procedure SysResetFPU;
var
  localfpucw: word;
begin
  localfpucw:=fpucw;
  asm
    fninit
    fwait
    fldcw   localfpucw
  end;
end;

var i: integer;
    a: single;
begin
  SysResetFPU;
   a := 1e-40;
  for i := 1 to 40 do begin
    a := a/2;
    writeln (a);
  end;
end.


Addendum: SSE supports DAZ and FTZ without raising an exception, which is probably what you want. However it only works with SSE, and not on all processors. You need to set SSE bits but read this before:

https://software.intel.com/en-us/articles/x87-and-sse-floating-point-assists-in-ia-32-flush-to-zero-ftz-and-denormals-are-zero-daz/

The code may look something like this

Code: [Select]

const

  MM_MaskDAZMode   = %0000000000100000;
  MM_MaskInvalidOp = %0000000010000000;
  MM_MaskDenorm    = %0000000100000000;
  MM_MaskDivZero   = %0000001000000000;
  MM_MaskOverflow  = %0000010000000000;
  MM_MaskUnderflow = %0000100000000000;
  MM_MaskPrecision = %0001000000000000;
  MM_MaskFlushZero = %0100000000000000;

  mxcsr : dword = MM_MaskUnderflow or MM_MaskPrecision
               or MM_MaskDenorm or MM_MaskDAZMode;


Procedure SysResetFPU;
var
  localmxcsr: dword;
begin
  if has_sse_support then
    begin
      localmxcsr:=mxcsr;
      asm
        ldmxcsr localmxcsr
      end;
    end;
end;


[veselic]

Thanks for so much  consideration!
Sorry to have been imprecise. Yes, the problem is possible
poor efficiency, if there are plenty of denormalized numbers around
This seems to be a general phenomenon.  Let me give an example.
The program

program axb;
var a,c: single;
begin
a := 1.175E-38;
c := a/100;
writeln('denormalized = ',c:20);
end.

returns

denormalized =      1.175002775E-40

and I would like it to return pure zero. This is what I called
"flush to zero" and I know from contributing to some LAPACK
routines (not as the FORTRAN programmer, though)
 that this setting had decisive impact on the efficiency.

I also understood that not all processors would be able to comply
with that requirement

Thanks again, K. Veselic

[rvk]

Code: [Select]

const
  MM_MaskDAZMode   = %0000000000100000;
  MM_MaskInvalidOp = %0000000010000000;
  MM_MaskDenorm    = %0000000100000000;
  MM_MaskDivZero   = %0000001000000000;
  MM_MaskOverflow  = %0000010000000000;
  MM_MaskUnderflow = %0000100000000000;
  MM_MaskPrecision = %0001000000000000;
  MM_MaskFlushZero = %0100000000000000;


Are these correct? Shouldn't they start with %0000000001000000 (bit 6 counting from 0) and %1000000000000000 (bit 15) ??
So:

Code: [Select]

const
  MM_MaskDAZMode   = %0000000001000000;
  MM_MaskInvalidOp = %0000000100000000;
  MM_MaskDenorm    = %0000001000000000;
  MM_MaskDivZero   = %0000010000000000;
  MM_MaskOverflow  = %0000100000000000;
  MM_MaskUnderflow = %0001000000000000;
  MM_MaskPrecision = %0010000000000000;
  MM_MaskFlushZero = %1000000000000000;


[veselic]

flush-to-zero mode

I have understood the message, so I leave it alone.
Pity. With FORTARN this is just a compiler option
and -- as I could understand -- without unexpected side
effects.

Here are just some respectful general remarks for the
 pascal community:

In fact, I am still designing numerical algorithms in pascal
(they are later rewritten in FORTRAN by others, if needed).
The advantages
of pascal are for me (i) the transparence and the readability
of the source code and (ii) easy way to compute with different
precisions: single, double, extended and equally easy possibilty of their
mixing e.g. in accumulating scalar products. Tools like matlab or octave
just cannot do.

All these features might be also useful for teaching numerical methods
such that pascal is used as a pseudocode easily convertible to
a real program.

Thanking again to all of you for kind attention and help.

K. Veselic.

[Nitorami]

@rvk

Are these correct? Shouldn't they start with %0000000001000000 (bit 6 counting from 0) and %1000000000000000 (bit 15) ??

I think they are correct. I copied the keys from the FPC source and only added the first and the last.

@veselic
As you see this is not a problem in pascal, but has not been made available as compiler option. The programmer must find out himself, and this is reasonable I think.
Fortran was made for efficient mathematical calculations (Formula Translating System), and is mainly used by specialists, while Pascal is a generic all-purpose language. If the complier would give easy access to features like FTZ or DAZ, programmers who do not understand it would misuse this and then complain about crashes and wrong numerical results.

Yes, denormalised numbers process slower, but should not normally occur. The programmer should better check his code rather than turning of features of the FPU which exist for a good reason.
But I agree that in your numerical algorithms, denormals may by unavoidable, and you may have reason to use FTZ. But you should be aware that you will lose precision and that certain functions like comparison may yield different results. Before using FTZ, I would definitely take some time to study fundamental papers such as Goldberg, "What Every Computer Scientist Should Know About Floating-Point Arithmetic".

[veselic]

Thanks Nitorami,

I think numerical experiments also a legal way to learn how an algorithm
behaves -- including such possibly nefarious options
as flush-to zero-mode. It is another matter to offer such
experimental work to the public use.

By the way, I put your code segment in the declaration part of my program,
called Procedure SysResetFPU, it compiled OK but the output did not change.


Anyhow, thanks again, K. Veselic

[Nitorami]

Veselic

Sorry, it does not work because I indeed confused a few bits, as rvk already suspected. If you use the correct mask as below, it will work (for SSE), and smaller numbers than approximately 7e-39 will be flushed to zero.

Code: [Select]

const

  MM_MaskDAZMode   = %0000000001000000;
  MM_MaskInvalidOp = %0000000010000000;
  MM_MaskDenorm    = %0000000100000000;
  MM_MaskDivZero   = %0000001000000000;
  MM_MaskOverflow  = %0000010000000000;
  MM_MaskUnderflow = %0000100000000000;
  MM_MaskPrecision = %0001000000000000;
  MM_RoundingCtrl  = %0110000000000000;
  MM_MaskFTZ       = %1000000000000000;

  mxcsr : dword = MM_MaskUnderflow
               or MM_MaskPrecision
               or MM_MaskDenorm
//               or MM_MaskFTZ
               or MM_MaskDAZMode;

Regards
Martin