From 57e712a0062aafc55b6278f93c93ad4b52b87248 Mon Sep 17 00:00:00 2001
From: Krad <kradchen123@163.com>
Date: Sun, 23 Apr 2023 11:12:00 +0800
Subject: [PATCH] Add polyval to Function1D

---
 src/Function.cpp         | 18 ------------------
 src/Function.h           |  1 -
 src/Function1D.cpp       | 19 +++++++++++++++++++
 src/Function1D.h         | 10 ++++++++++
 test/Function1D_Test.cpp | 22 ++++++++++++++--------
 5 files changed, 43 insertions(+), 27 deletions(-)

diff --git a/src/Function.cpp b/src/Function.cpp
index 9b2fe1d..8fa3845 100644
--- a/src/Function.cpp
+++ b/src/Function.cpp
@@ -20,24 +20,6 @@ namespace Aurora {
         return result;
     }
 
-    double *polyval(double *pX, double *pP, int size) {
-        auto result = new double[size];
-        auto powArg = new double[size];
-        for (int j = size, i = 0; j > 0; --j, ++i) {
-            powArg[i] = (double) (j - 1);
-        }
-        auto temp = new double[size];
-        for (int i = 0; i < size; ++i) {
-            vdPowI(size, pX + i, 0, powArg, 1, temp, 1);
-            vdMul(size, pP, temp, temp);
-            result[i] = cblas_dasum(3, temp, 1);
-        }
-        delete[] powArg;
-        delete[] temp;
-
-        return result;
-    }
-
     double *std(int rows, int cols, double *input) {
         auto std = new double[cols];
 //    #pragma omp parallel for num_threads(4)
diff --git a/src/Function.h b/src/Function.h
index ced7354..04ee914 100644
--- a/src/Function.h
+++ b/src/Function.h
@@ -14,7 +14,6 @@ namespace Aurora{
     double* mul( double* inputA, double * inputB, int size);
     double* mulz( std::complex<double> *inputA, std::complex<double> * inputB, int size);
     double immse(double * dataA, double * dataB, int size);
-    double* polyval(double *pX, double *pP, int size);
     double* std(int rows, int cols, double * input);
 
     double* inv(int cols,double *pMatrix);
diff --git a/src/Function1D.cpp b/src/Function1D.cpp
index 985939b..a3bbcb2 100644
--- a/src/Function1D.cpp
+++ b/src/Function1D.cpp
@@ -303,3 +303,22 @@ Matrix Aurora::repmat(const Matrix& aMatrix,int aRowTimes, int aColumnTimes, int
 
     return Matrix(std::shared_ptr<double>(resultData, Aurora::free), resultInfo, aMatrix.getValueType());
 }
+
+
+Matrix Aurora::polyval(const Matrix &aP, const Matrix &aX) {
+    auto result = malloc(aX.getDataSize());
+    auto powArg = new double[aP.getDataSize()];
+    for (int j = aP.getDataSize(), i = 0; j > 0; --j, ++i) {
+        powArg[i] = (double) (j - 1);
+    }
+    auto temp = new double[aP.getDataSize()];
+    for (int i = 0; i < aX.getDataSize(); ++i) {
+        vdPowI(aP.getDataSize(), aX.getData() + i, 0, powArg, 1, temp, 1);
+        vdMul(aP.getDataSize(), aP.getData(), temp, temp);
+        result[i] = cblas_dasum(3, temp, 1);
+    }
+    delete[] powArg;
+    delete[] temp;
+
+    return Matrix::New(result,aX);
+}
diff --git a/src/Function1D.h b/src/Function1D.h
index 79a5ed0..58dd1cd 100644
--- a/src/Function1D.h
+++ b/src/Function1D.h
@@ -48,6 +48,16 @@ namespace Aurora {
     Matrix repmat(const Matrix& aMatrix,int aRowTimes, int aColumnTimes, int aSliceTimes);
 
     Matrix log(const Matrix& aMatrix, int aBaseNum = -1);
+
+    /**
+     * 多项式计算
+     * @brief 例如p[1 0 1],x[3 2 5],代表对多项式 y = x^2 + 1 求(x=3, x=2, x=5)时所有的y
+     * @attention 本函数的数据一律视为向量，在计算时会完全忽视维度，值考虑作为数组输入
+     * @param aP 多项式系数，指定为向量
+     * @param aX 查询点，指定为向量
+     * @return 查询结果
+     */
+    Matrix polyval(const Matrix& aP, const Matrix& aX);
 };
 
 
diff --git a/test/Function1D_Test.cpp b/test/Function1D_Test.cpp
index 282152a..2651769 100644
--- a/test/Function1D_Test.cpp
+++ b/test/Function1D_Test.cpp
@@ -38,6 +38,11 @@ TEST_F(Function1D_Test, sign) {
     EXPECT_EQ(-1, result[0]);
     EXPECT_EQ(0, result[3]);
     EXPECT_EQ(1, result[4]);
+    C = sign(A);
+    result = C.getData();
+    EXPECT_EQ(-1, result[0]);
+    EXPECT_EQ(0, result[3]);
+    EXPECT_EQ(1, result[4]);
 }
 
 TEST_F(Function1D_Test, repmat) {
@@ -102,13 +107,14 @@ TEST_F(Function1D_Test, interp1) {
 }
 
 TEST_F(Function1D_Test, polyval){
-    double dataP[3]={3,2,1};
-    double dataX[3]={5,7,9};
-    double*resultP = Aurora::polyval(dataX,dataP,3);
-    EXPECT_DOUBLE_EQ(86., resultP[0])<<" polyval error;";
-    EXPECT_DOUBLE_EQ(162., resultP[1])<<" polyval error;";
-    EXPECT_DOUBLE_EQ(262., resultP[2])<<" polyval error;";
-    delete [] resultP;
+    double *dataP = new double[3]{3,2,1};
+    double *dataX = new double[3]{5,7,9};
+    auto matrixP = Aurora::Matrix::fromRawData(dataP,3);
+    auto matrixX = Aurora::Matrix::fromRawData(dataX,3);
+    auto result = Aurora::polyval(matrixP,matrixX);
+    EXPECT_DOUBLE_EQ(86., result.getData()[0]);
+    EXPECT_DOUBLE_EQ(162., result.getData()[1]);
+    EXPECT_DOUBLE_EQ(262., result.getData()[2]);
 }
 
 TEST_F(Function1D_Test, complexAndEtc){
@@ -200,7 +206,7 @@ TEST_F(Function1D_Test, absAndSqrt) {
     EXPECT_EQ(Aurora::Normal, matrixB.getValueType());
     EXPECT_DOUBLE_EQ(1, matrixB.getData()[0]);
     EXPECT_DOUBLE_EQ(2, matrixB.getData()[1]);
-    EXPECT_TRUE(isnanf(matrixB.getData()[2]));
+    EXPECT_TRUE(std::isnan(matrixB.getData()[2]));
 }