dino.celebic/scf_celebic
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

many changes

Browse source
This commit is contained in:
dino.celebic 2025-11-03 22:35:52 +01:00
commit 7e2626266e
35 changed files with 276 additions and 6417 deletions
Download patch file Download diff file Expand all files Collapse all files

30
ex1/code/Makefile Normal file
View file

@ -0,0 +1,30 @@
PROGRAM = main
SOURCES = main.cpp mylib.cpp
OBJECTS = ${SOURCES:.cpp=.o}
CXX = g++
LINKER = g++
WARNINGS = -Wall -pedantic -Wextra -Weffc++ -Woverloaded-virtual -Wfloat-equal -Wshadow \
-Wredundant-decls -fmax-errors=1
CXXFLAGS = -g -flto -O3 -ffast-math -march=native ${WARNINGS}
LINKFLAGS = -g -flto -O3
all: ${PROGRAM}
%.o: %.cpp
${CXX} ${CXXFLAGS} -c $< -o $@
${PROGRAM}: ${OBJECTS}
$(LINKER) ${OBJECTS} ${LINKFLAGS} -o ${PROGRAM}
clean:
rm -f ${OBJECTS} ${PROGRAM}
run: ${PROGRAM}
# run: clean ${PROGRAM}
./${PROGRAM}

500
ex1/code/data_1.txt Normal file
View file

@ -0,0 +1,500 @@
141
261
87
430
258
298
425
120
496
707
244
786
75
394
4
221
2
190
143
269
175
139
599
902
940
222
483
377
524
265
69
437
174
27
955
431
962
763
8
681
706
646
553
219
773
229
371
891
857
403
319
609
911
910
592
333
854
443
905
34
533
717
180
337
188
322
404
549
49
553
275
242
244
155
957
936
819
729
176
361
189
2
317
700
626
544
440
288
502
762
763
577
748
646
124
505
348
93
148
199
673
432
695
257
10
533
280
947
907
393
25
672
838
972
57
451
583
687
720
651
727
374
582
117
58
980
285
595
963
186
194
342
933
391
274
152
398
375
132
436
92
615
11
574
790
236
449
570
62
497
643
222
838
972
847
506
279
747
237
958
621
601
173
91
256
859
912
700
726
230
577
811
404
989
90
321
512
61
726
557
530
830
859
790
318
453
753
110
110
270
525
973
711
312
292
851
912
640
256
89
839
585
949
62
585
286
828
191
443
394
827
677
208
319
134
672
571
170
148
477
909
553
33
54
806
452
383
790
365
533
712
872
329
651
975
76
588
414
310
264
759
996
187
782
196
993
803
425
729
499
809
357
74
591
911
194
433
750
40
947
764
559
184
498
518
995
855
963
679
404
935
480
232
397
706
559
757
996
963
536
964
116
52
305
581
531
902
541
432
543
713
17
801
143
479
257
370
662
170
279
199
196
327
881
472
404
180
969
408
845
616
377
878
785
465
814
899
430
335
597
902
703
378
735
955
543
541
312
72
182
93
464
10
916
643
2
31
209
455
128
9
728
355
781
437
437
50
50
92
595
242
842
858
964
489
221
227
537
763
348
462
640
918
162
716
578
434
885
394
179
634
625
328
803
1000
981
128
233
24
608
111
408
885
549
370
209
441
957
125
471
857
44
692
979
284
134
686
910
611
900
194
755
347
419
156
820
625
739
806
68
951
498
756
743
832
157
458
619
933
836
896
583
583
855
35
886
408
37
747
155
144
606
255
325
402
407
387
610
167
189
95
324
770
235
741
693
825
828
294
310
524
326
832
811
557
263
681
234
457
385
539
992
756
981
235
529
52
757
602
858
989
930
410
1
541
208
220
326
96
748
749
544
339
833
553
958
893
357
547
347
623
797
746
126
823
26
415
732
782
368

386
ex1/code/main.cpp Normal file
View file

@ -0,0 +1,386 @@
// g++ *.cpp -o main
// g++ -g -ffast-math -O3 -march=native -Wall -pedantic -Wextra -Weffc++ -Woverloaded-virtual -Wfloat-equal -Wshadow -Wredundant-decls -fmax-errors=1 *.cpp -o main
#include "mylib.h"
#include "timing.h"
#include <cassert> // assert
#include <vector>
#include <iostream>
#include <cmath>
#include <tuple>
#include <string>
#include <algorithm>
#include <fstream>
#include <list>
#include <stdexcept>
using namespace std;
using namespace std::chrono; // timing
static void task_a() {
printf("\n\n-------------- Task A --------------\n\n");
auto [a,b,c] = means0(1,4,16);
auto [d,e,f] = means0(2,3,5);
auto [g,h,i] = means0(1000,4000,16000);
printf("means(1,4,16) = (%f, %f, %f)\n", a, b, c);
printf("means(2,3,5) = (%f, %f, %f)\n", d, e, f);
printf("means(1000,4000,16000) = (%f, %f, %f)\n", g, h, i);
vector<double> v = {4,8,15,16,23,42};
auto [j,k,l] = means(v);
printf("means(4,8,15,16,23,42) = (%f, %f, %f)\n", j, k, l);
}
static void task_b() {
printf("\n\n-------------- Task B --------------\n\n");
// Read vector
vector<double> a;
read_vector_from_file("data_1.txt", a);
// Print numbers
// for (unsigned int k=0; k<a.size(); ++k)
// {
// cout << " " << a.at(k);
// }
// cout << endl;
// min and max
auto min = min_element(a.begin(), a.end());
auto max = max_element(a.begin(), a.end());
printf("Minimum: %f\n", *min);
printf("Maximum: %f\n", *max);
// means
auto [x,y,z] = means(a);
printf("Arithmetic: %f\n", x);
printf("Geometric: %f\n", y);
printf("Harmonic: %f\n", z);
// deviation
double deviation(0.0);
for (long unsigned int i=0; i<a.size(); i++){
deviation += pow(x - a.at(i),2);
}
deviation = sqrt(deviation/static_cast<double>(a.size()));
printf("Deviation: %f\n", deviation);
// write results to file
vector<double> b = {*min,*max,x,y,z,deviation};
write_vector_to_file("out_1.txt", b);
}
static void task_c() {
printf("\n\n-------------- Task C --------------\n\n");
vector<int> n_values = {15, 1001, 1432987};
for (int n : n_values) {
printf("n = %d\n", n);
long int sum = 0;
double loops = 1000;
// Timing first function
tic();
for (int i=0; i<loops; i++){
sum = sum_of_spec(n);
}
double sec1 = toc();
printf("For-loop funtion: result = %ld | time = %f milliseconds\n", sum, sec1*1000);
// Timing second function
tic();
for (int i=0; i<loops; i++){
sum = static_cast<long int>(formula(n));
}
double sec2 = toc();
printf("Formula funtion: result = %ld | time = %f milliseconds\n", sum, sec2*1000);
}
}
#ifdef __GNUC__
#pragma GCC push_options
#pragma GCC optimize("O1")
#endif
static void task_d() {
printf("\n\n-------------- Task D --------------\n\n");
int const NLOOPS = 25; // chose a value such that the benchmark runs at least 10 sec.
unsigned int N = 50000001;
//##########################################################################
// Memory allocation
cout << "Memory allocation\n";
vector<double> x(N), y(N);
cout.precision(2);
cout << 2.0 * N *sizeof(x[0]) / 1024 / 1024 / 1024 << " GByte Memory allocated\n";
cout.precision(6);
//##########################################################################
// Data initialization
// Special: x_i = i+1; y_i = 1/x_i ==> <x,y> == N
for (unsigned int i = 0; i < N; ++i)
{
x[i] = i + 1;
y[i] = 1.0 / pow(x[i], 2);
}
//##########################################################################
cout << "\nStart Benchmarking Normal sum\n";
// Do calculation
auto t1 = system_clock::now(); // start timer
double sk1(0.0),ss(0.0);
for (int i = 0; i < NLOOPS; ++i)
{
sk1 = normal_sum(y);
ss += sk1; // prevents the optimizer from removing unused calculation results.
}
auto t2 = system_clock::now(); // stop timer
auto duration = duration_cast<microseconds>(t2 - t1); // duration in microseconds
double t_diff = static_cast<double>(duration.count()) / 1e6; // overall duration in seconds
t_diff = t_diff/NLOOPS;
// Print result
printf("\nSum = %.16f\n", sk1);
//##########################################################################
// Timings and Performance
cout << endl;
cout.precision(2);
cout << "Timing in sec. : " << t_diff << endl;
cout << "GFLOPS : " << 2.0 * N / t_diff / 1024 / 1024 / 1024 << endl;
cout << "GiByte/s : " << 2.0 * N / t_diff / 1024 / 1024 / 1024 * sizeof(x[0]) << endl;
//##########################################################################
cout << "\nStart Benchmarking Kahan summation\n";
// Do calculation
t1 = system_clock::now(); // start timer
double sk2(0.0),sss(0.0);
for (int i = 0; i < NLOOPS; ++i)
{
sk2 = Kahan_skalar(y);
sss += sk2; // prevents the optimizer from removing unused calculation results.
}
t2 = system_clock::now(); // stop timer
duration = duration_cast<microseconds>(t2 - t1); // duration in microseconds
t_diff = static_cast<double>(duration.count()) / 1e6; // overall duration in seconds
t_diff = t_diff/NLOOPS; // duration per loop seconds
// duration per loop seconds
// Print result
printf("\nSum = %.16f\n", sk2);
//##########################################################################
// Timings and Performance
cout << endl;
cout.precision(2);
cout << "Timing in sec. : " << t_diff << endl;
cout << "GFLOPS : " << 2.0 * N / t_diff / 1024 / 1024 / 1024 << endl;
cout << "GiByte/s : " << 2.0 * N / t_diff / 1024 / 1024 / 1024 * sizeof(x[0]) << endl;
//##########################################################################
// Print limit
printf("\nLimit = %.16f\n\n", pow(M_PI,2) / 6.0f);
}
#ifdef __GNUC__
#pragma GCC pop_options
#endif
static void task_e() {
printf("\n\n-------------- Task E --------------\n\n");
for (int n : {100, 1000, 10000}) {
vector<int> vec(n);
list<int> lst(n);
// Initialize
for (int i = 1; i < n+1; ++i) {
vec.push_back(i);
lst.push_back(i);
}
// Insert into vector
tic();
insert_into_vector(vec, n);
double sec1 = toc();
printf("Vector insertion time for n = %d: %.f microseconds.\n", n, sec1*1000*1000);
// Insert into list
tic();
insert_into_list(lst, n);
double sec2 = toc();
printf("List insertion time for n = %d: %.f microseconds.\n", n, sec2*1000*1000);
}
}
static void task_f() {
printf("\n\n-------------- Task F --------------\n\n");
// single_goldbach(k)
int k = 694;
printf("single_goldbach(k = %d) = %d\n", k, single_goldbach(k));
// Prints decompositions
print_decomps(k);
// count_goldbach(n)
printf("\nNOTE: For n=2'000'000 it will take ~30 seconds.\n");
for (int n : {10'000, 100'000, 400'000, 1'000'000}) { //, 2'000'000}) {
tic();
vector<int> counts = count_goldbach(n);
double sec = toc();
auto max = max_element(counts.begin(), counts.end());
printf("count_goldbach(n = %d): k = %ld, decompositions = %d, time elapsed: %f milliseconds\n", n, max-counts.begin(), *max, sec*1000);
}
// Results
// count_goldbach(n = 10'000): k = 9240, decompositions = 329, time elapsed: 1.235096 milliseconds
// count_goldbach(n = 100'000): k = 99330, decompositions = 2168, time elapsed: 39.003922 milliseconds
// count_goldbach(n = 400'000): k = 390390, decompositions = 7094, time elapsed: 497.282572 milliseconds
// count_goldbach(n = 1'000'000): k = 990990, decompositions = 15594, time elapsed: 3236.044944 milliseconds
// count_goldbach(n = 2'000'000): k = 1981980, decompositions = 27988, time elapsed: 29864.384370 milliseconds
// count_goldbach(n = 10'000'000): k = 9699690, decompositions = 124180, time elapsed: 825392.110981 milliseconds
}
static void task_g() {
printf("\n\n-------------- Task G --------------\n\n");
DenseMatrix const M(5,3);
vector<double> const u{{1,2,3}};
vector<double> f1 = M.Mult(u);
vector<double> const v{{-1,2,-3,4,-5}};
vector<double> f2 = M.MultT(v);
cout << "M = " << endl;
M.print();
cout << endl << "u = ";
for (size_t i=0; i<u.size(); i++){ cout << u[i] << " ";}
cout << endl << "M * u = ";
for (size_t i=0; i<f1.size(); i++){cout << f1[i] << " ";}
cout << endl << "v = ";
for (size_t i=0; i<v.size(); i++){cout << v[i] << " ";}
cout << endl << "M^T * v = ";
for (size_t i=0; i<f2.size(); i++){cout << f2[i] << " ";}
cout << endl << endl;
// #######################################################
int const NLOOPS = 100;
int const n = 3000;
// Time initialization
tic();
DenseMatrix const M2(n,n);
vector<double> w(n,0);
for (int i=0; i<n; i++){w[i]=i+1;}
double t1 = toc();
// Time Mult
tic();
vector<double> f3 = M2.Mult(w);
for (size_t k=1; k<NLOOPS; ++k){
f3 = M2.Mult(w);
}
double t2 = toc();
// Time MultT
tic();
vector<double> f4 = M2.MultT(w);
for (size_t k=1; k<NLOOPS; ++k){
f4 = M2.MultT(w);
}
double t3 = toc();
// Print results
printf("Results for %dx%d matrix vector multiplication doing %d loops\n", n, n, NLOOPS);
printf("Time for initialization: %f seconds.\n", t1);
printf("Time for Mult : %f seconds, %f per loop.\n", t2, t2/NLOOPS);
printf("Time for MultT : %f seconds, %f per loop.\n", t3, t3/NLOOPS);
// Check if resulting vectors are equal
for (size_t i = 0; i < n; i++)
{
double err = f3[i] - f4[i];
if(abs(err) > 1e-4)
{
cout << "Resulting vectors are not equal" << endl;
}
}
// ################################################
// Time initialization
tic();
vector<double> x(n,0);
for (int i=0; i<n; i++){x[i] = sigmoid( (10.0*static_cast<double>(i))/(n-1) - 5 );}
DenseMatrix2 M3(x,x);
double t4 = toc();
// Time Mult
tic();
vector<double> f5 = M3.Mult(w);
for (size_t k=1; k<NLOOPS; ++k){
f5 = M3.Mult(w);
}
double t5 = toc();
// Time MultT
tic();
vector<double> f6 = M3.MultT(w);
for (size_t k=1; k<NLOOPS; ++k){
f6 = M3.MultT(w);
}
double t6 = toc();
// Print results
printf("\nResults for %dx%d matrix vector multiplication doing %d loops taking advantage of tensor product structure of the matrix\n", n, n, NLOOPS);
printf("Time for initialization: %f seconds.\n", t4);
printf("Time for Mult : %f seconds, %f per loop.\n", t5, t5/NLOOPS);
printf("Time for MultT : %f seconds, %f per loop.\n", t6, t6/NLOOPS);
// Check if resulting vectors are equal
for (size_t i = 0; i < n; i++)
{
double err = f5[i] - f6[i];
if(abs(err) > 1e-4)
{
cout << "Resulting vectors are not equal" << endl;
}
}
}
int main(){
task_a();
task_b();
task_c();
task_d();
task_e();
task_f();
task_g();
return 0;
}

73
ex1/code/mayer_primes.h Normal file
View file

@ -0,0 +1,73 @@
#pragma once
#include <cstring> //memset
#include <vector>
//using namespace std;
/** \brief Determines all prime numbers in interval [2, @p max].
*
* The sieve of Eratosthenes is used.
*
* The implementation originates from <a href="http://code.activestate.com/recipes/576559-fast-prime-generator/">Florian Mayer</a>.
*
* \param[in] max end of interval for the prime number search.
* \return vector of prime numbers @f$2,3,5, ..., p<=max @f$.
*
* \copyright
* Copyright (c) 2008 Florian Mayer (adapted by Gundolf Haase 2018)
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to deal
* in the Software without restriction, including without limitation the rights
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in
* all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
* INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*
*/
template <class T>
std::vector<T> get_primes(T max)
{
std::vector<T> primes;
char *sieve;
sieve = new char[max / 8 + 1];
// Fill sieve with 1
memset(sieve, 0xFF, (max / 8 + 1) * sizeof(char));
for (T x = 2; x <= max; x++)
{
if (sieve[x / 8] & (0x01 << (x % 8))) {
primes.push_back(x);
// Is prime. Mark multiplicates.
for (T j = 2 * x; j <= max; j += x)
{
sieve[j / 8] &= ~(0x01 << (j % 8));
}
}
}
delete[] sieve;
return primes;
}
//---------------------------------------------------------------
//int main() // by Florian Mayer
//{g++ -O3 -std=c++14 -fopenmp main.cpp && ./a.out
// vector<unsigned long> primes;
// primes = get_primes(10000000);
// // return 0;
// // Print out result.
// vector<unsigned long>::iterator it;
// for(it=primes.begin(); it < primes.end(); it++)
// cout << *it << " ";
//
// cout << endl;
// return 0;
//}

354
ex1/code/mylib.cpp Normal file
View file

@ -0,0 +1,354 @@
#include "mylib.h"
#include "mayer_primes.h"
#include <cassert> // assert
#include <vector>
#include <iostream>
#include <cmath>
#include <tuple>
#include <string>
#include <algorithm>
#include <fstream>
#include <stdexcept>
#include <random>
#include <list>
using namespace std;
// -------------- Task A --------------
tuple<double, double, double> means0(double a, double b, double c){
double arith = (a+b+c) / 3;
double geo = pow((a*b*c), 1.0/3);
double harm = 3 / ((1.0/a) + (1.0/b) + (1.0/c));
return make_tuple(arith, geo, harm);
}
tuple<double, double, double> means(const vector<double>& v){
size_t n = v.size();
double sum = 0;
double logsum = 0;
double invsum = 0;
for (size_t i = 0; i<n; ++i){
sum += v[i];
logsum += log(v[i]);
invsum += 1.0/v[i];
}
double arith = sum / static_cast<double>(n);
double geo = exp(1.0/static_cast<double>(n) * logsum);
double harm = static_cast<double>(n) / invsum;
return make_tuple(arith, geo, harm);
}
// -------------- Task B --------------
void fill_vector(istream& istr, vector<double>& v)
{
double d=0;
while ( istr >> d) v.push_back(d); // Einlesen
if (!istr.eof())
{ // Fehlerbehandlung
cout << " Error handling \n";
if ( istr.bad() ) throw runtime_error("Schwerer Fehler in istr");
if ( istr.fail() ) // Versuch des Aufraeumens
{
cout << " Failed in reading all data.\n";
istr.clear();
}
}
v.shrink_to_fit(); // C++11
return;
}
void read_vector_from_file(const string& file_name, vector<double>& v)
{
ifstream fin(file_name); // Oeffne das File im ASCII-Modus
if( fin.is_open() ) // File gefunden:
{
v.clear(); // Vektor leeren
fill_vector(fin, v);
}
else // File nicht gefunden:
{
cout << "\nFile " << file_name << " has not been found.\n\n" ;
assert( fin.is_open() && "File not found." ); // exeption handling for the poor programmer
}
return;
}
void write_vector_to_file(const string& file_name, const vector<double>& v)
{
ofstream fout(file_name); // Oeffne das File im ASCII-Modus
if( fout.is_open() )
{
for (size_t k=0; k<v.size(); ++k)
{
fout << v.at(k) << endl;
}
}
else
{
cout << "\nFile " << file_name << " has not been opened.\n\n" ;
assert( fout.is_open() && "File not opened." ); // exeption handling for the poor programmer
}
return;
}
// -------------- Task C --------------
long int sum_of_spec(int n)
{
long int sum = 0;
for (long int i=1; i<n+1; i++){
if (i % 3 == 0 || i % 5 == 0){
sum += i;
}
}
return sum;
}
double formula(int n)
{
double div_by_3 = floor( n / 3.0 );
double div_by_5 = floor( n / 5.0 );
double div_by_15 = floor( n / 15.0 );
double S_3 = 3.0 * (div_by_3*((div_by_3+1)/2.0));
double S_5 = 5.0 * (div_by_5*((div_by_5+1)/2.0));
double S_15 = 15.0 * (div_by_15*((div_by_15+1)/2.0));
double sum = S_3 + S_5 - S_15;
return sum;
}
// -------------- Task D --------------
#ifdef __GNUC__
#pragma GCC push_options
#pragma GCC optimize("O1")
#endif
long double Kahan_skalar(vector<double> const &input)
{
long double sum = 0.0;
long double c = 0.0;
for (long unsigned int i=0; i<input.size(); i++){
long double y = input[i] - c;
long double t = sum + y;
c = (t-sum) - y;
sum = t;
}
return sum;
}
long double normal_sum(vector<double> const &input)
{
long double sum = 0.0;
for (long unsigned int i=0; i<input.size(); i++){
sum += input[i];
}
return sum;
}
#ifdef __GNUC__
#pragma GCC pop_options
#endif
// -------------- Task E --------------
void insert_into_vector(vector<int>& vec, int n) {
random_device rd; // random device
mt19937 gen(rd()); // seed
uniform_int_distribution<> dist(1, n); // define range
for (int i = 0; i < n; ++i) {
int rand_num = dist(gen);
auto pos = lower_bound(vec.begin(), vec.end(), rand_num);
vec.insert(pos, rand_num);
}
assert(is_sorted(vec.begin(), vec.end()));
}
void insert_into_list(list<int>& lst, int n) {
random_device rd;
mt19937 gen(rd());
uniform_int_distribution<> dist(1, n);
for (int i = 0; i < n; ++i) {
int rand_num = dist(gen);
auto pos = lower_bound(lst.begin(), lst.end(), rand_num);
lst.insert(pos, rand_num);
}
assert(is_sorted(lst.begin(), lst.end()));
}
// -------------- Task F --------------
int single_goldbach(int k) {
const vector<int> primes = get_primes(k);
int count = 0;
for (size_t i = 0; i < primes.size(); i++) {
for (size_t j = i; j < primes.size(); j++) {
if (primes[i] + primes[j] == k) {
count++;
}
}
}
return count;
}
vector<int> count_goldbach(int n) {
const vector<int> primes = get_primes(n);
vector<int> counts(n+1);
for (size_t i = 1; i < primes.size(); i++) {
for (size_t j = i; j < primes.size(); j++) {
int sum = primes[i] + primes[j];
if (sum <= n) {
counts[sum]++;
}
}
}
return counts;
}
void print_decomps(int k) {
const vector<int> primes = get_primes(k);
cout << "\nDecompositions for k = " << k << ": ";
for (size_t i = 0; i < primes.size(); i++) {
for (size_t j = i; j < primes.size(); j++) {
if (primes[i] + primes[j] == k) {
cout << primes[i] << " + " << primes[j] << ", ";
}
}
}
cout << endl;
}
// -------------- Task G --------------
double sigmoid(double x)
{
return 1.0 / (1.0 + std::exp(-x));
}
double DenseMatrix::sigmoid(double x)
{
return 1.0 / (1.0 + exp(-x));
}
DenseMatrix::DenseMatrix(size_t n, size_t m) : rows(n), cols(m), matrix(n*m)
{
size_t nm = max(n, m);
for (size_t i = 0; i < rows; i++)
{
for (size_t j = 0; j < cols; j++)
{
double x_i = 10.0*static_cast<double>(i)/(nm-1) - 5.0;
double x_j = 10.0*static_cast<double>(j)/(nm-1) - 5.0;
matrix[i*cols + j] = sigmoid(x_i) * sigmoid(x_j);
}
}
}
vector<double> DenseMatrix::Mult(const vector<double>& vec) const
{
assert(vec.size() == cols);
vector<double> result(rows, 0.0);
for (size_t i = 0; i < rows; ++i)
{
for (size_t j = 0; j < cols; ++j)
{
result[i] += matrix[i*cols + j] * vec[j];
}
}
return result;
}
vector<double> DenseMatrix::MultT(const vector<double>& vec) const
{
assert(vec.size() == rows);
vector<double> result(cols, 0.0);
for (size_t i = 0; i < cols; ++i)
{
for (size_t j = 0; j < rows; ++j)
{
result[i] += matrix[j*cols + i] * vec[j];
}
}
return result;
}
void DenseMatrix::print() const {
size_t count(0);
cout.precision(5);
for (double val : matrix) {
printf("%.6f ", val);
count++;
if(count%cols == 0){cout << endl;}
}
}
// #################################################
DenseMatrix2::DenseMatrix2(vector<double> const u, vector<double> const v) : rows(u.size()), cols(v.size()), u_(u), v_(v) {}
vector<double> DenseMatrix2::Mult(const vector<double>& vec) const
{
assert(vec.size() == cols);
vector<double> result(rows, 0.0);
double scalar(0.0);
for (size_t i = 0; i < vec.size(); i++)
{
scalar += v_[i] * vec[i];
}
for (size_t i = 0; i < u_.size(); i++)
{
result[i] = scalar*u_[i];
}
return result;
}
vector<double> DenseMatrix2::MultT(const vector<double>& vec) const
{
assert(vec.size() == rows);
vector<double> result(cols, 0.0);
double scalar(0.0);
for (size_t i = 0; i < u_.size(); i++)
{
scalar += u_[i] * vec[i];
}
for (size_t i = 0; i < v_.size(); i++)
{
result[i] = scalar*v_[i];
}
return result;
}
void task_a();
void task_b();
void task_c();
void task_d();
void task_e();
void task_f();
void task_g();

118
ex1/code/mylib.h Normal file
View file

@ -0,0 +1,118 @@
#pragma once
#include <cassert> // assert
#include <vector>
#include <iostream>
#include <cmath>
#include <tuple>
#include <string>
#include <algorithm>
#include <fstream>
#include <stdexcept>
#include <list>
using namespace std;
// -------------- Task A --------------
// Returns arithmetic, geometric and harmonic mean for 3 values a,b,c.
tuple<double, double, double> means0(double a, double b, double c);
// Returns arithmetic, geometric and harmonic mean for a vector.
tuple<double, double, double> means(const vector<double>& v);
// -------------- Task B --------------
/**
This function opens the ASCII-file named @p file_name and reads the
double data into the C++ vector @p v.
If the file @p file_name does not exist then the code stops with an appropriate message.
@param[in] file_name name of the ASCII-file
@param[out] v C++ vector with double values
*/
void read_vector_from_file(const string& file_name, vector<double>& v);
/**
This function opens the ASCII-file named @p file_name and rewrites its with the
double data from the C++ vector @p v.
If there are problems in opening/generating file @p file_name
then the code stops with an appropriate message.
@param[in] file_name name of the ASCII-file
@param[in] v C++ vector with double values
*/
void write_vector_to_file(const string& file_name, const vector<double>& v);
/**
Fills the double-vector @p v with data from an input stream @p istr until this input stream
ends regularily. The vector is cleared and its memory is automatically allocated.
@param[in] istr input stream
@param[out] v C++ vector with double values
@warning An exception is thrown in case of wrong data format or corrupted data.
*/
void fill_vector(istream& istr, vector<double>& v);
// -------------- Task C --------------
// Sums up all positive integers less or equal n which are multiples of 3 or of 5 (including or!) by brute force.
long int sum_of_spec(int n);
// Sums up all positive integers less or equal n which are multiples of 3 or of 5 (including or!) by inclusion-exclusion principle.
double formula(int n);
// -------------- Task D --------------
long double Kahan_skalar(std::vector<double> const &input);
long double normal_sum(std::vector<double> const &input);
// -------------- Task E --------------
// Inserts n random numbers into sorted vector v such that v remains sorted.
void insert_into_vector(vector<int>& vec, int n);
// Inserts n random numbers into sorted list such that the list remains sorted.
void insert_into_list(list<int>& lst, int n);
// -------------- Task F --------------
// Counts number of possible decompositions with 2 primes that sum up to k.
int single_goldbach(int k);
// Counts number of possible decompositions with 2 primes that sum up to k for all even numbers k \in {4,...,n}.
vector<int> count_goldbach(int n);
// Prints all decompositions of k.
void print_decomps(int k);
// -------------- Task G --------------
// Sigmoid function 1/(1+exp(-x))
double sigmoid(double x);
class DenseMatrix {
private:
double sigmoid(double x);
size_t rows;
size_t cols;
vector<double> matrix;
public:
DenseMatrix(size_t n, size_t m); // Constructor
vector<double> Mult(const vector<double>& vec) const;
vector<double> MultT(const vector<double>& vec) const;
void print() const;
};
class DenseMatrix2 {
private:
size_t rows;
size_t cols;
vector<double> u_;
vector<double> v_;
public:
DenseMatrix2(vector<double> const u, vector<double> const v); // Constructor
vector<double> Mult(const vector<double>& vec) const;
vector<double> MultT(const vector<double>& vec) const;
};

51
ex1/code/timing.h Normal file
View file

@ -0,0 +1,51 @@
//
// Gundolf Haase, Oct 18 2024
//
#pragma once
#include <chrono> // timing
#include <stack>
//using Clock = std::chrono::system_clock; //!< The wall clock timer chosen
using Clock = std::chrono::high_resolution_clock;
using TPoint= std::chrono::time_point<Clock>;
// [Galowicz, C++17 STL Cookbook, p. 29]
inline
std::stack<TPoint> MyStopWatch; //!< starting time of stopwatch
/** Starts stopwatch timer.
* Use as @code tic(); myfunction(...) ; double tsec = toc(); @endcode
*
* The timining can be nested and the recent time point is stored on top of the stack.
*
* @return recent time point
* @see toc
*/
inline auto tic()
{
MyStopWatch.push(Clock::now());
return MyStopWatch.top();
}
/** Returns the elapsed time from stopwatch.
*
* The time point from top of the stack is used
* if time point @p t_b is not passed as input parameter.
* Use as @code tic(); myfunction(...) ; double tsec = toc(); @endcode
* or as @code auto t_b = tic(); myfunction(...) ; double tsec = toc(t_b); @endcode
* The last option is to be used in the case of
* non-nested but overlapping time measurements.
*
* @param[in] t_b start time of some stop watch
* @return elapsed time in seconds.
*
*/
inline double toc(TPoint const &t_b = MyStopWatch.top())
{
// https://en.cppreference.com/w/cpp/chrono/treat_as_floating_point
using Unit = std::chrono::seconds;
using FpSeconds = std::chrono::duration<double, Unit::period>;
auto t_e = Clock::now();
MyStopWatch.pop();
return FpSeconds(t_e-t_b).count();
}