09-13-2022, 06:59 PM
bplus, the log and exp are done in double precision as part of the first approximation, I don't see a problem with that, also what if you want to take the root of a huge number?
for example
31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989 ^ (1/3)
I have split the first approximation to allow for larger exponents than double allows
instead of t2 = Log(t) / p + Log(10) * ex / p, do
t2 = Log(t) / p
t2 = Exp(t2)
str2fp ry, Str$(t2) 'convert the double t2 to decfloat in ry
t = Log(10) * ex / p
t2 = Exp(t - Fix(t))
str2fp tmp, Str$(t2) 'convert the double t2 to decfloat in tmp
fpmul ry, ry, tmp, prec
str2fp tmp, "2.7182818284590452353602874713527"
fpipow tmp, tmp, Fix(t), 24
fpmul ry, ry, tmp, prec
basically it does the following
here the value of x is split into a fractional part and a base 10 exponent, so if x = 3.14159e2000 then frac_x=3.14159 and ex=2000
first approximation =
exp(log(frac_x) / p) * exp(fix(ex/p)) * exp(frac(ex/p)) 'here frac is the fractional part of expression
for example
31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989 ^ (1/3)
I have split the first approximation to allow for larger exponents than double allows
instead of t2 = Log(t) / p + Log(10) * ex / p, do
t2 = Log(t) / p
t2 = Exp(t2)
str2fp ry, Str$(t2) 'convert the double t2 to decfloat in ry
t = Log(10) * ex / p
t2 = Exp(t - Fix(t))
str2fp tmp, Str$(t2) 'convert the double t2 to decfloat in tmp
fpmul ry, ry, tmp, prec
str2fp tmp, "2.7182818284590452353602874713527"
fpipow tmp, tmp, Fix(t), 24
fpmul ry, ry, tmp, prec
basically it does the following
here the value of x is split into a fractional part and a base 10 exponent, so if x = 3.14159e2000 then frac_x=3.14159 and ex=2000
first approximation =
exp(log(frac_x) / p) * exp(fix(ex/p)) * exp(frac(ex/p)) 'here frac is the fractional part of expression