1 // Copyright (c) 2015-2016, Massachusetts Institute of Technology
2 // Copyright (c) 2016-2017 Sandia Corporation
3 // Copyright (c) 2017 NTESS, LLC.
4
5 // This file is part of the Compressed Continuous Computation (C3) Library
6 // Author: Alex A. Gorodetsky
7 // Contact: alex@alexgorodetsky.com
8
9 // All rights reserved.
10
11 // Redistribution and use in source and binary forms, with or without modification,
12 // are permitted provided that the following conditions are met:
13
14 // 1. Redistributions of source code must retain the above copyright notice,
15 // this list of conditions and the following disclaimer.
16
17 // 2. Redistributions in binary form must reproduce the above copyright notice,
18 // this list of conditions and the following disclaimer in the documentation
19 // and/or other materials provided with the distribution.
20
21 // 3. Neither the name of the copyright holder nor the names of its contributors
22 // may be used to endorse or promote products derived from this software
23 // without specific prior written permission.
24
25 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
26 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
28 // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
29 // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
30 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
31 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
32 // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
33 // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
34 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
35
36 //Code
37
38
39
40
41
42 /** \file hpoly.c
43 * Provides routines for manipulating hermite polynomials
44 */
45
46 #include <stdlib.h>
47 #include <stdio.h>
48 #include <math.h>
49 #include <assert.h>
50
51 #include "hpoly.h"
52
53 //for orthogonal
zero_seq11(size_t n)54 inline static double zero_seq11(size_t n){ (void) n; return 0.0; }
55 /* inline static double one_seq11(size_t n) { (void) n; return 1.0; } */
56 /* inline static double hermc_seq(size_t n) */
57 /* { */
58 /* return -((double)n - 1.0); */
59 /* } */
60
61 // for orthonormal
an(size_t n)62 inline static double an(size_t n){ return 1.0 / sqrt((double)n);}
cn(size_t n)63 inline static double cn(size_t n){ return - sqrt(((double)n-1)/(double)n);}
64
65
hermortho(size_t n)66 static inline double hermortho(size_t n)
67 {
68 (void)n;
69 return 1.0;
70 /* double val = 1.0; */
71 /* for (size_t ii = 1; ii < n; ii++){ */
72 /* val *= (ii+1); */
73 /* } */
74 /* //return sqrt(2*M_PI)*val; */
75 /* return val; */
76 }
77
78 /********************************************************//**
79 * Initialize a Hermite polynomial
80 *
81 * \return p - polynomial
82 *************************************************************/
init_hermite_poly()83 struct OrthPoly * init_hermite_poly(){
84
85 struct OrthPoly * p;
86 if ( NULL == (p = malloc(sizeof(struct OrthPoly)))){
87 fprintf(stderr, "failed to allocate memory for poly exp.\n");
88 exit(1);
89 }
90 p->ptype = HERMITE;
91 /* p->an = &one_seq11; */
92 p->an = an;
93 p->bn = &zero_seq11;
94 /* p->cn = &hermc_seq; */
95 p->cn = cn;
96
97 p->lower = -DBL_MAX;
98 p->upper = DBL_MAX;
99
100 p->const_term = 1.0;
101 p->lin_coeff = 1.0;
102 p->lin_const = 0.0;
103
104 p->norm = hermortho;
105
106 return p;
107 }
108
109 /********************************************************//**
110 * Integrate a Hermite approximation
111 *
112 * \param[in] poly - polynomial to integrate
113 *
114 * \return out - integral of approximation
115 *************************************************************/
hermite_integrate(const struct OrthPolyExpansion * poly)116 double hermite_integrate(const struct OrthPolyExpansion * poly)
117 {
118 double out = poly->coeff[0];//sqrt(2.0*M_PI);
119 return out;
120 }
121