Remainder Computation (Modulo) on Floating Numbers, Delphi and Python

The remainder computation on two integers a and b can be normally simply expressed as %, or mod in most of the current modern programming languages. The process of obtaining this remainder is fairly simple: where a / b refers the integer part of the division. Based on this, it can be extended to floating number using the following.

def fr(a, b):
    c = int(a / b)
    return a - c * b;

if __name__ == "__main__":
    print fr(-1, 5.234) # -1 

It is also already implemented in package math by the fmod method.

import math
print math.fmod(-1, 5.234) # -1

However, the % operator on floating numbers in Python directly yield different results when one of the floating numbers is negative.

print (-1) % 5.234 # 4.234

Using Inline Assembly (IA-32), the delphi implementation can be expressed to the following.

function asmFmod(const x: single; const y: single): single; assembler; register;
asm
  fld dword ptr[y]
  fld dword ptr[x]
@r:
  fprem
  fstsw ax
  sahf
  jp @r
  fstp st(1)
end

The IA-32 floating assembly operates on a stack size of 8. By loading instructions (e.g. fld) the floating numbers are pushed into the stack. fprem computes the partial remainder i.e. divide stack element 0 by stack element 1 and return the (partial) remainder to stack element 0.

The final result is returned using store and pop instruction (fstp). The single float type takes 4 byte which needs to be loaded using dword. If double type is used (8 bytes), the qword should do the tricks. To adopt both types, you can define using switches in Delphi.

{$DEFINE S_FLOAT}
{$IFDEF S_FLOAT}
type
  float = single;
{$ELSE}
type
  float = double;
{$ENDIF}

function asmFmod(const x: float; const y: float): float; assembler; register;
asm
{$IFDEF S_FLOAT}
  fld dword ptr[y]
  fld dword ptr[x]
@r:
  fprem
  fstsw ax
  sahf
  jp @r
  fstp st(1)
{$ELSE}
  fld qword ptr[y]
  fld qword ptr[x]
@r:
  fprem
  fstsw ax
  sahf
  jp @r
  fstp st(1)
{$ENDIF}
end;

The y can be replaced by which yields a new function.

function asmFmod2PI(const x: float): float; assembler; register;
asm
{$IFDEF S_FLOAT}
  fldpi
  fadd st, st
  fld dword ptr[x]
@r:
  fprem
  fstsw ax
  sahf
  jp @r
  fstp st(1)
{$ELSE}
  fldpi
  fadd st, st
  fld qword ptr[x]
@r:
  fprem
  fstsw ax
  sahf
  jp @r
  fstp st(1)
{$ENDIF}
end;

The fldpi loads the mathmatics constant into the floating-point stack. fadd st, st makes it .

–EOF (The Ultimate Computing & Technology Blog) —