Issues with machine number consistency in 7
- To: mathgroup at smc.vnet.net
- Subject: [mg94768] Issues with machine number consistency in 7
- From: budgemano <mdalpiaz at email.arizona.edu>
- Date: Tue, 23 Dec 2008 07:01:07 -0500 (EST)
We have code that is designed to perform a ray trace calculation and compare the results to a look up table. This worked great in 6, but once we switched to 7 things stopped working. Running the code numerous times produces slightly different results in a seemingly random fashion, though on the machine level it seems. A snippet of the code that does this is below (and it looks messy here I know), and it is only linear algebra type calculations, so nothing iterative I would think. The question is: how does this produce slightly different results with each run, and even more so what happened between 6 and 7 that makes this happen?
And yes the accuracy is most likely fine, in that a few parts per quadrillion is not a big deal, but I would like to understand the reason before just slapping a bandage on the problem.
CreatSP[n_, k_] := Module[{nhat, khat, s, shat, p, phat},
nhat = n/Norm[n];
If[k == {0, 0, 0}, khat = k, khat = k/Norm[k]];
If[Norm[Cross[nhat, khat]] == 0,(*if rays||normal,
pick arbitrarily*)
If[khat == {0, 0, 1} || khat == {0, 0, -1},
s = Cross[khat + {0, 1*^-15, 0}, nhat],
s = Cross[khat + {0, 0, 1*^-15}, nhat]],
s = Cross[khat, nhat]];
shat = s/Norm[s];
p = Cross[khat, shat];
If[p == {0, 0, 0}, phat = p, phat = p/Norm[p]];
{shat, phat}];
JonesTo3D[
j_] := {{j[[1, 1]], j[[1, 2]], 0}, {j[[2, 1]], j[[2, 2]], 0}, {0,
0, 0}};
GlobalToLocal[n_, k_] := Module[{shat, phat, nhat, khat},
nhat = n/Norm[n];
If[k == {0, 0, 0}, khat = k, khat = k/Norm[k]];
{shat, phat} = CreatSP[nhat, khat];
Return[{shat, phat, khat}]];
KReflect[n_, k_] := Module[{vs, vpl, nhat, khat},
If[n == {0, 0, 0}, nhat = n, nhat = n/Norm[n]];
If[k == {0, 0, 0}, khat = k, khat = k/Norm[k]];
If[nhat == khat || nhat == -khat, -k,
(*else*){vs, vpl} =
CreatSP[khat, nhat];(*feed the vectors in reverse order,
to get s& p projected on the surface*)
Dot[k, vpl]*vpl -
Dot[k, n]*n (*output k vector*)] (*end if*)
]
Reflect3D[n_, kin_, j_] := Module[{kout},
kout = KReflect[n, kin];
Transpose[GlobalToLocal[n, kout]].JonesTo3D[j].GlobalToLocal[n,
kin]
]
Surf3DReflect[n_, kin_, j_] := {KReflect[n, kin],
Reflect3D[n, kin, j]};
SetArithmetic[4, RoundingRule -> RoundToEven] ;
nhat = {0.186658373151350784576457645764576, \
-0.908774534497317234675674568458,
0.3732124022345396456745768874687} // N;
kIn = {0.5018839897852984567457648, 0.544293499215831476457645674576,
0.67220312964800834576456745674567};
J = {{-4.094099334937029435673467356734 +
1.474574404164915678487487 I, -1.10844925387366011246358345623 -
5.4374686620520058634643181 I}, \
{8.179637040999307235787945456132 + 2.51644639657543619834750187643 I,
2.331233134964963234512765301876094 -
7.360368412868619230947509283762 I}};
For[ii = 1, ii <= 20, ii++,
If[
Surf3DReflect[nhat, kIn, J] =!= Surf3DReflect[nhat, kIn, J]
,
Print["Iteration " <> ToString[ii] <> " fails comparison "]
];
]