analitza/analitza/polynomial.cpp

0001 /*************************************************************************************
0002  *  Copyright (C) 2011 by Aleix Pol <aleixpol@kde.org>                               *
0003  *                                                                                   *
0004  *  This program is free software; you can redistribute it and/or                    *
0005  *  modify it under the terms of the GNU General Public License                      *
0006  *  as published by the Free Software Foundation; either version 2                   *
0007  *  of the License, or (at your option) any later version.                           *
0008  *                                                                                   *
0009  *  This program is distributed in the hope that it will be useful,                  *
0010  *  but WITHOUT ANY WARRANTY; without even the implied warranty of                   *
0011  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                    *
0012  *  GNU General Public License for more details.                                     *
0013  *                                                                                   *
0014  *  You should have received a copy of the GNU General Public License                *
0015  *  along with this program; if not, write to the Free Software                      *
0016  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA   *
0017  *************************************************************************************/
0018
0019 #include "polynomial.h"
0020 #include "apply.h"
0021 #include "value.h"
0022 #include "analitzautils.h"
0023 #include "operations.h"
0024 #include "analyzer.h"
0025
0026 using namespace Analitza;
0027
0028 // static QDebug operator<<(QDebug dbg, const Object* c)
0029 // {
0030 //     dbg.nospace() << (c ? c->toString() : "<null>");
0031 //
0032 //     return dbg.space();
0033 // }
0034 //
0035 // static QDebug operator<<(QDebug dbg, const Monomial &c)
0036 // {
0037 //     dbg.nospace() << "(" << c.first << ", " << c.second << ")";
0038 //
0039 //     return dbg.space();
0040 // }
0041
0042 static Object* negateObject(Object* o)
0043 {
0044     if(o->type()==Operator::value) {
0045         Cn* v = static_cast<Cn*>(o);
0046         v->rvalue() *= -1;
0047         return v;
0048     } else {
0049         Apply* a = new Apply;
0050         a->appendBranch(new Operator(Operator::minus));
0051         a->appendBranch(o);
0052         return a;
0053     }
0054 }
0055
0056 static Object* simpExpression(Object* o)
0057 {
0058     Analyzer a;
0059     a.setExpression(Expression(o));
0060     a.simplify();
0061     Expression ret = a.expression();
0062     return ret.takeTree();
0063 }
0064
0065 bool Monomial::isScalar(const Object* o)
0066 {
0067     return o->type()==Object::value || (o->type()==Object::vector && !AnalitzaUtils::hasVars(o));
0068 }
0069
0070 Object* Monomial::createMono(const Operator& o) const
0071 {
0072     Operator::OperatorType mult = o.multiplicityOperator();
0073
0074     Object* toAdd=nullptr;
0075     if(first==0.)
0076         delete second;
0077     else if(first==1.)
0078         toAdd=second;
0079     else if(first==-1. && mult==Operator::times)
0080         toAdd=negateObject(second);
0081     else if(mult==Operator::times && second->isApply() && static_cast<Apply*>(second)->firstOperator()==Operator::times) {
0082         Apply* res = static_cast<Apply*>(second);
0083         res->prependBranch(new Cn(first));
0084         toAdd=res;
0085     } else {
0086         Apply *cint = new Apply;
0087         cint->appendBranch(new Operator(mult));
0088         if(mult==Operator::times) {
0089             cint->appendBranch(new Cn(first));
0090             cint->appendBranch(second);
0091         } else {
0092             cint->appendBranch(second);
0093             cint->appendBranch(new Cn(first));
0094         }
0095         toAdd=cint;
0096     }
0097     return toAdd;
0098 }
0099
0100 Monomial::Monomial(const Operator& o, Object* o2, bool& sign)
0101 {
0102     bool ismono=false;
0103     Operator::OperatorType mult = o.multiplicityOperator();
0104
0105     if(o2->isApply()) {
0106         Apply *cx = (Apply*) o2;
0107         if(cx->firstOperator()==mult) {
0108             if(cx->countValues()==2) {
0109                 bool valid=false;
0110                 int scalar, var;
0111
0112                 if(mult!=Operator::power && cx->m_params[0]->type()==Object::value) {
0113                     scalar=0;
0114                     var=1;
0115                     valid=true;
0116                 } else if(cx->m_params[1]->type()==Object::value) {
0117                     scalar=1;
0118                     var=0;
0119                     valid=true;
0120                 }
0121
0122                 if(valid) {
0123 //                     qDebug() << ";;;;;;;" << o2->toString();
0124                     Cn* sc= (Cn*) cx->m_params[scalar];
0125                     first = sc->value();
0126                     second = cx->m_params[var];
0127
0128                     cx->m_params[var]=nullptr;
0129                     delete cx;
0130
0131                     ismono=true;
0132 //                     qDebug() << ":::::::" << *this;
0133                 }
0134             } else if(mult==Operator::times) {
0135                 first=1;
0136                 Apply::iterator it=cx->firstValue(), itEnd=cx->end();
0137                 QVector<Object*> vars;
0138                 QVector<Object*> values;
0139
0140                 for(; it!=itEnd; ++it) {
0141                     if((*it)->type()==Object::value) {
0142                         first *= static_cast<Cn*>(*it)->value();
0143                         values += *it;
0144                         ismono=true;
0145                     } else {
0146                         vars += *it;
0147                     }
0148                 }
0149
0150                 if(ismono) {
0151                     cx->m_params = vars;
0152                     second = cx;
0153                     qDeleteAll(values);
0154                 }
0155             }
0156         } else if(cx->firstOperator()==Operator::minus && cx->isUnary()) {
0157             *this = Monomial(o, *cx->firstValue(), sign);
0158             *cx->firstValue()=nullptr;
0159             delete cx;
0160             ismono=true;
0161
0162             if(o==Operator::times)
0163                 sign = !sign;
0164             else if(o==Operator::plus || o==Operator::minus)
0165                 first *= -1;
0166         }
0167     } else if(o2->type()==Object::value && (o==Operator::plus || o==Operator::minus)) {
0168         Cn* v=static_cast<Cn*>(o2);
0169         if(v->value()<0) {
0170             v->rvalue()*=-1;
0171             first = -1;
0172             second = o2;
0173             ismono = true;
0174         }
0175     }
0176
0177     if(!ismono) {
0178         first = 1.;
0179         second = o2;
0180     }
0181 }
0182
0183 bool Monomial::isValue() const
0184 {
0185     return isScalar(second);
0186 }
0187
0188 Polynomial::Polynomial(Apply* c)
0189     : m_operator(c->firstOperator())
0190     , m_sign(true)
0191 {
0192     bool first = true, firstValue = false;
0193     QList<Monomial> monos;
0194 //     qDebug() << "1.........." << c->toString();
0195     for(Apply::const_iterator it=c->constBegin(), itEnd=c->constEnd(); it!=itEnd; ++it, first=false) {
0196         Monomial imono(m_operator, *it, m_sign);
0197
0198         bool added=false;
0199         if(imono.second->isApply() && imono.first == 1) {
0200             Apply* a = static_cast<Apply*>(imono.second);
0201             Operator op=a->firstOperator();
0202             if(a->firstOperator()==m_operator
0203                 || (m_operator==Operator::minus && op==Operator::plus)
0204                 || (m_operator==Operator::plus && op==Operator::minus)
0205             ) {
0206                 Polynomial p(a);
0207
0208                 a->m_params.clear();
0209                 delete a;
0210
0211                 if((!first && m_operator==Operator::plus && op==Operator::minus)
0212                     || (first && m_operator==Operator::minus && op==Operator::plus)
0213                     || (first && m_operator==Operator::minus && op==Operator::minus)
0214                 )
0215                 {
0216                     //if it's a minus, the first is already positive
0217                     p.negate(1);
0218                 }
0219
0220                 monos.append(p);
0221                 added=true;
0222             }
0223         }
0224
0225         if(!added)
0226             monos.append(imono);
0227     }
0228
0229     first = true;
0230     for(auto it=monos.begin(), itEnd=monos.end(); it!=itEnd; ++it, first=false) {
0231         if(m_operator==Operator::minus && !first)
0232             it->first*=-1;
0233     }
0234
0235     for(auto it=monos.begin(), itEnd=monos.end(); it!=itEnd; ++it, first=false)
0236         addMonomial(*it);
0237
0238     simpScalars(firstValue);
0239 }
0240
0241 void Polynomial::addMonomial(const Monomial& m)
0242 {
0243     if(m.isValue()) {
0244         m_scalars += m.createMono(m_operator);
0245         return;
0246     }
0247
0248     bool found = false;
0249     auto it1(begin());
0250     for(; it1!=end(); ++it1) {
0251         Object *o1=it1->second, *o2=m.second;
0252         found = AnalitzaUtils::equalTree(o1, o2);
0253         if(found)
0254             break;
0255     }
0256
0257     if(found) {
0258         it1->first += m.first;
0259         delete m.second;
0260
0261         if(it1->first==0.) {
0262             delete it1->second;
0263             erase(it1);
0264         }
0265     } else {
0266         append(m);
0267     }
0268 }
0269
0270 void Polynomial::simpScalars(bool firstValue)
0271 {
0272     Object *value=nullptr;
0273     if(!firstValue && m_operator==Operator::minus && !m_scalars.isEmpty())
0274         m_scalars.first() = negateObject(m_scalars.first());
0275
0276     for(auto i=m_scalars.constBegin(); i!=m_scalars.constEnd(); ++i) {
0277         bool d=false;
0278
0279         Object* aux = simpExpression(*i);
0280         if(value) {
0281             QString* err=nullptr;
0282             value=Operations::reduce(m_operator.operatorType(), value, aux, &err);
0283             d=err;
0284             delete err;
0285         } else
0286             value=aux;
0287
0288         if(d) {
0289             addValue(aux);
0290             value=nullptr;
0291         }
0292     }
0293
0294     addValue(value);
0295     m_scalars.clear();
0296 }
0297
0298 void Polynomial::addValue(Object* value)
0299 {
0300     bool sign=false;
0301     if(!value)
0302         return;
0303
0304     if(value->isZero())
0305         delete value;
0306     else {
0307         Monomial imono(m_operator, value, sign);
0308
0309         if(m_operator==Operator::plus)
0310             append(imono);
0311         else if(m_operator==Operator::minus) {
0312             imono.first *= -1;
0313             append(imono);
0314         } else
0315             prepend(imono);
0316     }
0317 }
0318
0319 Object* Polynomial::toObject()
0320 {
0321     Object* root = nullptr;
0322     if(count()==1) {
0323         root = first().createMono(m_operator);
0324     } else if(count()>1) {
0325         Apply* c = new Apply;
0326         c->appendBranch(new Operator(m_operator));
0327
0328         auto i=begin();
0329         bool first=true;
0330         for(; i!=end(); ++i, first=false) {
0331             if(!first && m_operator==Operator::minus)
0332                 i->first *= -1;
0333
0334             Object* toAdd=i->createMono(m_operator);
0335             if(toAdd)
0336                 c->appendBranch(toAdd);
0337         }
0338         root=c;
0339     }
0340
0341     if(!root) {
0342         root=new Cn(0.);
0343     } else if(!m_sign) {
0344         root = negateObject(root);
0345     }
0346
0347     return root;
0348 }
0349
0350 void Polynomial::negate(int i)
0351 {
0352     for(iterator it=begin(); it!=end(); ++it, --i) {
0353         if(i<=0)
0354             it->first *= -1;
0355     }
0356 }