UDFSoftBlog

Конфигурация CDLC 1.1 позволяет работать с вещественными числами, поддерживая тип double. Однако стандартная библиотека Math включает в себя очень скудный набор математических функций: sin, cos, tan, sqrt. Существует несколько сторонних математических классов (например, Real.java, которую можно скачать по адресу http://sourceforge.net/projects/real-java). В этой статье я предлагаю написать собственную реализацию нескольких популярных математических функций.

В основе компьютерной математики лежит использование разложений функций в математические ряды. Подробно этот вопрос рассматривается в курсе Математического анализа любого Вуза. Если вкратце, то математическую функцию можно представить в виде бесконечной суммы слагаемых, причем каждое следующее слагаемое по модулю меньше предыдущего. Поэтому для вычисления функции с заданной точностью нужно выполнять сложение до тех пор, пока следующее слагаемое не станет меньше, чем требуемая точность вычисления.



Вычисление Экспоненты

Из курса математического анализа известно, что экспоненту можно представить в виде бесконечного ряда

.

При этом, чем больше аргумент x, тем больше слагаемых требуется брать для удовлетворения требуемой точности. При x>706 exp(x) не умещается в переменную типа double, поэтому перед вычислением экспоненты целесообразно проверить значение x на превышение порогового предела. Для реализации эффективного алгоритма вычисления экспоненты, нужно воспользоваться известным из курса школьной алгебры фактом:



Очевидно, число x<706 можно представить в виде суммы:

где коэффициенты a могут принимать значения 0 или 1, а b<1.



Величины можно вычислить заранее и записать в массив. Затем вычленить из x целую(b) и дробную часть. Для дробной части вычислить экспоненту как сумму ряда Тейлора. Найти ai можно, переведя целую часть числа x в двоичную систему. Тогда самому правому символу в двоичном представлении числа соответствует a0, второму справа числу – a0, и так далее. Ниже приведен код функции, вычисляющей экспоненту числа:

 

private double MExp(double x0){ 
    double x=x0;
    if (x0<0){x=-x0;} 
   //Массив со степенями экспоненты. 
   double ExpConst[]={
        2.718281828459045, //e^1
        7.389056098930649, //e^2
        54.59815003314422, //e^4
        2980.957987041726, //e^8
        8886110.520507860, //e^16
        78962960182680.61, //e^32
        6.235149080811582e27, //e^64
        3.887708405994552e55, //e^128
        1.511427665004070e111, //e^256
        2.284413586539655e222 //e^512 
    };
   int x1=(int)x; //Выделяем целую часть числа
   //Вычисляем экспоненту как сумму ряда Тейлора
   int long n=1;
   double b=1;
   double sn=1; 
   while (sn>1E-16){
        sn=sn(x-x1)/n;
        b=b+sn;
        n=n++;
   }
   //Переводим показатель экспоненты в двоичную систему.
   StringBuffer s1=new StringBuffer(10);
   s1.append(Integer.toBinaryString(x1)); 
   int len=s1.length();
   for (int i=s1.length(); i>0;i--)
   {
       if (s1.charAt(i-1)==\'1\'){b=bExpConst[len-i];}
   }
    if (x0<0){b=1/b;}
    return b;
}

Гиперболические функции

Имея функцию для вычисления экспоненты, легко найти значения гиперболических функций.

, , .

Имеет смысл предварительно записать значение exp(x) во вспомогательную переменную, которую затем использовать для вычисления по формулам.
Вычисление натурального логарифма

Натуральный логарифм можно представить в виде суммы:



Вычислять значение логарифма непосредственно по этой формуле не очень эффективно. Для оптимизации алгоритма вычисления можно воспользоваться известными соотношениями

Представим число x в виде



где b<1, a-целое, тогда

.

Для того чтобы представить число x в требуемом виде, нужно найти позицию символа “E” в строковом представлении числа x, тогда a равно числу, записанному правее символа “E” плюс 1, a b – числу записанному левее “E”, деленному на 10.

Ниже представлен код, реализующий алгоритм вычисления натурального логарифма:

 

private double MLn(double x0){
    double x=x0;
    double y=0;
        //Получаем показатель степени.
        String s0=""+x;
        int i=s0.indexOf("E");
        String s1=s0.substring(i+1, s0.length());//Правее E 
        String s2=s0.substring(0, i);//Левее E
        double a=0,b=0; 
        a=Double.parseDouble(s1)+1;
        b=Double.parseDouble(s2)/10; 
  //Вычисление Логарифма b как суммы ряда Тейлора
       int n=1;
       double sn=1;
       while (sn>(1E-16)n){
          sn=-sn(b-1);
          y=y+sn/n;
          n=n++;
      }
    y=y+a2.302585092994046;
    return y;
}

Вычисление корней, степеней и логарифмов

Имея функцию для вычисления натурального логарифма, легко вычислить следующие функции:


Вычисление арксинуса и арккосинуса

Вычисление обратных тригонометрических функций arcsin и arcos сводится к вычислению рядов.

Ниже приведен алгоритм вычисления арксинуса.

 

private double Marcsin(double x0){
    double x=x0;
    if (x0<0){x=-x0;}
    double y=x;
    int n=1;
    double sn=x; 
    while (sn>1E-16){
            sn=sn(2+1.0/n)0.5xx;
            y=y+sn/(2n+1)/(2n+1); 
            n=n+1;
    }
    if (x0<0){y=-y;}
    return y; 
}

Вычисление арктангенса 

Предложенная формула эффективна для вычисления арктангенса малого аргумента. Для построения алгоритма эффективного вычисления арктангенса от больших x, целесообразно воспользоваться формулой:


Ниже предлагается переписанный для java алгоритма, который написал товарищ Nikitin V.F. в 2000 году.
Вначале нужно проверить знак x, если x<0, нужно изменить знак, сделав аргумент неотрицательным.
Затем если x>1, обратить его: x=1/x.
Затем сокращаем область определения (запоминая число сделанных шагов), используя формулу

до тех пор, пока x не окажется в интервале [0,pi/12]. 
После этого можно вычислять y=arctg(x) по формуле Тейлора.
Сместим результат y на p/6 необходимое число раз
Если x было больше 1, то результат 
Если x было отрицательным, то результат 

 

private double MArctg(double x0) {
  int sp=0;
  double x,x2,y;
  x=x0;
  if(x<0) {x=-x;}
  if(x>1) {x=1.0/x;}
  //Уменьшаем интервал области аргумента
  while(x>0.2617993877991495) {
    sp++; //Вспомогательный счетчик шагов 
    x=(x1.732050807569-1)/(x+1.732050807569);
  }
  //Вычисляем ряд Тейлора
    y=x;
    int n=1;
    double sn=x; 
    while (sn>1E-16){
            sn=sn(2+1.0/n)0.5xx;
            y=y+sn/(2n+1)/(2n+1); 
            n=n+1;
    }
        
  //Смещаем все на pi/6 необходимое число раз 
  y=y+sp0.523598775598 
 
  if(x0>1) a=0.2617993877991495-a;
  if(x0<0) y=-y;
 
  return y;
}

08 мая 2008 г.