-
Notifications
You must be signed in to change notification settings - Fork 0
/
IntHZF.cpp
124 lines (95 loc) · 2.09 KB
/
IntHZF.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
#include<iostream>
#include<stdlib.h> // atof is in this library...
#include<fstream>
#include<cmath>
#include<ctime>
//#include<tgmath.h>
using namespace std;
double pi = 3.141592653589793;
double s;
double q;
double myfunc(double *x);
double myfunc2(double x);
double newfunc(double x);
double mean(double *x,int len);
double stdev(double *x,double xm, int len);
double mean(double *x,int len){
double sum = 0.0;
for(int i=0;i<len;i++){
sum += x[i];
}
sum /= (double)len;
return sum;
}
double stdev(double *x, double xm, int len){
double sum = 0.0;
for(int i=0;i<len;i++){
sum += pow(x[i]-xm,2);
}
sum /= (double)len;
return sqrt(sum);
}
double myfunc2(double x){
return exp(-q*x)/(1 - exp(-x));
}
double newfunc(double x){
return pow(x,(s - 1));
}
int main(int argc, char *argv[]){
int seed;
double r;
double *x, *y;
double xmag = 0.0;
int N, d;
int M = 1000;
int h = 0;
double a,b;
double *areas;
double ave, err;
clock_t start, stop;
//int s;
//int q;
cout<<"Give me the value of s ";
cin>>s;
cout<<"Give me the value of q ";
cin>>q;
a = 0.;
b = 1e6;
if(argc<2){
N = 100000;
cout<<"N = "<<N<<endl;
}
else{
N = atoi(argv[1]);
}
double area;
seed = 100;//time(NULL);
srand(seed);// "Seed the random number generator"
x = new double [1];
y = new double [1];
areas = new double [M];
start = clock();
for(int k=0;k<M;k++){
areas[k] = 0.0;
for(int i=0;i<N;i++){
xmag = 0.0;
for(int j=0;j<1;j++){
r = ((double) rand())/((double)RAND_MAX);
x[j] = r*(b - a) + a; // 0 < x[j] < 1
y[j] = -log(1-x[j]/b); // 0 < x[j] < 1
}
areas[k] += newfunc(y[0]);
}// end of for i loop
areas[k] /= ((double) N);
} // end of for k loop
delete y;
delete x;
stop = clock();
double totaltime = (stop - start)/((double)CLOCKS_PER_SEC);
ave = mean(areas,M);
err = stdev(areas,ave,M)/sqrt((double)M - 1.);
delete areas;
cout<<"The integral is "<<ave<<"+/-"<<err<<endl;
// cout<<"It should be "<<pi<<endl;
return 0;
}