1// -*- mode: c++ -*- 2/* 3 Copyright (c) 2010-2011, Intel Corporation 4 All rights reserved. 5 6 Redistribution and use in source and binary forms, with or without 7 modification, are permitted provided that the following conditions are 8 met: 9 10 * Redistributions of source code must retain the above copyright 11 notice, this list of conditions and the following disclaimer. 12 13 * Redistributions in binary form must reproduce the above copyright 14 notice, this list of conditions and the following disclaimer in the 15 documentation and/or other materials provided with the distribution. 16 17 * Neither the name of Intel Corporation nor the names of its 18 contributors may be used to endorse or promote products derived from 19 this software without specific prior written permission. 20 21 22 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS 23 IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED 24 TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 25 PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER 26 OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 27 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 28 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 29 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 30 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 31 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 32 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 33*/ 34 35#include "options_defs.h" 36 37// Cumulative normal distribution function 38static inline float 39CND(float X) { 40 float L = abs(X); 41 42 float k = 1.0 / (1.0 + 0.2316419 * L); 43 float k2 = k*k; 44 float k3 = k2*k; 45 float k4 = k2*k2; 46 float k5 = k3*k2; 47 48 const float invSqrt2Pi = 0.39894228040f; 49 float w = (0.31938153f * k - 0.356563782f * k2 + 1.781477937f * k3 + 50 -1.821255978f * k4 + 1.330274429f * k5); 51 w *= invSqrt2Pi * exp(-L * L * .5f); 52 53 if (X > 0.f) 54 w = 1.0 - w; 55 return w; 56} 57 58task void 59bs_task(uniform float Sa[], uniform float Xa[], uniform float Ta[], 60 uniform float ra[], uniform float va[], 61 uniform float result[], uniform int count) { 62 uniform int first = taskIndex * (count/taskCount); 63 uniform int last = min(count, (uniform int)((taskIndex+1) * (count/taskCount))); 64 65 foreach (i = first ... last) { 66 float S = Sa[i], X = Xa[i], T = Ta[i], r = ra[i], v = va[i]; 67 68 float d1 = (log(S/X) + (r + v * v * .5f) * T) / (v * sqrt(T)); 69 float d2 = d1 - v * sqrt(T); 70 71 result[i] = S * CND(d1) - X * exp(-r * T) * CND(d2); 72 } 73} 74 75export void 76black_scholes_ispc_tasks(uniform float Sa[], uniform float Xa[], uniform float Ta[], 77 uniform float ra[], uniform float va[], 78 uniform float result[], uniform int count) { 79 uniform int nTasks = max((uniform int)64, (uniform int)count/16384); 80 launch[nTasks] bs_task(Sa, Xa, Ta, ra, va, result, count); 81} 82 83 84export void 85black_scholes_ispc(uniform float Sa[], uniform float Xa[], uniform float Ta[], 86 uniform float ra[], uniform float va[], 87 uniform float result[], uniform int count) { 88 foreach (i = 0 ... count) { 89 float S = Sa[i], X = Xa[i], T = Ta[i], r = ra[i], v = va[i]; 90 91 float d1 = (log(S/X) + (r + v * v * .5f) * T) / (v * sqrt(T)); 92 float d2 = d1 - v * sqrt(T); 93 94 result[i] = S * CND(d1) - X * exp(-r * T) * CND(d2); 95 } 96} 97 98 99static inline float 100binomial_put(float S, float X, float T, float r, float v) { 101 float V[BINOMIAL_NUM]; 102 103 float dt = T / BINOMIAL_NUM; 104 float u = exp(v * sqrt(dt)); 105 float d = 1. / u; 106 float disc = exp(r * dt); 107 float Pu = (disc - d) / (u - d); 108 109 for (uniform int j = 0; j < BINOMIAL_NUM; ++j) { 110 float upow = pow(u, (uniform float)(2*j-BINOMIAL_NUM)); 111 V[j] = max(0., X - S * upow); 112 } 113 114 for (uniform int j = BINOMIAL_NUM-1; j >= 0; --j) 115 for (uniform int k = 0; k < j; ++k) 116 V[k] = ((1 - Pu) * V[k] + Pu * V[k + 1]) / disc; 117 return V[0]; 118} 119 120 121export void 122binomial_put_ispc(uniform float Sa[], uniform float Xa[], uniform float Ta[], 123 uniform float ra[], uniform float va[], 124 uniform float result[], uniform int count) { 125 foreach (i = 0 ... count) { 126 float S = Sa[i], X = Xa[i], T = Ta[i], r = ra[i], v = va[i]; 127 result[i] = binomial_put(S, X, T, r, v); 128 } 129} 130 131 132task void 133binomial_task(uniform float Sa[], uniform float Xa[], 134 uniform float Ta[], uniform float ra[], 135 uniform float va[], uniform float result[], 136 uniform int count) { 137 uniform int first = taskIndex * (count/taskCount); 138 uniform int last = min(count, (uniform int)((taskIndex+1) * (count/taskCount))); 139 140 foreach (i = first ... last) { 141 float S = Sa[i], X = Xa[i], T = Ta[i], r = ra[i], v = va[i]; 142 result[i] = binomial_put(S, X, T, r, v); 143 } 144} 145 146 147export void 148binomial_put_ispc_tasks(uniform float Sa[], uniform float Xa[], 149 uniform float Ta[], uniform float ra[], 150 uniform float va[], uniform float result[], 151 uniform int count) { 152 uniform int nTasks = max((uniform int)64, (uniform int)count/16384); 153 launch[nTasks] binomial_task(Sa, Xa, Ta, ra, va, result, count); 154} 155