5

Define the following four rational numbers. $$ a = \frac{4243257079864535154162785598448178442416}{41016602865234375} \\ b = -\frac{143308384621912247542172258992236503771301}{1210966757832031250} \\ c = \frac{350687399375274064088342133116344593437371021}{4109863607096484375000} \\ d = -\frac{762284492856611655417326017768244278005511063}{12085448243163671875000} $$ Let $$ p = \frac{501}{10}, \qquad m = \frac{499}{10}. $$ Compute $$ \mathrm{Result} = a \cos [p] + b \cos [m] + c \sin [p] + d \sin [m]. $$ Each term in this sum is roughly $10^{23}$. There is a curious cancellation (using 40 digits) happening amongst these four terms; the correct answer is $\mathrm{Result}=7.32 \times 10^{-18}$.

My question: Where does this cancellation come from, analytically? Can you massage the terms into a form where the cancellation is manifest, and machine precision can evaluate the answer with a semblance of accuracy?

(If you're curious, the result came from the analytic integral of a highly oscillatory function.)

Tito Piezas III
  • 60,745
Jolyon
  • 242
  • Can you refer any web link on this topic? – MAN-MADE Jun 22 '17 at 04:10
  • @MANMAID No sorry. Best I can do is give you the integral it came from. Integrate[x^6 SphericalBesselJ[10, a x] SphericalBesselJ[10, b x], x]. I can't remember what values I used for the constants or limits of integration, but it was something like b/a ~ 500. – Jolyon Jun 22 '17 at 04:17
  • @Jolyon: I trimmed some extra detail so as not to distract from your result. – Tito Piezas III Jun 22 '17 at 05:19

1 Answers1

1

Here's two ideas that might be useful. They are in no way a solution.

Idea 1:

Let $r = 50, s=\frac1{10}$. Then $p = r+s, m = r-s$ so that $\cos(p) =\cos(r+s) =\cos(r)\cos(s)-\sin(r)\sin(s), \cos(m) =\cos(r-s) =\cos(r)\cos(s)+\sin(r)\sin(s), \sin(p) =\sin(r+s) =\sin(r)\cos(s)+\cos(r)\sin(s), \sin(m) =\sin(r-s) =\sin(r)\cos(s)-\cos(r)\sin(s) $.

Therefore

$\begin{array}\\ \mathrm{Result} &= a \cos [p] + b \cos [m] + c \sin [p] + d \sin [m]\\ &= a (\cos(r)\cos(s)-\sin(r)\sin(s)) + b (\cos(r)\cos(s)+\sin(r)\sin(s))\\ &\quad + c (\sin(r)\cos(s)+\cos(r)\sin(s)) + d (\sin(r)\cos(s)-\cos(r)\sin(s))\\ &=\cos(r)(a\cos(s)+b\cos(s)+c\sin(s)-d\sin(s))\\ &\quad+\sin(r)(-a\sin(s)+b\sin(s)+c\cos(s)+d\cos(s))\\ &=\cos(r)((a+b)\cos(s)+(c-d)\sin(s)) +\sin(r)((-a+b)\sin(s)+(c+d)\cos(s))\\ &=\cos(r)(w\cos(s+x)) +\sin(r)(y\sin(s+z))\\ \end{array} $

where $w^2 =(a+b)^2+(c-d)^2, y^2 =(-a+b)^2+(c+d)^2, x =\arctan(\frac{c-d}{a+b}), z =\arctan(\frac{-a+b}{c+d}) $.

Computing these values might be interesting.

Idea 2:

$\begin{array}\\ \mathrm{Result} &= a \cos [p] + b \cos [m] + c \sin [p] + d \sin [m]\\ &= a \cos [p] + c \sin [p]+ b \cos [m] + d \sin [m]\\ &= u \sin [p+v] + w \sin [m+x]\\ \end{array} $

where $u^2 = a^2+c^2, w^2 =b^2+d^2, v = \arctan(\frac{a}{c}), x = \arctan(\frac{b}{d})$.

I don't know where to go from here.

marty cohen
  • 110,450
  • So, w=1.4910^23 and y=2.2310^23. Cos(50)=0.96, Sin(50)=-0.26. Nothing interesting there. However, x=-1.47079632679489661923132169163975144209858446 to 45 digits, and z = -1.47079632679489661923132169163975144209871044 (also to 45 digits). Note that these vary only in the 41st digit! I can see writing z = x + delta, but still not sure where the cancellation results from after that. – Jolyon Jun 22 '17 at 04:34
  • In the second idea, u and w both equal 1.341016829893478479580822451711498255166*10^23 to 40 digits (the last digit is a 5 for u, 6 for w). The two arctans are less interesting; v=0.88, x=1.08. However, p+v = 50.9811102752095508810978862056203824512371 to 42 digits, and m+x is the same, except for the last two digits being 69. So, this gets really close to explaining where the cancellation is arising from. Not so close on being able to do anything with machine precision. – Jolyon Jun 22 '17 at 04:44
  • One extra note - on idea 2, the negative root is needed for w. The two sine terms are then essentially identical, so the problem boils down to taking the difference of u and w. There's probably a clever way of doing so... – Jolyon Jun 22 '17 at 06:25