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