Sdcb.Arithmetic.Mpfr 4.2.0

Sdcb.Arithmetic

Sdcb.Arithmetic is a modern .NET library that can PInvoke to gmp and mpfr, that enable both high performance and best .NET convenience.

Known classes in Sdcb.Arithmetic:

NuGet Packages

libgmp

Package Id	Version	License	Notes
Sdcb.Arithmetic.Gmp		MIT	.NET binding for libgmp
Sdcb.Arithmetic.Gmp.runtime.win64		LGPL	native lib in windows x64
Sdcb.Arithmetic.Gmp.runtime.win32		LGPL	native lib in windows x86
Sdcb.Arithmetic.Gmp.runtime.linux64		LGPL	native lib in Linux x64

mpfr

Package Id	Version	License	Notes
Sdcb.Arithmetic.Mpfr		MIT	.NET binding for libmpfr
Sdcb.Arithmetic.Mpfr.runtime.win64		LGPL	native lib in windows x64
Sdcb.Arithmetic.Mpfr.runtime.win32		LGPL	native lib in windows x86
Sdcb.Arithmetic.Mpfr.runtime.linux64		LGPL	native lib in linux x64

Question - Why not linux-64 can't run in my linux?

The linux-x64 package is compiled using vcpkg in Ubuntu 22.04, so it may not run in other linux distributions.

If you want to run in other linux distributions, you can compile it yourself, or install it using apt or yum or other package manager.

This is the all native dynamic library name in case you wondering(defined in GmpNativeLoader.cs and MpfrNativeLoader.cs):

Examples

Calculate 1,000,000 length of π using Sdcb.Arithmetic.Gmp:

// Install NuGet package: Sdcb.Arithmetic.Gmp
// Install NuGet package: Sdcb.Arithmetic.Gmp.runtime.win-x64(for windows)
using Sdcb.Arithmetic.Gmp;

Console.WriteLine(CalcPI().ToString("N1000000"));

GmpFloat CalcPI(int inputDigits = 1_000_000)
{
    const double DIGITS_PER_TERM = 14.1816474627254776555; // = log(53360^3) / log(10)
    int DIGITS = (int)Math.Max(inputDigits, Math.Ceiling(DIGITS_PER_TERM));
    uint PREC = (uint)(DIGITS * Math.Log2(10));
    int N = (int)(DIGITS / DIGITS_PER_TERM);
    const int A = 13591409;
    const int B = 545140134;
    const int C = 640320;
    const int D = 426880;
    const int E = 10005;
    const double E3_24 = (double)C * C * C / 24;

    using PQT pqt = ComputePQT(0, N);

    GmpFloat pi = new(precision: PREC);
    // pi = D * sqrt((mpf_class)E) * PQT.Q;
    pi.Assign(GmpFloat.From(D, PREC) * GmpFloat.Sqrt((GmpFloat)E, PREC) * (GmpFloat)pqt.Q);
    // pi /= (A * PQT.Q + PQT.T);
    GmpFloat.DivideInplace(pi, pi, GmpFloat.From(A * pqt.Q + pqt.T, PREC));
    return pi;

    PQT ComputePQT(int n1, int n2)
    {
        int m;

        if (n1 + 1 == n2)
        {
            PQT res = new()
            {
                P = GmpInteger.From(2 * n2 - 1)
            };
            GmpInteger.MultiplyInplace(res.P, res.P, 6 * n2 - 1);
            GmpInteger.MultiplyInplace(res.P, res.P, 6 * n2 - 5);

            GmpInteger q = GmpInteger.From(E3_24);
            GmpInteger.MultiplyInplace(q, q, n2);
            GmpInteger.MultiplyInplace(q, q, n2);
            GmpInteger.MultiplyInplace(q, q, n2);
            res.Q = q;

            GmpInteger t = GmpInteger.From(B);
            GmpInteger.MultiplyInplace(t, t, n2);
            GmpInteger.AddInplace(t, t, A);
            GmpInteger.MultiplyInplace(t, t, res.P);
            // res.T = (A + B * n2) * res.P;            
            if ((n2 & 1) == 1) GmpInteger.NegateInplace(t, t);
            res.T = t;

            return res;
        }
        else
        {
            m = (n1 + n2) / 2;
            PQT res1 = ComputePQT(n1, m);
            using PQT res2 = ComputePQT(m, n2);
            GmpInteger p = res1.P * res2.P;
            GmpInteger q = res1.Q * res2.Q;

            // t = res1.T * res2.Q + res1.P * res2.T
            GmpInteger.MultiplyInplace(res1.T, res1.T, res2.Q);
            GmpInteger.MultiplyInplace(res1.P, res1.P, res2.T);
            GmpInteger.AddInplace(res1.T, res1.T, res1.P);
            res1.P.Dispose();
            res1.Q.Dispose();
            return new PQT
            {
                P = p,
                Q = q,
                T = res1.T,
            };
        }
    }
}

public ref struct PQT
{
    public GmpInteger P;
    public GmpInteger Q;
    public GmpInteger T;

    public readonly void Dispose()
    {
        P?.Dispose();
        Q?.Dispose();
        T?.Dispose();
    }
}

Technical notes

Why choosing struct in class design instead of raw memory IntPtr design?

• (1) Struct in class design:

  class GmpInteger
  {
    public readonly Mpz_t Raw;
  
    public unsafe void DoWork()
    {
        fixed (Mpz_t* ptr = &Raw)
        {
            GmpLib.__dowork((IntPtr)ptr);
        }
    }
  }
  
  struct Mpz_t
  {
    public int A, B;
    public IntPtr Limbs;
  }
  

• (2) Raw memory IntPtr design:

  class GmpInteger : IDisposable
  {
    public readonly IntPtr Raw;
  
    public unsafe GmpInteger()
    {
        Raw = Marshal.AllocHGlobal(sizeof(Mpz_t));
    }
  
    public void DoWork()
    {
        GmpLib.__dowork(Raw);
    }
  
    public void Dispose()
    {
        Marshal.FreeHGlobal(Raw);
    }
  }
  

  Here is some benchmark I tested for both DoWork senario and initialize-dispose senario:

  Details:

  • init & dispose combines following actions:
    • allocating struct memory
    • calling mpz_init
    • calling mpz_free
    • free the memory
    • Measure the operations-per-seconds, higher ops is better
  • dowork contains following actions:
    • create a GmpFloat to 1.5 with precision=1000(v = 1, a = 1.5)
    • Calling MultiplyInplace(v *= a) 10*1024*1024 times
    • Measure the duration, lower is better

  Here is the tested results in my laptop:

  case/senario	init & dispose	dowork
  Struct in class	82,055,792 ops	1237ms
  Raw memory IntPtr	15,543,619 ops	1134ms

  As you can see, raw memory IntPtr design will benifits ~8.33% faster in dowork senario above, but struct in class design will be 5.2x faster in init & dispose senario.

  Finally I choosed the struct in class design, here is some existing Raw memory IntPtr design work if you also wants to check or test:

  • branch: https://github.com/sdcb/Sdcb.Arithmetic/tree/feature/gmp-raw-ptr
  • nuget-package: https://www.nuget.org/packages/Sdcb.Arithmetic.Gmp/1.0.10-preview.12

  Struct in class performance results environment:

  • commit: https://github.com/sdcb/Sdcb.Arithmetic/tree/976d4271c487a554be936cf56ed55f2b1314042e
  • nuget-package: https://www.nuget.org/packages/Sdcb.Arithmetic.Gmp/1.0.10-preview.11

  In the future, Raw memory IntPtr design can be pick-up if a handy, good performance memory allocator was found.