| 3039 | } |
| 3040 | |
| 3041 | void ecall_sample(float* input1,float* input2, int N,int nums, float* output){ |
| 3042 | float *x = (float*)malloc(sizeof(float)*N); |
| 3043 | memset(x, 0, sizeof(float)*N); |
| 3044 | for(int j=0;j<N;j++){ |
| 3045 | x[j] = 0.0f; |
| 3046 | float sum = input1[insert*N+j]; |
| 3047 | float c = 0.0f; |
| 3048 | for(int i=0;i<Nt;i++){ |
| 3049 | float y = input1[indexs[i]*N+j] - c; |
| 3050 | float t = sum + y; |
| 3051 | c = (t - sum) - y; |
| 3052 | sum = t; |
| 3053 | } |
| 3054 | x[j] = sum; |
| 3055 | } |
| 3056 | |
| 3057 | float *y = (float*)malloc(sizeof(float)*N); |
| 3058 | memset(y, 0, sizeof(float)*N); |
| 3059 | for(int j=0;j<N;j++){ |
| 3060 | y[j] = 0.0f; |
| 3061 | float sum = input2[insert*N+j]; |
| 3062 | float c = 0.0f; |
| 3063 | for(int i=0;i<Nt;i++){ |
| 3064 | float z = input2[indexs[i]*N+j] - c; |
| 3065 | float t = sum + z; |
| 3066 | c = (t - sum) - z; |
| 3067 | sum = t; |
| 3068 | } |
| 3069 | y[j] = sum; |
| 3070 | } |
| 3071 | |
| 3072 | float *z = (float*)malloc(sizeof(float)*N); |
| 3073 | memset(z, 0, sizeof(float)*N); |
| 3074 | for(int num = 0;num<nums;num++){ |
| 3075 | for(int i=0;i<N;i++){ |
| 3076 | z[i] = gaussrand(x[i], y[i]); |
| 3077 | } |
| 3078 | float mean = 0.0f, std, mean_sqr = 0.0f; |
| 3079 | float c_mean = 0.0f, c_mean_sqr = 0.0f; |
| 3080 | for (int j = 0; j < N; j++) { |
| 3081 | float y = (z[j]/static_cast<float>(N)) - c_mean; |
| 3082 | float t = mean + y; |
| 3083 | c_mean = (t - mean) - y; |
| 3084 | mean = t; |
| 3085 | |
| 3086 | float y2 = (z[j]*z[j]/static_cast<float>(N)) - c_mean_sqr; |
| 3087 | float t2 = mean_sqr + y2; |
| 3088 | c_mean_sqr = (t2 - mean_sqr) - y2; |
| 3089 | mean_sqr = t2; |
| 3090 | } |
| 3091 | std = sqrtf(mean_sqr - mean * mean)/static_cast<float>(Nt); |
| 3092 | mean = mean/static_cast<float>(Nt); |
| 3093 | for(int i=0;i<N*(Ne-1);i++){ |
| 3094 | output[num*Ne*N+i] =(z[i%N] - gaussrand(mean, std))/static_cast<float>(Nt+1); |
| 3095 | } |
| 3096 | for(int j=N*(Ne-1);j<N*Ne;j++){ |
| 3097 | output[num*Ne*N+j] = 0.0f; |
| 3098 | float sum = 0.0f; |