Discussion:
Determine that last digit of a number is 5
(too old to reply)
Robert AH Prins
2015-03-15 19:49:44 UTC
Permalink
I can divide by 10 or multiply with some magic, shift, multiply and subtract,
but is there a shorter way of determining that a number ends in the digit 5 -
I'm not interested in divisibility by 5, I just need to know that it ends in 5.

Why? For some input a program (or rather the RTL added by the compiler) has
trouble correctly rounding, no matter what I use for the FPU control word, so I
don't have any other choice than doing it myself.

Robert
--
Robert AH Prins
robert(a)prino(d)org
Terje Mathisen
2015-03-15 21:05:42 UTC
Permalink
Post by Robert AH Prins
I can divide by 10 or multiply with some magic, shift, multiply and
subtract, but is there a shorter way of determining that a number ends
in the digit 5 - I'm not interested in divisibility by 5, I just need to
know that it ends in 5.
Why? For some input a program (or rather the RTL added by the compiler)
has trouble correctly rounding, no matter what I use for the FPU control
word, so I don't have any other choice than doing it myself.
What is your platform? x87 stack variables? double/single/extended
precision?

Can you give an example?

FP or always integer values?

Terje
--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"
Robert Prins
2015-03-19 21:56:27 UTC
Permalink
Post by Robert AH Prins
I can divide by 10 or multiply with some magic, shift, multiply and
subtract, but is there a shorter way of determining that a number ends
in the digit 5 - I'm not interested in divisibility by 5, I just need to
know that it ends in 5.
Why? For some input a program (or rather the RTL added by the compiler)
has trouble correctly rounding, no matter what I use for the FPU control
word, so I don't have any other choice than doing it myself.
What is your platform? x87 stack variables? double/single/extended precision?
Can you give an example?
FP or always integer values?
The platform is Virtual Pascal, and the variables are double, but the FPU
control word is set to extended. One of the values for which I've been able to
figure out it's sequence of calculations is:

32-bit longint (37989) is multiplied by double (0.1) and stored in double (X)
32-bit longint (6) is stored in double (Y)

x is divided by Y, resulting in 633.15 and that passed to the routine that
should print it in a six character field with a single decimal position. The
receiving routine reads, in in-line assembler:

procedure add_double_2_line(var _double: double; _i, _j: longint); assembler;
{&uses none} {&frame-}
asm
mov eax, _double
mov edx, _i

fld qword ptr [eax]
mov eax, _j
mov ecx, eax
fld st

dec eax
jnz @01

fmul d_10v0
jmp @02

@01:
dec eax
jnz @02

fmul d_100v0

@02:
call System._Frac

push 6
push 4
push offset print
push type print - 1
call System._StrFlt

mov eax, ecx

cmp byte ptr print[3], "5"
jb @05

dec ecx
jnz @03

fadd d_0v05
jmp @05

@03:
dec ecx
jnz @05

fadd d_0v005
jmp @05

@04:
fadd d_0v5

@05:
movzx ecx, byte ptr _line
push edx // _i
add edx, 3
add byte ptr _line, dl
add ecx, offset _line

push eax // _j
push ecx
push type _line - 1
call System._StrFlt

mov byte ptr [ecx], " "
mov dword ptr [ecx + edx - 2], " | "
end; {add_double_2_line}

which was converted from the following Pascal code:

procedure add_double_2_line(var _double: double; _i, _j: longint);
var _corr: double;
var _d : double;

begin
_d:= _double;

case _j of
1: _corr:= d_10v0;
2: _corr:= d_100v0;
0: _corr:= d_1v0;
end;

{}write(_i:5, ' / ', _j:5, ' / ', _double:20:13, ' / ', round(_double * _corr *
10.0):10, ' / ', _corr:20:10, ' / ');

_corr:= frac(_double * _corr);

str(_corr:6:4, print);

{}write(_corr:20:10, ' / ', print:7, ' / ');

if print[3] >= '5' then
begin
case _j of
1: _corr:= d_0v05;
2: _corr:= d_0v005;
0: _corr:= d_0v5;
end;

_d:= _double + _corr;

{}write(_d:20:13, ' / ', _corr:20:10, ' / ');
end;

str(_double:_i:_j, print);
{}write(print:15, ' / ');
str(_d:_i:_j, print);
{}writeln(print:15);
add_print_2_line;
end; {add_double_2_line}

The output for the above values is, wrapping line:

6 / 1 / 633.1500000000000 / 63315 / 10.0000000000 / 0.5000000000 / 0.5000 /
633.1999999999999 / 0.0500000000 / 633.1 / 633.2

and the penultimate 633.1 is wrong. The 63315 value is the value where the last
digit is a 5 and that's the 5 I want to detect.

A case where things are OK?

6 / 1 / 975.8500000000000 / 97585 / 10.0000000000 / 0.5000000000 / 0.5000 /
975.9000000000000 / 0.0500000000 / 975.9 / 975.9

Here the original longints are 19517, giving X=1951.7 and 2, Y=2.

In both cases the round(_double * _corr * 10.0) gives the correct value with a
"5" as final digit.

I presume the whole problem occurs because of the various conversions and the
fact that the intermediate values (in some cases) are too small by epsilon...

Robert
--
Robert AH Prins
robert(a)prino(d)org
Terje Mathisen
2015-03-20 10:04:01 UTC
Permalink
Post by Robert Prins
Post by Robert AH Prins
I can divide by 10 or multiply with some magic, shift, multiply and
subtract, but is there a shorter way of determining that a number ends
in the digit 5 - I'm not interested in divisibility by 5, I just need to
know that it ends in 5.
Why? For some input a program (or rather the RTL added by the compiler)
has trouble correctly rounding, no matter what I use for the FPU control
word, so I don't have any other choice than doing it myself.
What is your platform? x87 stack variables? double/single/extended precision?
Can you give an example?
FP or always integer values?
The platform is Virtual Pascal, and the variables are double, but the
FPU control word is set to extended. One of the values for which I've
32-bit longint (37989) is multiplied by double (0.1) and stored in double (X)
32-bit longint (6) is stored in double (Y)
x is divided by Y, resulting in 633.15 and that passed to the routine
that should print it in a six character field with a single decimal
I am very sorry to say so, but you are living in a state of (floating
point) sin. :-(

It is in fact possible to use fp variables to store decimal data which
needs exact rounding witha given number of decimals, but you can only do
that if you know _exactly_ what you are doing, and spend quite a bit of
time on numerical analysis.

The easiest method is to work with scaled integers, i.e. if you know
that you need exact dollar results with two decimals, then you work in
cents instead.

Since Euro requires all financial calculations to be carried to 5
decimals, most library code has switched from the old 4 decimal format
to 5 in order to support Euro.

In your example above you could start with a scale factor of 10, so that
the calculation becomes:

37989 * 10 = 379890 (internal format)

Multiply by 0.1 to get back to 37989

Dividing that number by 6 gives 6331 with a remainder of 3, and since
3*2 is exactly equal to the divisor we know that the calculation should
be rounded up if the last result digit is odd (which it is), so the
proper result is 6332.

When you display this you have to re-insert the decimal point before the
last digit.

Terje
--
- <Terje.Mathisen at tmsw.no>
"almost all programming can be viewed as an exercise in caching"
Robert Prins
2015-03-20 21:41:31 UTC
Permalink
Post by Robert Prins
Post by Robert AH Prins
I can divide by 10 or multiply with some magic, shift, multiply and
subtract, but is there a shorter way of determining that a number ends
in the digit 5 - I'm not interested in divisibility by 5, I just need to
know that it ends in 5.
Why? For some input a program (or rather the RTL added by the compiler)
has trouble correctly rounding, no matter what I use for the FPU control
word, so I don't have any other choice than doing it myself.
What is your platform? x87 stack variables? double/single/extended precision?
Can you give an example?
FP or always integer values?
The platform is Virtual Pascal, and the variables are double, but the
FPU control word is set to extended. One of the values for which I've
32-bit longint (37989) is multiplied by double (0.1) and stored in double (X)
32-bit longint (6) is stored in double (Y)
x is divided by Y, resulting in 633.15 and that passed to the routine
that should print it in a six character field with a single decimal
I am very sorry to say so, but you are living in a state of (floating point)
sin. :-(
It is in fact possible to use fp variables to store decimal data which needs
exact rounding witha given number of decimals, but you can only do that if you
know _exactly_ what you are doing, and spend quite a bit of time on numerical
analysis.
I'm too used to PL/I where "fixed dec" variables are the norm, and where the
language very explicitly tells you what happens if you multiply a 7.1 (+/- 0.1
to 999999.9 + 000000.0) variable with a 3.2 (+/- 0.01 to 9.99 + 0.00) variable.
The easiest method is to work with scaled integers, i.e. if you know that you
need exact dollar results with two decimals, then you work in cents instead.
I'm already doing that for distances (hm rather than km) and times (minutes),
but doing it for velocities would require a lot of changes as there are
borderline cases, one of them in the input where the velocity is 100.0, rounded
from from 99.99529...
Since Euro requires all financial calculations to be carried to 5 decimals, most
library code has switched from the old 4 decimal format to 5 in order to support
Euro.
Been there done that, around AD 1999, using PL/I :)
In your example above you could start with a scale factor of 10, so that the
37989 * 10 = 379890 (internal format)
Multiply by 0.1 to get back to 37989
Dividing that number by 6 gives 6331 with a remainder of 3, and since 3*2 is
exactly equal to the divisor we know that the calculation should be rounded up
if the last result digit is odd (which it is), so the proper result is 6332.
When you display this you have to re-insert the decimal point before the last
digit.
This is probably the way to go, but having once done a conversion from the
ancient Turbo Pascal 6-byte reals to doubles, it's not something I'm really
looking forward to, especially given the fact that this problem only occurs in
less than 0.6% (221 out of 39658) of all printed doubles.

Robert
--
Robert AH Prins
robert(a)prino(d)org
David W Noon
2015-03-15 19:31:01 UTC
Permalink
On Sun, 15 Mar 2015 19:49:44 +0000, in message
<***@mid.individual.net>, Robert Ah Prins
(***@prino.org) wrote:

[snip]
Post by Robert AH Prins
Why? For some input a program (or rather the RTL added by the
compiler) has trouble correctly rounding, no matter what I use for
the FPU control word, so I don't have any other choice than doing
it myself.
I would infer that you mean a floating point number, since you refer
to the FPU control word.

Do you want to determine if the integer part ends in 5 or if the last
decimal digit of the fraction is 5? If it's the former, do you want
it rounded, forced or truncated to obtain the whole number?

E.g. 4.99999 ends in 5 if you round or force the fraction away, but
does not if you truncate it. Alternatively, 4.5 ends in 5 if you are
talking about the fraction, but it does not if you truncate. Real
problems arise when the decimal fraction cannot be represented in the
binary mantissa. e.g. 4.05.

So, what is the real question you are trying to ask?

As a first guess, I would suggest that PICTURE variables could be
useful. E.g.

last_digit_5:
PROC(some_number) RETURNS(BIT(1) ALIGNED) REORDER;
DCL some_number BIN FLOAT(53) BYVALUE;

/* Round the number to 7 decimal places. */
DCL my_number PIC '(8)9V(7)9' INIT(some_number + 5D-7),
last_digit CHAR(1) DEF my_number POS(15);

RETURN(last_digit = '5');
END;

If your real need is rather simpler than this, simplify the above as
convenient.
- --
Regards,

Dave [RLU #314465]
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
***@spamtrap.ntlworld.com (David W Noon)
Remove spam trap to reply by e-mail.
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
David W Noon
2015-03-16 01:11:43 UTC
Permalink
On Sun, 15 Mar 2015 19:31:01 +0000, in message
<me4mk8$mqt$***@dont-email.me>, David W Noon
(***@nospicedham.spamtrap.ntlworld.com) wrote:

On 15/03/15 19:31, David W Noon wrote:
[snip]
Post by David W Noon
E.g. 4.99999 ends in 5 if you round or force the fraction away,
but does not if you truncate it. Alternatively, 4.5 ends in 5 if
you are talking about the fraction, but it does not if you
truncate. Real problems arise when the decimal fraction cannot be
represented in the binary mantissa. e.g. 4.05.
So, what is the real question you are trying to ask?
Sorry about the PL/I code -- wrong newsgroup. ... :-(

Anyhow, my question remains: what does "ending in 5" actually mean?

[As an aside, that is the first message I have had get through to
c.l.a.x. for some few years, ever since I have been using
Eternal-September as my Usenet newsserver. I hope the connection
holds up.]
- --
Regards,

Dave [RLU #314465]
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
***@spamtrap.ntlworld.com (David W Noon)
Remove spam trap to reply by e-mail.
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Robert Prins
2015-03-19 18:02:36 UTC
Permalink
On 2015-03-15 19:31, David W Noon wrote:> -----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
On Sun, 15 Mar 2015 19:49:44 +0000, in message
[snip]
Post by Robert AH Prins
Why? For some input a program (or rather the RTL added by the
compiler) has trouble correctly rounding, no matter what I use for
the FPU control word, so I don't have any other choice than doing
it myself.
I would infer that you mean a floating point number, since you refer
to the FPU control word.
Do you want to determine if the integer part ends in 5 or if the last
decimal digit of the fraction is 5? If it's the former, do you want
it rounded, forced or truncated to obtain the whole number?
E.g. 4.99999 ends in 5 if you round or force the fraction away, but
does not if you truncate it. Alternatively, 4.5 ends in 5 if you are
talking about the fraction, but it does not if you truncate. Real
problems arise when the decimal fraction cannot be represented in the
binary mantissa. e.g. 4.05.
So, what is the real question you are trying to ask?
As a first guess, I would suggest that PICTURE variables could be
useful. E.g.
PROC(some_number) RETURNS(BIT(1) ALIGNED) REORDER;
DCL some_number BIN FLOAT(53) BYVALUE;
/* Round the number to 7 decimal places. */
DCL my_number PIC '(8)9V(7)9' INIT(some_number + 5D-7),
last_digit CHAR(1) DEF my_number POS(15);
RETURN(last_digit = '5');
END;
If your real need is rather simpler than this, simplify the above as
convenient.
You mentioning PL/I is quite amazing, or you must have followed me for longer,
as the problem I have is to match the output from PL/I and Pascal versions
working on the same data.

I don't need the test in the PL/I program, but if I did, given that the data is
held in fixed dec format, a simple

if substr(unspec(field), n, 4) = '0101'b then

would do the trick, in the current Enterprise PL/I compilers, not tested, this
may very well compile into a single instruction.

The PL/I ROUND builtin nicely rounds 0.5 up, but the Virtual Pascal RTL routine
to convert from double to character doesn't seem to honour the FPU control word,
and on some data I have to jump through the hoops of using code like:

procedure add_double_2_line(var _double: double; _i, _j: longint);
var _corr: double;
var _d : double;

begin
_d:= _double;

case _j of
1: _corr:= d_10v0;
2: _corr:= d_100v0;
0: _corr:= d_1v0;
end;

_corr:= frac(_double * _corr);

str(_corr:6:4, print);

if print[3] >= '5' then
begin
case _j of
1: _corr:= d_0v05;
2: _corr:= d_0v005;
0: _corr:= d_0v5;
end;

_d:= _double + _corr;
end;

str(_d:_i:_j, print);
add_print_2_line;
end; {add_double_2_line}

Which requires a double conversion of the floating point to string as the RTL
routine behind the Str procedure occasionally, no matter what the FPU control
word is, rounds the wrong way, probably for values that cannot be held exactly
into the double format.

Robert
--
Robert AH Prins
robert(a)prino(d)org
David W Noon
2015-03-21 21:13:56 UTC
Permalink
On Thu, 19 Mar 2015 18:02:36 +0000, in message
Post by Robert Prins
You mentioning PL/I is quite amazing, or you must have followed me
for longer, as the problem I have is to match the output from PL/I
and Pascal versions working on the same data.
So it's really a Pascal programming problem. [And Virtual Pascal at
that! That's a name I remember from 20+ years ago on OS/2.]

Try the following:
=======================================================================
function round_val(value : double; width, places : integer) : string;
var
scaled, fraction : double;
truncated : longint;
temp : string;
begin
case places of
0: scaled := value;
1: scaled := value*10.0;
2: scaled := value*100.0
end;

{ Split the floater into its whole and fractional parts. }
fraction := frac(scaled);
truncated := trunc(scaled);

{ Do our rounding. }
if fraction >= 0.5 then
truncated := truncated + 1;

{ Make the scaled, truncated value a string. }
Str(truncated:width, temp);

{ Insert a decimal point, if required. }
if places > 0 then
Insert('.', temp, Length(temp)-places+1);

{ Return the formatted value. }
round_val := temp
end;
=======================================================================
You can manipulate the rounding algorithm as you need.

HTH.
- --
Regards,

Dave [RLU #314465]
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
***@spamtrap.ntlworld.com (David W Noon)
Remove spam trap to reply by e-mail.
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Robert Prins
2015-03-21 23:38:45 UTC
Permalink
Post by Robert Prins
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
On Thu, 19 Mar 2015 18:02:36 +0000, in message
Post by Robert Prins
You mentioning PL/I is quite amazing, or you must have followed me
for longer, as the problem I have is to match the output from PL/I
and Pascal versions working on the same data.
So it's really a Pascal programming problem. [And Virtual Pascal at
that! That's a name I remember from 20+ years ago on OS/2.]
=======================================================================
function round_val(value : double; width, places : integer) : string;
var
scaled, fraction : double;
truncated : longint;
temp : string;
begin
case places of
0: scaled := value;
1: scaled := value*10.0;
2: scaled := value*100.0
end;
{ Split the floater into its whole and fractional parts. }
fraction := frac(scaled);
truncated := trunc(scaled);
{ Do our rounding. }
if fraction >= 0.5 then
truncated := truncated + 1;
{ Make the scaled, truncated value a string. }
Str(truncated:width, temp);
{ Insert a decimal point, if required. }
if places > 0 then
Insert('.', temp, Length(temp)-places+1);
{ Return the formatted value. }
round_val := temp
end;
=======================================================================
You can manipulate the rounding algorithm as you need.
Your code is a variation on what I'm using now, but rather more efficient, my
code is:

procedure add_double_2_line(var _dub: double; _i, _j: longint);
var _corr: double;
var _d : double;

begin
_d:= _dub;

case _j of
1: _corr:= d_10v0;
2: _corr:= d_100v0;
0: _corr:= d_1v0;
end;

_corr:= frac(_dub * _corr);

{
Obviously, directly testing _corr >= 0.5 is far more efficient, ouch!
}
str(_corr:6:4, print);
if print[3] >= '5' then
begin
case _j of
1: _corr:= d_0v05;
2: _corr:= d_0v005;
0: _corr:= d_0v5;
end;

_d:= _dub + _corr;
end;

str(_d:_i:_j, print);
add_print_2_line;
end; {add_double_2_line}

Using your code should remove the double float-to-string conversion. Still don't
like the frac(), as that involves a change of the FPU control word, but there's
not a lot that can be done about that.

Robert
--
Robert AH Prins
robert(a)prino(d)org
Robert Prins
2015-03-22 00:30:52 UTC
Permalink
Post by Robert Prins
Post by Robert Prins
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
On Thu, 19 Mar 2015 18:02:36 +0000, in message
Post by Robert Prins
You mentioning PL/I is quite amazing, or you must have followed me
for longer, as the problem I have is to match the output from PL/I
and Pascal versions working on the same data.
So it's really a Pascal programming problem. [And Virtual Pascal at
that! That's a name I remember from 20+ years ago on OS/2.]
=======================================================================
function round_val(value : double; width, places : integer) : string;
var
scaled, fraction : double;
truncated : longint;
temp : string;
begin
case places of
0: scaled := value;
1: scaled := value*10.0;
2: scaled := value*100.0
end;
{ Split the floater into its whole and fractional parts. }
fraction := frac(scaled);
truncated := trunc(scaled);
{ Do our rounding. }
if fraction >= 0.5 then
truncated := truncated + 1;
{ Make the scaled, truncated value a string. }
Str(truncated:width, temp);
{ Insert a decimal point, if required. }
if places > 0 then
Insert('.', temp, Length(temp)-places+1);
{ Return the formatted value. }
round_val := temp
end;
=======================================================================
You can manipulate the rounding algorithm as you need.
Your code is a variation on what I'm using now, but rather more efficient, my
procedure add_double_2_line(var _dub: double; _i, _j: longint);
var _corr: double;
var _d : double;
begin
_d:= _dub;
case _j of
1: _corr:= d_10v0;
2: _corr:= d_100v0;
0: _corr:= d_1v0;
end;
_corr:= frac(_dub * _corr);
{
Obviously, directly testing _corr >= 0.5 is far more efficient, ouch!
}
str(_corr:6:4, print);
if print[3] >= '5' then
begin
case _j of
1: _corr:= d_0v05;
2: _corr:= d_0v005;
0: _corr:= d_0v5;
end;
_d:= _dub + _corr;
end;
str(_d:_i:_j, print);
add_print_2_line;
end; {add_double_2_line}
Using your code should remove the double float-to-string conversion. Still don't
like the frac(), as that involves a change of the FPU control word, but there's
not a lot that can be done about that.
Except for the fact that I tried this before, and it does not work. 192.1 / 2
gives 96.0 where I need 96.1, fraction of the double is < 0.5, and 43 / 5 gives
8.61, and there are hundreds(!) more.

Robert
--
Robert AH Prins
robert(a)prino(d)org
David W Noon
2015-03-22 00:02:54 UTC
Permalink
On Sun, 22 Mar 2015 00:30:52 +0000, in message
<mekuqi$bf9$***@dont-email.me>, Robert Prins
(***@nospicedham.prino.org) wrote:

[snip]
Post by Robert Prins
Except for the fact that I tried this before, and it does not
work. 192.1 / 2 gives 96.0 where I need 96.1, fraction of the
double is < 0.5, and 43 / 5 gives 8.61, and there are hundreds(!)
more.
You are facing representational errors there.

You can deal with this if you come up with a fuzzy comparison against
0.5 for the fraction. Probably the first kluge you can try is to use
0.4999999999999999 [or the like] for the comparison.

It helps if you know the characteristics of your input data. If an
input value never carries more than, say, 2 places, you can likely
calculate a lower bound from 0.5 for representational errors. So it's
now a numerical analysis problem.
- --
Regards,

Dave [RLU #314465]
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
***@spamtrap.ntlworld.com (David W Noon)
Remove spam trap to reply by e-mail.
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Robert Prins
2015-03-22 11:50:31 UTC
Permalink
Post by Robert Prins
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
On Sun, 22 Mar 2015 00:30:52 +0000, in message
[snip]
Post by Robert Prins
Except for the fact that I tried this before, and it does not
work. 192.1 / 2 gives 96.0 where I need 96.1, fraction of the
double is < 0.5, and 43 / 5 gives 8.61, and there are hundreds(!)
more.
You are facing representational errors there.
Yes, I know all too well. IEEE-754 doubles and IBM's PL/I FIXED DEC are sadly
less than compatible, as are the rules for rounding.
Post by Robert Prins
You can deal with this if you come up with a fuzzy comparison against
0.5 for the fraction. Probably the first kluge you can try is to use
0.4999999999999999 [or the like] for the comparison.
That's why I now, in essence, use double rounding by converting the fraction
into a string. Ugly, but it works.
Post by Robert Prins
It helps if you know the characteristics of your input data. If an
input value never carries more than, say, 2 places, you can likely
calculate a lower bound from 0.5 for representational errors. So it's
now a numerical analysis problem.
The data passed to the routine ranges from 0.003 (0.003396546844041891) to
159822 (159822.225229205826) and the input leading to these values (don't know
where the 0.003... comes from) comes from simple divisions like 0.8 * 60 / 1 up
to 430404.0 * 60 / (4169 * 60 + 35) to more complicated expressions like
430404.0^2 *(636 + 3444) / (4169 + 35 / 60) / 636 / 3444, but the number of
decimal places is limited to 0, 1, or 2.

There's probably a way, by adding some magic numbers to the double values of
shifting the data into a 32-bit longint format, I remember having once seen this
technique, but sofar I've been unable to retrieve it, to speed up the code, but
given that using frac() alters the FPU controlword, a solution to change that to
round up rather than round to even may be possible.

Robert
--
Robert AH Prins
robert(a)prino(d)org
a***@nospicedham.gmail.com
2015-04-22 16:27:45 UTC
Permalink
But why round?
If it is in calculation, it is better the right value...
For store, what one would store
the right value etc

Round is need after the conversion
for show one result in a monitor
so it is a easy
routine that round array of
decimal char digits just for show
in a monitor or print in a paper
right digits.

But the machine number stored
for next calculation is the one
not round.

But if someone want to find one
round routine I not know/remember
for speak. I know it was a problem
a little difficult..
Robert Prins
2015-03-22 12:11:28 UTC
Permalink
Post by Robert Prins
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
On Sun, 22 Mar 2015 00:30:52 +0000, in message
[snip]
Post by Robert Prins
Except for the fact that I tried this before, and it does not
work. 192.1 / 2 gives 96.0 where I need 96.1, fraction of the
double is < 0.5, and 43 / 5 gives 8.61, and there are hundreds(!)
more.
You are facing representational errors there.
You can deal with this if you come up with a fuzzy comparison against
0.5 for the fraction. Probably the first kluge you can try is to use
0.4999999999999999 [or the like] for the comparison.
It helps if you know the characteristics of your input data. If an
input value never carries more than, say, 2 places, you can likely
calculate a lower bound from 0.5 for representational errors. So it's
now a numerical analysis problem.
I've added a few writeln()'s and the problem only occurs if the fractional part
of the input is of the format 0.x500000 and it needs to be printed with one
decimal place. In those cases the VP RTL consistently rounds down.

Robert
--
Robert AH Prins
robert(a)prino(d)org
Robert Prins
2015-03-22 17:14:17 UTC
Permalink
Post by Robert Prins
Post by Robert Prins
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
On Sun, 22 Mar 2015 00:30:52 +0000, in message
[snip]
Post by Robert Prins
Except for the fact that I tried this before, and it does not
work. 192.1 / 2 gives 96.0 where I need 96.1, fraction of the
double is < 0.5, and 43 / 5 gives 8.61, and there are hundreds(!)
more.
You are facing representational errors there.
You can deal with this if you come up with a fuzzy comparison against
0.5 for the fraction. Probably the first kluge you can try is to use
0.4999999999999999 [or the like] for the comparison.
It helps if you know the characteristics of your input data. If an
input value never carries more than, say, 2 places, you can likely
calculate a lower bound from 0.5 for representational errors. So it's
now a numerical analysis problem.
I've added a few writeln()'s and the problem only occurs if the fractional part
of the input is of the format 0.x500000 and it needs to be printed with one
decimal place. In those cases the VP RTL consistently rounds down.
Having concentrated on one of the problems, 3798.9 / 6, I've found that storing
this as a double, and subsequently loading the double, with the FPU
automatically converting it to extended leads to a value that is fractionally
too low, and that causes the problem. No clue why this only seems to happen for
values where the fraction is of the format 0.n500000, and other than using
extended throughout, which I don't even want to contemplate as if would require
replacing shl's with mul's all over the place, or use padding, I may stick with
what I have, despite this resulting in a double conversion of float to text. :(

Robert
--
Robert AH Prins
robert(a)prino(d)org
David W Noon
2015-03-22 18:18:10 UTC
Permalink
On Sun, 22 Mar 2015 17:14:17 +0000, in message
<mempjr$of0$***@dont-email.me>, Robert Prins
(***@nospicedham.prino.org) wrote:

[snip]
Post by Robert Prins
Having concentrated on one of the problems, 3798.9 / 6, I've found
that storing this as a double, and subsequently loading the double,
with the FPU automatically converting it to extended leads to a
value that is fractionally too low, and that causes the problem.
How attached are you to Virtual Pascal?

I have been using 32-bit Free Pascal for my testing and that value
works perfectly. I now suspect it is your compiler or RTL causing the
problem. Is it possible for you to switch to Free Pascal?

Anyhow, my suggested kluge of using 0.4999999 as the rounding
comparison could well fix it too.

[We should really use another newsgroup for this, unless you decide to
use assembler. However, the Pascal newsgroups are either empty in
recent years or full of New Age spam. Alternatively, follow the
instructions in my signature block to use e-mail.]
- --
Regards,

Dave [RLU #314465]
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
***@spamtrap.ntlworld.com (David W Noon)
Remove spam trap to reply by e-mail.
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Robert AH Prins
2015-03-22 19:12:37 UTC
Permalink
Post by Robert Prins
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
On Sun, 22 Mar 2015 17:14:17 +0000, in message
[snip]
Post by Robert Prins
Having concentrated on one of the problems, 3798.9 / 6, I've found
that storing this as a double, and subsequently loading the double,
with the FPU automatically converting it to extended leads to a
value that is fractionally too low, and that causes the problem.
How attached are you to Virtual Pascal?
I have been using 32-bit Free Pascal for my testing and that value
works perfectly. I now suspect it is your compiler or RTL causing the
problem. Is it possible for you to switch to Free Pascal?
I'm very attached to VP, as it allowed me to convert the original BP code almost
without changes. VP also generates really small executables. I've tried FPC on
several occasions, but it does weird things with alignment of variables and even
my pure Pascal programs abend and the executables are ridiculously big when
compared to VP. Converting VP to assembler is fun, but I have no plans to switch
to VP, as the IDE and debugger are, for me, still not in a state that I feel
comfortable with.
Post by Robert Prins
Anyhow, my suggested kluge of using 0.4999999 as the rounding
comparison could well fix it too.
My current setup also works, and using in-line assembler there may actually be a
more optimal way as I would be able to use 63 bit intermediate results (edx:eax)
that are impossible to use in VP - I don't want to use COMP's. Anyway, it
doesn't seem to be a VP or VP RTL issue, the issue seems to be that storing an
FPU value as double followed by reloading it changes the value, not by much, but
by enough to make some values round wrongly.
Post by Robert Prins
[We should really use another newsgroup for this, unless you decide to
use assembler. However, the Pascal newsgroups are either empty in
recent years or full of New Age spam. Alternatively, follow the
instructions in my signature block to use e-mail.]
I've stopped posting the c.l.p.borland FAQ years ago, I am still the maintainer,
but the group is dead. c.l.p.misc has been taken over by Amine Moulay Ramadine
and c.l.p.ansi-iso has not seen any postings since September 2013. Maybe, just
maybe it would be a good idea to merge all of the pascal groups back into
c.l.pascal...

For what it's worth, I'm using VP to do the heavy lifting, once the code works,
I take the resulting assembler listing and start hacking around. Probably 95% of
my code is nowadays in-line assembler, and that makes the program(s) both
smaller and a lot faster. The VP generated code is of the same quality of the
code generated by BP7, not very optimal, to be very polite, if you want I can
send you the programs I work on, or I can put them on my ftp site.

When I'm stuck with something I usually post here, which has led to some fairly
long threads in this group. :)

Robert
--
Robert AH Prins
robert(a)prino(d)org
David W Noon
2015-03-24 21:48:24 UTC
Permalink
On Sun, 22 Mar 2015 19:12:37 +0000, in message
Post by Robert AH Prins
My current setup also works, and using in-line assembler there may
actually be a more optimal way as I would be able to use 63 bit
intermediate results (edx:eax) that are impossible to use in VP -
I don't want to use COMP's. Anyway, it doesn't seem to be a VP or
VP RTL issue, the issue seems to be that storing an FPU value as
double followed by reloading it changes the value, not by much, but
by enough to make some values round wrongly.
I have been able to reproduce your problem, by munging the FPU control
word.

Here is an implementation with inline assembler that works on my
system here, and should be robust under Virtual Pascal too.

=========================================================================
{$ASMMODE Intel}

function round_val(value : double; width, places : integer) : string;
const
scaling_factor : array [0..2] of double = (1.0, 10.0, 100.0);
var
truncated : longint;
fpu_cw : shortint;
temp : string;
begin
fpu_cw := 0;
truncated := 0;

asm
FCLEX

FSTCW fpu_cw
MOV AX, WORD PTR fpu_cw
AND AX, 0F0FFH
OR AX, 0200H
MOV WORD PTR fpu_cw, AX
FLDCW fpu_cw

MOVZX EAX, places
FLD QWORD PTR scaling_factor[EAX*8]
FMUL value
FISTP DWORD PTR truncated
end ['EAX', 'ST(0)'];

{ Make the scaled, truncated value a string. }
Str(truncated, temp);

{ Insert a decimal point, if required. }
if places > 0 then
Insert('.', temp, Length(temp)-places+1);

if Length(temp) < width then
Insert(StringOfChar(' ', width-Length(temp)), temp, 1);

{ Return the formatted value. }
round_val := temp
end;
=========================================================================

Again, I used FPC to test this. You might need to massage the
assembler syntax if VP uses different syntax from FPC.

HTH.
- --
Regards,

Dave [RLU #314465]
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
***@spamtrap.ntlworld.com (David W Noon)
Remove spam trap to reply by e-mail.
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Robert Prins
2015-03-25 22:55:08 UTC
Permalink
Post by Robert Prins
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
On Sun, 22 Mar 2015 19:12:37 +0000, in message
Post by Robert AH Prins
My current setup also works, and using in-line assembler there may
actually be a more optimal way as I would be able to use 63 bit
intermediate results (edx:eax) that are impossible to use in VP -
I don't want to use COMP's. Anyway, it doesn't seem to be a VP or
VP RTL issue, the issue seems to be that storing an FPU value as
double followed by reloading it changes the value, not by much, but
by enough to make some values round wrongly.
I have been able to reproduce your problem, by munging the FPU control
word.
Here is an implementation with inline assembler that works on my
system here, and should be robust under Virtual Pascal too.
=========================================================================
{$ASMMODE Intel}
function round_val(value : double; width, places : integer) : string;
const
scaling_factor : array [0..2] of double = (1.0, 10.0, 100.0);
var
truncated : longint;
fpu_cw : shortint;
temp : string;
begin
fpu_cw := 0;
truncated := 0;
asm
FCLEX
FSTCW fpu_cw
MOV AX, WORD PTR fpu_cw
AND AX, 0F0FFH
OR AX, 0200H
MOV WORD PTR fpu_cw, AX
FLDCW fpu_cw
MOVZX EAX, places
FLD QWORD PTR scaling_factor[EAX*8]
FMUL value
FISTP DWORD PTR truncated
end ['EAX', 'ST(0)'];
{ Make the scaled, truncated value a string. }
Str(truncated, temp);
{ Insert a decimal point, if required. }
if places > 0 then
Insert('.', temp, Length(temp)-places+1);
if Length(temp) < width then
Insert(StringOfChar(' ', width-Length(temp)), temp, 1);
{ Return the formatted value. }
round_val := temp
end;
=========================================================================
Again, I used FPC to test this. You might need to massage the
assembler syntax if VP uses different syntax from FPC.
To round down, I use

push eax
fisttp [esp]
pop eax

where fisttp [esp] is coded as

db $db,$0c,$24 // fisttp dword [esp]

given the lack of this instruction in VP. I mercilessly assume that most current
CPUs will have this instruction.

Your solution always truncates, and that's what I do not want. As I wrote, some
correct values are munged by storing the extended and subsequently loading the
double. Your solution also has the problem that you don't restore the FPU
control word. I do like your loading of the scaling factor, and I will
incorporate that in my logic.

Robert
--
Robert AH Prins
robert(a)prino(d)org
David W Noon
2015-03-25 22:48:07 UTC
Permalink
On Wed, 25 Mar 2015 22:55:08 +0000, in message
<mevamn$ie$***@dont-email.me>, Robert Prins
(***@nospicedham.prino.org) wrote:

[snip]
Post by Robert Prins
To round down, I use
push eax fisttp [esp] pop eax
where fisttp [esp] is coded as
db $db,$0c,$24 // fisttp dword [esp]
given the lack of this instruction in VP. I mercilessly assume that
most current CPUs will have this instruction.
That instruction was not on the Pentium III and is not on the Athlon
MP that I am using here. [Yes, those are both museum pieces.]
Post by Robert Prins
Your solution always truncates,
No it doesn't. It always rounds to nearest.
Post by Robert Prins
and that's what I do not want. As I wrote, some correct values are
munged by storing the extended and subsequently loading the double.
Your solution also has the problem that you don't restore the FPU
control word.
That's easily remedied -- and largely irrelevant since VP's library
functions munge the control word anyway.
Post by Robert Prins
I do like your loading of the scaling factor, and I will
incorporate that in my logic.
On modern CPUs it's often faster (and rather simpler to code) to avoid
the branching logic of a case block and just do the FMUL anyway.
- --
Regards,

Dave [RLU #314465]
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
***@spamtrap.ntlworld.com (David W Noon)
Remove spam trap to reply by e-mail.
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Robert Prins
2015-03-29 21:52:27 UTC
Permalink
Post by Robert Prins
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
On Wed, 25 Mar 2015 22:55:08 +0000, in message
[snip]
Post by Robert Prins
To round down, I use
push eax fisttp [esp] pop eax
where fisttp [esp] is coded as
db $db,$0c,$24 // fisttp dword [esp]
given the lack of this instruction in VP. I mercilessly assume that
most current CPUs will have this instruction.
That instruction was not on the Pentium III and is not on the Athlon
MP that I am using here. [Yes, those are both museum pieces.]
Post by Robert Prins
Your solution always truncates,
No it doesn't. It always rounds to nearest.
Post by Robert Prins
and that's what I do not want. As I wrote, some correct values are
munged by storing the extended and subsequently loading the double.
Your solution also has the problem that you don't restore the FPU
control word.
That's easily remedied -- and largely irrelevant since VP's library
functions munge the control word anyway.
Post by Robert Prins
I do like your loading of the scaling factor, and I will
incorporate that in my logic.
On modern CPUs it's often faster (and rather simpler to code) to avoid
the branching logic of a case block and just do the FMUL anyway.
I've managed to remove part of the code that deals with doubles, the original
code would calculate

km * 6.0 / time

and print the resulting double. I've now changed that into

using the upper 32 bits of

((km * 600 / time + 5) * 3435973837) shr 3,

which, through reciprocal multiplication (3435973837 = 2^35 / 10), gives me 10
times the velocity, properly truncated and that is passed to a routine that
prints integers.

Robert
--
Robert AH Prins
robert(a)prino(d)org
David W Noon
2015-03-29 21:16:11 UTC
Permalink
On Sun, 29 Mar 2015 21:52:27 +0000, in message
<mf9ogn$ele$***@dont-email.me>, Robert Prins
(***@nospicedham.prino.org) wrote:

[snip]
Post by Robert Prins
I've managed to remove part of the code that deals with doubles,
the original code would calculate
km * 6.0 / time
and print the resulting double. I've now changed that into using
the upper 32 bits of
((km * 600 / time + 5) * 3435973837) shr 3,
which, through reciprocal multiplication (3435973837 = 2^35 / 10),
gives me 10 times the velocity, properly truncated and that is
passed to a routine that prints integers.
Ah! So you have bypassed the problem completely. ... :-)
- --
Regards,

Dave [RLU #314465]
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
***@spamtrap.ntlworld.com (David W Noon)
Remove spam trap to reply by e-mail.
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Robert Prins
2015-03-30 09:59:03 UTC
Permalink
Post by Robert Prins
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
On Sun, 29 Mar 2015 21:52:27 +0000, in message
[snip]
Post by Robert Prins
I've managed to remove part of the code that deals with doubles,
the original code would calculate
km * 6.0 / time
and print the resulting double. I've now changed that into using
the upper 32 bits of
((km * 600 / time + 5) * 3435973837) shr 3,
which, through reciprocal multiplication (3435973837 = 2^35 / 10),
gives me 10 times the velocity, properly truncated and that is
passed to a routine that prints integers.
Ah! So you have bypassed the problem completely. ... :-)
At least for these values. I've not yet looked into those where the one-decimal
place value is not a velocity, but may be the result of a more complicated
expression. I'm also not very happy with the final 'shr 3', but despite my best
efforts and going through Agner Fog's manuals, I've not been able to get rid of
it. Still, printing the integer will be substantially faster than printing the
double. ;)

Robert
--
Robert AH Prins
robert(a)prino(d)org
Bob Masta
2015-03-30 13:53:53 UTC
Permalink
On Sun, 29 Mar 2015 21:52:27 +0000, Robert Prins
Post by Robert Prins
I've managed to remove part of the code that deals with doubles, the original
code would calculate
km * 6.0 / time
and print the resulting double. I've now changed that into
using the upper 32 bits of
((km * 600 / time + 5) * 3435973837) shr 3,
which, through reciprocal multiplication (3435973837 = 2^35 / 10), gives me 10
times the velocity, properly truncated and that is passed to a routine that
prints integers.
I am missing the point of the reciprocal multiplication.

If you are looking for 10 times the velocity, rounded to the
nearest integer (ie nearest 10th of the raw velocity),
what's wrong with:

int(10 * km * 600 / time + 0.5)

where int() truncates without rounding?

What am I missing?

Best regards,


Bob Masta

DAQARTA v7.60
Data AcQuisition And Real-Time Analysis
www.daqarta.com
Scope, Spectrum, Spectrogram, Sound Level Meter
Frequency Counter, Pitch Track, Pitch-to-MIDI
FREE Signal Generator, DaqMusiq generator
Science with your sound card!
Robert Prins
2015-03-30 21:26:20 UTC
Permalink
Post by Bob Masta
On Sun, 29 Mar 2015 21:52:27 +0000, Robert Prins
Post by Robert Prins
I've managed to remove part of the code that deals with doubles, the original
code would calculate
km * 6.0 / time
and print the resulting double. I've now changed that into
using the upper 32 bits of
((km * 600 / time + 5) * 3435973837) shr 3,
which, through reciprocal multiplication (3435973837 = 2^35 / 10), gives me 10
times the velocity, properly truncated and that is passed to a routine that
prints integers.
I am missing the point of the reciprocal multiplication.
If you are looking for 10 times the velocity, rounded to the
nearest integer (ie nearest 10th of the raw velocity),
int(10 * km * 600 / time + 0.5)
where int() truncates without rounding?
What am I missing?
Nothing, well except for the fact that you cannot add 0.5 when you're working
with integers only, and that's what the code below is doing:

mov eax, dword ptr [esi + offset dtv.km] // km are already * 10
mov ecx, 600
mul ecx

mov ecx, dword ptr [esi + offset dtv.time] // minutes
div ecx

add eax, 5
mov ecx, $cccccccd
mul ecx
shr edx, 3

Robert
--
Robert AH Prins
robert(a)prino(d)org
Robert Prins
2015-03-30 23:04:27 UTC
Permalink
Post by Bob Masta
On Sun, 29 Mar 2015 21:52:27 +0000, Robert Prins
Post by Robert Prins
I've managed to remove part of the code that deals with doubles, the original
code would calculate
km * 6.0 / time
and print the resulting double. I've now changed that into
using the upper 32 bits of
((km * 600 / time + 5) * 3435973837) shr 3,
which, through reciprocal multiplication (3435973837 = 2^35 / 10), gives me 10
times the velocity, properly truncated and that is passed to a routine that
prints integers.
I am missing the point of the reciprocal multiplication.
If you are looking for 10 times the velocity, rounded to the
nearest integer (ie nearest 10th of the raw velocity),
int(10 * km * 600 / time + 0.5)
where int() truncates without rounding?
What am I missing?
As I already said...

However, some further thought about your simple solution led me to change my code to

push 0 // placeholder for dtv.v as longint

cmp dword ptr [esi + offset dtv.time], 0
je @01

fild dword ptr [esi + offset dtv.km]
fmul d_6v0
fidiv dword ptr [esi + offset dtv.time]
fst qword ptr [esi + offset dtv.v] // I (still) need the double
fmul d_10v0
fadd d_0v5

{$ifdef sse3}
db $db,$0c,$24 // fisttp dword [esp] VP doesn't have fisttp
{$else}
call System._Trunc
mov [esp], eax
{$endif}

@01:

And that is even better and shorter than what I had. Happy bunny!

Robert
--
Robert AH Prins
robert(a)prino(d)org
Robert AH Prins
2015-04-28 21:54:33 UTC
Permalink
Post by Robert Prins
Post by Bob Masta
On Sun, 29 Mar 2015 21:52:27 +0000, Robert Prins
Post by Robert Prins
I've managed to remove part of the code that deals with doubles, the original
code would calculate
km * 6.0 / time
and print the resulting double. I've now changed that into
using the upper 32 bits of
((km * 600 / time + 5) * 3435973837) shr 3,
which, through reciprocal multiplication (3435973837 = 2^35 / 10), gives me 10
times the velocity, properly truncated and that is passed to a routine that
prints integers.
I am missing the point of the reciprocal multiplication.
If you are looking for 10 times the velocity, rounded to the
nearest integer (ie nearest 10th of the raw velocity),
int(10 * km * 600 / time + 0.5)
where int() truncates without rounding?
What am I missing?
As I already said...
However, some further thought about your simple solution led me to change my code to
push 0 // placeholder for dtv.v as longint
cmp dword ptr [esi + offset dtv.time], 0
fild dword ptr [esi + offset dtv.km]
fmul d_6v0
fidiv dword ptr [esi + offset dtv.time]
fst qword ptr [esi + offset dtv.v] // I (still) need the double
fmul d_10v0
fadd d_0v5
{$ifdef sse3}
db $db,$0c,$24 // fisttp dword [esp] VP doesn't have fisttp
{$else}
call System._Trunc
mov [esp], eax
{$endif}
And that is even better and shorter than what I had. Happy bunny!
I finally solved the problem, the Pascal output files still the ones coming from
PL/I, and I've managed to get rid of the (slow) float-to-string conversion
routine from the RTL, replacing it with the integer-to-string version. For those
interested:

Common:

const d_rmul: array[0..2] of double = (1.0, 10.0, 100.0);
const d_radd: array[0..2] of double = (0.5, 0.05, 0.005);
var _line : string[255];
var print : string[19];

The "print_2_line" procedure just concatenates the contents of "print" to
"_line" and adds ' | '.

Original Pascal:

procedure add_double_2_line(var _double: double; _i, _j: longint);
var _d: double;

begin
_d:= _double;

str(frac(_double * d_rmul[_j]):6:4, print);

if print[3] >= '5' then
_d:= _double + d_radd[_j];

str(_d:_i:_j, print);
add_print_2_line;
end; {add_double_2_line}

Original assembler:

procedure add_double_2_line(var _double: double; _i, _j: longint); assembler;
{&uses none} {&frame-}
asm
//a-in add_double_2_line
mov eax, _double
mov edx, _i

fld qword ptr [eax]
mov eax, _j

fld st

fmul qword ptr d_rmul[eax * 8]

fld st
push ecx
fisttp dword [esp]
fild dword ptr [esp]
fsubp st(1), st
pop ecx

push 6
push 4
push offset print
push type print - 1
call System._StrFlt

cmp byte ptr print[3], "5"
jb @01

fadd qword ptr d_radd[eax * 8]

@01:
movzx ecx, byte ptr _line
push edx // _i
add edx, 3
add byte ptr _line, dl
add ecx, offset _line

push eax // _j
push ecx
push type _line - 1
call System._StrFlt

mov byte ptr [ecx], " "
mov dword ptr [ecx + edx - 2], " ³ "
//a-out
end; {add_double_2_line}

Current routine:

procedure add_double_2_line(var _double: double; _i, _j: longint); assembler;
{&uses none} {&frame+}
asm
//a-in add_double_2_line

mov eax, _double
fld qword ptr [eax]
fld st

// if (trunc(frac(_d * d_rmul[_j]) * 10000 + 0.5) >= 5000) and
// (trunc(frac(_d * d_rmul[_j]) * 10000 + 0.5) < 10000) then

mov ecx, _j
fmul qword ptr d_rmul[ecx * 8]

fld st
push eax
db $db,$0c,$24 // fisttp dword [esp]
fild dword ptr [esp]
fsubp st(1),st
fmul d_10000v0
fadd d_0v5

db $db,$0c,$24 // fisttp dword [esp]
pop eax
cmp eax, 5000
jl @01

cmp eax, 10000
jge @01

// _d:= _d + d_radd[_j]

fadd qword ptr d_radd[ecx * 8]

@01:
fmul qword ptr d_rmul[ecx * 8]
fadd d_0v5

push eax // value to be printed
db $db,$0c,$24 // fisttp dword [esp]

mov edx, _i
push edx // width

movzx eax, byte ptr _line
add edx, 3
add byte ptr _line, dl
add eax, offset _line

push eax // position in _line
push type _line - 1 // length of receiving field, big enough
call System._StrInt

mov byte ptr [eax], " " // blank out length set by StrInt
mov dword ptr [edx + eax - 2], " | "
test ecx, ecx // no further processing if rounded to 0 decimals
jz @03

// The following code moves the value one position to the left and
// adds a decimal point, and converts blanks into significant zeros

neg ecx
lea ecx, [ecx + edx - 3]

@02:
mov dl, [eax + 1]
mov [eax], dl
inc eax
dec ecx
jnz @02

mov byte ptr [eax], "."
or dword ptr [eax - 1], $00300030 // "#0 0 #0 0"

@03:
//a-out
end; {add_double_2_line}

Still not entirely happy about it, especially the small loop at the end, but
given that at most 5 characters have to be moved, soit...

Robert
--
Robert AH Prins
robert(a)prino(d)org
No programming @ http://hitchwiki.org/prino/
Waldek Hebisch
2015-03-18 01:27:39 UTC
Permalink
Post by Robert AH Prins
I can divide by 10 or multiply with some magic, shift, multiply
and subtract, but is there a shorter way of determining that a
number ends in the digit 5 - I'm not interested in divisibility
by 5, I just need to know that it ends in 5.
Number ends in 4 when it is divisible by 5 and not divisible by 2.
Divisibility by 2 requirs testing single bit, so it is enough
to handle divisibility by 5. For this you can multiply by
modular inverse of 5 and test if result is in aproprate range.
For 32 bit numbers the inverse is a = 3435973837. When dealing
with unsigned numbers you need just to check that a*x <= 2^32/5,
if yes number is divisible. For signed you need to test
-2^31/5 <= a*x <= (2^31 - 1)/5. Note: this uses 32-bit
unsigned multiplication (truncating higher bits).

This requres single multiplication by magic number, comparison
(or two in signed case) and a bit test (for divisibility by 2).
Seem to be a bit shorter than division with remainder done
via multiplication by the inverse.
--
Waldek Hebisch
***@math.uni.wroc.pl
Catalin George Festila
2015-03-20 17:53:57 UTC
Permalink
I don't have time to give you one good replay , but try to use 64 bit.
Hardware make sense ...
Robert Prins
2015-03-20 20:51:20 UTC
Permalink
Post by Catalin George Festila
I don't have time to give you one good replay , but try to use 64 bit.
Hardware make sense ...
Impossible, Virtual Pascal is 32-bit.

Robert
--
Robert AH Prins
robert(a)prino(d)org
Bob Masta
2015-03-21 13:15:17 UTC
Permalink
On Sun, 15 Mar 2015 19:49:44 +0000, Robert AH Prins
Post by Robert AH Prins
I can divide by 10 or multiply with some magic, shift, multiply and subtract,
but is there a shorter way of determining that a number ends in the digit 5 -
I'm not interested in divisibility by 5, I just need to know that it ends in 5.
Why? For some input a program (or rather the RTL added by the compiler) has
trouble correctly rounding, no matter what I use for the FPU control word, so I
don't have any other choice than doing it myself.
Instead of examining the last digit, I wonder if there is a
simpler way based on the old "add 1/2 LSD and truncate"
method of rounding. I assume that's not adequate for your
needs here, but could it be the basis for something that
would work for you? Just trying to explore all the
options... <g>

Best regards,


Bob Masta

DAQARTA v7.60
Data AcQuisition And Real-Time Analysis
www.daqarta.com
Scope, Spectrum, Spectrogram, Sound Level Meter
Frequency Counter, Pitch Track, Pitch-to-MIDI
FREE Signal Generator, DaqMusiq generator
Science with your sound card!
Continue reading on narkive:
Loading...