[iiv 16]

Добрый день и с наступающим Вас Новым Годом!

 

У меня есть template (кордик синус и косинус), который в зависимости от того float это или double должен считать соответсвующие пары синусов и косинусов. Для этого мне надо в зависимости от того float это или double

1. правильно проинициализировать целочисленную переменную BL,
2. сделать int IntA или long long IntA и вычислить ее,
3. предварительно заполнить BinarySin, BinaryCos предвычисленными числами.

Я как-то это на С++ запрограммировал, но не все мне нравиться, а именно

1. Пункт: хочется реально сравнивать только на float и double и в остальных случаях давать ошибку,
2. Пункт: я все привел к типу long long, и не понимаю как сделать int для float, и long long для double,
3. Пункт: я хочу чтобы инициализация этих массивов происходила на стадии компиляции, а сейчас у меня это происходит в runtime.

Пожалуйста, посоветуйте, что и как тут изменить, чтобы это все исправить!

 

Спасибо!

 

Код тут:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>


template<typename T> void CalcSinCos(T a, T &s, T &c) // a должно быть больше 0. и меньше 1.
{ const int BL = ((sizeof(T)==4))?24:53;
  long long IntA;

  if constexpr (sizeof(T)==4)
    IntA=(int)(a*(1024.*1024*16.));
  else
    IntA=(long long)(a*(1024.*1024.*1024.*1024.*1024.*8.));

static T BinarySin[BL];
static T BinaryCos[BL];
static int NotYetInizialized=0;

  if(NotYetInizialized==0)
  { NotYetInizialized=1;
    float a=0.5;

    for(int i=0; i<BL; i++)
    { BinarySin[i]=sin(a);
      BinaryCos[i]=cos(a);
      a*=0.5; // a/=2.;
    }
  }

  s=0.;
  c=1.;

  for(int i=0; i<BL; i++)
    if((IntA&(1<<(BL-1-i)))!=0)
    { // s, c на данный момент - это sin(a), cos(a), а BinarySin[i] и BinaryCos[i] - это sin(b), cos(b)
      // результат sin(a+b) и cos(a+b) нам надо снова получить в s, c
      T temp_s = s*BinaryCos[i] + BinarySin[i]*c; // sin(a+b) = sin(a)*cos(b) + sin(b)*cos(a)
      T temp_c = c*BinaryCos[i] - s*BinarySin[i]; // cos(a+b) = cos(a)*cos(b) - sin(a)*sin(b)
      s = temp_s;
      c = temp_c;
    }
  return;
}


int main()
{

  for(int i=0; i<20; i++)
  { float a = ((float)rand()) / (float)RAND_MAX;
    float s, c;
    CalcSinCos(a, s, c);
    printf("a=%g, s=%g (%g), c=%g (%g)\n", a, s, sin(a), c, cos(a));
  }

  return 0;
}

 

[x893 34]

Это шаблон. вызываете с float будет сделан класс для float.

Для double будет другой класс.

 

[Forger 17]

2 hours ago, iiv said:

посоветуйте, что и как тут изменить, чтобы это все исправить!

Шаблон тут видится как "собаке - пятая нога", поскольку у вас тут предполагается всего ДВЕ инстанции шаблона. Сделайте проще, так сказать "вручную", в стиле if else. 

 

 

[krux 8]

не помню где я это обсуждал, но повторю еще раз здесь.

вся тригонометрия, все-все sin и cos которые заложены в x86 стандартных библиотеках вычисляют его по полиномам(многочленам) Чебышева с коэффициентами, определяющими необходимую точностиь (они тупо зашиты в константы в шестнадцатиреричной форме чтобы не переполнять IEEE754). А дальше оптимизация до уровня ассемблера ведется стандартными методами, не учитывающми буквально ничего.

Думаю необходимая и достаточная информация тут есть.

 

Короче, кордик, работающий бьстрее стандартных функций и с необходимой точностью можно сделать намного проще. ИМХО

[iiv 16]

Огромное спасибо за советы!

 

Формально - да, только две инстанции, но я-то хотел туда же еще добавить свои типы, которых у меня суммарно 6.

 

Пока сделал как ниже, что порешило мой пункт 3, в ассемблере теперь видно, что синус и косинус вызываются один раз и есть предвычисленные константы.

 

Но вот с первыми моими двумя пунктами пока еще не придумал как сделать.

 

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <array>

template<typename T>constexpr auto generate()
{ const int BL = ((sizeof(T)==4))?24:53;
  std::array<std::array<T, BL>, 2> res{};

  T a=0.5;

  for(int i=0; i<BL; i++)
  { res[0][i]=sin(a);
    res[1][i]=cos(a);
    a*=0.5; // a/=2.;
  }
  return res;
}


template<typename T> void CalcSinCos(T a, T &s, T &c) // a должно быть больше 0 и меньше 1
{ const int BL = ((sizeof(T)==4))?24:53;
  long long IntA;

  if constexpr (sizeof(T)==4)
    IntA=(int)(a*(1024.*1024*16.));
  else
    IntA=(long long)(a*(1024.*1024.*1024.*1024.*1024.*8.));

  constexpr auto BinarySinCos = generate<T>();

  s=0.;
  c=1.;

  for(int i=0; i<BL; i++)
    if((IntA&(1<<(BL-1-i)))!=0)
    { // s, c на данный момент - это sin(a), cos(a), а BinarySin[i] и BinaryCos[i] - это sin(b), cos(b)
      // результат sin(a+b) и cos(a+b) нам надо снова получить в s, c
      T temp_s = s*BinarySinCos[1][i] + c*BinarySinCos[0][i]; // sin(a+b) = sin(a)*cos(b) + sin(b)*cos(a)
      T temp_c = c*BinarySinCos[1][i] - s*BinarySinCos[0][i]; // cos(a+b) = cos(a)*cos(b) - sin(a)*sin(b)
      s = temp_s;
      c = temp_c;
    }
  return;
}


int main()
{

  for(int i=0; i<20; i++)
  { float a = ((float)rand()) / (float)RAND_MAX;
    float s, c;
    CalcSinCos(a, s, c);
    printf("a=%g, s=%g (%g), c=%g (%g)\n", a, s, sin(a), c, cos(a));
  }

  return 0;
}

 

20 minutes ago, krux said:

вся тригонометрия, все-все sin и cos которые заложены в x86 стандартных библиотеках вычисляют его по полиномам(многочленам) Чебышева с коэффициентами, определяющими необходимую точностиь (они тупо зашиты в константы в шестнадцатиреричной форме чтобы не переполнять IEEE754). А дальше оптимизация до уровня ассемблера ведется стандартными методами, не учитывающми буквально ничего.

Спасибо! Это понятно, что Чебышевым получается обычно не плохо. Другое дело, что кордик по основанию 16 или 256 (не по основанию 2, как у меня в примере) с Тейлором для маленьких составных кордика реально в 2-3 раза быстрее, чем стандартный вызов системной sincosf. А на некоторых МК ой как плохо в стандартных библиотеках все реализовано, что кордик даже по основанию 2 начинает обыгрывать системные функции. А у меня еще есть float-float тип, для которого стандартных синусов-то и нет.

[krux 8]

On 12/31/2022 at 5:11 PM, iiv said:

 

Современные штуки, разработчики которых не понимают что они стали настолько ненадежными что их лучше заменить таблицами Брадиса.

Да брат, так и есть.

[xvr 12]

Пункты 1 и 2:

template<class Type> struct TypeTraits {};
template<> struct TypeTraits<float> {using IntAType = int;};
template<> struct TypeTraits<double> {using IntAType = long long;};
  
template<typename T> void CalcSinCos(T a, T &s, T &c) // a должно быть больше 0 и меньше 1
{ const int BL = ((sizeof(T)==4))?24:53;
  typename TypeTraits<T>::IntAType IntA = 0;
 ...
 }

  

 

[krux 8]

кстати еще момент. CORDIC разрабатывался изначально как аналогоое устройство с известным временем задержки от входа до выхода для одновременного получения результатов sin и  cos. На современных CPU время от входа до выхода может быть произвольным, особенно если там операционка. Если вам нужен реалтайм, возможно потребуется Tin Tout в качестве дополнительных параметров CORDIC-а и потом фильтрация, хоть по среднему, хоть по интерполяции, эксраполяции, а может и по Калману.

[iiv 16]

1 hour ago, xvr said:

  typename TypeTraits<T>::IntAType IntA = 0;

Супер, спасибо большое, xvr!

1 hour ago, krux said:

кстати еще момент. CORDIC разрабатывался изначально как аналогоое устройство

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

[Edit2007 3]

Цитата

Меня, кстати очень удивляет, почему этот алгоритм называют кордиком

WIKI: CORDIC (Метод CORDIC от англ. COordinate Rotation DIgital Computer — цифровой вычислитель поворота системы координат; метод «цифра за цифрой», алгоритм Волдера) — итерационный метод сведения прямых вычислений сложных функций к выполнению простых операций сложения и сдвига.

[iiv 16]

Продолжаю тормозить, теперь совсем не могу понять почему я не могу проинициализировать константным выражением массив из 7 элементов

 

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <array>

//                             0    1    2    -    4    5 остальные более дальние точки не будем принимать во внимание
const double WCoefficients[]={1., 0.8, 0.7, 0.0, 0.5, 0.4};


template<typename T> constexpr void InvSpec3Mat(T *x, T *z)
{ T q0=(x[0]+x[1])*x[3]-2.*x[2]*x[2];
  if(q0) q0=1./q0;
 
  T q=x[3]*q0;
 
  T p=x[0]-x[1];
  if(p) p=1./p;
 
  z[0] = (q+p)*0.5;
  z[1] = (q-p)*0.5;
 
  z[2] = -x[2]*q0;
  z[3] = (x[3])?(1-2.*x[2]*z[2])/x[3]:0.;
}


template<typename T>constexpr auto MatCoefficientInit()
{ T A[6][6]={{0,0,0,0,0,0},{0,0,0,0,0,0},{0,0,0,0,0,0},{0,0,0,0,0,0},{0,0,0,0,0,0},{0,0,0,0,0,0}};
  std::array<T, 7> RES{};

  for(int i=-2; i<=2; i++)
    for(int j=-2; j<=2; j++)
    { int q=i*i+j*j;
      if(q<6)
      { T v[6];
        v[0]=(T)(i*i);
        v[1]=(T)(j*j);
        v[2]=1.;
        v[3]=(T)(i*j);
        v[4]=(T)i;
        v[5]=(T)j;
        for(int k1=0; k1<6; k1++)
          for(int k2=0; k2<6; k2++)
            A[k1][k2]+=v[k1]*v[k2]*WCoefficients[q]*WCoefficients[q];
      }
    }
  for(int i=3; i<6; i++)
    RES[i+1]=A[i][i];
  T X[4];
  X[0]=A[0][0];
  X[1]=A[0][1];
  X[2]=A[0][2];
  X[3]=A[2][2];
  InvSpec3Mat(X, &(RES[0]));
  return RES;
}


constexpr auto MatCoefficients=MatCoefficientInit<double>(); // здесь хочется привязать этот тип к типу WCoefficients

То есть у меня есть WCoefficients, я этим кодом хочу вычислить массив MatCoefficients из 7ми элементов того же типа, причем так, чтобы все эти темплейты были посчитаны во время компиляции и не занимали бы вообще ни одной ассемблерной инструкции во время работы. Понятно, что можно написать отдельную программу, и результат вставлять в код, но ведь в С++ наделали столько, что это все должно быть автоматически!

 

Ругается так:

Quote

gg.cpp:59:59: error: ‘constexpr auto MatCoefficientInit() [with T = double]’ called in a constant expression
   59 | constexpr auto MatCoefficients=MatCoefficientInit<double>();
      |                                                           ^
gg.cpp:27:36: note: ‘constexpr auto MatCoefficientInit() [with T = double]’ is not usable as a ‘constexpr’ function because:
   27 | template<typename T>constexpr auto MatCoefficientInit()
 

помогите, пожалуйста, понять что я делаю не так!

 

Спасибо!

Ой... Оказывается я g++-9 использовал, на g++-10 все работает. Вопрос закрыт!