Lazarus
Free Pascal => General => Topic started by: HappyLarry on January 18, 2015, 05:42:31 pm
-
I have been using the gmp unit and have a reasonable understanding of it.
I have used the floating point 'sqrt' function and it has worked correctly up until now. However when I try to find root of
1247257156019619631959228850474915324796444794921068388944443347008749568000000000000000
I get the answer
35316528085580830944581390125742259962343219.93...
The correct answer is
35316528085580830944581390125749276564383499.14....
I have set the precision as 4*(number of digits in the original number). I have noticed that the error occurs on the 32nd digit and 32*4 = 128 which may be important, but I don't know why nor how to get around it.
Any ideas would be welcome.
-
Here is the procedure which might help).
procedure sqrt2(var strAnswer:string; var Expon:integer; strNum1: string; NumPlaces:integer);
var
a,b: MPFloat;
n:integer;
p:mp_exp_t;
begin
f_init (a);
f_init (b);
f_set_str (a, '1', 10);
f_set_str (b, strNum1, 10);
p:=4*NumPlaces;
f_set_prec (a, p);
f_sqrt (a,b);
strAnswer:=f_get_str(n, 10, p, a);
Expon:=n;
f_clear (a);
f_clear (b);
end;
-
Precision maybe binary based (number of bits) instead of decimal ?
-
It is binary-based but I have assumed that the number of bits needed to represent a number in decimal will be 4*(number of decimal digits).
e.g
any 1 digit number can be represented in 4 bits since 24 = 16
any 2 digit number can be represented in 8 bits since 28 = 256
any 3 digit number can be represented in 12 bits since 212 = 4096
any 4 digit number can be represented in 16 bits since 216 = 65536
etc.
That explains why I have set the precision as 4*(number of digits).
-
It is binary-based but I have assumed that the number of bits needed to represent a number in decimal will be 4*(number of decimal digits).
Ok. Then I don't know :-)
(it is 3.322 btw, log 10 / log 2, which comes down to slightly over 106 bits for 32 digits.
-
I assumed that too little precision was a problem rather than too much!
I have just tried the procedure with
p:=trunc(3.323 * NumPlaces)+1;
but . . . the answer is the same (wrong).
-
The correct answer is
35316528085580830944581390125749276564383499.14....
35316528085580830944581390125749276564383499.172 is closer to the right answer than the "correct" answer you gave?
35316528085580830944581390125749276564383499.172 ^ 2 = 1247257156019619631959228850474915324796444794943759140060606358369027292979933081630924.685584
35316528085580830944581390125749276564383499.14 ^ 2 = 1247257156019619631959228850474915324796444792683501342583433177915818324931979381510380.7396
Bart
-
@Bart
You are right. I had checked for the correct answer using bc but copied it out wrongly (bc is easy to use but I can't copy-and-paste off the screen).
The correct answer is (as you say)
35316528085580830944581390125749276564383499.1716787516731376083479
Either way, the answer I get using gmp is not even close. I wondered what I was doing wrong.
-
The correct answer is (as you say)
35316528085580830944581390125749276564383499.1716787516731376083479
Needs more digits:
35316528085580830944581390125749276564383499.171678751673137608347918603494756890875859839579687481657183436994157146340168806365847166873813522
-
I wondered what I was doing wrong.
Two small typos:
1-Perc is number of "bits".
2-In your call to f_get_str you do not pass the number of "digits"
Change your code to:
procedure sqrt2(var strAnswer:string; var Expon:integer; strNum1: string; NumPlaces:integer);
var
a,b: MPFloat;
n:integer;
p:mp_exp_t;
prec: ValUInt;
digits: ValUInt;
begin
digits := 200;
prec := 8*digits;
f_set_default_prec(prec);
f_init (a);
f_init (b);
f_set_str (a, '1', 10);
f_set_str (b, strNum1, 10);
//p:=4*NumPlaces;
//f_set_prec (a, p);
f_sqrt (a,b);
strAnswer:=f_get_str(n, 10, digits, a);
Expon:=n;
f_clear (a);
f_clear (b);
end;
You should get 200 digits:
35316528085580830944581390125749276564383499.17167875167313760834791860349475689087585983957968748165718343699415714634016880636584716687381352199364066613445449266509253059227517946683412221654516511
-
This works. Thank you.
Conclusion:
Precision is measured in bits and each digit required needs 8 bits of precision.
A number is stored digit-by-digit like the characters in a string rather than as a binary number like an integer.
Effectively I had half the precision needed.
Thanks to all contributors.