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

Option Explicit

Public Const eps As Double = 1E-17
Public Const NHX  As Double = 16

' Эта программа реализует алгоритм 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 As Double
    Dim s1 As Double
    Dim s2 As Double
    Dim s3 As Double
    Dim s4 As Double
    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
  
    Public Function ihex(ByVal x As Double, ByVal NHX As Integer, ByRef chx)
    
    'Возвращает, в СНХ, первые NHX шестигранные цифры фракции х.
    
    Dim i As Integer
    Dim y As Double
    Dim hx
        
        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

    
    
    Public Function series(ByVal m As Integer, ByVal 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

Option Explicit

Public Const eps As Double = 1E-17
Public Const NHX  As Double = 16

' Код программы помещают в модуль VBA Excel, печать результата производится на активный лист,
' подключение библиотек, форм и т.п. не требуется.

    Public Sub main()
    Dim i
    
    Dim pid As Double
    Dim s1 As Double
    Dim s2 As Double
    Dim s3 As Double
    Dim s4 As Double
    'Dim series#(m%, N%)
    Dim id 'as integer
    '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)
    i = i + 1
    ActiveSheet.Cells(i, 1) = id 'выводит результат на активный лист
    ActiveSheet.Cells(i, 2) = pid
    ActiveSheet.Cells(i, 3) = "chx"
    End Sub
  
    Public Function ihex(ByVal x As Double, ByVal NHX As Integer, ByRef chx)
    
    'Возвращает, в СНХ, первые NHX шестигранные цифры фракции х.
    
    Dim i As Integer
    Dim y As Double
    Dim hx
        
        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

    
    
    Public Function series(ByVal m As Integer, ByVal id) '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



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

/*  
    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
 
    Public Const eps As Variant = 1E-17
    Public Const NHX  As Integer = 16

' Код программы помещают в модуль VBA Excel, печать результата производится на активный лист,
' подключение библиотек, форм и т.п. не требуется.

    Public Sub main()
    
    Dim pid As Variant
    Dim s1 As Variant
    Dim s2 As Variant
    Dim s3 As Variant
    Dim s4 As Variant
    Dim id As Integer
    'void ihex(x%, m%, char c())
    
    For id = 1 To 1000 Step 1
    'id = 1000
    '#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 - Int(pid) + 1
    'ihex (pid, NHX, chx)
    'printf (" position = %i\n fraction = %.15f \n hex digits = %10.10s\n",id, pid, chx)

    ActiveSheet.Cells(Int(id), 1) = id 'выводит результат на активный лист
    ActiveSheet.Cells(Int(id), 2) = pid
   ' ActiveSheet.Cells(i, 3) = "chx"
    Next id
    
    End Sub
 
    Public Function ihex(ByVal x As Variant, ByVal NHX As Integer, ByVal chx)
    
    'Возвращает, в СНХ, первые NHX шестигранные цифры фракции х.
    
    Dim i As Integer
    Dim y As Variant
    Dim hx
              
        hx = Array("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "A", "B", "C", "D", "E", "F")
    
        y = Abs(x)
    
    For i = 1 To NHX - 1
        y = 16 * (y - Int(y))
        chx(i) = hx(Int(y))
    Next i
    
    End Function
 
    Public Function series(ByVal m As Integer, ByVal id As Integer) 'id as integer
    'Эта процедура оценивает серию sum_k 16 ^ (ID-K) / (8 * K + M)
    'используя модульную технику возведения в степень.
    Dim k As Integer
    Dim ak As Variant
   
    Dim p As Variant
    Dim s As Variant
    Dim t As Variant
    
    s = 0

    'Суммирует ряд до идентификатора
 
    For k = 1 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 Variant, ByVal ak As Variant)
    'expm = 16 ^ р мод ak. Эта процедура использует бинарную схему
    'возведения в степень слева-направо.

    Dim i As Variant
    Dim j As Variant
    Dim p1 As Variant
    Dim pt As Variant
    Dim r As Variant
    Dim ntp As Variant
    
    Static tp(25) As Variant
    Static tp1 As Variant
    
    ntp = 25 '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 = 1 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