From 01b1847bdd455bf3ad075b27b225f585cdaa5f4c Mon Sep 17 00:00:00 2001 From: rdbyk Date: Mon, 26 Feb 2018 12:05:03 +0100 Subject: [PATCH] refactored complex division; use C99/C++11 complex arithmetic --- .../operations/matrix_right_division.h | 40 +-- .../ast/src/c/operations/matrix_division.c | 243 ++++++------------ .../ast/src/cpp/operations/types_divide.cpp | 50 ++-- .../ast/tests/unit_tests/dotdivide.dia.ref | 8 +- .../ast/tests/unit_tests/dotdivide.tst | 9 +- .../modules/elementary_functions/src/c/expm.c | 13 +- 6 files changed, 113 insertions(+), 250 deletions(-) diff --git a/scilab/modules/ast/includes/operations/matrix_right_division.h b/scilab/modules/ast/includes/operations/matrix_right_division.h index 7c20d911029..5678cad3867 100644 --- a/scilab/modules/ast/includes/operations/matrix_right_division.h +++ b/scilab/modules/ast/includes/operations/matrix_right_division.h @@ -2,7 +2,7 @@ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab * Copyright (C) 2008-2008 - DIGITEO - Antoine ELIAS * Copyright (C) 2012 - 2016 - Scilab Enterprises - * Copyright (C) 2017 - Dirk Reusch, Kybernetik Dr. Reusch + * Copyright (C) 2017 - 2018 Dirk Reusch, Kybernetik Dr. Reusch * * This file is hereby licensed under the terms of the GNU GPL v2.0, * pursuant to article 5.3.4 of the CeCILL v.2.1. @@ -26,39 +26,9 @@ int iRightDivisionOfComplexMatrix( double *_pdblReal2, double *_pdblImg2, int _iRows2, int _iCols2, double *_pdblRealOut, double *_pdblImgOut, int _iRowsOut, int _iColsOut, double *_pdblRcond); -int iRightDivisionRealMatrixByRealMatrix( - double *_pdblReal1, int _iInc1, - double *_pdblReal2, int _iInc2, - double *_pdblRealOut, int _iIncOut, int _iSize); - -int iRightDivisionComplexByReal( - double _dblReal1, double _dblImg1, - double _dblReal2, - double *_pdblRealOut, double *_pdblImgOut); - -int iRightDivisionComplexMatrixByRealMatrix( - double *_pdblReal1, double *_pdblImg1, int _iInc1, - double *_pdblReal2, int _iInc2, - double *_pdblRealOut, double *_pdblImgOut, int _iIncOut, int _iSize); - -int iRightDivisionRealByComplex( - double _dblReal1, - double _dblReal2, double _dblImg2, - double *_pdblRealOut, double *_pdblImgOut); - -int iRightDivisionRealMatrixByComplexMatrix( - double *_pdblReal1, int _iInc1, - double *_pdblReal2, double *_pdblImg2, int _iInc2, - double *_pdblRealOut, double *_pdblImgOut, int _iIncOut, int _iSize); - -int iRightDivisionComplexByComplex( - double _dblReal1, double _dblImg1, - double _dblReal2, double _dblImg2, - double *_pdblRealOut, double *_pdblImgOut); - -int iRightDivisionComplexMatrixByComplexMatrix( - double *_pdblReal1, double *_pdblImg1, int _iInc1, - double *_pdblReal2, double *_pdblImg2, int _iInc2, - double *_pdblRealOut, double *_pdblImgOut, int _iIncOut, int _iSize); +int iRightDivisionRealMatrixByReal(double* A, double b, double* X, int n); +int iRightDivisionComplexMatrixByReal(double* A, double* B, double c, double* X, double* Y, int n); +int iRightDivisionRealMatrixByComplex(double* A, double c, double d, double* X, double* Y, int n); +int iRightDivisionComplexMatrixByComplex(double* A, double* B, double c, double d, double* X, double* Y, int n); #endif /* __MATRIX_RDIV__ */ diff --git a/scilab/modules/ast/src/c/operations/matrix_division.c b/scilab/modules/ast/src/c/operations/matrix_division.c index 43a303b07ec..8ea07f50bf9 100644 --- a/scilab/modules/ast/src/c/operations/matrix_division.c +++ b/scilab/modules/ast/src/c/operations/matrix_division.c @@ -2,7 +2,7 @@ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab * Copyright (C) 2008-2008 - INRIA - Antoine ELIAS * Copyright (C) 2012 - 2016 - Scilab Enterprises - * Copyright (C) 2017 - Dirk Reusch, Kybernetik Dr. Reusch + * Copyright (C) 2017 - 2018 Dirk Reusch, Kybernetik Dr. Reusch * * This file is hereby licensed under the terms of the GNU GPL v2.0, * pursuant to article 5.3.4 of the CeCILL v.2.1. @@ -16,234 +16,135 @@ #include "configvariable_interface.h" #include "matrix_division.h" #include +#include -int iRightDivisionComplexMatrixByComplexMatrix( - double *_pdblReal1, double *_pdblImg1, int _iInc1, - double *_pdblReal2, double *_pdblImg2, int _iInc2, - double *_pdblRealOut, double *_pdblImgOut, int _iIncOut, int _iSize) +int iRightDivisionComplexMatrixByComplex(double* A, double* B, double c, double d, + double* X, double* Y, int n) { - int iErr = 0; - int iIndex = 0; //Main loop index - int iIndex1 = 0; //Loop index on left operand - int iIndex2 = 0; //Loop index on right operand - int iIndexOut = 0; //Lopp index on result matrix + // X + Y*i = (A + B*i) / (c + d*i) - if (_iInc2 == 0) - { - if ((getieee() == 0) && (fabs(_pdblReal2[iIndex2]) + fabs(_pdblImg2[iIndex2]) == 0)) - { - return 3; - } - } + int iErr = 0; - for (iIndex = 0 ; iIndex < _iSize ; iIndex++) + if (c == 0 && d == 0) { - iErr = iRightDivisionComplexByComplex(_pdblReal1[iIndex1], _pdblImg1[iIndex1], _pdblReal2[iIndex2], _pdblImg2[iIndex2], &_pdblRealOut[iIndexOut], &_pdblImgOut[iIndexOut]); - iIndexOut += _iIncOut; - iIndex1 += _iInc1; - iIndex2 += _iInc2; - } - - return iErr; -} + int mode = getieee(); -int iRightDivisionComplexByComplex( - double _dblReal1, double _dblImg1, - double _dblReal2, double _dblImg2, - double *_pdblRealOut, double *_pdblImgOut) -{ - int iErr = 0; - if (_dblImg2 == 0) - { - if (_dblReal2 == 0) + if (mode == 0) { - //got NaN + i NaN - iErr = 10; - *_pdblRealOut = _dblImg2 / _dblReal2; - *_pdblImgOut = *_pdblRealOut; + iErr = 3; // error + return iErr; } - else + else if (mode == 1) { - double dblReal2Inv = 1.0 / _dblReal2; - *_pdblRealOut = _dblReal1 * dblReal2Inv; - *_pdblImgOut = _dblImg1 * dblReal2Inv; + iErr = 4; // warning } } - else if (_dblReal2 == 0) + + int i; + for (i = 0; i < n; ++i) { - double dblImg2Inv = 1.0 / _dblImg2; - *_pdblRealOut = _dblImg1 * dblImg2Inv; - *_pdblImgOut = (-_dblReal1) * dblImg2Inv; + double complex z = (A[i] + B[i] * I) / (c + d * I); + X[i] = creal(z); + Y[i] = cimag(z); } - else - { - //Generic division algorithm - if (fabs(_dblReal2) >= fabs(_dblImg2)) - { - double dblRatio2 = _dblImg2 / _dblReal2; - double dblSumInv = 1.0 / (_dblReal2 + dblRatio2 * _dblImg2); - *_pdblRealOut = (_dblReal1 + _dblImg1 * dblRatio2) * dblSumInv; - *_pdblImgOut = (_dblImg1 - _dblReal1 * dblRatio2) * dblSumInv; - } - else - { - double dblRatio2 = _dblReal2 / _dblImg2; - double dblSumInv = 1.0 / (_dblImg2 + dblRatio2 * _dblReal2); - *_pdblRealOut = (_dblReal1 * dblRatio2 + _dblImg1) * dblSumInv; - *_pdblImgOut = (_dblImg1 * dblRatio2 - _dblReal1) * dblSumInv; - } - } return iErr; } -int iRightDivisionRealMatrixByComplexMatrix( - double *_pdblReal1, int _iInc1, - double *_pdblReal2, double *_pdblImg2, int _iInc2, - double *_pdblRealOut, double *_pdblImgOut, int _iIncOut, int _iSize) +int iRightDivisionRealMatrixByComplex(double* A, double c, double d, double* X, double* Y, int n) { - int iErr = 0; - int iIndex = 0; //Main loop index - int iIndex1 = 0; //Loop index on left operand - int iIndex2 = 0; //Loop index on right operand - int iIndexOut = 0; //Lopp index on result matrix + // X + Y*i = A / (c + d*i) - if (_iInc2 == 0) + int iErr = 0; + + if (c == 0 && d == 0) { - if ((getieee() == 0) && (fabs(_pdblReal2[iIndex2]) + fabs(_pdblImg2[iIndex2]) == 0)) + int mode = getieee(); + + if (mode == 0) + { + iErr = 3; // error + return iErr; + } + else if (mode == 1) { - return 3; + iErr = 4; // warning } } - for (iIndex = 0 ; iIndex < _iSize ; iIndex++) + int i; + for (i = 0; i < n; ++i) { - iErr = iRightDivisionRealByComplex(_pdblReal1[iIndex1], _pdblReal2[iIndex2], _pdblImg2[iIndex2], &_pdblRealOut[iIndexOut], &_pdblImgOut[iIndexOut]); - iIndexOut += _iIncOut; - iIndex1 += _iInc1; - iIndex2 += _iInc2; + double complex z = A[i] / (c + d * I); + X[i] = creal(z); + Y[i] = cimag(z); } return iErr; } -int iRightDivisionRealByComplex( - double _dblReal1, - double _dblReal2, double _dblImg2, - double *_pdblRealOut, double *_pdblImgOut) +int iRightDivisionComplexMatrixByReal(double* A, double* B, double c, double* X, double* Y, int n) { + // X + Y*i = (A + B*i) / c + int iErr = 0; - if (_dblImg2 == 0) - { - *_pdblRealOut = _dblReal1 / _dblReal2; - *_pdblImgOut = 0; - } - else if (_dblReal2 == 0) - { - *_pdblRealOut = 0; - *_pdblImgOut = -_dblReal1 / _dblImg2; - } - else + + if (c == 0) { - double dblAbsSum = fabs(_dblReal2) + fabs(_dblImg2); + int mode = getieee(); - if (dblAbsSum == 0) + if (mode == 0) { - iErr = 10; - *_pdblRealOut = _dblReal1 / dblAbsSum; - *_pdblImgOut = 0; + iErr = 3; // error + return iErr; } - else + else if (mode == 1) { - double dblAbsSumInv = 1.0 / dblAbsSum; - double dblReal1Sum = _dblReal1 * dblAbsSumInv; - double dblReal2Sum = _dblReal2 * dblAbsSumInv; - double dblImg2Sum = _dblImg2 * dblAbsSumInv; - double dblSumInv = 1.0 / (pow(dblReal2Sum, 2) + pow(dblImg2Sum, 2)); - *_pdblRealOut = (dblReal1Sum * dblReal2Sum) * dblSumInv; - *_pdblImgOut = (-dblReal1Sum * dblImg2Sum) * dblSumInv; + iErr = 4; // warning } } - return iErr; -} -int iRightDivisionComplexMatrixByRealMatrix( - double *_pdblReal1, double *_pdblImg1, int _iInc1, - double *_pdblReal2, int _iInc2, - double *_pdblRealOut, double *_pdblImgOut, int _iIncOut, int _iSize) -{ - int iErr = 0; - int iIndex = 0; //Main loop index - int iIndex1 = 0; //Loop index on left operand - int iIndex2 = 0; //Loop index on right operand - int iIndexOut = 0; //Lopp index on result matrix - for (iIndex = 0 ; iIndex < _iSize ; iIndex++) + double invc = 1.0 / c; + + int i; + for (i = 0; i < n; ++i) { - iErr = iRightDivisionComplexByReal(_pdblReal1[iIndex1], _pdblImg1[iIndex1], _pdblReal2[iIndex2], &_pdblRealOut[iIndexOut], &_pdblImgOut[iIndexOut]); - iIndexOut += _iIncOut; - iIndex1 += _iInc1; - iIndex2 += _iInc2; + X[i] = A[i] * invc; + Y[i] = B[i] * invc; } + return iErr; } -int iRightDivisionComplexByReal( - double _dblReal1, double _dblImg1, - double _dblReal2, - double *_pdblRealOut, double *_pdblImgOut) +int iRightDivisionRealMatrixByReal(double* A, double b, double* X, int n) { + // X = A / b + int iErr = 0; - if (_dblReal2 == 0) + + if (b == 0) { - if (getieee() == 0) + int mode = getieee(); + + if (mode == 0) { - iErr = 3; + iErr = 3; // error return iErr; } - else if (getieee() == 1) + else if (mode == 1) { - //Warning - iErr = 4; + iErr = 4; // warning } } - double dblReal2Inv = 1.0 / _dblReal2; - *_pdblRealOut = _dblReal1 * dblReal2Inv; - *_pdblImgOut = _dblImg1 * dblReal2Inv; - return iErr; -} + double invb = 1 / b; -int iRightDivisionRealMatrixByRealMatrix( - double *_pdblReal1, int _iInc1, - double *_pdblReal2, int _iInc2, - double *_pdblRealOut, int _iIncOut, int _iSize) -{ - int iIndex = 0; //Main loop index - int iIndex1 = 0; //Loop index on left operand - int iIndex2 = 0; //Loop index on right operand - int iIndexOut = 0; //Lopp index on result matrix - int iErr = 0; - - for (iIndex = 0 ; iIndex < _iSize ; iIndex++) + int i; + for (i = 0; i < n; ++i) { - if (_pdblReal2[iIndex2] == 0) - { - if (getieee() == 0) - { - iErr = 3; - return iErr; - } - else if (getieee() == 1) - { - iErr = 4; - } - } - - _pdblRealOut[iIndexOut] = _pdblReal1[iIndex1] / _pdblReal2[iIndex2]; - iIndexOut += _iIncOut; - iIndex1 += _iInc1; - iIndex2 += _iInc2; + X[i] = A[i] * invb; } + return iErr; } diff --git a/scilab/modules/ast/src/cpp/operations/types_divide.cpp b/scilab/modules/ast/src/cpp/operations/types_divide.cpp index 5fafd332207..fb9843c3e25 100644 --- a/scilab/modules/ast/src/cpp/operations/types_divide.cpp +++ b/scilab/modules/ast/src/cpp/operations/types_divide.cpp @@ -2,7 +2,7 @@ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab * Copyright (C) 2008-2008 - DIGITEO - Antoine ELIAS * Copyright (C) 2012 - 2016 - Scilab Enterprises - * Copyright (C) 2017 - Dirk Reusch, Kybernetik Dr. Reusch + * Copyright (C) 2017 - 2018 Dirk Reusch, Kybernetik Dr. Reusch * * This file is hereby licensed under the terms of the GNU GPL v2.0, * pursuant to article 5.3.4 of the CeCILL v.2.1. @@ -103,7 +103,6 @@ InternalType *GenericRDivide(InternalType *_pLeftOperand, InternalType *_pRightO sciprint(_("Warning : Division by zero...\n")); } break; - // default : throw ast::InternalError(_W("Operator / : Error %d not yet managed.\n"), iResult); default : sciprint(_("Operator / : Error %d not yet managed.\n"), iResult); } @@ -131,45 +130,42 @@ int RDivideDoubleByDouble(Double *_pDouble1, Double *_pDouble2, Double **_pDoubl if (_pDouble2->isScalar()) { - //Y / x - int iInc1 = 1; - int iInc2 = 0; + // Y / x bool bComplex1 = _pDouble1->isComplex(); bool bComplex2 = _pDouble2->isComplex(); - *_pDoubleOut = new Double(_pDouble1->getDims(), _pDouble1->getDimsArray(), bComplex1 || bComplex2); + *_pDoubleOut = new Double(_pDouble1->getDims(), _pDouble1->getDimsArray(), bComplex1 || bComplex2); if (bComplex1 == false && bComplex2 == false) { // Real1 \ Real2 -> Real2 / Real1 - iErr = iRightDivisionRealMatrixByRealMatrix( - _pDouble1->get(), iInc1, - _pDouble2->get(), iInc2, - (*_pDoubleOut)->get(), 1, _pDouble1->getSize()); + iErr = iRightDivisionRealMatrixByReal( + _pDouble1->get(), + _pDouble2->getFirst(), + (*_pDoubleOut)->get(), _pDouble1->getSize()); } else if (bComplex1 == false && bComplex2 == true) { // Real \ Complex -> Complex / Real - iErr = iRightDivisionRealMatrixByComplexMatrix( - _pDouble1->get(), iInc1, - _pDouble2->get(), _pDouble2->getImg(), iInc2, - (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), 1, _pDouble1->getSize()); + iErr = iRightDivisionRealMatrixByComplex( + _pDouble1->get(), + _pDouble2->getFirst(), _pDouble2->getImgFirst(), + (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), _pDouble1->getSize()); } else if (bComplex1 == true && bComplex2 == false) { // Complex \ Real -> Real / Complex - iErr = iRightDivisionComplexMatrixByRealMatrix( - _pDouble1->get(), _pDouble1->getImg(), iInc1, - _pDouble2->get(), iInc2, - (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), 1, _pDouble1->getSize()); + iErr = iRightDivisionComplexMatrixByReal( + _pDouble1->get(), _pDouble1->getImg(), _pDouble2->getFirst(), + (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), _pDouble1->getSize()); } else if (bComplex1 == true && bComplex2 == true) { // Complex \ Complex - iErr = iRightDivisionComplexMatrixByComplexMatrix( - _pDouble1->get(), _pDouble1->getImg(), iInc1, - _pDouble2->get(), _pDouble2->getImg(), iInc2, - (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), 1, _pDouble1->getSize()); + iErr = iRightDivisionComplexMatrixByComplex( + _pDouble1->get(), _pDouble1->getImg(), + _pDouble2->getFirst(), _pDouble2->getImgFirst(), + (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg(), _pDouble1->getSize()); } return iErr; @@ -186,7 +182,7 @@ int RDivideDoubleByDouble(Double *_pDouble1, Double *_pDouble2, Double **_pDoubl // x / eye() = x if (_pDouble2->isIdentity() ) { - *_pDoubleOut = new Double(*_pDouble1); + *_pDoubleOut = new Double(*_pDouble1); return 0; } @@ -319,19 +315,19 @@ int RDividePolyByDouble(Polynom* _pPoly, Double* _pDouble, Polynom** _pPolyOut) if (bComplex1 == false && bComplex2 == false) { - iRightDivisionRealMatrixByRealMatrix(pC->get(), 1, &dblDivR, 0, pC->get(), 1, pC->getSize()); + iRightDivisionRealMatrixByReal(pC->get(), dblDivR, pC->get(), pC->getSize()); } else if (bComplex1 == true && bComplex2 == false) { - iRightDivisionComplexMatrixByRealMatrix(pC->get(), pC->getImg(), 1, &dblDivR, 0, pC->get(), pC->getImg(), 1, pC->getSize()); + iRightDivisionComplexMatrixByReal(pC->get(), pC->getImg(), dblDivR, pC->get(), pC->getImg(), pC->getSize()); } else if (bComplex1 == false && bComplex2 == true) { - iRightDivisionRealMatrixByComplexMatrix(pC->get(), 1, &dblDivR, &dblDivI, 0, pC->get(), pC->getImg(), 1, pC->getSize()); + iRightDivisionRealMatrixByComplex(pC->get(), dblDivR, dblDivI, pC->get(), pC->getImg(), pC->getSize()); } else if (bComplex1 == true && bComplex2 == true) { - iRightDivisionComplexMatrixByComplexMatrix(pC->get(), pC->getImg(), 1, &dblDivR, &dblDivI, 0, pC->get(), pC->getImg(), 1, pC->getSize()); + iRightDivisionComplexMatrixByComplex(pC->get(), pC->getImg(), dblDivR, dblDivI, pC->get(), pC->getImg(), pC->getSize()); } } return 0; diff --git a/scilab/modules/ast/tests/unit_tests/dotdivide.dia.ref b/scilab/modules/ast/tests/unit_tests/dotdivide.dia.ref index 146d284a53e..5d48c07ccd6 100644 --- a/scilab/modules/ast/tests/unit_tests/dotdivide.dia.ref +++ b/scilab/modules/ast/tests/unit_tests/dotdivide.dia.ref @@ -141,7 +141,7 @@ assert_checkequal(e .\ i32, int32((- 32)*eye())); assert_checkequal(e .\ ui32, uint32((32)*eye())); assert_checkequal(ec .\ empty, []); assert_checkalmostequal(ec .\ r, (10/26-%i* 1/13)*eye()); -assert_checkequal(ec .\ c, (7/26+%i* 9/26)*eye()); +assert_checkalmostequal(ec .\ c, (7/26+%i* 9/26)*eye()); assert_checkalmostequal(ec .\ e, (5/26-%i* 1/26)*eye()); assert_checkequal(ec .\ ec, (1+%i*0)*eye()); assert_checkequal(p .\ empty, []); @@ -253,7 +253,7 @@ assert_checkequal(SP .\ SP1, sparse([1, 2; 3, 10; 4, 5],[10; (3+1/3); 5],[4, 10] assert_checkequal(SP .\ SPC1, sparse([1, 2; 3, 10; 4, 5],[2+%i* 6; (2/3)+%i* 2; 1+%i* 3],[4, 10])); assert_checkequal(SPC .\ empty, []); assert_checkequal(SPC .\ r, sparse([1, 2; 3, 10; 4, 5],[2/17-%i* 8/17; 2/51-%i* 8/51; 1/17-%i* 4/17],[4, 10])); -assert_checkequal(SPC .\ c, sparse([1, 2; 3, 10; 4, 5],[9/17-%i* 2/17; 3/17-%i* 2/51; 9/34-%i* 1/17],[4, 10])); +assert_checkalmostequal(SPC .\ c, sparse([1, 2; 3, 10; 4, 5],[9/17-%i* 2/17; 3/17-%i* 2/51; 9/34-%i* 1/17],[4, 10])); assert_checkequal(SPC .\ SP, sparse([1, 2; 3, 10; 4, 5],[1/17-%i* 4/17; 1/17-%i* 4/17; 1/17-%i* 4/17],[4, 10])); assert_checkequal(SPC .\ SPC, sparse([1, 2; 3, 10; 4, 5],[1; 1; 1],[4, 10])); assert_checkequal(SPC .\ SP1, sparse([1, 2; 3, 10; 4, 5],[10/17-%i* 40/17; 10/51-%i* 40/51; 10/34-%i* 20/17],[4, 10])); @@ -500,7 +500,7 @@ assert_checkequal(c ./ c, 1+%i*0); assert_checkequal(c ./ R, [ 1+%i*2, 0.5+%i; (1/3)+%i*(2/3), 0.25+%i*0.5]); assert_checkequal(c ./ C, [1, 0.5; (1/3), 0.25]+%i*0); assert_checkequal(c ./ e, (1+%i* 2)*eye()); -assert_checkequal(c ./ ec, (7/26+%i* 9/26)*eye()); +assert_checkalmostequal(c ./ ec, (7/26+%i* 9/26)*eye()); test = c ./ p; assert_checkequal(test(2), 1+%i*2+0*%s); assert_checkequal(test(3), 1+(1)*s+(- 1)*s^ 2); @@ -514,7 +514,7 @@ test = c ./ PC; assert_checkequal(test(2), [1,1;1,1]*(1+%i*2)+0*%s); assert_checkequal(test(3),[2+%i* 4+(2-%i* 6)*s+(- 2+%i* 8)*s^ 2,- 3-%i* 6+(- 3+%i* 9)*s+(3-%i* 12)*s^ 2;4+%i* 8+(4-%i* 12)*s+(- 4+%i* 16)*s^ 2,- 5-%i* 10+(- 5+%i* 15)*s+(5-%i* 20)*s^ 2]); assert_checkequal(c ./ SP, sparse([1, 2; 3, 10; 4, 5],[1+%i* 2; (1/3)+%i* (2/3); 0.5+%i],[4, 10])); -assert_checkequal(c ./ SPC, sparse([1, 2; 3, 10; 4, 5],[9/17-%i* 2/17; 3/17-%i* 2/51; 9/34-%i* 1/17],[4, 10])); +assert_checkalmostequal(c ./ SPC, sparse([1, 2; 3, 10; 4, 5],[9/17-%i* 2/17; 3/17-%i* 2/51; 9/34-%i* 1/17],[4, 10])); assert_checkequal(c ./ SP1, sparse([1, 1], 0.1+%i* 0.2,[1, 1])); assert_checkalmostequal(c ./ SPC1, sparse([1, 1], 0.35-%i* 0.05,[1, 1])); assert_checkequal(R ./ empty, []); diff --git a/scilab/modules/ast/tests/unit_tests/dotdivide.tst b/scilab/modules/ast/tests/unit_tests/dotdivide.tst index 541ee38eae8..3a0a4618c64 100644 --- a/scilab/modules/ast/tests/unit_tests/dotdivide.tst +++ b/scilab/modules/ast/tests/unit_tests/dotdivide.tst @@ -1,6 +1,7 @@ // ============================================================================ // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab // Copyright (C) 2014 - Scilab Enterprises - Sylvain GENIN +// Copyright (C) 2018 - Dirk Reusch, Kybernetik Dr. Reusch // // This file is distributed under the same license as the Scilab package. // ============================================================================ @@ -153,7 +154,7 @@ assert_checkequal(e .\ ui32, uint32((32)*eye())); assert_checkequal(ec .\ empty, []); assert_checkalmostequal(ec .\ r, (10/26-%i* 1/13)*eye()); -assert_checkequal(ec .\ c, (7/26+%i* 9/26)*eye()); +assert_checkalmostequal(ec .\ c, (7/26+%i* 9/26)*eye()); assert_checkalmostequal(ec .\ e, (5/26-%i* 1/26)*eye()); assert_checkequal(ec .\ ec, (1+%i*0)*eye()); @@ -272,7 +273,7 @@ assert_checkequal(SP .\ SPC1, sparse([1, 2; 3, 10; 4, 5],[2+%i* 6; (2/3)+%i* 2; assert_checkequal(SPC .\ empty, []); assert_checkequal(SPC .\ r, sparse([1, 2; 3, 10; 4, 5],[2/17-%i* 8/17; 2/51-%i* 8/51; 1/17-%i* 4/17],[4, 10])); -assert_checkequal(SPC .\ c, sparse([1, 2; 3, 10; 4, 5],[9/17-%i* 2/17; 3/17-%i* 2/51; 9/34-%i* 1/17],[4, 10])); +assert_checkalmostequal(SPC .\ c, sparse([1, 2; 3, 10; 4, 5],[9/17-%i* 2/17; 3/17-%i* 2/51; 9/34-%i* 1/17],[4, 10])); assert_checkequal(SPC .\ SP, sparse([1, 2; 3, 10; 4, 5],[1/17-%i* 4/17; 1/17-%i* 4/17; 1/17-%i* 4/17],[4, 10])); assert_checkequal(SPC .\ SPC, sparse([1, 2; 3, 10; 4, 5],[1; 1; 1],[4, 10])); assert_checkequal(SPC .\ SP1, sparse([1, 2; 3, 10; 4, 5],[10/17-%i* 40/17; 10/51-%i* 40/51; 10/34-%i* 20/17],[4, 10])); @@ -537,7 +538,7 @@ assert_checkequal(c ./ c, 1+%i*0); assert_checkequal(c ./ R, [ 1+%i*2, 0.5+%i; (1/3)+%i*(2/3), 0.25+%i*0.5]); assert_checkequal(c ./ C, [1, 0.5; (1/3), 0.25]+%i*0); assert_checkequal(c ./ e, (1+%i* 2)*eye()); -assert_checkequal(c ./ ec, (7/26+%i* 9/26)*eye()); +assert_checkalmostequal(c ./ ec, (7/26+%i* 9/26)*eye()); test = c ./ p; assert_checkequal(test(2), 1+%i*2+0*%s); assert_checkequal(test(3), 1+(1)*s+(- 1)*s^ 2); @@ -551,7 +552,7 @@ test = c ./ PC; assert_checkequal(test(2), [1,1;1,1]*(1+%i*2)+0*%s); assert_checkequal(test(3),[2+%i* 4+(2-%i* 6)*s+(- 2+%i* 8)*s^ 2,- 3-%i* 6+(- 3+%i* 9)*s+(3-%i* 12)*s^ 2;4+%i* 8+(4-%i* 12)*s+(- 4+%i* 16)*s^ 2,- 5-%i* 10+(- 5+%i* 15)*s+(5-%i* 20)*s^ 2]); assert_checkequal(c ./ SP, sparse([1, 2; 3, 10; 4, 5],[1+%i* 2; (1/3)+%i* (2/3); 0.5+%i],[4, 10])); -assert_checkequal(c ./ SPC, sparse([1, 2; 3, 10; 4, 5],[9/17-%i* 2/17; 3/17-%i* 2/51; 9/34-%i* 1/17],[4, 10])); +assert_checkalmostequal(c ./ SPC, sparse([1, 2; 3, 10; 4, 5],[9/17-%i* 2/17; 3/17-%i* 2/51; 9/34-%i* 1/17],[4, 10])); assert_checkequal(c ./ SP1, sparse([1, 1], 0.1+%i* 0.2,[1, 1])); assert_checkalmostequal(c ./ SPC1, sparse([1, 1], 0.35-%i* 0.05,[1, 1])); diff --git a/scilab/modules/elementary_functions/src/c/expm.c b/scilab/modules/elementary_functions/src/c/expm.c index 586c94a10d7..876a4370162 100644 --- a/scilab/modules/elementary_functions/src/c/expm.c +++ b/scilab/modules/elementary_functions/src/c/expm.c @@ -3,7 +3,7 @@ * Copyright (C) 2006 - INRIA - Allan CORNET * Copyright (C) 2012 - Digiteo - Cedric Delamarre * Copyright (C) 2012 - 2016 - Scilab Enterprises - * Copyright (C) 2017 - Dirk Reusch, Kybernetik Dr. Reusch + * Copyright (C) 2017 - 2018 Dirk Reusch, Kybernetik Dr. Reusch * * This file is hereby licensed under the terms of the GNU GPL v2.0, * pursuant to article 5.3.4 of the CeCILL v.2.1. @@ -106,17 +106,12 @@ int zexpms2(double *_pdblReal, double *_pdblImg, double *_pdblReturnReal, double //A = A./2^s if (iComplex) { - iRightDivisionComplexMatrixByRealMatrix( - _pdblReal, _pdblImg , 1, - &dblS2, 0, - pdblMatrixRealA, pdblMatrixImgA, 1, iSquare); + iRightDivisionComplexMatrixByReal(_pdblReal,_pdblImg, dblS2, + pdblMatrixRealA, pdblMatrixImgA, iSquare); } else { - iRightDivisionRealMatrixByRealMatrix( - _pdblReal, 1, - &dblS2, 0, - pdblMatrixRealA, 1, iSquare); + iRightDivisionRealMatrixByReal(_pdblReal, dblS2, pdblMatrixRealA, iSquare); } // Pade approximation for exp(A)