News:

Masm32 SDK description, downloads and other helpful links
Message to All Guests

Main Menu

REAL16 for statistics and Gamma function

Started by bigbadbob, June 18, 2018, 02:44:37 PM

Previous topic - Next topic

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.

        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.

        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

#2
deleted

nidud

#3
deleted

jj2007

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


No need for ffree.

nidud

#5
deleted

felipe

Quote from: nidud on June 19, 2018, 04:14:13 AM
So what happens then?
fmul clear out st(1) and fstp st(0) ?

This is from raymond's tutorial on the fpu:
QuoteNote: 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

Quote from: bigbadbob on June 18, 2018, 02:44:37 PM
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

Quote from: nidud on June 19, 2018, 04:14:13 AM
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)".

    ; 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

#9
deleted

jj2007

Quote from: nidud on June 19, 2018, 09:50:13 AMYou have to link with the LIBC library and use the stdio functions.
I downloaded the required inc files, it builds fine but
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

#11
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

#12
deleted

jj2007

Quote from: nidud on June 20, 2018, 12:50:19 AMLIBC 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"?