Free42 with IEEE 754-2008 decimal floating-point -- interested in a sneak preview?

[Paul Dale]

(03-16-2014 01:45 PM)Werner Wrote:  in DEG mode:
SIN(180°+x) = -SIN(x) fails for x=10^-31, 10^-30, 10^-16 ... and many more (these fail on the WP34 as well btw)

I'm not surprised here. 180 + 10^-31 is pushing the number of digits in double precision whereas -10^-31 isn't stretching anything. In single precision there doesn't seem to be an issue thankfully.

I know the 34S trigonometric functions don't do the usual full reduction to the first octant before evaluating the functions. I brute force things with a longer series -- it is smaller codewise. There is also the conversion to radians involved to disturb the results a bit.

In summary:

* not unexpected,
* doesn't impact single precision,
* will take space to fix.

Is it worthwhile trying? I'll see if I can come up with something nice and simple.


- Pauli

[Thomas Okken]

(03-17-2014 11:08 AM)J-F Garnier Wrote:  I noted that there is also a similar issue in the reverse situation:
1E5 LOG 5 - --> 1E-33 (also for 1E40 and other)
This didn't happen in versions 1.4.7x

OK, I'll add code to handle powers of 10 in LOG. (I'm amazed that code isn't there already, since testing for exact powers of 10 is so easy in decimal. Continuity issues maybe?)

[Thomas Klemm]

(03-16-2014 03:06 PM)Thomas Okken Wrote:  By the way, where can I find out what formatting is supported on this forum? I thought it was strictly plain text here until I saw your [b] tags. :-)
(UPDATE: Duh, Help -> Posting -> MyCode. Never mind!)

Besides the obvious buttons of the editor there are some hidden features:

• In addition to the aforementioned tags sub and sup have been added: a² or x_i
• But of course you can achieve the same result using LaTeX: \(a^2\) or \(x_i\)
• And then we have Katie's compilation of smilies.
  

Cheers
Thomas

[Thomas Okken]

(03-17-2014 01:29 PM)Thomas Okken Wrote:  

(03-17-2014 11:08 AM)J-F Garnier Wrote:  I noted that there is also a similar issue in the reverse situation:
1E5 LOG 5 - --> 1E-33 (also for 1E40 and other)
This didn't happen in versions 1.4.7x

OK, I'll add code to handle powers of 10 in LOG. (I'm amazed that code isn't there already, since testing for exact powers of 10 is so easy in decimal. Continuity issues maybe?)

In order to make sure LOG doesn't just return exact results for exact powers of 10, but also satisfies

  10ⁿ ≤ x < 10ⁿ⁺¹ → n ≤ LOG(x) < n+1

which is necessary for code like FWIW to work reliably, I'm thinking of this argument reduction code:

Code:

function log10(x) {

  // These cases aren't relevant in Free42, since it never calls log10 with x ≤ 0, but still.
  if (x < 0)
    return NaN;
  if (x == 0)
    return -Infinity;

  // Extract normalized decimal mantissa and exponent
  // (make sure 1 ≤ mantissa < 10)
  m = decimal mantissa of x;
  e = decimal exponent of x;
  while (m ≥ 10) {
    m = m / 10;
    e = e + 1;
  }
  while (m < 1) {
    m = m * 10;
    e = e - 1;
  }

  // Handle the trivial case
  if (m == 1)
    return e;

  // Handle the general case, while enforcing
  // 10^n ≤ x < 10^(n+1) → n ≤ LOG(x) < n+1
  r = e + bid128_log10(m);
  if (r < e)
    return e;
  else if (r ≥ e + 1)
    return highest representable number less than e + 1;
  else
    return r;
}

[Thomas Okken]

(03-17-2014 10:03 PM)Thomas Klemm Wrote:  Besides the obvious buttons of the editor there are some hidden features: [...]

Thanks! :-)
BTW, the "obvious buttons" weren't obvious to me because I was using Firefox with NoScript at home, and NoScript nixes those buttons, unless you tell it not to. My mistake!

- Thomas

[Werner]

(03-17-2014 01:29 PM)Thomas Okken Wrote:  In order to make sure LOG doesn't just return exact results for exact powers of 10, but also satisfies

  10ⁿ ≤ x < 10ⁿ⁺¹ → n ≤ LOG(x) < n+1

This is not possible.
On a real 42S, LOG(9.99999999999) = 1, exactly. Meaning it is the closest number to the exact result. The same holds for any precision. On Free42 Decimal-34, LOG(1e34-1) is 34 and LOG(1e33-1) is 33, exactly, and that's what the result should be.
Extracting the exponent should rather be done dividing by the mantissa, that can be determined as follows (thanks to Dieter):

Code:


00 { 25-Byte Prgm }
01*LBL "MANT"
02 ABS
03 1        avoid upper exponent limit. ONLY needed for 9.99999999999e499
04 X<Y?
05 10^X
06 /
07 ENTER
08 X#0?
09 LOG
10 INT
11 10^X
12 /
13 1        multiply by 10 if needed
14 X>Y?
15 10^X
16 *
17 END

Werner

[J-F Garnier]

(03-18-2014 07:16 AM)Werner Wrote:  

(03-17-2014 01:29 PM)Thomas Okken Wrote:    10ⁿ ≤ x < 10ⁿ⁺¹ → n ≤ LOG(x) < n+1

This is not possible.
On a real 42S, LOG(9.99999999999) = 1, exactly. Meaning it is the closest number to the exact result. The same holds for any precision. On Free42 Decimal-34, LOG(1e34-1) is 34 and LOG(1e33-1) is 33, exactly, and that's what the result should be.

I tend to agree with Werner, and the condition above is not needed. Strictly speaking, all we need is that 10^x is exact for x an integer, and it will be nice that LOG(x) is an integer for x exactly a power of ten.

J-F

[Werner]

(03-17-2014 11:40 AM)Paul Dale Wrote:  I know the 34S trigonometric functions don't do the usual full reduction to the first octant before evaluating the functions. I brute force things with a longer series -- it is smaller codewise. There is also the conversion to radians involved to disturb the results a bit.

In summary:

* not unexpected,
* doesn't impact single precision,
* will take space to fix.

Is it worthwhile trying? I'll see if I can come up with something nice and simple

- Pauli

Strange that singe precision works, if you don't do the range reduction to the first octant. For me personally, these things are much more important than having 34 digits most of the time. Accuracy over precision, anytime. And of course the venerable 42S gets it right. And Free42, latest version.

range reduction to 45°

Werner

[Paul Dale]

It isn't that strange that single precision works properly here: internally I always use a higher precision regardless of the single/double precision mode setting. A bit of brute force here saves a lot of code space and thinking time.

I agree accuracy over precision too, that has been my goal all along. I'm not sure we can justify the space to do the octant reduction separately in degrees, radians and gradians and I suspect doing this after the conversion to radians has occurred will result in this problem.

Fortunately, it is documented that double precision results can be incorrect in the latter digits for some (most?) functions Double precision was a bit of an after thought which allowed keystroke programs to produce reasonably accurate single precision results -- at least that was its intention.


- Pauli

[Werner]

In that case, might I humbly request that you do consider it for the 43S?
Werner

[Paul Dale]

(03-18-2014 09:22 AM)Werner Wrote:  In that case, might I humbly request that you do consider it for the 43S?

That sounds more than reasonable. We'll have space there. We will have space on the 31S as well so I guess I'll have to implement it.

I've had a bit of an idea about how to do this without using lots of space. Not sure if it will work or not yet...


- Pauli

[Thomas Okken]

(03-18-2014 07:16 AM)Werner Wrote:  On a real 42S, LOG(9.99999999999) = 1, exactly. Meaning it is the closest number to the exact result. The same holds for any precision. On Free42 Decimal-34, LOG(1e34-1) is 34 and LOG(1e33-1) is 33, exactly, and that's what the result should be.

Good point. I don't want to return incorrect results just for the convenience of a flawed algorithm!
OK, I'll fix LOG like this: use the bid128_ilogb and bid128_scalbn functions to pull out the exponent and normalize the mantissa to the range [1, 10), then use bid128_log10 on the normalized mantissa, and finally add the exponent to the result. That will return exact results for powers of ten, without introducing numerical inconsistencies.

[Paul Dale]

The 31S will have appropriate range reduction code in degrees mode but not in radians (it doesn't support gradians). The code space for this is likely more than we want to commit to on the 34S but I coded the gradians modulo reduction too for future use.


- Pauli

(03-16-2014 01:45 PM)Werner Wrote:  in DEG mode:
SIN(180°+x) = -SIN(x) fails for x=10^-31, 10^-30, 10^-16 ... and many more (these fail on the WP34 as well btw)
UPDATE (Free42, but not WP34): this is probably due to 1e-31 - 10^-31 not being zero.. in fact 10^x does not produce the correct result for x= integer; even -1..)
SIN(45°) ^= COS(45°)
I can't find the range reduction code straight away, to see what would be wrong with it.. can you point me to it?
Cheers, Werner