last day (17 days later) » 

2:00 PM
-4
Q: BBP - formula for PI calculation

Sebi2020The Problem I want to implement a program that calculates the n-th number of pi. But something is wrong with my implementation. I can't find the error. Any Ideas? My Thoughts I think I have a wrong understanding of the BBP-algorithm. Here is my Source: BBP formula digit extraction algorithm ...

 
the debugger can't say if it's a mathematical problem... so?
 
Better formatting attracts more attention. Also, I am pretty sure you can simplify your egcd function. Finally, do not use tags that tell you not to use them.
 
Which part of my posting needs to be better formatted? Yes the ecgd function could be simpler. But for now this function works.
 
@Sebi2020 Basic indentation, consistant spacing (and probably a bit more spacing)... Oh, I forgot : consistant lines, too : sometimes you have one-line blocks, sometimes not. And "this function works" does not mean it is readable.
 
You can't understand a mathematical function if you don't know the algorithm. I've added indents. Thank you for downvoting!!
A can fix the problem by myself, but I can't see it. Because of this I'm posting this question, more eyes see more. If it would be so easy i won't post such a thing.
 
2:00 PM
@Sebi2020 I did not downvote, I was (and am still) waiting for you to improve your formatting. You have done something, you are on the right track ; but you can still do way more. Just an exemple to show why formatting is not to be underrated : before your edit, I had trouble locating the main function.
 
Sorry, but I'm looking for help, if everyone has to get put in front of the nose things, so you can guess it yourself, you probably need no help.
I have more than once explained that I do not know where the problem lies and because of this I'm asking for help. I think I had a wrong understanding of the algorithm, so this implementation don't work, the functions in the math.c part work, but the BBP function and calcSingleSum() is the critical part.
 
@Sebi2020 Then why do you not specify it in your question instead of writing a comment ? I was looking at the egcd function because it was so unreadable (and it often causes problems). If you had guaranteed at the very beginning that it works and that the problem lies elsewhere, and if you had done that for the entire code, we could have solved your problem by now.
 
I have added some details, but my question is already put off-topic by @Floris -.-
 
@Sebi2020 That is not surprising. By the way, did you even try to compile this code ? Because right now, what you provide is not compilable. So, once again, improve your code, and only then ask for help.
 
Yes it's compilable! Have done this for several times! If you want to compile this you need the header definitions of math.c for main.c One error is a forgotten ")" in main. fixed that.
 
2:00 PM
@Sebi2020 Now turn on your compiler warnings (at least -Wall) and correct them. You might get surprises.
 
@Nelxiost Yes two warnings: erg is not used in main.c and I haven't returned a value in main.c but that's not the problem.
 
I have significantly more warnings, and the return value of main is not one of them.
 
2:11 PM
btw: Nice page ^^
 
You need to share before linking (bottom-right)
 
I want to write some unit testing for this program to test all the mathemtical functions. But I think really it's not a code problem. Really a mathematical problem.
 
Well yes, this is better. What means PiN ? Do you really need a function for doing "*16" ?
Unit testing is a good idea, and you do not have many functions to test, so it can be easy.
Is BBP meant to work well with 7 as second argument (samples) ?
BBP(0,13) == BBP(0,1000)
 
2:33 PM
Thats only for the sum to infty part. for the first numbers of pi, a few samples suffice. But not if you want to get 100000th's digit. with BBP(10000,10000) != BBP(10000,10)
PiN is useless but no problem, have you looked at the wiki's page for BBP-formula? Have seen something that a do wrong in BBP or calcSingleSum?
*Have you seen something that I do wrong
 
Maybe you should use a standard function instead of a1 -= (long) a1;. Not sure this works as expected.
 
It works, I've tried that before and I'd found it on stack overflow.
If you cast some time to integer types the mantissa is cut off.
*type
 
2:51 PM
Ok then. Are you sure modExp works too ?
 
Yes, I've tried that with big values and compared it with the results of the windows calculator:
W-Calc: 200^22 Mod 123 = 4,194304e+50 Mod 123 = 25
modExp: 200^22 Mod 123 = 25
 
Do you have a way to test calcSingleSum?
 
Okay found one error for the case e==0 and m==1, should return 0. Have fixed this, but nothing changes. for 0 I should get Pi with 3.14... I try this in geogebra to test calcSingleSum
 
Link the last code
 
3:13 PM
Okay, tested calcSingleSum, works as expected. One moment for the new code
Because the calculator and my program get the same results something must be wrong with the formula
 
I get 0 for the first sum.
It looks like an integer division error
I get different results if I cast the division
erg += (double)modExp(16,n-k,8*k+b)/(8*k+b);
But it still does not give pi.
 
But geogebra also gets 0 for the first sum
 
3:29 PM
For the 4 sums ?
It is an error anyway, because it will always yield 0
Whatever the integers, (n % m)/m is 0
 
I'm wondering they there is a need for a cast, thought it must return double because of the division
 
int/int = int
However, (double)int/int = double
 
I'm wondering because my calculator return 0 for every of the foiur sums
 
Well it should not
Take the second sum with n = 0 and b = 4
 
second with samples ?
 
3:37 PM
modExp(16,n-k,8*k+b) is modExp(16, 0, 4)
Nah, the second 'a'
modExp(16, 0, 4) is 16^0 % 4 = 1%4 = 1
(and the result is divided by b = 4, so it yields 0.25)
 
I found the problem with my calculator if you use sum for 0 with upper limit 0 no summerization should be executed
so the result is 0
But it's also an bug in the calculator. the Sum function doesn't work properly
But have you looked yet in the Wikipedia article about the Pi digit extraction?
 
Yes, but I do not see where the problem is
 
Do you think my function does the right steps?
 
Although I would consider not using the sum separation to simplify
When you write code, you need to make sure everything is fine step by step
And the "mod 8k + 1" trick is another step
For exemple, I would try that :
double calcSingleSum(int n, int b, int samples) {
    int k;
double sum = 0;
for(k = 0; k <= samples; ++k) {
	sum += pow(16, n-k) / (8*k+b);
}
printf("Sum: %4.10f\n", sum);
return sum;
}
Wth I cannot format it properly
Anyway you get the idea
Also, I just found out something
 
3:58 PM
But the separation is needed (say the wikipedia article) because this are different operations. one modExp operation and the other is pow
 
It is not
modExp is the same as pow but with a modulo to get rid of the integer part
pow(16, n-k) is meant to be 16^(0-k) with k between 0 and someBigNumber, right ?
 
yes but the need the modulo operation and your version has no part like this. please read the article in wikipedia:
Since we only care about the fractional part of the sum, we look at our two terms and realise that only the first sum is able to produce whole numbers; conversely, the second sum cannot produce whole numbers since the numerator can never be larger than the denominator for k > n. Therefore, we need a trick to remove the whole numbers for the first sum. That trick is mod 8k + 1. Our sum for the first fractional part then becomes:
**Notice how the modulo operator always guarantees **
 
> Since we only care about the fractional part of the sum
Just cut it out after computing the sum
 
But wikipedia says the need both sums
But however, some other problem. The CPU's mod operation doesn't support doubles. So if we do 0.25 % 1123 you get an error, thats one problem.
So the modExp has some issues
 
Huh ?
You do not have any modulo on double
 
4:05 PM
But every calculator says 0.25 % x = 0.25
have to go away for a hour
there is also an modulo exponentation for negativ exponents e.g. 2^-2 % 4 = 0.25
 
2^(-2) = 0.25 anyway...
 
4:21 PM
By the way this algorithm will give you the n'th digit of pi written in hexadecimal. I doubt that you will get 3.14159 with that.
 
4:36 PM
yes I know. But its also not right in hexadezimal, the first digits are 3,24
yes 2^(-2) = 0.25 ... 0.25 % x = 0.25 (says the calculator), but the modular exponentation says something other: Modular exponentiation can be performed with a negative exponent e by finding the modular multiplicative inverse d of b modulo m using the extended Euclidean algorithm
 
@Sebi2020 Then do not use %f in printf, use %a
 
have f only for testing purpourses
 
But how do you know what you should get with %f ?
 
because 3_10 = 3_16
it's only an test, because %x or %a can't deal with doubles
 
4:52 PM
%a deals with doubles.
 
okay I've changed that, but it's not pi yet
omg complicated...
 
5:10 PM
@Sebi2020 Note that it will not work for some other pi digit
 

  last day (17 days later) »