Re: yet another question on Mathematica numerical, how to force it to use single precision?
- From: Daniel Lichtblau <danl@xxxxxxxxxxx>
- Date: Mon, 17 Nov 2008 14:44:34 -0800 (PST)
On Nov 17, 3:07 pm, "Nasser Abbasi" <n...@xxxxxxxxx> wrote:
How can I make Mathemtica give me the same digits to this little sum below
as I get from Fortran or Matlab?
I tried using just N[], which is should use machine precision. I want to
tell it to do the computation in "single precision", i.e. 32 bit as compared
to double precision.
Could some one remind me how to do this?
This is the Mathematica code
---------------------
In[95]:= s = Sum[N[1/k^2], {k, 1, 20000}]
Precision[s]
Accuracy[s]
Out[95]= 1.644884068098209
Out[96]= MachinePrecision
Out[97]= 15.738454476026273
----------------------
This is what I get from Fortran
PROGRAM main
IMPLICIT NONE
REAL :: s
INTEGER :: n,MAX
s = 0.0;
MAX = 20000;
DO n = 1,MAX
s = s + (1./n**2);
END DO
WRITE(*,*), 'sum', s
END PROGRAM main
$ f77 sum_problem.f
$ ./a.exe
sum 1.64472532
------------------------------
This is what I ge from Matlab (by forcing matlab to use single precision
format long g;
s=single(0);
MAX = 20000;
for i=1:MAX
s=s+(1/i^2);
end
s
s =
1.644725
So, Fortran and Matlab generate the same first 6 digits beyond the decimal
point. I want Mathematica to generate those same digits.
Any ideas?
thanks,
Nasser
I do not know how to hit this value dead on, but the code below comes
close to emulating this. I use a certain binary precision, and
truncate sums when summands are smaller than that precision times the
running sum. I also force, in effect, fixed precision arithmetic (this
alone does not quite suffice because it seems to keep too many guard
digits).
myAdd[t1_,t2_,_] /; t1==0 :=t2
myAdd[t1_,t2_,_] /; t2==0 :=t1
myAdd[t1_,t2_,bprec_] := With[{ratio=Abs[t1/t2]},
If[ratio>2^(bprec),t1,
If[ratio<2^(-bprec),t2,t1+t2]]]
sumAtPrecision[prec_,term_,index_,lo_,hi_] := Block[
{s=0,l10prec=Log[10,2]*prec},
Block[{$MinPrecision=l10prec,$MaxPrecision=l10prec},
Do[s=N[myAdd[s,term,prec],l10prec], {index,lo,hi}]
];
s]
If I specify 25 or 26 bits, I get results close to and bracketing what
you indicate.
In[62]:= sumAtPrecision[25,1/k^2,k,1,20000]
Out[62]= 1.6447127
In[63]:= sumAtPrecision[26,1/k^2,k,1,20000]
Out[63]= 1.6447775
Daniel Lichtblau
Wolfram Research
.
- References:
- Prev by Date: Re: yet another question on Mathematica numerical, how to force it to use single precision?
- Next by Date: Solving solvable Polynomials of Degree 11
- Previous by thread: Re: yet another question on Mathematica numerical, how to force it to use single precision?
- Next by thread: Re: yet another question on Mathematica numerical, how to force it to use single precision?
- Index(es):
Relevant Pages
|
