Add BPP implementation from David Bailey.
This commit is contained in:
parent
c7539cd353
commit
5a728d7718
|
@ -2,6 +2,6 @@ INCLUDES = $(FUSE_CFLAGS) $(EXTRA_CFLAGS)
|
|||
|
||||
bin_PROGRAMS=pifs
|
||||
|
||||
pifs_SOURCES=pifs.c
|
||||
pifs_SOURCES=pifs.c piqpr8.c
|
||||
|
||||
pifs_LDADD = $(FUSE_LIBS)
|
||||
pifs_LDADD = $(FUSE_LIBS) -lm
|
||||
|
|
|
@ -0,0 +1,119 @@
|
|||
/*
|
||||
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>
|
||||
|
||||
static 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;
|
||||
}
|
||||
|
||||
static 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, p, s, t;
|
||||
#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;
|
||||
}
|
||||
|
||||
unsigned char get_byte(int id)
|
||||
{
|
||||
double s1 = series (1, id);
|
||||
double s2 = series (4, id);
|
||||
double s3 = series (5, id);
|
||||
double s4 = series (6, id);
|
||||
double pid = 4. * s1 - 2. * s2 - s3 - s4;
|
||||
pid = pid - (int) pid + 1.;
|
||||
|
||||
double y = fabs(pid);
|
||||
y = 16. * (y - floor (y));
|
||||
unsigned char first = y;
|
||||
y = 16. * (y - floor (y));
|
||||
unsigned char second = y;
|
||||
return (first << 4) | second;
|
||||
}
|
||||
|
||||
|
Loading…
Reference in New Issue