Reputation: 183
Good night, I'm having a issue with writeln on Embarcadero 10.2, I'm trying to write the sin() function as a Taylor expansion infinite serie.
Everything it's running fine, but the output is being in Scientific Expression like: 3.60448486921676E-0158 when the correct was 0.912945250727627654376099983845.
I need a 30 digits precision. My code below, It's a console program.
program Project1;
{$APPTYPE CONSOLE}
uses
System.SysUtils,
Windows,
Math;
function fact(n: LongInt): extended;
begin
if (n = 0) then
fact := 1
else
fact := n * fact(n - 1);
end;
var
one:integer;
fin:extended;
number:LongInt;
i:integer;
cons:LongInt;
cons1:extended;
begin
one := -1;
writeln('Digite o angulo em radianos para o Seno: ');
readln(number);
for i := 1 to 120 do
begin
cons := (2*i)+1;
if(i mod 2) = 0 then
one := 1
else
one := -1;
cons1 := fact(cons);
fin := (one/cons1)*(power(number,cons));
cons := 0;
end;
writeln(fin);
readln;
end.
Upvotes: 3
Views: 232
Reputation: 310
program Project1;
{$APPTYPE CONSOLE}
uses
System.SysUtils;
function sine(n: real): real;
asm
fld n
fsin
fstp Result
end;
var
n:real;
begin
write('Write the number: ');
readln(n);
writeln(sine(n):10:11);
readln;
end.
Upvotes: -1
Reputation: 80197
To perform calculation with very long precision, you need some library for arbitrary precision arithmetics - use Delphi wrapper/interface to GMP or some Delphi library. Example. Another one.
Note that for libraries with only long integers support (not floats) you still can calculate series like this Python code (it has internal support for big integers. **
is power, //
is div
).
The only division (slow operation) is performed once.
d30 = 10**30
d60 = d30 * d30
x = 3141592653589793238462643383279 // 4 #Pi/4 * d30
xx = - x * x
denom = d30
nom = x
sum = nom
i = 2
while abs(nom) * d30 > denom:
nom = nom * xx
mul = i * (i + 1) * d60
denom = denom * mul
sum = sum * mul + nom
i += 2
print('0.' + str(d30 * sum // denom))
>> 0.707106780551955811469384356540
Approach is based on this:
x - x^3/6 + x^5/120 = (x*120 - x^3 * 20 + x^5) / 120
Upvotes: 3
Reputation: 1141
Computing sin(x)
with Taylor expansion is stable only for small arguments, say less than Pi/2. For larger arguments you suffer from catastrophic cancelation, because intermediate terms are very large, e.g. for x=100
the largest terms is the 49'th with a value
0.107151028812546692318354675952E43
. All these large terms are added together and give a sin() value with magnitude <= 1
. You would need about 70 decimal digits, do compute the result with about 30 digits.
The solutions is: Use range reduction and multiprecision floating-point routines
(see mpf_sin
in my opensource Pascal MPArith package). Using the interactive calculator you get
sin(100) = -0.5063656411097587936565576104597854321
in 0.2 ms.
Upvotes: 3