# Nth digit of Pi

#### Juggler_IN

##### Active Member
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
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:
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

Last edited by a moderator:

### Excel Facts

Format cells as time
Select range and press Ctrl+Shift+2 to format cells as time. (Shift 2 is the @ sign).
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:

VBA Code:
``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.

@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.

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.
VBA Code:
``````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:
Code:
``````namespace System.Linq
{
class P
{
static void Main()
{
for (int i = 0; i <= 50; ++i)
{
Console.Write(NthPi(i));
}

Console.WriteLine();
}

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;
}
}
}``````

Replies
0
Views
523
Replies
21
Views
336
Replies
6
Views
116
Replies
1
Views
116
Replies
4
Views
101

1,196,510
Messages
6,015,619
Members
441,908
Latest member
Guitar1st

### We've detected that you are using an adblocker.

We have a great community of people providing Excel help here, but the hosting costs are enormous. You can help keep this site running by allowing ads on MrExcel.com.

### Which adblocker are you using?

1)Click on the icon in the browser’s toolbar.
2)Click on the icon in the browser’s toolbar.
2)Click on the "Pause on this site" option.
Go back

1)Click on the icon in the browser’s toolbar.
2)Click on the toggle to disable it for "mrexcel.com".
Go back

### Disable uBlock Origin

Follow these easy steps to disable uBlock Origin

1)Click on the icon in the browser’s toolbar.
2)Click on the "Power" button.
3)Click on the "Refresh" button.
Go back

### Disable uBlock

Follow these easy steps to disable uBlock

1)Click on the icon in the browser’s toolbar.
2)Click on the "Power" button.
3)Click on the "Refresh" button.
Go back