/*  
    This program implements the BBP algorithm to generate a few hexadecimal
    digits beginning immediately after a given position id, or in other words
    beginning at position id + 1.  On most systems using IEEE 64-bit floating-
    point arithmetic, this code works correctly so long as d is less than
    approximately 1.18 x 10^7.  If 80-bit arithmetic can be employed, this limit
    is significantly higher.  Whatever arithmetic is used, results for a given
    position id can be checked by repeating with id-1 or id+1, and verifying 
    that the hex digits perfectly overlap with an offset of one, except possibly
    for a few trailing digits.  The resulting fractions are typically accurate 
    to at least 11 decimal digits, and to at least 9 hex digits.  
*/

/*  David H. Bailey     2006-09-08 */

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

int main()
{
  double pid, s1, s2, s3, s4;
  double series (int m, int n);
  void ihex (double x, int m, char c[]);
  int id = 1000000;
#define NHX 16
  char chx[NHX];

/*  id is the digit position.  Digits generated follow immediately after id. */

  s1 = series (1, id);
  s2 = series (4, id);
  s3 = series (5, id);
  s4 = series (6, id);
  pid = 4. * s1 - 2. * s2 - s3 - s4;
  pid = pid - (int) pid + 1.;
  ihex (pid, NHX, chx);
  printf (" position = %i\n fraction = %.15f \n hex digits =  %10.10s\n",
  id, pid, chx);
}

void ihex (double x, int nhx, char chx[])

/*  This returns, in chx, the first nhx hex digits of the fraction of x. */

{
  int i;
  double y;
  char hx[] = "0123456789ABCDEF";

  y = fabs (x);

  for (i = 0; i < nhx; i++){
    y = 16. * (y - floor (y));
    chx[i] = hx[(int) y];
  }
}

double series (int m, int id)

/*  This routine evaluates the series  sum_k 16^(id-k)/(8*k+m) 
    using the modular exponentiation technique. */

{
  int k;
  double ak, eps, p, s, t;
  double expm (double x, double y);
#define eps 1e-17

  s = 0.;

/*  Sum the series up to id. */

  for (k = 0; k < id; k++){
    ak = 8 * k + m;
    p = id - k;
    t = expm (p, ak);
    s = s + t / ak;
    s = s - (int) s;
  }

/*  Compute a few terms where k >= id. */

  for (k = id; k <= id + 100; k++){
    ak = 8 * k + m;
    t = pow (16., (double) (id - k)) / ak;
    if (t < eps) break;
    s = s + t;
    s = s - (int) s;
  }
  return s;
}

double expm (double p, double ak)

/*  expm = 16^p mod ak.  This routine uses the left-to-right binary 
    exponentiation scheme. */

{
  int i, j;
  double p1, pt, r;
#define ntp 25
  static double tp[ntp];
  static int tp1 = 0;

/*  If this is the first call to expm, fill the power of two table tp. */

  if (tp1 == 0) {
    tp1 = 1;
    tp[0] = 1.;

    for (i = 1; i < ntp; i++) tp[i] = 2. * tp[i-1];
  }

  if (ak == 1.) return 0.;

/*  Find the greatest power of two less than or equal to p. */

  for (i = 0; i < ntp; i++) if (tp[i] > p) break;

  pt = tp[i-1];
  p1 = p;
  r = 1.;

/*  Perform binary exponentiation algorithm modulo ak. */

  for (j = 1; j <= i; j++){
    if (p1 >= pt){
      r = 16. * r;
      r = r - (int) (r / ak) * ak;
      p1 = p1 - pt;
    }
    pt = 0.5 * pt;
    if (pt >= 1.){
      r = r * r;
      r = r - (int) (r / ak) * ak;
    }
  }

  return r;
}

Option Explicit

' Эта программа реализует алгоритм BBP генерировать несколько
' шестнадцатеричных цифр, начиная сразу после заданной позиции ID,
' или, другими словами, начиная с позиции ID + 1. На большинстве систем,
' использующих 64-разрядную арифметику IEEE с плавающей запятой,
' этот код работает правильно, пока г меньше примерно 1,18 х 10 ^ 7.
' Если 80-битный арифметическое могут быть использованы,
' этот предел значительно выше. Независимо арифметика используется,
' результаты для данной позиции ID можно проверить путем повторения с ID-1 или ID 1,
' и проверки, что шестнадцатеричные цифры прекрасно перекрывают со смещением одного,
' за исключением, возможно в течение нескольких хвостовых цифр. Полученные фракции,
' как правило, с точностью не менее 11 десятичных цифр, и, по крайней мере 9 шестнадцатеричных цифр.


Private Function main() As Integer

Dim pid$, s1$, s2$, s3$, s4$
Dim series$(m%, n%)
Dim id%
'void ihex(x%, m%, char c())
id% = 1000000
'$define NHX 16
' char chx[NHX];
'идентификатор цифра положение. Цифры генерируется следовать сразу после идентификатором.

 s1$ = series$(1, id%)
 s2$ = series$(4, id%)
 s3$ = series$(5, id%)
 s4$ = series$(6, id%)
' pid$ = 4$ * s1$ - 2$ * s2$ - s3$ - s4$
' pid$ = pid$ - pid$ + 1$
 'ihex (pid, NHX, chx);
 'printf (" position = %i\n fraction = %.15f \n hex digits = %10.10s\n",
 'id, pid, chx);

End Function

Function ihex(x$, nhx%, ByRef chx As String)
Dim i%
Dim y$
Dim hx$
hx = "0123456789ABCDEF"
y = Abs(x)
For i% = 0 To nhx - 1
'y = 16$ * (y - Floor(y))
End Function

Function series$(m%, id%)
'Эта процедура оценивает сериала sum_k 16 ^ (ID-K) / (8 * K + M)
'используя модульную технику возведения в степень.
'{
 Dim k%
 Dim ak$
 Dim eps$
 Dim p$
 Dim s$
 Dim t$
 Dim expm$(x$, y$)
 Dim eps$
 
 eps$ = 1E-17
 s$ = 0

'Просуммировать ряд до идентификатором
 
 For k% = 0 To id% Step 1
 ak$ = 8 * k% + m
 p$ = id% - k%
 t$ = expm(p$, ak$)
 s$ = s$ + t$ / ak$
 s$ = s$ - Int(s$)
 Next k%

'Вычислить несколько терминов, где к> = ID.

 For k = id To id + 100 Step 1
    ak$ = 8 * k% + m
    't = pow (16., (double) (id - k)) / ak;
    If t < eps Then Exit For
        s = s + t
        s = s - Int(s)
 Next k
 
 series$() = s

End Function

Public Function expm$(ByVal p$, ByVal ak$)

    'expm = 16 ^ р мод ак. Эта процедура использует бинарный схему
    'возведения в степень слева-направо.

    Dim i%
    Dim j%
    Dim p1$
    Dim pt$
    Dim r$
    Dim ntp%
    
    Dim tp$(25)
    Dim tp1%
    
    tp1% = 0
    ntp% = 25
    
'Если это первый вызов expm, заполнить мощь двух таблицы тп.
    If tp1% = 0 Then
        tp1% = 1
        tp$(0) = 1

            For i% = 1 To ntp% Step 1
            tp$(i%) = 2 * tp$(i% - 1)
            Next i%
    End If

 If ak$ = "1" Then expm$ = "0"

'Найти наибольшую степень двойки меньше или равен р.
    For i% = 0 To ntp Step 1
        If tp(i%) > p Then Exit For
            pt$ = tp$(i% - 1)
            p1$ = p$
            r$ = 1
    Next i%

'Выполните бинарный алгоритм возведения в степень по модулю ак.
     For j% = 1 To i% Step 1
        If p1$ >= pt$ Then
          r$ = 16 * r$
           r$ = r$ - Int((r$ / ak$) * ak$)
            p1$ = p1$ - pt$
        End If
 
    pt$ = 0.5 * pt$
  
        If pt$ >= 1 Then
            r$ = r$ * r$
            r$ = r$ - Int((r$ / ak$) * ak$)
        End If
    Next j%

    expm$ = r$

End Function

Private Sub CommandButton1_Click()
         
      Dim n As Long, v As Double
      For n = 1 To 10001 Step 4 'For n = 1 To 100000001 Step 4
        v = v + 4 / n - 4 / (n + 2)
      Next n
      Cells(1, 1) = v
End Sub

Option Explicit

    Public Function expm(ByVal p As Double, ByVal ak As Double)
    'expm = 16 ^ р мод ak. Эта процедура использует бинарную схему
    'возведения в степень слева-направо.

    Dim i As Integer
    Dim j As Integer
    Dim p1 As Double
    Dim pt As Double
    Dim r As Double
    Dim ntp As Integer
    
    Static tp(25) As Double
    Static tp1 As Integer
    
    tp1% = 0
    ntp% = 25
    
    'Если это первый вызов expm, то заполнить массив tp.
    If tp1 = 0 Then
        tp1 = 1
        tp(0) = 1

            For i = 1 To ntp Step 1
            tp(i) = 2 * tp(i - 1)
            Next i
    End If

    If ak = 1 Then expm = 0

    'Найти наибольшую степень двойки меньше или равной р.
    For i = 0 To ntp Step 1
        If tp(i) > p Then Exit For
            pt = tp(i - 1)
            p1 = p
            r = 1
    Next i%

    'Выполните бинарный алгоритм возведения в степень по модулю ак.
    For j% = 1 To i% Step 1
        If p1 >= pt Then
            r = 16 * r
            r = r - Int((r / ak) * ak)
            p1 = p1 - pt
        End If
 
    pt = 0.5 * pt
  
        If pt >= 1 Then
            r = r * r
            r = r - Int((r / ak) * ak)
        End If
    Next j

    expm = r

End Function

double expm (double p, double ak)

/*  expm = 16^p mod ak.  This routine uses the left-to-right binary 
    exponentiation scheme. */

{
  int i, j;
  double p1, pt, r;
#define ntp 25
  static double tp[ntp];
  static int tp1 = 0;

/*  If this is the first call to expm, fill the power of two table tp. */

  if (tp1 == 0) {
    tp1 = 1;
    tp[0] = 1.;

    for (i = 1; i < ntp; i++) tp[i] = 2. * tp[i-1];
  }

  if (ak == 1.) return 0.;

/*  Find the greatest power of two less than or equal to p. */

  for (i = 0; i < ntp; i++) if (tp[i] > p) break;

  pt = tp[i-1];
  p1 = p;
  r = 1.;

/*  Perform binary exponentiation algorithm modulo ak. */

  for (j = 1; j <= i; j++){
    if (p1 >= pt){
      r = 16. * r;
      r = r - (int) (r / ak) * ak;
      p1 = p1 - pt;
    }
    pt = 0.5 * pt;
    if (pt >= 1.){
      r = r * r;
      r = r - (int) (r / ak) * ak;
    }
  }

  return r;
}

double expm (double p, double ak)

/*  expm = 16^p mod ak.  This routine uses the left-to-right binary 
    exponentiation scheme. */

{
  int i, j;
  double p1, pt, r;
#define ntp 25
  static double tp[ntp];
  static int tp1 = 0;

/*  If this is the first call to expm, fill the power of two table tp. */

  if (tp1 == 0) {
    tp1 = 1;
    tp[0] = 1.;

    for (i = 1; i < ntp; i++) tp[i] = 2. * tp[i-1];
  }

  if (ak == 1.) return 0.;

/*  Find the greatest power of two less than or equal to p. */

  for (i = 0; i < ntp; i++) if (tp[i] > p) break;

  pt = tp[i-1];
  p1 = p;
  r = 1.;

/*  Perform binary exponentiation algorithm modulo ak. */

  for (j = 1; j <= i; j++){
    if (p1 >= pt){
      r = 16. * r;
      r = r - (int) (r / ak) * ak;
      p1 = p1 - pt;
    }
    pt = 0.5 * pt;
    if (pt >= 1.){
      r = r * r;
      r = r - (int) (r / ak) * ak;
    }
  }

  return r;
}

Option Explicit

    Public Function expm(ByVal p As Double, ByVal ak As Double)
    'expm = 16 ^ р мод ak. Эта процедура использует бинарную схему
    'возведения в степень слева-направо.

    Dim i As Integer
    Dim j As Integer
    Dim p1 As Double
    Dim pt As Double
    Dim r As Double
    Dim ntp As Integer
    
    Static tp(25) As Double
    Static tp1 As Integer
    
    ntp = 25
    
    'Если это первый вызов expm, то заполнить массив tp.
    If tp1 <> 1 Then
        tp1 = 1
        tp(0) = 1

            For i = 1 To ntp - 1 Step 1
            tp(i) = 2 * tp(i - 1)
            Next i
    End If

    If ak = 1 Then expm = 0

    'Найти наибольшую степень двойки меньше или равной р.
        For i = 0 To ntp - 1 Step 1
            If tp(i) > p Then Exit For
        Next i
            
            pt = tp(i - 1)
            p1 = p
            r = 1
    
    'Выполняет бинарный алгоритм возведения в степень по модулю ак.
    For j = 1 To i - 1 Step 1
        If p1 >= pt Then
            r = 16 * r
            r = r - Int((r / ak) * ak)
            p1 = p1 - pt
        End If
 
    pt = 0.5 * pt
  
        If pt >= 1 Then
            r = r * r
            r = r - Int((r / ak) * ak)
        End If
    Next j

    expm = r

End Function

double series (int m, int id)

/*  This routine evaluates the series  sum_k 16^(id-k)/(8*k+m) 
    using the modular exponentiation technique. */

{
  int k;
  double ak, eps, p, s, t;
  double expm (double x, double y);
#define eps 1e-17

  s = 0.;

/*  Sum the series up to id. */

  for (k = 0; k < id; k++){
    ak = 8 * k + m;
    p = id - k;
    t = expm (p, ak);
    s = s + t / ak;
    s = s - (int) s;
  }

/*  Compute a few terms where k >= id. */

  for (k = id; k <= id + 100; k++){
    ak = 8 * k + m;
    t = pow (16., (double) (id - k)) / ak;
    if (t < eps) break;
    s = s + t;
    s = s - (int) s;
  }
  return s;
}

 Function series(m As Integer, id As Integer)
    'Эта процедура оценивает серию sum_k 16 ^ (ID-K) / (8 * K + M)
    'используя модульную технику возведения в степень.
    Dim k As Integer
    Dim ak As Double
    Dim eps As Double
    Dim p As Double
    Dim s As Double
    Dim t As Double
    Dim expm(x, y) As Double
    Dim eps
 
    eps = 1E-17
    s = 0

    'Суммирует ряд до идентификатора
 
    For k = 0 To id Step 1
        ak = 8 * k + m
        p = id - k
        t = expm(p, ak)
        s = s + t / ak
        s = s - Int(s)
    Next k%

    'Вычисляет несколько терминов, где к> = ID.

    For k = id To id + 100 Step 1
            ak = 8 * k + m
            
           't = pow (16., (double) (id - k)) / ak;
            t = 16 * (id - k) / ak
            
        If t < eps Then Exit For
        
            s = s + t
            s = s - Int(s)
    Next k
 
    series = s

    End Function

    Function series(m As Integer, id As Integer)
    'Эта процедура оценивает серию sum_k 16 ^ (ID-K) / (8 * K + M)
    'используя модульную технику возведения в степень.
    Dim k As Integer
    Dim ak As Double
    Dim eps As Double
    Dim p As Double
    Dim s As Double
    Dim t As Double
    Dim expm(x, y) As Double
    Const eps As Double = 1E-17

    s = 0

    'Суммирует ряд до идентификатора
 
    For k = 0 To id - 1 Step 1
        ak = 8 * k + m
        p = id - k
        t = expm(p, ak)
        s = s + t / ak
        s = s - Int(s)
    Next k%

    'Вычисляет несколько терминов, где к> = ID.

    For k = id To id - 1 + 100 Step 1
            ak = 8 * k + m
            
           't = pow (16., (double) (id - k)) / ak;
            t = 16 ^ ((id - k)) / ak
            
        If t < eps Then Exit For
        
            s = s + t
            s = s - Int(s)
    Next k
 
    series = s

    End Function

Option Explicit

Public Const eps As Double = 1E-17

    Function series(m As Integer, id As Integer)
    'Эта процедура оценивает серию sum_k 16 ^ (ID-K) / (8 * K + M)
    'используя модульную технику возведения в степень.
    Dim k As Integer
    Dim ak As Double
    Dim eps As Double
    Dim p As Double
    Dim s As Double
    Dim t As Double
    
    s = 0

    'Суммирует ряд до идентификатора
 
    For k = 0 To id - 1 Step 1
        ak = 8 * k + m
        p = id - k
        t = expm(p, ak)
        s = s + t / ak
        s = s - Int(s)
    Next k%

    'Вычисляет несколько терминов, где к> = ID.

    For k = id To id - 1 + 100 Step 1
            ak = 8 * k + m
            
           't = pow (16., (double) (id - k)) / ak;
            t = 16 ^ ((id - k)) / ak
            
        If t < eps Then Exit For
        
            s = s + t
            s = s - Int(s)
    Next k
 
    series = s

    End Function


    Public Function expm(ByVal p As Double, ByVal ak As Double)
    'expm = 16 ^ р мод ak. Эта процедура использует бинарную схему
    'возведения в степень слева-направо.

    Dim i As Integer
    Dim j As Integer
    Dim p1 As Double
    Dim pt As Double
    Dim r As Double
    Dim ntp As Integer
    
    Static tp(25) As Double
    Static tp1 As Integer
    
    ntp = 25
    
    'Если это первый вызов expm, то заполнить массив tp.
    If tp1 <> 1 Then
        tp1 = 1
        tp(0) = 1

            For i = 1 To ntp - 1 Step 1
            tp(i) = 2 * tp(i - 1)
            Next i
    End If

    If ak = 1 Then expm = 0

    'Найти наибольшую степень двойки меньше или равной р.
        For i = 0 To ntp - 1 Step 1
            If tp(i) > p Then Exit For
        Next i
            
            pt = tp(i - 1)
            p1 = p
            r = 1
    
    'Выполняет бинарный алгоритм возведения в степень по модулю ак.
    For j = 1 To i - 1 Step 1
        If p1 >= pt Then
            r = 16 * r
            r = r - Int((r / ak) * ak)
            p1 = p1 - pt
        End If
 
    pt = 0.5 * pt
  
        If pt >= 1 Then
            r = r * r
            r = r - Int((r / ak) * ak)
        End If
    Next j

    expm = r

End Function

void ihex (double x, int nhx, char chx[])

/*  This returns, in chx, the first nhx hex digits of the fraction of x. */

{
  int i;
  double y;
  char hx[] = "0123456789ABCDEF";

  y = fabs (x);

  for (i = 0; i < nhx; i++){
    y = 16. * (y - floor (y));
    chx[i] = hx[(int) y];
  }
}

  Function ihex(x As Double, nhx As Integer, ByRef chx As String)
    Dim i As Integer
    Dim y As Double
    Dim hx As Double
        hx = "0123456789ABCDEF"
    
        y = Abs(x)
    
    For i = 0 To nhx - 1
        'y = 16# * (y - Floor(y))
        'chx(i) = hx(int( y))
    Next nhx
    
    End Function

void ihex (double x, int nhx, char chx[])

/*  This returns, in chx, the first nhx hex digits of the fraction of x. */

{
  int i;
  double y;
  char hx[] = "0123456789ABCDEF";

  y = fabs (x);

  for (i = 0; i < nhx; i++){
    y = 16. * (y - floor (y));
    chx[i] = hx[(int) y];
  }
}

   Function ihex(x As Double, nhx As Integer, ByRef chx As String)
    
    'Возвращает, в СНХ, первые NHX шестигранные цифры фракции х.
    
    Dim i As Integer
    Dim y As Double
    Dim hx As Double
        
        hx = Array("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "A", "B", "C", "D", "E", "F")
    
        y = Abs(x)
    
    For i = 0 To nhx - 1
        y = 16 * (y - Int(y))
        chx(i) = hx(Int(y))
    Next nhx
    
    End Function