В этом параграфе описаны алгоритмы, позволяющие, в принципе, вычислять с произвольной точностью: ошибки всегда возникают из-за конечной разрядности. Алгоритмы основаны на разложении в ряды Тейлора и цепные дроби функций, связанных с простыми соотношениями. Структура всех этих алгоритмов описывается в приложении А; коды на Си содержатся в приложениях В и Г. Вероятно, наиболее известными соотношениями, позволяющими вычислять с произвольной точностью, являются = 0.5 + I(x), (1) . Это разложение легко получить, интегрируя почленно разложение Тейлора в нуле функции exp(-x2/2). Хотя при этом получается знакопеременный ряд, первый отброшенный его член не является, вообще говоря, верхней границей для остатка, т.к. его члены монотонно убывают, лишь начиная с некоторого. Тем не менее, довольно легко убедиться в том, что при =0.3989... это приятное свойство имеет место. А так как в любой практической задаче <<0.1, можно считать, что если выполнено неравенство (2) , то сумма n первых членов ряда доставляет I(x) с точностью (если, конечно, не учитывать накопление ошибки). АлгоритмI. Этот алгоритм основывается на разложении (1). Обратите внимание: при n-м выполнении шага I2 переменная m равна 2n, так что t здесь делится все-таки на 2nn!. Разложение (1) было использовано в алгоритме 272 [12], который, однако, значительно менее экономен, чем нижеследующая версия. I1.Положить t=sum=x; x2=x*x; m = 2. I2.Положить t=-x2*t/m; s=sum; sum=s+t/(m+1); m=m+2. I3.Если |s-sum| > e, то перейти к I2. I4.Положить F = 0.5 + sum/Ц2p. Менее известными, но столь же просто получаемыми, как (1), являются соотношения (3) , где . Это разложение для можно получить почленным интегрированием ряда Тейлора функции j(x)/j(t) в нуле. Для остатка ряда (3) справедливы оценки , поэтому оценка Fnошибки ограничения для F(x) удовлетворяет неравенству (4) . Неравенства (2) и (4) показывают, что скорость сходимости разложений (1) и (3) падает с ростом x. Алгоритм T. Этот алгоритм основан на разложении (3). В работах [1, 11] в аналогичном алгоритме на шаге T1 производится присвоение t = sum = |x|*exp(-x2/2)/Ц2pс очевидным изменением шага T4. Так как, однако, шаг T3 там не изменен, условием останова оказывается |sn+1 - sn| < e*exp(x2/2)/Ц2p, что, согласно (4), не гарантирует необходимую точность вычисления F(x). Тем не менее, если вычисление F(x)производится с машинной точностью (и машинный нуль достаточно мал) и ошибка при вычислении экспоненты не слишком велика, алгоритмы в [1, 11] дают в пределах представления чисел с плавающей запятой правильный результат. T1. Положить t = sum = |x|; x2 = x*x; n = 3. T2. Положить t = x2*t/n; s = sum; sum = s+t; n = n+2. T3. Если |sum - s| > e , перейти к T2. T4. Положить F = 0.5 + sign(x)*sum*exp(-x2/2)/Ц2p . Вот еще одно разложение функции , которое можно применять для вычисления F(x)при небольших x: (5)  Эту цепную дробь можно получить разными способами (см., например, [2,3]). Члены ее звеньев альтернируют и потому характер сходимости ее подходящих дробей довольно сложен. Однако, для подпоследовательности подходящих дробей с четными номерами можно показать, что ее четные члены сходятся к монотонно возрастая, а нечетные - монотонно убывая. Поэтому модуль разности |w2n-2-w2n| двух соседних подходящих дробей этой подпоследовательности оценивает сверху уклонение каждой из них от . Алгоритм S. В данном алгоритме, основанном на разложении (5), для получения подпоследовательности подходящих дробей с четными номерами мы производим на каждом витке цикла вычисление двух очередных подходящих дробей; соответствующее рекурсивное соотношение для сжатой дроби слишком сложно и потому не использовано. S1. Положить x2 = x*x; S2. /* Вычисление подходящих дробей с номерами 0 и 2:         предпервое состояние цикла */      rho = 1/(3-x2); sum = |x|*rho;      term = x2*sum; rho=3*rho; s = sum = 3*sum;      n = 1; S3. /* Вычисление четной (когда x2 < 0)         или нечетной (x2 > 0) подходящей дроби */      Положить n = n+1; r = x2*m/(4*n*n-1);      rho=1/(1+r*rho); term=(rho-1)*term;      t = s; s=sum; sum=sum+term;      x2=-x2; S4. Если x2 < 0 (т.е. на шаге S3вычислялась нечетная     подходящая дробь), то перейти к шагу S3. S5. /* Проверка точности */    Если |s-t| > e или |s-sum| > e, перейти к S3. S6. Положить F = 0.5 + sign(x)*sum*exp(-x2/2)/Ц2p;      на этом работа алгоритма заканчивается. Соотношения [13] (6') wn-1-wn = (-1)n(4n-1)(2n-2)! x4n-1/(q2n-2q2n), (6")  показывают, что необходимое число подходящих дробей растет с ростом x. Поэтому разложение следует применять лишь при небольших значениях x; сам Л.Шентон [13] рекомендует применять его только при x <Ц3. До сих пор мы имели дело с разложениями, "качество" которых ухудшалось с ростом x. Рассмотрим теперь соотношения, которые с ростом x позволяют вычислять F(x) все быстрее: (7') F(x) = 0.5 + j(x)R(x), (7")  Функция R(x) называется отношением Миллса. Разложение в цепную дробь для него было предложено еще Лапласом в 1805 г. Члены звеньев этой цепной дроби при x > 0 положительны и потому (см., например, [6]) подходящие дроби с четными номерами сходятся к R(x) монотонно возрастая, а с нечетными - монотонно убывая. Поэтому модуль разности двух соседних подходящих дробей этой подпоследовательности оценивает сверху уклонение каждой из них от R(x). Можно показать, что имеет место оценка (8) , из которой следует, что скорость сходимости разложения растет с ростом x. Алгоритм L. В этом алгоритме F(x)вычисляется с помощью разложения отношения Миллса R(x) в цепную дробь Лапласа (7). L1. Положить x2=x*x; y=exp(-x2/2)/Ц2 p; x2=1/x2;       sum=term=y/|x|; r=s=0; rho=1. L2. /* Вычисление очередной подходящей дроби */     Положить r=r+x2; rho=1/(1+r*rho); term=(rho-1)*term;       t=s; s=sum; sum=sum+term. L3. /* Проверка точности */     Если |s-t| > e или |s-sum| > e , перейти к L2. L4. Если x > 0, положить F = 1-sum и остановить алгоритм. L5. Положить F = sum. Для вычисления F(x) при больших x используют также асимптотический ряд (9) , который легко получить, проинтегрировав отношение Миллса R(x) по частям. По заданному e можно найти xe, такой, что при x ? xe, ряд (9) позволяет вычислять F(x) с точностью e. Однако, с помощью этого ряда невозможно, вообще говоря, достичь в любой заданной точке x произвольной точности, поэтому мы и не исследуем алгоритм, основанный на этом разложении. Для проведения численных экспериментов алгоритмы I, T, S и L были реализованы на Си. В нижеследующей таблице для каждого алгоритма в столбце N приводится количество "витков" цикла, потребовавшихся для вычисления значения с e = 0 для нескольких x; поскольку все алгоритмы дали одинаковые 8 знаков после запятой, сами вычисленные значения опущены. В столбце O приводится число потребовавшихся операций. Для алгоритма I оно вычислялось по формуле 8+11?N, для алгоритма T - по формуле 14+8?N, для S - по формуле 23+25?N, наконец, для L - по формуле 16+10?N. Обратите внимание - всем операциям, в том числе, таким, как вызов функции fabs или exp, придается равный вес. Для иллюстративных целей такая точность признана достаточной. x I T S L   N O N O N O N O 0.25 8 96 8 78 8 223 6027 60286 0.50 10 118 10 94 10 273 1507 15086 0.75 12 140 13 118 12 323 694 6956 1.00 14 162 15 134 14 373 391 3926 1.25 16 184 17 150 14 373 259 2606 1.50 18 206 19 166 16 423 180 1816 1.75 20 228 21 182 18 473 137 1386 2.00 23 261 23 198 20 523 108 1096 2.25 25 283 25 214 22 573 88 896 2.50 27 305 27 230 22 573 74 756 2.75 30 338 29 246 24 623 63 646 3.00 33 371 32 270 26 673 54 556 3.25 35 393 34 286 28 723 48 496 3.50 38 426 36 302 30 773 44 456 3.75 41 459 39 326 32 823 39 406 4.00 44 492 41 342 34 873 35 366 4.25 48 536 44 366 36 923 33 346 4.50 51 569 47 390 38 973 31 326 4.75 55 613 49 406 40 1023 28 296 5.00 58 646 52 430 42 1073 26 276 5.25 62 690 54 446 44 1123 25 266 5.50 66 734 57 470 48 1223 24 256 5.75 70 778 60 494 50 1273 22 236 6.00 74 822 63 518 52 1323 21 226 6.25 79 877 66 542 54 1373 21 226 6.50 83 921 69 566 58 1473 20 216 6.75 88 976 72 590 60 1523 19 206 7.00 93 1031 76 622 64 1623 18 196 7.25 98 1086 79 646 66 1673 18 196 7.50 103 1141 82 670 70 1773 17 186 7.75 108 1196 86 702 72 1823 17 186 Построим графики, отражающие данные из этой таблицы, - их проще анализировать. На нижеследующих рисунках слева - столбцы N, справа - O. а) число витков цикла б) количество операций Мы видим, что по числу N витков цикла алгоритмы при не слишком больших x практически не различаются до x'3.5, после чего лидирует алгоритм S. Начиная с x' 4, бесспорным лидером становится алгоритм L. По числу операций, конечно, более важной характеристике, до x' 4 лидером является алгоритм T, после - снова алгоритм L. Помимо превосходства в скорости алгоритм T выгодно отличается от I и S устойчивостью по отношению к накоплению ошибки. Для доказательства высказанного утверждения оценим распространяемую ошибку округления. Пусть -относительная ошибка n-й частичной суммы в алгоритме I, а - относительная ошибка tn. Для простоты предполагаем, что представление целых чисел и действия над ними осуществляются без ошибок, а ошибки округления всех арифметических операций над действительными числами не превосходят r. Тогда для имеем следующее соотношение (10) , причем (11)  где - относительная ошибка представления x2.Будем считать также, что , , . Из (11) с легкостью следует , а из (10) получаем , откуда (12) . Верхняя оценка ошибки (12) растет с пугающей скоростью: так, если r=10-10, то , а .При этом следует учесть, что при выводе мы игнорировали знакопеременность ряда (1), а это сильно смягчает грубость оценки (12). Для алгоритма T относительная ошибка частичной суммы sn также удовлетворяет соотношению (10), однако теперь . Действуя, как раньше, получаем оценку (13) . Эта оценка возрастает лишь немного медленнее, чем оценка (12), Однако, члены ряда (6) положительны, поэтому истинная относительная ошибка значительно меньше полученной оценки (13). Хотя общий член ряда (2) убывает быстрее, чем общий член ряда (3), машинные эксперименты показывают, что число исполнений шага T2 меньше числа 'витков цикла' в алгоритме I; это также приводит к меньшей распространяемой ошибке. Тем не менее, на точность вычисления F(x) сильно влияет точность вычисления exp(-x2/2): возникающая здесь ошибка может 'съесть' все преимущества разложения (3). Полученные данные позволяют рекомендовать для вычисления F(x)с произвольной точностью сочетание алгоритмов T и L. Нужно только выбрать точку p, разделяющую области действия этих алгоритмов. Итак, при возрастающих x время счета по алгоритму T возрастает, а по алгоритму L убывает. Поэтому в качестве p следовало бы выбирать точку, в которой эти времена становятся равными. К сожалению такая точка зависит от точности счета e , от разрядности используемого компьютера, от его скорости. Поэтому выбрать навеки оптимальную' точку p не удастся; я рекомендую выбирать ее при e , равном машинному нулю.  8  Комментарии к статье  8 8  Обсудить в чате