Nth digit of Pi

[Juggler_IN]

I am not getting the desired outcome after the translation of a c code to VBA for finding n'th decimal digit of pi. Attached is the VBA and the reference C code.

Any suggestion as to where it is failing? (Few of the support codes are not an exact translation of c.)

VBA Code

Function fmod(a, b)
' /* return a mod b */

fmod = a - (b * (a \ b))

End Function
Function inv_mod(x, y)
' /* return the inverse of x mod y */

    Dim q, u, v, a, c, t

    u = x
    v = y
    c = 1
    a = 0

    Do
        q = v / u
        t = c
        c = a - q * c
        a = t
        t = u
        u = v - q * u
        v = t
    Loop While (u <> 0)

    a = a Mod y
    If (a < 0) Then
        a = y + a
    End If
   
    inv_mod = a

End Function
Function pow_mod(a, b, m)
' /* return (a^b) mod m */
   
    Dim y

    y = 1

    Do While b > 0
        y = (y * a) Mod m
        b = b - 1
    Loop

    pow_mod = y

End Function
Public Function mul_mod(a, b, m)
' /* return (a*b) mod m */
   
    Dim u

    u = 0
    a = a Mod m

    While (b > 0)
        If (b Mod 2 = 1) Then u = (u + a) Mod m
        a = (a * 2) Mod m
        b = b / 2
    Wend

    mul_mod = u Mod m

End Function
Public Function is_prime(x)
' /* return true if n is prime */
       
        Dim i

        If x = 2 Then
            is_prime = True
        ElseIf x <= 1 Or x Mod 2 = 0 Then
            is_prime = False
        Else
            is_prime = True
            For i = 3 To Int(Sqr(x)) Step 2
                If x Mod i = 0 Then
                    is_prime = False
                    Exit For
                End If
            Next i
        End If

End Function
Function next_prime(n)
' /* return the prime number immediatly after n */

    Do
        n = n + 1
    Loop While Not is_prime(n)
   
    next_prime = n

End Function
Sub main()
' Note: I have modified capital N in c code to m; lowercase n stays as is.
    Dim av, a, vmax, m, n, num, den, k, kq, kq2, t, v, s, i, sum, h

    n = 10

    m = Int(((n + 20) * Log(10) / Log(2)))

    sum = 0

    For a = 3 To (2 * m) Step next_prime(a)

        vmax = Int((Log(2 * m) / Log(a)))
        av = 1
        For i = 0 To vmax - 1
            av = av * a
        Next

        s = 0
        num = 1
        den = 1
        v = 0
        kq = 1
        kq2 = 1

        For k = 1 To m - 1

            t = k
            If (kq >= a) Then
                Do
                    t = t / a
                    v = v - 1
                Loop While ((t Mod a) = 0)
                kq = 0
            End If
            kq = kq + 1
            num = mul_mod(num, t, av)

            t = (2 * k - 1)
            If (kq2 >= a) Then
                If (kq2 = a) Then
                    Do
                        t = t / a
                        v = v + 1
                    Loop While ((t Mod a) = 0)
                End If
                kq2 = kq2 - a
            End If
            den = mul_mod(den, t, av)
            kq2 = kq2 + 2

            If (v > 0) Then
                t = inv_mod(den, av)
                t = mul_mod(t, num, av)
                t = mul_mod(t, k, av)
                For i = v To vmax - 1
                    t = mul_mod(t, a, av)
                Next i
                s = s + t
                If (s >= av) Then
                    s = s - av
                End If
            End If

        Next

        t = pow_mod(10, n - 1, av)
        s = mul_mod(s, t, av)
        sum = fmod(sum + s / av, 1#)

    Next

    Debug.Print Int(sum * 1000000000#)

End Sub

And, reference C Code:

/*
 * Computation of the n'th decimal digit of pi with very little memory.
 * Written by Fabrice Bellard on February 26, 1997.
 *
 * We use a slightly modified version of the method described by Simon
 * Plouffe in "On the Computation of the n'th decimal digit of various
 * transcendental numbers" (November 1996). We have modified the algorithm
 * to get a running time of O(n^2) instead of O(n^3log(n)^3).
 *
 * This program uses a variation of the formula found by Gosper in 1974 :
 *
 * pi = sum( (25*n-3)/(binomial(3*n,n)*2^(n-1)), n=0..infinity);
 *
 * This program uses mostly integer arithmetic. It may be slow on some
 * hardwares where integer multiplications and divisons must be done by
 * software. We have supposed that 'int' has a size of at least 32 bits. If
 * your compiler supports 'long long' integers of 64 bits, you may use the
 * integer version of 'mul_mod' (see HAS_LONG_LONG).
 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/* uncomment the following line to use 'long long' integers */
/* #define HAS_LONG_LONG */

#ifdef HAS_LONG_LONG
#define mul_mod(a,b,m) (( (long long) (a) * (long long) (b) ) % (m))
#else
#define mul_mod(a,b,m) fmod( (double) a * (double) b, m)
#endif

/* return the inverse of x mod y */
int inv_mod(int x, int y)
{
    int q, u, v, a, c, t;

    u = x;
    v = y;
    c = 1;
    a = 0;
    do {
    q = v / u;

    t = c;
    c = a - q * c;
    a = t;

    t = u;
    u = v - q * u;
    v = t;
    } while (u != 0);
    a = a % y;
    if (a < 0)
    a = y + a;
    return a;
}

/* return the inverse of u mod v, if v is odd */
int inv_mod2(int u, int v)
{
    int u1, u3, v1, v3, t1, t3;

    u1 = 1;
    u3 = u;

    v1 = v;
    v3 = v;

    if ((u & 1) != 0) {
    t1 = 0;
    t3 = -v;
    goto Y4;
    } else {
    t1 = 1;
    t3 = u;
    }

    do {

    do {
        if ((t1 & 1) == 0) {
        t1 = t1 >> 1;
        t3 = t3 >> 1;
        } else {
        t1 = (t1 + v) >> 1;
        t3 = t3 >> 1;
        }
      Y4:;
    } while ((t3 & 1) == 0);

    if (t3 >= 0) {
        u1 = t1;
        u3 = t3;
    } else {
        v1 = v - t1;
        v3 = -t3;
    }
    t1 = u1 - v1;
    t3 = u3 - v3;
    if (t1 < 0) {
        t1 = t1 + v;
    }
    } while (t3 != 0);
    return u1;
}

/* return (a^b) mod m */
int pow_mod(int a, int b, int m)
{
    int r, aa;

    r = 1;
    aa = a;
    while (1) {
    if (b & 1)
        r = mul_mod(r, aa, m);
    b = b >> 1;
    if (b == 0)
        break;
    aa = mul_mod(aa, aa, m);
    }
    return r;
}

/* return true if n is prime */
int is_prime(int n)
{
    int r, i;
    if ((n % 2) == 0)
    return 0;

    r = (int) (sqrt(n));
    for (i = 3; i <= r; i += 2)
    if ((n % i) == 0)
        return 0;
    return 1;
}

/* return the prime number immediatly after n */
int next_prime(int n)
{
    do {
    n++;
    } while (!is_prime(n));
    return n;
}

#define DIVN(t,a,v,vinc,kq,kqinc)        \
{                        \
  kq+=kqinc;                    \
  if (kq >= a) {                \
    do { kq-=a; } while (kq>=a);        \
    if (kq == 0) {                \
      do {                    \
    t=t/a;                    \
    v+=vinc;                \
      } while ((t % a) == 0);            \
    }                        \
  }                        \
}

int main(int argc, char *argv[])
{
    int av, a, vmax, N, n, num, den, k, kq1, kq2, kq3, kq4, t, v, s, i, t1;
    double sum;

    if (argc < 2 || (n = atoi(argv[1])) <= 0) {
    printf("This program computes the n'th decimal digit of pi\n"
           "usage: pi n , where n is the digit you want\n");
    exit(1);
    }

    N = (int) ((n + 20) * log(10) / log(13.5));
    sum = 0;

    for (a = 2; a <= (3 * N); a = next_prime(a)) {
    vmax = (int) (log(3 * N) / log(a));
    if (a == 2) {
        vmax = vmax + (N - n);
        if (vmax <= 0)
        continue;
    }
    av = 1;
    for (i = 0; i < vmax; i++)
        av = av * a;

    s = 0;
    den = 1;
    kq1 = 0;
    kq2 = -1;
    kq3 = -3;
    kq4 = -2;
    if (a == 2) {
        num = 1;
        v = -n;
    } else {
        num = pow_mod(2, n, av);
        v = 0;
    }

    for (k = 1; k <= N; k++) {

        t = 2 * k;
        DIVN(t, a, v, -1, kq1, 2);
        num = mul_mod(num, t, av);

        t = 2 * k - 1;
        DIVN(t, a, v, -1, kq2, 2);
        num = mul_mod(num, t, av);

        t = 3 * (3 * k - 1);
        DIVN(t, a, v, 1, kq3, 9);
        den = mul_mod(den, t, av);

        t = (3 * k - 2);
        DIVN(t, a, v, 1, kq4, 3);
        if (a != 2)
        t = t * 2;
        else
        v++;
        den = mul_mod(den, t, av);

        if (v > 0) {
        if (a != 2)
            t = inv_mod2(den, av);
        else
            t = inv_mod(den, av);
        t = mul_mod(t, num, av);
        for (i = v; i < vmax; i++)
            t = mul_mod(t, a, av);
        t1 = (25 * k - 3);
        t = mul_mod(t, t1, av);
        s += t;
        if (s >= av)
            s -= av;
        }
    }

    t = pow_mod(5, n - 1, av);
    s = mul_mod(s, t, av);
    sum = fmod(sum + (double) s / (double) av, 1.0);
    }
    printf("Decimal digits of pi at position %d: %09d\n", n,
       (int) (sum * 1e9));
    return 0;
}

Cross-posted at: ExcelForum.com

 

[rlv01]

Hard to tell without good reference calculation to compare results against. Having done this kind of thing a number of times, I can only tell you that if it were me, I would first spend my time insuring that the sub functions (inv_mod, pow_mod, etc) are tested and bulletproof. And I would probably not rely on using the variant data type as you are doing:

Function pow_mod(a, b, m)

for the function parameters. I've always found that it is best to explicitly type variables to help uncover any issues. I realize this is not too specific, but I hope that helps.

 

[Juggler_IN]

@rlv01; Sorry for the late reply. The support codes are tested. I am not pursuing this code and have done a new code posting: NthPi.

 

[Juggler_IN]

I need help in fixing a C to VBA conversion code.

The code outputs the specified n-th digit of pi in base-10 (Decimal System) on the lines of the Bailey–Borwein–Plouffe formula which outputs the specified n-th digit of pi in base-16 (Hexadecimal System).

The converted, but not functional, VBA code and the reference C code are enclosed.

A fully-functional execution of the C code can be viewed at tio.run which is an online C compiler.

Sub Main0()
' Expected output for i = 1 to 50 is 1 4 1 5 9 2 6 5 3 5 8 9 7 9 3 2 3 8 4 6 2 6 4 3 3 8 3 2 7 9 4 0 2 8 8 4 1 9 7 1 6 9 3 9 9 3 7 5 1 0 5.

    For i = 0 To 50
        Debug.Print NthPi(i) & " ";
    Next i

End Sub
Function NthPi(ByVal d As Integer) As Long

    Dim ll, j, i

    ll = Int((d + 2) * 10 / 3 + 2)
    j = 0
    i = 0

    Dim x, r

    ReDim x(0 To ll)
    ReDim r(0 To ll)

    Do
        x(j) = 20
        j = j + 1
    Loop While j < ll

    Dim c, n, e, p

    p = 0

    Do While i < d

        j = 0
        c = 0

        Do
            n = ll - j - 1
            e = n * 2 + 1
            r(j) = (x(j) + c) Mod e
            c = (x(j) / e) * n
            j = j + 1
        Loop While j < ll

        ll = ll - 1
        p = x(ll) / 10
        r(ll) = x(ll + 1) Mod 10

        j = 0
        
        Do
            x(j) = r(j) * 10
            j = j + 1
        Loop While j < ll

        i = i + 1

    Loop

    NthPi = p Mod 10

End Function

C code:

namespace System.Linq
{
    class P
    {
        static void Main()
        {
            for (int i = 0; i <= 50; ++i)
            {
                Console.Write(NthPi(i));
            }

            Console.WriteLine();
            Console.ReadLine();
        }

        static long NthPi(int d)
        {
            int l = (d+=2) * 10 / 3 + 2, j = 0, i = 0;

            long[] x = new long[l], r = new long[l];

            for (; j < l;)
                x[j++] = 20;

            long c, n, e, p = 0;

            for (; i < d; ++i)
            {
                for (j = 0, c = 0; j < l; c = (x[j++] / e) * n)
                {
                    n = l - j - 1;
                    e = n * 2 + 1;
                    
                    r[j] = (x[j] += c) % e;
                }
                
                p = x[--l] / 10;
                
                r[l] = x[l++] % 10;

                for (j = 0; j < l;)
                    x[j] = r[j++] * 10;
            }
            
            return p % 10;
        }
    }
}