types.h Определяет составные типы данных, используемые в данной задаче, и объявляет прототипы функций.
#if !defined(types_h)
#define types_h
#include <vector>
#include <list>
#include <string>
using namespace std;
//typedef struct
class j_k_powers
{
public:
double val; // factor;
int jp; // power of J^2;
int kp; // power of J_z;
bool operator==(const j_k_powers& ) const;
};
typedef vector<j_k_powers> Term; /* regular expression for Tlm */
class Tlm {
public:
Term clm;
list<int> a; /* Простые множители числителя под корнем */
list<int> b; /* Простые множители знаменателя под корнем */
/* list используется, чтобы было легче упрощать a/b */
};
typedef list<int> lstarr;
void PrimeFactors(int ,vector<int>& );
Term Terms_multiplication(Term ,Term );
void EqualAdds(Term& );
bool compj_k_powers (j_k_powers , j_k_powers );
void division(list<int>& , list<int>& );
bool comp_int (int , int );
list<int> lists_devide(list<int>& , list<int>& );
list<int> lists_mult(list<int> , vector<int> );
void print_vect(const string,list<int> );
void clear(Tlm& );
string IntToString (int );
int SqrtInt(list<int> );
int fullnumber(list<int>);
void list_erase(void);
#endif //types_h
main.cpp Основная и вспомогательные фйнкции.
// main.cpp : Defines the entry point for the console application.
// Program to evaluate expressions for irreducible tensor operators.
// Next version after tlm.cpp. Using rearanged formulas.
#include "types.h"
#include <fstream>
#include <algorithm>
#include <string>
#include <sstream>
#define wmn 1000 // to sort vector
using namespace std;
ofstream ff;
int main(int argc, char* argv[])
{
ff.open("out.txt",ios::out);
int M = 4; /* Значение M цепочки TLM */
int QuantityOfTlm = 12;
ff << "Рассчитываются TLM для M=" << M << endl;
vector<Tlm> TlmArray; /* Массив для всех найденных Tlm. Первую позицию занимает TMM */
TlmArray.reserve(QuantityOfTlm);
string s1 = IntToString(M);
string s2;
int Lmax = M + QuantityOfTlm;
// Рабочие переменные
Tlm tmw;
j_k_powers w;
vector<int> f; /* prime factors */
f.resize(40);
// == T^(M)_(M) =========================
// T^(M)_(M) has just one term. T^(M)_(M) = 1*(T^(M)_(M)).
ff << endl << "=========== " << "T^"+s1+"_"+s1 << " ===========" << endl;
clear(tmw);
w.val = 1;
w.jp = 0;
w.kp = 0;
tmw.clm.push_back(w);
tmw.a.push_back(1);
tmw.b.push_back(1);
TlmArray.push_back(tmw);
ff << "val " << "jp " << "kp" << endl;
ff << w.val << " " << w.jp << " " << w.kp << endl;
// --------------------------------------
// == T^(M+1)_(M) =======================
s2 = IntToString(M+1);
ff << endl << "=========== " << "T^"+s2+"_"+s1 << " ===========" << endl;
clear(tmw);
PrimeFactors(2*M+1,f); // Раскладывает на простые множители.
// Вставляем обязательно 1, для того, чтобы при сокращении множителей
// не получить нулевого результата.
// Числитель
tmw.a.push_back(1);
for(int i=0; i < f.size(); i++)
{
tmw.a.push_back(f[i]);
}
// Знаменатель
tmw.b.push_back(1);
tmw.b.push_back(2);
tmw.b.push_back(2);
ff << "val " << "jp " << "kp" << endl;
w.val = M;
w.jp = 0;
w.kp = 0;
if(M != 0){
tmw.clm.push_back(w);
ff << w.val << " " << w.jp << " " << w.kp << endl;
}
w.val = 2;
w.jp = 0;
w.kp = 1;
tmw.clm.push_back(w);
ff << w.val << " " << w.jp << " " << w.kp << endl;
TlmArray.push_back(tmw);
// ------------------------------------------
// == all Tlms с использованием рекурентной формулы ==
// Соответсвия: i=0 -> T^(M)_(M), i=1 -> T^(M+1)_(M), i=2 -> T^(M+2)_(M),....
Tlm tmw1, tmw2; /* две составляющие нового Tlm */
Term r, g1, g2;
for(i = 2; i < QuantityOfTlm; i++)
{
clear(tmw1);
clear(tmw2);
r.clear();
int L=M+i; /* очередное L, для которого строим выражение */
s2 = IntToString(L);
ff << endl << "=========== " << "T^"+s2+"_"+s1 << " ===========" << endl;
int s = i-1; /* номер последнего вычисленного Tlm */
// Включаем в выражение член с T^{L-1}_{M} -------------------
w.val = M;
w.jp = 0;
w.kp = 0;
r.push_back(w);
w.val = 2;
w.jp = 0;
w.kp = 1;
r.push_back(w);
/*
* Умножение P_{L-1}*(2J_z+M) и приведение подобных.
*/
/* g1 содержит полином от степеней J^2 и J_z в первой части формулы,
то есть той, которая содержит T^{L-1}_{M} */
g1 = Terms_multiplication(TlmArray[s].clm,r);
EqualAdds(g1);
/* tmw1.a, tmw1.b - числитель и знаменатель подкоренного выражения в первой части формулы,
то есть той, которая содержит T^{L-1}_{M} */
// Знаменатель добавляемого подкоренного выражения
PrimeFactors((L+M)*(L-M)*16,f);
tmw1.b = lists_mult(TlmArray[s].b, f);
// Числитель добавляемого подкоренного выражения
PrimeFactors((2*L-1)*(2*L-1)*4,f);
tmw1.a = lists_mult(TlmArray[s].a, f);
int j1;
// Член в выражении с T^{L-2}_{M} -------------------
r.clear();
w.val = L*(L-2);
w.jp = 0;
w.kp = 0;
r.push_back(w);
w.val = -4;
w.jp = 1;
w.kp = 0;
r.push_back(w);
/* g2 содержит полином от степеней J^2 и J_z во второй части формулы,
то есть той, которая содержит T^{L-2}_{M} */
g2 = Terms_multiplication(TlmArray[s-1].clm,r);
EqualAdds(g2);
/* tmw2.a, tmw2.b - числитель и знаменатель подкоренного выражения во второй части формулы,
то есть той, которая содержит T^{L-2}_{M} */
// Знаменатель добавляемого подкоренного выражения
PrimeFactors((L+M)*(L-M)*16,f);
tmw2.b = lists_mult(TlmArray[s-1].b, f);
PrimeFactors((L-1+M)*(L-1-M),f);
tmw2.a = lists_mult(TlmArray[s-1].a, f);
list<int> bl = lists_devide(tmw1.b, tmw2.b);
/* bl - общий множитель подкоренного выражения.
В tmw1.b, tmw2.b остаются индивидуальные множители */
/* Анализируем и приводим к общему знаменателю, используя bl, tmw1.b и tmw2.b*/
for(list<int>::iterator k = tmw1.b.begin(); k != tmw1.b.end(); k++)
bl.push_back(*k);
for(k = tmw2.b.begin(); k != tmw2.b.end(); k++)
if(*k != 1)
bl.push_back(*k);
bl.sort();
/* bl содержит общий знаменатель */
/* Числитель подкоренного выражения при T^{L-1}_{M}*/
for(k = tmw2.b.begin(); k != tmw2.b.end(); k++)
if(*k != 1)
tmw1.a.push_back(*k);
tmw1.a.sort();
/* Числитель подкоренного выражения при T^{L-2}_{M}*/
for(k = tmw1.b.begin(); k != tmw1.b.end(); k++)
if(*k != 1)
tmw2.a.push_back(*k);
tmw2.a.sort();
list<int> a1=tmw1.a, a2=tmw2.a;
// bc - общий множитель числителей под корнем, который нужно вынести за скобки
// и попытаться сократить.
list<int> bc = lists_devide(a1, a2);
int sqa1 = SqrtInt(a1);
int sqa2 = SqrtInt(a2);
if(sqa1*sqa1 != fullnumber(a1))
{
ff << "Какие-то проблемы. 1" << endl;
exit(1);
}
if(sqa2*sqa2 != fullnumber(a2))
{
ff << "Какие-то проблемы. 2" << endl;
exit(2);
}
//Корректировка левого полинома
for(int j=0; j < g1.size(); j++)
g1[j].val*=sqa1;
//Корректировка правого полинома
// ff << "g2s " << g2.size() << endl;
for(j=0; j < g2.size(); j++)
g2[j].val*=sqa2;
// Сложение левого и правого полиномов с результатом в g1.
for(j=0; j < g2.size(); j++)
g1.push_back(g2[j]);
// Сбор подобных членов в полиноме.
EqualAdds(g1);
// Определение общего множителя в полиноме
vector<lstarr> mutual;
mutual.clear();
list<int> wlist;
for(j=0; j < g1.size(); j++)
{
f.clear();
wlist.clear();
PrimeFactors(g1[j].val,f);
for(vector<int>::iterator k = f.begin(); k != f.end(); k++)
wlist.push_back(*k);
mutual.push_back(wlist);
}
wlist=mutual[0];
for(j=1; j < mutual.size(); j++)
wlist=lists_devide(wlist, mutual[j]);
int icf = fullnumber(wlist);
tmw.clm = g1;
ff << "Результирующий полином" << endl;
for(j=0; j < tmw.clm.size(); j++)
{
tmw.clm[j].val/=icf;
ff << tmw.clm[j].val << " " << tmw.clm[j].jp << " " << tmw.clm[j].kp << endl;
}
tmw.a = bc;
tmw.b = bl;
// Переносим общий множитель полинома в числитель множителя всего выражения под корень
for(list<int>::iterator j2=wlist.begin(); j2 != wlist.end(); j2++)
{
if(*j2 != 1)
{
tmw.a.push_back(*j2);
tmw.a.push_back(*j2);
}
}
lists_devide(tmw.a, tmw.b);
print_vect("Числитель подкоренного выражение",tmw.a);
print_vect("Знаменатель подкоренного выражение",tmw.b);
TlmArray.push_back(tmw);
}
return 0;
}
//======================================================================
/* ===== Разложение числа на простые множители ===== */
//----------------------------------------------------------------------
void PrimeFactors(int s,vector<int>& f)
{
f.clear();
int t, i=2;
// Число, раскладываемое на простые множители должно быть положительным.
t = abs(s);
while(i<=t)
{
if(t%i==0)
{
f.push_back(i);
t=t/i;
}
else
i=i+1;
}
}
//======================================================================
// Умножение двух полиномов от степеней J^2 и J_z
//----------------------------------------------------------------------
Term Terms_multiplication(Term A,Term B)
{
Term C;
j_k_powers w;
C.clear();
for(int k=0; k<A.size(); k++)
for(int i=0; i<B.size(); i++)
{
w.jp=A[k].jp+B[i].jp;
w.kp=A[k].kp+B[i].kp;
w.val=A[k].val*B[i].val;
C.push_back(w);
}
return C;
}
//======================================================================
// Сбор подобных членов
//----------------------------------------------------------------------
void EqualAdds(Term& A)
{
Term B=A,C;
j_k_powers w;
std::sort(B.begin(),B.end(),compj_k_powers);
//print_Term(B,"Sorted path");
C.clear();
int i=0, k=0;
while(i<B.size())
{
w=B[i];
//ff << i << endl;
C.push_back(w);
i++;
while((i<B.size()) && (w==B[i]))
{
C[k].val+=B[i].val;
//ff << i << " " << k << endl;
i++;
}
k++;
}
// print_Term(C,"Ultimate one path");
// Deleting terms with zero coefficients/
k=0;
B.clear();
for(i=0; i<C.size(); i++)
if(C[i].val != 0)
B.push_back(C[i]);
A=B;
}
//======================================================================
// Логическая функция сравнения для сортировки вектора Term.
//----------------------------------------------------------------------
bool compj_k_powers (const j_k_powers A, const j_k_powers B)
{
int a = wmn*A.jp+A.kp;
int b = wmn*B.jp+B.kp;
return (a > b);
// return (a >= b);
// return (wmn*A.jp+A.kp >= wmn*B.jp+B.kp);
}
//======================================================================
// Оператор равенства степеней.
//----------------------------------------------------------------------
bool j_k_powers::operator==(const j_k_powers& A) const
{
return(this->jp == A.jp && this->kp == A.kp);
}
//======================================================================
// Деление двух целых чисел представленных в виде произведения простых чисел.
//----------------------------------------------------------------------
void division(std::list<int>& a, std::list<int>& b)
{
// sort(a.begin(),a.end(),comp_int);
// sort(b.begin(),b.end(),comp_int);
a.sort();
b.sort();
}
//======================================================================
// Логическая функция сравнения для сортировки вектора std::vector<int>.
//----------------------------------------------------------------------
bool comp_int (int A, int B)
{
return (A >= B);
}
//======================================================================
// Исключение общих простых множителей из двух списков простых множителей.
// a - числитель.
// b - знаменатель.
// Общие множители - возвращаемое значение.
//------------------------------------------------------------------
list<int> lists_devide(list<int>& a, list<int>& b)
{
list<int> o;
for( list<int>::iterator k = b.begin(); k != b.end();)
{
int c = 0;
for( list<int>::iterator i = a.begin(); i != a.end();)
{
if(*k == *i && *k != 1)
{
i=a.erase(i);
c = 1;
break;
}
else
i++;
}
if(c == 1)
{
o.push_back(*k);
k=b.erase(k);
/*
ff<< "a" << endl;
for( i = a.begin(); i != a.end();i++)
ff<<*i<<" ";
ff<<endl;
ff<< "b" << endl;
for( i = b.begin(); i != b.end();i++)
ff<<*i<<" ";
ff<<endl;
*/
}
else
k++;
}
return o;
}
//======================================================================
// Умножение (объединение) двух строк простых множителей.
//------------------------------------------------------------------
list<int> lists_mult(list<int> a, vector<int> b)
{
list<int> c;
c.clear();
c = a;
for(vector<int>::iterator i = b.begin(); i != b.end(); i++)
c.push_back(*i);
c.sort();
return c;
}
//======================================================================
// Перемножение всех членов list<int>.
//------------------------------------------------------------------
int fullnumber(list<int> a)
{
int m=1;
for( list<int>::iterator k = a.begin(); k != a.end();k++)
m =m* (*k);
return m;
}
/* ======================================================================
* Извлечение корня из целого числа представленного в виде произведения
* простых множителей и являющегося квадратом целого числа.
* В функции выполняется, что число действительно есть квадрат целого числаю
* Проверка осуществляется путем определения количества простых множителей.
------------------------------------------------------------------ */
int SqrtInt(list<int> f)
{
list<int> fin = f;
fin.sort();
int c=0, o=1;
for(list<int>::iterator i=fin.begin(); i != fin.end();)
{
int m;
int k; // количество одинаковых
m=*i; k=0;
while(m == *i && i != fin.end())
{
k++; i++;
}
if(k%2 != 0 && m != 1)
{
ff << "Корень извлечь нельзя" << endl;
return -1;
}
k=k/2;
for(int l=0; l<k; l++)
o*=m;
}
return o;
}
//======================================================================
// Очистка полей Tlm.
//------------------------------------------------------------------
void clear(Tlm& tmw)
{
tmw.a.clear();
tmw.b.clear();
tmw.clm.clear();
}
//======================================================================
// Преобразование int to string.
//------------------------------------------------------------------
string IntToString (int a)
{
string str;
ostringstream temp;
temp<<a;
return temp.str();
}
//======================================================================
// Вспомогательный код для отработки разных алгоритмов.
//------------------------------------------------------------------
void list_erase(void)
{
list<int> x;
for ( int i =0;i<10; i++)
x.push_back(i);
for( list<int>::iterator k = x.begin(); k != x.end();k++)
ff<<*k<<" ";
ff<<endl;
for( k = x.begin(); k != x.end();)
if( !((*k)%2) )
k=x.erase(k);
else
++k;
for( k = x.begin(); k != x.end();k++)
ff<<*k<<" ";
ff<<endl;
}
//======================================================================
// Вывод в поток ff компонент list<int> f.
// h = 1,2 - два разных формата вывода.
//------------------------------------------------------------------
void print_vect(const string inf,list<int> f)
{
int h=1;
// ff << inf << endl;
ff << inf+" = " << fullnumber(f) << " -> ";
if(h==1)
{
for(list<int>::iterator i=f.begin(); i != f.end(); i++)
{
ff << *i << " ";
}
ff << endl;
}
h=2;
if(h==2)
{
// Извлечение корня, а затем вывод в поток
list<int> o, fin = f;
fin.sort();
int c=0;
for(list<int>::iterator i=fin.begin(); i != fin.end();)
{
int m=*i;
int k=0; // количество одинаковых
while(m == *i && i != fin.end())
{
k++; i++;
}
ff << m << ", ко-во =" << k << endl;
}
}
}