REAL16 for statistics and Gamma function

[bigbadbob]

I was wondering how to do this in MASM.  It is hard to do this in any language because you have to mask it and do a lot of manipulation.
I also think that it is even harder to convert from BINARY FLOATING POINT to decimal notation as it means converting it to a string.

I have an algorithm for calculating Gamma that is currently written in a high level language.  It does not have an accurate answer in REAL8 (double) because of rounding error.  I think that calculating it in REAL16 will reduce the round off error, because I could later convert the result to REAL8 dropping the inaccurate bits.  Not really sure but I know that this will significantly slow down the calculation.

This Gamma function keeps multiplying until the iteration stopped changing the delta or 100 iterations.

Code Select Expand


        public static double Gamma(double z)
        {
            double factor = 1.0;
            double delta = 1.0;
            for (int n = 1; n <= 100 && delta > 0.0; n++)
            {
                double previous = factor;

                factor *= (1.0 / (1.0 + (z / n))) * Math.Pow((1.0 + (1.0 / n)), z);

                delta = Math.Abs(factor - previous);
            }

            return factor / z;
        }


This Gamma function starts from 0 increments by some small number and stops at a fake value of infinity called approxInfinity.

Code Select Expand


        public static double GammaIntegral(double z, double approxInfinity, double increment)
        {
            double series = 0.0;
            double delta = 1.0;

            for (double x = 0; x <= approxInfinity; x += increment)
            {
                series += (Math.Pow(x, z - 1) * Math.Exp(-x))*increment;
            }

            return series;
        }


I've passed a few numbers for a single value of z.  Up to the first few decimal points they match, but different parameters yield different results.
I believe that the round off error because of not enough bits in the mantissa is part of the problem.

[jj2007]

If REAL8 is not enough, try REAL10 - that is the native precision of the FPU. Ten bits of precision more.

We have a dozen threads about REAL16 precision, though. Forum search is here. The available REAL16 libraries are much slower than the FPU, of course.

[nidud]

deleted

[nidud]

deleted

[jj2007]

Check again:

Code Select Expand

    fld   tbyte ptr [rdx]
    fld   tbyte ptr [r8]
    fmul  ; NO: st(0),st(1)
    fstp  tbyte ptr [rcx]

No need for ffree.

[nidud]

deleted

[felipe]

So what happens then?
fmul clear out st(1) and fstp st(0) ?

This is from raymond's tutorial on the fpu:

Note: FMUL without operands can also be used with the MASM assembler but such an instruction is coded as FMULP ST(1),ST(0).

:idea:

[felipe]

I also think that it is even harder to convert from BINARY FLOATING POINT to decimal notation as it means converting it to a string.

Don't know if you are going to use the real16 or the real10 types, but for real10 this convertion isn't so difficult and there are already functions (binary fp to ascii and viceversa) in the fpulib for this,already available. You can use it as is or you can check the source and make an custom adaption.  :idea:

[jj2007]

So what happens then?
fmul clear out st(1) and fstp st(0) ?

Yes, more or less. Raymond's tutorial is a great source for that.
In short: fmul stands for "multiply ST(0) with ST(1), discard ST(1) and put the result in ST(0)".

Code Select Expand

    ; The args are REAL10, so this is converted to REAL16 below.
    ;
    ; void test( result:ptr, a:ptr REAL10, b:ptr REAL10 );

result:ptr is what? ptr to REAL16? How do you print that? No criticism, just curious...

[nidud]

deleted

[jj2007]

You have to link with the LIBC library and use the stdio functions.

I downloaded the required inc files, it builds fine but

Code Select Expand

POLINK: error: Unresolved external symbol 'printf'.
POLINK: fatal error: 1 unresolved external(s).

Is there an installer for AsmC that takes care of such problems?

[aw27]

Of course, people can have a lot of fun arguing whether REAL10 is better than REAL8 or whether REAL16 from software is fast enough.

But for serious stuff, people should use an arbitrary precision library, in this case MPIR is the best alternative for Windows and can be easily used from ASM as I have shown here http://masm32.com/board/index.php?topic=6438.0
All right, it was done for UASM but can be easily converted to MASM.

I made a little change to that code and put it calculating Gamma(3000) and Gamma (1/2)
For people that has no idea, Gamma (3000) is the same as 2999!, a really huge number. And Gamma(1/2) is the same as square root of PI.

Let's see a picture, yes does not fit it all in the capture, we are using 10000 decimal digits precision for the calculation of Gamma (3000) . I selected only 100 decimal digits precision for Gamma (1/2).

MPIR has a function to directly calculate Gamma, that's what I have used, alternatively we could have made the calculation using one of the many existing formulas.

[nidud]

deleted

[jj2007]

LIBC is a static library supplied by a specific tool chain

Sure. Is there a download link for libc.lib, or do I have to fight some hours to build a source with another "specific tool chain"?