restructurization

_snippets/B-trees.c

@@ -0,0 +1,317 @@
+/*
+ * B-trees.c - my realisation of binary serch trees
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <err.h>
+#include <assert.h>
+#include <stdint.h>
+
+#include "B-trees.h"
+
+#ifndef Malloc
+#define Malloc my_alloc
+/**
+ * Memory allocation with control
+ * @param N - number of elements
+ * @param S - size of one element
+ * @return allocated memory or die
+ */
+void *my_alloc(size_t N, size_t S){
+    void *p = calloc(N, S);
+    if(!p) err(1, "malloc");
+    return p;
+}
+#endif
+#ifndef FREE
+#define FREE(arg) do{free(arg); arg = NULL;}while(0)
+#endif
+
+#ifdef EBUG
+#ifndef DBG
+	#define DBG(...) do{fprintf(stderr, __VA_ARGS__);}while(0)
+#endif
+#else
+#ifndef DBG
+	#define DBG(...)
+#endif
+#endif
+
+static inline uint32_t log_2(const uint32_t x) {
+    if(x == 0) return 0;
+    return (31 - __builtin_clz (x));
+}
+
+// later this function maybe mode difficult
+void bt_destroy(BTree *node){
+	FREE(node);
+}
+
+// get tree leaf nearest to v
+BTree *bt_get(BTree *root, int v){
+	if(!root) return NULL;
+	int D = root->data;
+	if(D == v) return root;
+	if(root->data > v) return root->left;
+	else return root->right;
+}
+
+// find leaf with v or nearest to it (prev)
+// prev is nearest parent or NULL for root
+BTree *bt_search(BTree *root, int v, BTree **prev){
+	if(!root){
+		if(prev) *prev = NULL;
+		return NULL;
+	}
+	BTree *p;
+	do{
+		p = root;
+		root = bt_get(root, v);
+	}while(root && root->data != v);
+	if(prev){
+		if(p != root) *prev = p;
+		else *prev = NULL;
+	}
+	return root;
+}
+
+BTree *bt_create_leaf(int v){
+	BTree *node = Malloc(1, sizeof(BTree));
+	if(!node) return NULL;
+	node->data = v;
+	return node;
+}
+
+// insert new leaf (or leave it as is, if exists)
+// return new node
+BTree *bt_insert(BTree **root, int v){
+	BTree *closest = NULL, *node = NULL;
+	if(root && *root){
+		node = bt_search(*root, v, &closest);
+		if(node) return node; // value exists
+	}
+	node = bt_create_leaf(v);
+	if(!node) return NULL;
+	if(!closest){ // there's no values at all!
+		if(root) *root = node; // modify root node
+		return node;
+	}
+	// we have a closest node to our data, so create node and insert it:
+	int D = closest->data;
+	if(D > v) // there's no left leaf of closest node, add our node there
+		closest->left = node;
+	else
+		closest->right = node;
+	return node;
+}
+
+// find mostly right element or NULL if no right children
+BTree *max_leaf(BTree *root, BTree **prev){
+	if(!root)
+		return NULL;
+	BTree *ret = root->left;
+	if(prev) *prev = root;
+	while(ret->right){
+		if(prev) *prev = ret;
+		ret = ret->right;
+	}
+	return ret;
+}
+// find mostly left element or NULL if no left children
+BTree *min_leaf(BTree *root, BTree **prev){
+	if(!root) return NULL;
+	BTree *ret = root->right;
+	if(prev) *prev = root;
+	while(ret->left){
+		if(prev) *prev = ret;
+		ret = ret->left;
+	}
+	return ret;
+}
+
+void print_tree(BTree *leaf){
+	if(!leaf){
+		DBG("null");
+		return;
+	}
+	DBG("(%d: ", leaf->data);
+	print_tree(leaf->left);
+	DBG(", ");
+	print_tree(leaf->right);
+	DBG(")");
+}
+
+// remove element, correct root if it is removed
+void bt_remove(BTree **root, int v){
+	assert(root);
+	void conn(BTree *p, BTree *n, BTree *node){
+	DBG("connect: %d and %d through %d\n", p ? p->data : -1,
+			n? n->data : -1, node->data);
+		if(p->left == node) p->left = n;
+		else p->right = n;
+	}
+	void conn_node(BTree *prev, BTree *next, BTree *node){
+		if(!prev) *root = next; // this is a root
+		else conn(prev, next, node);
+		bt_destroy(node);
+	}
+	void node_subst(BTree *nold, BTree *nnew, BTree *pold, BTree *pnew){
+		DBG("Substitute: %d by %d\n", nold-> data, nnew->data);
+		if(!pold) *root = nnew; // this is a root
+		else conn(pold, nnew, nold); // connect parent of old element to new
+		conn(pnew, NULL, nnew); // remove node from it's parent
+		nnew->left = nold->left;
+		nnew->right = nold->right;
+		bt_destroy(nold);
+	}
+	BTree *node, *prev, *Lmax, *Rmin, *Lmp, *Rmp;
+	int L, R;
+	node = bt_search(*root, v, &prev);
+	if(!node) return; // there's no node v
+	if(!node->left){ // there's only right child or none at all
+		conn_node(prev, node->right, node);
+		return;
+	}
+	if(!node->right){ // the same but like mirror
+		conn_node(prev, node->left, node);
+		return;
+	}
+	// we have situation that element have both children
+	// 1) find max from left child and min from right
+	Lmax = max_leaf(node, &Lmp);
+	L = Lmax->data;
+	Rmin = min_leaf(node, &Rmp);
+	R = Rmin->data;
+	// 2) now L is left max, R - right min
+	// we select the nearest to v
+	if(v - L > R - v){ // R is closest value, substitute node by Rmin
+		node_subst(node, Rmin, prev, Rmp);
+	}else{ // substitute node by Lmax
+		node_subst(node, Lmax, prev, Lmp);
+	}
+}
+
+// return node need rotation
+BTree *rotate_ccw(BTree *node, BTree *prev){
+	if(!node) return NULL;
+	#ifdef EBUG
+	DBG("\n\ntry to rotate near %d", node->data);
+	if(prev) DBG(", right element: %d", prev->data);
+	if(!node->right) DBG("\n");
+	#endif
+	if(!node->right) return node->left; // no right child -> return left
+	BTree *chld = node->right;
+	DBG(", child: %d\n", chld->data);
+	if(prev) prev->left = chld;
+	node->right = chld->left;
+	chld->left = node;
+	if(chld->right){DBG("R:%d\n",chld->right->data); return chld;} // need another rotation
+	else return node;
+}
+
+void tree_to_vine(BTree **root){
+	assert(root); assert(*root);
+	BTree *node = *root, *prev;
+	while((*root)->right) *root = (*root)->right; // max element will be a new root
+	while(node && node->right)
+		node = rotate_ccw(node, NULL); // go to new root
+	prev = node; node = prev->left;
+	while(rotate_ccw(node, prev)){
+		node = prev->left;
+		if(!node->right){
+			prev = node;
+			node = node->left;
+		}
+	}
+}
+
+// return node->left
+BTree *rotate_cw(BTree *node, BTree *parent){ // rotate cw near node->left
+	if(!node) return NULL;
+	DBG("\n\ntry to rotate cw near %d\n", node->data);
+	#ifdef EBUG
+	if(parent)DBG("\t\tparent: %d\n", parent->data);
+	#endif
+	if(!node->left) return NULL;
+	BTree *chld = node->left, *rght = chld->right;
+	if(parent) parent->left = chld;
+	DBG("\t\t, child: %d\n", chld->data);
+	chld->right = node;
+	node->left = rght;
+	return chld;
+}
+
+// make balanced tree
+void vine_to_tree(BTree **root){
+	uint32_t depth, i;
+	BTree *node = *root, *prev;
+	for(depth = 0; node; depth++, node = node->left);
+	DBG("depth: %u; ", depth);
+	depth = log_2(depth);
+	DBG("log: %d\n", depth);
+	for(i = 0; i < depth; i++){
+		DBG("\nIteration %d\n", i);
+		node = *root;
+		prev = NULL;
+		if((*root)->left) *root = (*root)->left;
+		while(node && node->left){
+			node = rotate_cw(node, prev);
+			if(!node) break;
+			prev = node;
+			node = node->left;
+		}
+	}
+}
+
+#ifdef STANDALONE
+#define INS(X)  do{if(!bt_insert(&root, X)) err(1, "alloc error");}while(0)
+int main(void){
+	BTree *tp, *root = NULL;
+	int i, ins[] = {7,3,12,1,5,9,20,0,4,6,13,21};
+	void search_bt(){
+		printf("\nRoot: %d\nTree", root->data);
+		print_tree(root); printf("\n\n");
+		for(i = 0; i < 22; i++){
+			tp = bt_search(root, i, NULL);
+			printf("%d ", i);
+			if(!tp) printf("not ");
+			printf("found\n");
+		}
+	}
+	for(i = 0; i < 12; i++)
+		INS(ins[i]);
+
+	search_bt();
+	printf("now we delete leafs 4, 12 and 7 (old root!)\n");
+	bt_remove(&root, 4); bt_remove(&root, 12); bt_remove(&root, 7);
+	search_bt();
+	printf("add items 2, 4, 19\n");
+	INS(2); INS(4); INS(19);
+	search_bt();
+	printf("\nnow make a vine\n");
+	tree_to_vine(&root);
+	search_bt();
+	printf("\nAnd balance tree\n");
+	vine_to_tree(&root);
+	search_bt();
+	return 0;
+}
+#endif

_snippets/B-trees.h

@@ -0,0 +1,41 @@
+/*
+ * B-trees.h
+ *
+ * Copyright 2013 Edward V. Emelianoff <eddy@sao.ru>
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#pragma once
+#ifndef __B_TREES_H__
+#define __B_TREES_H__
+
+typedef struct btree_node{
+	int data;
+	struct btree_node *left, *right;
+} BTree;
+
+BTree *bt_get(BTree *root, int v);
+BTree *bt_search(BTree *root, int v, BTree **prev);
+BTree *bt_create_leaf(int v);
+BTree *bt_insert(BTree **root, int v);
+BTree *max_leaf(BTree *root, BTree **prev);
+BTree *min_leaf(BTree *root, BTree **prev);
+void bt_remove(BTree **root, int v);
+
+void print_tree(BTree *leaf);
+
+#endif // __B_TREES_H__

_snippets/Bresenham_circle.c

@@ -0,0 +1,60 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+char *drawCircle(int R, int geom){
+	if(R > 200 || R < 1) return NULL;
+	int S, i, Y = 2 * R + 2;
+	if(geom)
+		S = Y * (R + 1);
+	else
+		S = Y * (Y - 1);
+	char *buf = malloc(S+1);
+	if(!buf) return NULL;
+	memset(buf, ' ', S);
+	buf[S] = 0;
+	for(i = Y-1; i < S; i+=Y) buf[i] = '\n';
+	inline void DrawPixel(int x, int y){
+		if(geom){
+			if(y%2==0) buf[Y * y/2 + x] = '*';
+		}else{
+			buf[Y * y + x] = '*';
+		}
+	}
+	// Bresenham's circle algorithm
+	int x = R;
+	int y = 0;
+	int radiusError = 1-x;
+	while(x >= y){
+		DrawPixel(x + R,   y + R);
+		DrawPixel(y + R,   x + R);
+		DrawPixel(-x + R,  y + R);
+		DrawPixel(-y + R,  x + R);
+		DrawPixel(-x + R, -y + R);
+		DrawPixel(-y + R, -x + R);
+		DrawPixel(x + R,  -y + R);
+		DrawPixel(y + R,  -x + R);
+		y++;
+		if (radiusError < 0){
+			radiusError += 2 * y + 1;
+		}else{
+			x--;
+			radiusError += 2 * (y - x) + 1;
+		}
+	}
+	return buf;
+}
+
+int main(int argc, char **argv){
+	int i, R;
+	char *buf;
+	for(i = 1; i < argc; i++){
+		if(!(buf = drawCircle(R = atoi(argv[i]), 1))){
+			printf("Wrong parameter %s\n", argv[i]);
+			continue;
+		}
+		printf("\nCircle with R = %d:\n%s\n", R, buf);
+		free(buf); buf = NULL;
+	}
+	return 0;
+}

_snippets/C-functions_overload.c

@@ -0,0 +1,20 @@
+#include <stdio.h>
+#include <stdint.h>
+
+#define FUNC(arg) _Generic(arg, uint16_t: funcu, int32_t: funci)(arg)
+
+void funcu(uint16_t arg){
+    printf("uint16_t: %u\n", arg);
+}
+
+void funci(int32_t arg){
+    printf("int32_t: %d\n", arg);
+}
+
+int main(){
+    uint16_t u = 32;
+    int32_t i = -50333;
+    FUNC(u);
+    FUNC(i);
+    return 0;
+}

_snippets/USB_reset.c

@@ -0,0 +1,39 @@
+// https://askubuntu.com/a/661/501233
+#include <stdio.h>
+#include <unistd.h>
+#include <fcntl.h>
+#include <errno.h>
+#include <sys/ioctl.h>
+
+#include <linux/usbdevice_fs.h>
+
+
+int main(int argc, char **argv)
+{
+    const char *filename;
+    int fd;
+    int rc;
+
+    if (argc != 2) {
+        fprintf(stderr, "Usage: usbreset device-filename\n");
+        return 1;
+    }
+    filename = argv[1];
+
+    fd = open(filename, O_WRONLY);
+    if (fd < 0) {
+        perror("Error opening output file");
+        return 1;
+    }
+
+    printf("Resetting USB device %s\n", filename);
+    rc = ioctl(fd, USBDEVFS_RESET, 0);
+    if (rc < 0) {
+        perror("Error in ioctl");
+        return 1;
+    }
+    printf("Reset successful\n");
+
+    close(fd);
+    return 0;
+}

_snippets/avx/Makefile

@@ -0,0 +1,9 @@
+CC=gcc
+CFLAGS= -march=native -O3
+
+all: add dotproduct
+
+%: %.c
+	@echo -e "\t\tCC $<"
+	$(CC) $(CFLAGS) -o $@ $<
+

_snippets/avx/add.c

@@ -0,0 +1,32 @@
+/* Construct a 256-bit vector from 4 64-bit doubles. Add it to itself
+ * and print the result.
+ */
+
+#include <stdio.h>
+#include <immintrin.h>
+
+int main() {
+
+  __m256i hello;
+  // Construction from scalars or literals.
+  __m256d a = _mm256_set_pd(1.0, 2.0, 3.0, 4.0);
+
+  // Does GCC generate the correct mov, or (better yet) elide the copy
+  // and pass two of the same register into the add? Let's look at the assembly.
+  __m256d b = _mm256_set_pd(0.0, 0.0, 0.0, 0.0), c;
+for(int i = 0; i < 1000000000; ++i){
+  // Add the two vectors, interpreting the bits as 4 double-precision
+  // floats.
+  c = _mm256_add_pd(a, b);
+  b = c;
+  
+ }
+
+  // Do we ever touch DRAM or will these four be registers?
+  __attribute__ ((aligned (32))) double output[4];
+  _mm256_store_pd(output, c);
+
+  printf("%f %f %f %f\n",
+         output[0], output[1], output[2], output[3]);
+  return 0;
+}

_snippets/avx/dotproduct.c

@@ -0,0 +1,61 @@
+/* Compute the dot product of two (properly aligned) vectors. */
+#include <stdio.h>
+#include <immintrin.h>
+#include <math.h>
+
+const int N = 83;
+
+double slow_dot_product(const double *a, const double *b) {
+    double answer = 0.0;
+    for(int ii = 0; ii < N; ++ii)
+        answer += a[ii]*b[ii];
+    return answer;
+}
+
+/* Horizontal add works within 128-bit lanes. Use scalar ops to add
+ * across the boundary. */
+double reduce_vector1(__m256d input) {
+    __m256d temp = _mm256_hadd_pd(input, input);
+    return ((double*)&temp)[0] + ((double*)&temp)[2];
+}
+
+/* Another way to get around the 128-bit boundary: grab the first 128
+ * bits, grab the lower 128 bits and then add them together with a 128
+ * bit add instruction. */
+double reduce_vector2(__m256d input) {
+    __m256d temp = _mm256_hadd_pd(input, input);
+    __m128d sum_high = _mm256_extractf128_pd(temp, 1);
+    __m128d result = _mm_add_pd(sum_high, _mm256_castpd256_pd128(temp));
+    return ((double*)&result)[0];
+}
+
+double dot_product(const double *a, const double *b) {
+    __m256d sum_vec = _mm256_set_pd(0.0, 0.0, 0.0, 0.0);
+
+    /* Add up partial dot-products in blocks of 256 bits */
+    for(int ii = 0; ii < N/4; ++ii) {
+        __m256d x = _mm256_load_pd(a+4*ii);
+        __m256d y = _mm256_load_pd(b+4*ii);
+        __m256d z = _mm256_mul_pd(x,y);
+        sum_vec = _mm256_add_pd(sum_vec, z);
+    }
+
+    /* Find the partial dot-product for the remaining elements after
+    * dealing with all 256-bit blocks. */
+    double final = 0.0;
+    for(int ii = N-N%4; ii < N; ++ii)
+        final += a[ii] * b[ii];
+
+    return reduce_vector2(sum_vec) + final;
+}
+
+int main() {
+    __attribute__ ((aligned (32))) double a[N], b[N];
+
+    for(int ii = 0; ii < N; ++ii)
+        a[ii] = b[ii] = ii/sqrt(N);
+
+    double answer = dot_product(a, b);
+    printf("%f\n", answer);
+    printf("%f\n", slow_dot_product(a,b));
+}

_snippets/calcAP/Makefile

@@ -0,0 +1,14 @@
+LDFLAGS = -lm -lsla -lsofa_c
+SRCS = $(wildcard *.c)
+CC = gcc
+DEFINES = -D_XOPEN_SOURCE=1111
+CFLAGS = -Wall -Werror -Wextra $(DEFINES)
+TARGS = $(SRCS:.c=)
+all : $(TARGS)
+slalib_and_sofa : slalib_and_sofa.c
+	$(CC) $(CFLAGS) $(LDFLAGS) $< -o $@
+slalib_sofa_nova : slalib_sofa_nova.c
+	$(CC) $(CFLAGS) $(LDFLAGS) -lnova -lerfa $< -o $@
+
+clean:
+	/bin/rm -f *.o *~

_snippets/calcAP/Readme

@@ -0,0 +1 @@
+Comparison of apparent place calculation in different libraries

_snippets/calcAP/slalib_and_sofa.c

@@ -0,0 +1,189 @@
+/*
+ * slalib_and_sofa.c - calculate apparent place by slalib & libsofa
+ *
+ * Copyright 2016 Edward V. Emelianov <eddy@sao.ru, edward.emelianoff@gmail.com>
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+
+#define _GNU_SOURCE 1111  // strcasecmp
+
+#include "sofa.h"
+#include <stdio.h>
+#include <stdbool.h>
+#include <time.h>
+
+#define DBG(...) printf(__VA_ARGS__)
+
+extern void sla_caldj(int*, int*, int*, double*, int*);
+extern void sla_amp(double*, double*, double*, double*, double*, double*);
+extern void sla_map(double*, double*, double*, double*, double*,double*, double*, double*, double*, double*);
+void slacaldj(int y, int m, int d, double *djm, int *j){
+	int iy = y, im = m, id = d;
+	sla_caldj(&iy, &im, &id, djm, j);
+}
+void slaamp(double ra, double da, double date, double eq, double *rm, double *dm ){
+	double r = ra, d = da, mjd = date, equi = eq;
+	sla_amp(&r, &d, &mjd, &equi, rm, dm);
+}
+// rm,dm - mean RA,Dec (rad), pr,pd - RA,Dec changes per Julian year (dRA/dt, dDec/dt)
+// px - parallax (arcsec), rv - radial speed (km/sec,  +ve if receding)
+// eq - epoch and equinox of star data (Julian)
+// date - TDB for apparent place (JD-2400000.5)
+void slamap(double rm, double dm, double pr, double pd,
+             double px, double rv, double eq, double date,
+             double *ra, double *da){
+	double r = rm, d = dm, p1 = pr, p2 = pd, ppx = px, prv = rv, equi = eq, dd = date;
+	sla_map(&r, &d, &p1, &p2, &ppx, &prv, &equi, &dd, ra, da);
+}
+
+void reprd(char* s, double ra, double dc){
+	char pm;
+	int i[4];
+	printf ( "%30s", s );
+	iauA2tf ( 7, ra, &pm, i );
+	printf ( " %2.2d %2.2d %2.2d.%7.7d", i[0],i[1],i[2],i[3] );
+	iauA2af ( 6, dc, &pm, i );
+	printf ( " %c%2.2d %2.2d %2.2d.%6.6d\n", pm, i[0],i[1],i[2],i[3] );
+}
+
+void radtodeg(double r){
+	int i[4]; char pm;
+	int rem = (int)(r / D2PI);
+	if(rem) r -= D2PI * rem;
+	if(r > DPI) r -= D2PI;
+	else if(r < -DPI) r += D2PI;
+	iauA2af (2, r, &pm, i);
+	printf("%c%02d %02d %02d.%2.d", pm, i[0],i[1],i[2],i[3]);
+}
+
+double getta(char *str){
+	int a,b,s = 1; double c;
+	if(3 != sscanf(str, "%d:%d:%lf", &a,&b,&c)) return -1;
+	if(a < 0){ s = -1; a = -a;}
+	c /= 3600.;
+	c += a + b/60.;
+	c *= s;
+	return c;
+}
+
+int main (int argc, char **argv){
+	double rc, dc;
+	if(argc == 3){
+		rc = getta(argv[1]) * DPI / 12;
+		dc = getta(argv[2]) * DD2R;
+	}else{
+		/* Star ICRS RA,Dec (radians). */
+		if ( iauTf2a ( ' ', 19, 50, 47.6, &rc ) ) return -1;
+		if ( iauAf2a ( '+', 8, 52, 12.3, &dc ) ) return -1;
+	}
+	reprd ( "ICRS, epoch J2000.0:", rc, dc );
+
+	struct tm tms;
+	time_t t = time(NULL);
+	gmtime_r(&t, &tms);
+	int y, m, d, err;
+	y = 1900 + tms.tm_year;
+	m = tms.tm_mon + 1;
+	d = tms.tm_mday;
+	double mjd, add = ((double)tms.tm_hour + (double)tms.tm_min/60.0 + tms.tm_sec/3600.0) / 24.;
+	DBG("Date: (d/m/y +frac)  %d/%d/%d +%g\n", d, m, y, add);
+	slacaldj(y, m, d, &mjd, &err);
+	if(err){
+		fprintf(stderr, "slacaldj(): Wrong %s!", (err == 1) ? "year" :
+			(err == 2? "month" : "day"));
+		return -1;
+	}
+	mjd += add;
+	DBG("MJD by slalib: %g\n", mjd);
+	double utc1, utc2;
+	/* UTC date. */
+	if(iauDtf2d("UTC", y, m, d, tms.tm_hour, tms.tm_min, tms.tm_sec,
+		&utc1, &utc2)) return -1;
+	DBG("UTC by sofa: %g, %g\n", utc1 - 2400000.5, utc2);
+	double tai1, tai2, tt1, tt2;
+	/* TT date. */
+	if ( iauUtctai ( utc1, utc2, &tai1, &tai2 ) ) return -1;
+	if ( iauTaitt ( tai1, tai2, &tt1, &tt2 ) ) return -1;
+	DBG("date by sofa (utc/tt): %g/%g & %g/%g\n", tai1 - 2400000.5, tt1 - 2400000.5, tai2, tt2);
+
+	double pmra=0, pr=0, pd=0, px=0, rv=0;
+	/*
+	// Proper motion: RA/Dec derivatives, epoch J2000.0.
+	pmra = 536.23e-3 * DAS2R;
+	pr = atan2 ( pmra, cos(dc) );
+	pd = 385.29e-3 * DAS2R;
+	// Parallax (arcsec) and recession speed (km/s).
+	px = 0.19495;
+	rv = -26.1;*/
+	double ri, di, eo;
+	/* ICRS to CIRS (geocentric observer). */
+	iauAtci13 ( rc, dc, pr, pd, px, rv, tt1, tt2, &ri, &di, &eo );
+	reprd ( "catalog -> CIRS:", ri, di );
+	double rca, dca;
+	/* CIRS to ICRS (astrometric). */
+	iauAtic13 ( ri, di, tt1, tt2, &rca, &dca, &eo );
+	reprd ( "CIRS -> astrometric:", rca, dca );
+	/* ICRS (astrometric) to CIRS (geocentric observer). */
+	iauAtci13 ( rca, dca, 0.0, 0.0, 0.0, 0.0, tt1, tt2, &ri, &di, &eo );
+	reprd ( "astrometric -> CIRS:", ri, di );
+	double ra, da;
+	/* Apparent place. */
+	ra = iauAnp ( ri - eo );
+	da = di;
+	reprd ( "geocentric apparent:", ra, da );
+	slamap(rc, dc, pmra, pd, px, rv, 2000., mjd, &ra, &da);
+	reprd ( "geocentric apparent (sla):", ra, da );
+	double ra2000, decl2000;
+	slaamp(ra, da, mjd, 2000.0, &ra2000, &decl2000);
+	reprd ( "apparent -> astrometric (sla):", ra2000, decl2000);
+
+	double elong, phi, hm, phpa, tc, rh, wl, xp, yp, dut1;
+	/* Site longitude, latitude (radians) and height above the geoid (m). */
+	iauAf2a ( '+', 41, 26, 26.45, &elong );
+	iauAf2a ( '+', 43, 39, 12.69, &phi );
+	hm = 2070.0;
+	/* Ambient pressure (HPa), temperature (C) and rel. humidity (frac). */
+	phpa = 770.0; // milliBar or hectopascal
+	tc = -5.0;
+	rh = 0.7;
+	/* Effective wavelength (microns) */
+	wl = 0.55;
+	/* EOPs: polar motion in radians, UT1-UTC in seconds. */
+	xp = 0.1074 * DAS2R; //polarX
+	yp = 0.2538 * DAS2R;//polarY
+	dut1 = 0.13026 ; // DUT1
+	/* ICRS to observed. */
+	double aob, zob, hob, dob, rob;
+	if ( iauAtco13 ( rc, dc, pr,
+	 pd, px, rv, utc1, utc2, dut1, elong, phi,
+	 hm, xp, yp, phpa, tc, rh, wl, &aob, &zob,
+	 &hob, &dob, &rob, &eo ) ) return -1;
+	reprd ( "ICRS -> observed:",  rob, dob );
+	printf("A(bta)/Z: ");
+	radtodeg(aob);
+	printf("("); radtodeg(DPI-aob);
+	printf(")/"); radtodeg(zob);
+	printf("\n");
+	if( iauAtoc13 ( "R", rc, dc, utc1, utc2, dut1,
+	//if( iauAtoc13 ( "R", rob, dob, 2451545, 0, dut1,
+	 elong, phi, hm, xp, yp, phpa, tc, rh, wl, &rca, &dca )) return -1;
+	reprd ( "observed -> astrometric:", rca, dca );
+	return 0;
+}
+
+

_snippets/calcAP/slalib_sofa_nova.c

@@ -0,0 +1,319 @@
+/*
+ * slalib_and_sofa.c - calculate apparent place by slalib & libsofa
+ *
+ * Copyright 2016 Edward V. Emelianov <eddy@sao.ru, edward.emelianoff@gmail.com>
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+
+#define _GNU_SOURCE 1111  // strcasecmp
+
+/**
+  SOFA - NOVA - SLA
+  comparison
+**/
+
+#include <stdio.h>
+#include <stdbool.h>
+#include <sys/time.h>
+#include <time.h>
+#include <sofa.h>
+#include <erfa.h> //- the same data as SOFA
+#include <libnova/libnova.h>
+
+#define DBG(...) printf(__VA_ARGS__)
+
+extern void sla_caldj(int*, int*, int*, double*, int*);
+extern void sla_amp(double*, double*, double*, double*, double*, double*);
+extern void sla_map(double*, double*, double*, double*, double*,double*, double*, double*, double*, double*);
+void slacaldj(int y, int m, int d, double *djm, int *j){
+    int iy = y, im = m, id = d;
+    sla_caldj(&iy, &im, &id, djm, j);
+}
+void slaamp(double ra, double da, double date, double eq, double *rm, double *dm ){
+    double r = ra, d = da, mjd = date, equi = eq;
+    sla_amp(&r, &d, &mjd, &equi, rm, dm);
+}
+// apparent->observed
+extern void sla_aop(double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*);
+/*
+*  Given:
+*     RAP    d      geocentric apparent right ascension
+*     DAP    d      geocentric apparent declination
+*     DATE   d      UTC date/time (Modified Julian Date, JD-2400000.5)
+*     DUT    d      delta UT:  UT1-UTC (UTC seconds)
+*     ELONGM d      mean longitude of the observer (radians, east +ve)
+*     PHIM   d      mean geodetic latitude of the observer (radians)
+*     HM     d      observer's height above sea level (metres)
+*     XP     d      polar motion x-coordinate (radians)
+*     YP     d      polar motion y-coordinate (radians)
+*     TDK    d      local ambient temperature (K; std=273.15D0)
+*     PMB    d      local atmospheric pressure (mb; std=1013.25D0)
+*     RH     d      local relative humidity (in the range 0D0-1D0)
+*     WL     d      effective wavelength (micron, e.g. 0.55D0)
+*     TLR    d      tropospheric lapse rate (K/metre, e.g. 0.0065D0)
+*
+*  Returned:
+*     AOB    d      observed azimuth (radians: N=0,E=90)
+*     ZOB    d      observed zenith distance (radians)
+*     HOB    d      observed Hour Angle (radians)
+*     DOB    d      observed Declination (radians)
+*     ROB    d      observed Right Ascension (radians)
+*/
+void slaaop(double rap, double dap, double date, double dut, double elongm, double phim,
+            double hm, double xp, double yp, double tdk, double pmb, double rh, double wl,
+            double tlr,
+            double* aob, double* zob, double* hob, double* dob, double* rob){
+    double _rap=rap, _dap=dap, _date=date, _dut=dut, _elongm=elongm, _phim=phim,
+            _hm=hm, _xp=xp, _yp=yp, _tdk=tdk, _pmb=pmb, _rh=rh, _wl=wl, _tlr=tlr;
+    sla_aop(&_rap, &_dap, &_date, &_dut, &_elongm, &_phim, &_hm, &_xp, &_yp, &_tdk,
+            &_pmb, &_rh, &_wl, &_tlr, aob, zob, hob, dob, rob);
+}
+
+// rm,dm - mean RA,Dec (rad), pr,pd - RA,Dec changes per Julian year (dRA/dt, dDec/dt)
+// px - parallax (arcsec), rv - radial speed (km/sec,  +ve if receding)
+// eq - epoch and equinox of star data (Julian)
+// date - TDB for apparent place (JD-2400000.5)
+void slamap(double rm, double dm, double pr, double pd,
+             double px, double rv, double eq, double date,
+             double *ra, double *da){
+    double r = rm, d = dm, p1 = pr, p2 = pd, ppx = px, prv = rv, equi = eq, dd = date;
+    sla_map(&r, &d, &p1, &p2, &ppx, &prv, &equi, &dd, ra, da);
+}
+
+void reprd(char* s, double ra, double dc){
+    char pm;
+    int i[4];
+    printf ( "%30s", s );
+    iauA2tf ( 7, ra, &pm, i );
+    printf ( " %2.2d %2.2d %2.2d.%7.7d", i[0],i[1],i[2],i[3] );
+    iauA2af ( 6, dc, &pm, i );
+    printf ( " %c%2.2d %2.2d %2.2d.%6.6d\n", pm, i[0],i[1],i[2],i[3] );
+}
+
+void radtodeg(double r){
+    int i[4]; char pm;
+    int rem = (int)(r / D2PI);
+    if(rem) r -= D2PI * rem;
+    if(r > DPI) r -= D2PI;
+    else if(r < -DPI) r += D2PI;
+    iauA2af (2, r, &pm, i);
+    printf("%c%02d %02d %02d.%2.d", pm, i[0],i[1],i[2],i[3]);
+}
+
+double getta(char *str){
+    int a,b,s = 1; double c;
+    if(3 != sscanf(str, "%d:%d:%lf", &a,&b,&c)) return -1;
+    if(a < 0 || *str == '-'){ s = -1; a = -a;}
+    c /= 3600.;
+    c += a + b/60.;
+    c *= s;
+    return c;
+}
+
+int main (int argc, char **argv){
+    double rc, dc;
+    if(argc == 3){
+        rc = getta(argv[1]) * DPI / 12;
+        dc = getta(argv[2]) * DD2R;
+    }else{
+        /* Star ICRS RA,Dec (radians). */
+        if ( iauTf2a ( ' ', 19, 50, 47.6, &rc ) ) return -1;
+        if ( iauAf2a ( '+', 8, 52, 12.3, &dc ) ) return -1;
+    }
+    reprd ( "ICRS (catalog), epoch J2000.0:", rc, dc );
+
+    struct tm tms;
+    struct timeval currentTime;
+    gettimeofday(&currentTime, NULL);
+    gmtime_r(&currentTime.tv_sec, &tms);
+    double tSeconds = tms.tm_sec + ((double)currentTime.tv_usec)/1e6;
+    int y, m, d, err;
+    y = 1900 + tms.tm_year;
+    m = tms.tm_mon + 1;
+    d = tms.tm_mday;
+    double mjd, add = ((double)tms.tm_hour + (double)tms.tm_min/60.0 + tSeconds/3600.0) / 24.;
+    DBG("Date: (d/m/y +frac)  %d/%d/%d +%.8f\n", d, m, y, add);
+    slacaldj(y, m, d, &mjd, &err);
+    if(err){
+        fprintf(stderr, "slacaldj(): Wrong %s!", (err == 1) ? "year" :
+            (err == 2? "month" : "day"));
+        return -1;
+    }
+    mjd += add;
+    DBG("MJD by slalib: %.8f\n", mjd);
+    double utc1, utc2;
+    /* UTC date. */
+    if(iauDtf2d("UTC", y, m, d, tms.tm_hour, tms.tm_min, tSeconds,
+        &utc1, &utc2)) return -1;
+    DBG("UTC by sofa: %g, %.8f\n", utc1 - 2400000.5, utc2);
+    double tai1, tai2, tt1, tt2;
+    /* TT date. */
+    if ( iauUtctai ( utc1, utc2, &tai1, &tai2 ) ) return -1;
+    if ( iauTaitt ( tai1, tai2, &tt1, &tt2 ) ) return -1;
+    DBG("date by sofa (TAI/TT): %g/%g & %g/%g\n", tai1 - 2400000.5, tt1 - 2400000.5, tai2, tt2);
+
+    double pmra=0;
+    double pr = 0.0;     // RA proper motion (radians/year; Note 2)
+    double pd = 0.0;     // Dec proper motion (radians/year)
+    double px = 0.0;     // parallax (arcsec)
+    double rv = 0.0;     // radial velocity (km/s, positive if receding)
+    /*
+    // Proper motion: RA/Dec derivatives, epoch J2000.0.
+    pmra = 536.23e-3 * DAS2R;
+    pr = atan2 ( pmra, cos(dc) );
+    pd = 385.29e-3 * DAS2R;
+    // Parallax (arcsec) and recession speed (km/s).
+    px = 0.19495;
+    rv = -26.1;*/
+    double ri, di, eo;
+    /* ICRS to CIRS (geocentric observer). */
+    iauAtci13 ( rc, dc, pr, pd, px, rv, tt1, tt2, &ri, &di, &eo );
+    reprd ( "ICRS -> CIRS:", ri, di );
+    double rca, dca, eo1;
+    /* CIRS to ICRS (astrometric). */
+    iauAtic13 ( ri, di, tt1, tt2, &rca, &dca, &eo1);
+    reprd ( "CIRS -> ICRSc:", rca, dca );
+    /* ICRS to CIRS without PM
+    iauAtci13 ( rca, dca, 0.0, 0.0, 0.0, 0.0, tt1, tt2, &ri, &di, &eo );
+    reprd ( "ICRSc -> CIRS (without PM):", ri, di ); */
+    double ra, da;
+    /* Apparent place. */
+    ra = iauAnp ( ri - eo );
+    da = di;
+    reprd ( "geocentric apparent:", ra, da );
+    slamap(rc, dc, pmra, pd, px, rv, 2000., mjd, &ra, &da);
+    reprd ( "geocentric apparent (sla):", ra, da );
+    double ra2000, decl2000;
+    slaamp(ra, da, mjd, 2000.0, &ra2000, &decl2000);
+    reprd ( "apparent -> astrometric (sla):", ra2000, decl2000);
+
+    double elong, phi, hm, phpa, tc, rh, wl, xp, yp, dut1;
+    /* Site longitude, latitude (radians) and height above the geoid (m). */
+    iauAf2a ( '+', 41, 26, 26.45, &elong );
+    iauAf2a ( '+', 43, 39, 12.69, &phi );
+    hm = 2070.0;
+    /* Ambient pressure (HPa), temperature (C) and rel. humidity (frac). */
+    phpa = 780.0; tc = -5.0; rh = 0.7;
+    /* Effective wavelength (microns) */
+    wl = 0.55;
+    /* EOPs: polar motion in radians, UT1-UTC in seconds. */
+    xp = 0.1074 * DAS2R; //polarX
+    yp = 0.2538 * DAS2R;//polarY
+    dut1 = 0.13026 ; // DUT1
+    /* ICRS to observed. */
+    double aob, zob, hob, dob, rob;
+    printf("iauAtco13(%g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g, %g)\n",
+           rc, dc, pr, pd, px, rv, utc1, utc2, dut1, elong, phi, hm, xp, yp, phpa, tc, rh, wl);
+    if ( iauAtco13 ( rc, dc, pr,
+     pd, px, rv, utc1, utc2, dut1, elong, phi,
+     hm, xp, yp, phpa, tc, rh, wl, &aob, &zob,
+     &hob, &dob, &rob, &eo ) ) return -1;
+    reprd ( "ICRS -> observed:",  rob, dob );
+    printf("A(bta)/Z: ");
+    radtodeg(aob);
+    printf("("); radtodeg(DPI-aob);
+    printf(")/"); radtodeg(zob);
+    printf("\n");
+    if( iauAtoc13 ( "R", rob, dob, utc1, utc2, dut1,
+     elong, phi, hm, xp, yp, phpa, tc, rh, wl, &rca, &dca )) return -1;
+    reprd ( "observed -> ICRS:", rca, dca );
+
+    /*
+    *  Given:
+    *     RAP    d      geocentric apparent right ascension
+    *     DAP    d      geocentric apparent declination
+    *     DATE   d      UTC date/time (Modified Julian Date, JD-2400000.5)
+    *     DUT    d      delta UT:  UT1-UTC (UTC seconds)
+    *     ELONGM d      mean longitude of the observer (radians, east +ve)
+    *     PHIM   d      mean geodetic latitude of the observer (radians)
+    *     HM     d      observer's height above sea level (metres)
+    *     XP     d      polar motion x-coordinate (radians)
+    *     YP     d      polar motion y-coordinate (radians)
+    *     TDK    d      local ambient temperature (K; std=273.15D0)
+    *     PMB    d      local atmospheric pressure (mb; std=1013.25D0)
+    *     RH     d      local relative humidity (in the range 0D0-1D0)
+    *     WL     d      effective wavelength (micron, e.g. 0.55D0)
+    *     TLR    d      tropospheric lapse rate (K/metre, e.g. 0.0065D0)
+    *
+    *  Returned:
+    *     AOB    d      observed azimuth (radians: N=0,E=90)
+    *     ZOB    d      observed zenith distance (radians)
+    *     HOB    d      observed Hour Angle (radians)
+    *     DOB    d      observed Declination (radians)
+    *     ROB    d      observed Right Ascension (radians)
+    */
+    slaaop(ra, da, mjd, dut1, elong, phi, hm, xp, yp, tc+273., phpa, rh, wl, 0.0065,
+           &aob, &zob, &hob, &dob, &rob);
+    reprd ( "ICRS -> observed (sla):",  rob, dob );
+    printf("A(bta)/Z: ");
+    radtodeg(aob);
+    printf("("); radtodeg(DPI-aob);
+    printf(")/"); radtodeg(zob);
+    printf("\n");
+
+    ////////////////////////////////////////////////////////////////////////////////////////////////
+    // libNOVA
+    ////////////////////////////////////////////////////////////////////////////////////////////////
+    struct ln_equ_posn mean_position;
+    mean_position.ra  = rc * ERFA_DR2D; // radians to degrees
+    mean_position.dec = dc * ERFA_DR2D;
+/*
+    struct timeval tv;
+    struct timezone tz;
+    gettimeofday(&tv, &tz); // number of seconds since the Epoch, 1970-01-01 00:00:00 +0000 (UTC) with microsecond precision
+    struct tm *utc_tm;
+    utc_tm = gmtime(&tv.tv_sec);
+    struct ln_date date;
+    date.seconds = utc_tm->tm_sec + ((double)tv.tv_usec / 1000000);
+    date.minutes = utc_tm->tm_min;
+    date.hours   = utc_tm->tm_hour;
+    date.days    = utc_tm->tm_mday;
+    date.months  = utc_tm->tm_mon + 1;
+    date.years   = utc_tm->tm_year + 1900;
+    double JNow = ln_get_julian_day(&date);*/
+    double JNow = ln_get_julian_from_sys();
+    struct ln_equ_posn propm={0,0}, apppl;//, equprec;
+    ln_get_apparent_posn(&mean_position, &propm, JNow, &apppl);
+    reprd ("geocentric apparent (NOVA):",  apppl.ra*ERFA_DD2R, apppl.dec*ERFA_DD2R);
+    // Calculate the effects of precession on equatorial coordinates, between arbitary Jxxxx epochs.
+/*
+    ln_get_equ_prec2(&mean_position, JD2000, JNow, &equprec);
+    reprd ("ln_get_equ_prec2 (NOVA):",  equprec.ra*ERFA_DD2R, equprec.dec*ERFA_DD2R);
+    */
+/*
+    // ERFA
+    struct tm *ts;
+    ts = gmtime(&t);
+    int result = eraDtf2d ( "UTC", ts->tm_year+1900, ts->tm_mon+1, ts->tm_mday, ts->tm_hour, ts->tm_min, ts->tm_sec, &utc1, &utc2 );
+    if (result != 0) {
+        printf("eraDtf2d call failed\n");
+        return 1;
+    }
+    // Make TT julian date for Atci13 call
+    result = eraUtctai( utc1, utc2, &tai1, &tai2 );
+    if (result != 0) {
+        printf("eraUtctai call failed\n");
+        return 1;
+    }
+    eraTaitt(  tai1, tai2, &tt1,  &tt2  );
+    eraAtci13 ( rc, dc, pr, pd, px, rv, tt1, tt2, &ri, &di, &eo );
+    reprd ( "geocentric apparent (ERFA):",  eraAnp(ri - eo), di);
+*/
+    return 0;
+}
+

_snippets/clr.c

@@ -0,0 +1,21 @@
+// make a simple "CLS"
+
+#include <stdio.h>
+#include <unistd.h>
+
+void printstrings(const char *add){
+    for(int i = 0; i < 40; ++i)
+        printf("String %d - %s\n", i, add);
+}
+
+const char *x[] = {"first", "second", "third"};
+
+int main(){
+    for(int i = 0; i < 3; ++i){
+        printf("\033c");
+        printstrings(x[i]);
+        sleep(1);
+    }
+    printf("\033c");
+    return 0;
+}

_snippets/comments_before.c

@@ -0,0 +1,97 @@
+/*
+ * comments_before.c - put multiline comments before text in string they follow
+ *
+ *
+ * Copyright 2015 Edward V. Emelianov <eddy@sao.ru, edward.emelianoff@gmail.com>
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#include <stdio.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <sys/mman.h>
+#include <unistd.h>
+#include <stdlib.h>
+#include <string.h>
+
+typedef struct{
+	char *data;
+	size_t len;
+} mmapbuf;
+
+mmapbuf *My_mmap(char *filename){
+	int fd;
+	char *ptr;
+	size_t Mlen;
+	struct stat statbuf;
+	if(!filename) exit(-1);
+	if((fd = open(filename, O_RDONLY)) < 0) exit(-2);
+	if(fstat (fd, &statbuf) < 0) exit(-3);
+	Mlen = statbuf.st_size;
+	if((ptr = mmap (0, Mlen, PROT_READ, MAP_PRIVATE, fd, 0)) == MAP_FAILED)
+		exit(-4);
+	if(close(fd)) exit(-5);
+	mmapbuf *ret = calloc(1, sizeof(mmapbuf));
+	ret->data = ptr;
+	ret->len = Mlen;
+	return  ret;
+}
+
+
+int main(int argc, char **argv){
+	mmapbuf *buf = My_mmap(argv[1]);
+	char *afile = strdup(buf->data), *eptr = afile + buf->len;
+	char *beforecom = NULL;
+	while(afile < eptr){
+		char *nl = strchr(afile, '\n');
+		if(!nl){
+			printf("%s\n", afile);
+			return 0;
+		}
+		*nl++ = 0;
+		char *commstart = strstr(afile, "/*");
+		if(commstart){
+			*commstart = 0;
+			free(beforecom);
+			beforecom = strdup(afile);
+			commstart += 2;
+			char *commend = NULL;
+			if(*commstart) commend = strstr(commstart, "*/");
+			if(commend){ // single comment
+				*commend = 0;
+				afile = commend + 2;
+				printf("/* %s */\n", commstart);
+				printf("%s%s\n", beforecom, afile);
+				afile = nl;
+			}else{ // multiline comment
+				printf("/*%s\n", commstart);
+				commend = strstr(nl, "*/");
+				if(!commend){
+					printf("%s*/\n", nl);
+					return -1;
+				}
+				*commend = 0;
+				afile = commend + 2;
+				if(*afile == '\n') ++afile;
+				printf("%s*/\n", nl);
+			}
+		}else{
+			printf("%s\n", afile);
+			afile = nl;
+		}
+	}
+}

_snippets/count_bits.c

@@ -0,0 +1 @@
+__builtin_popcount(x) - counts all non-zero bits in number

_snippets/cryptpass.c

@@ -0,0 +1,34 @@
+#define _GNU_SOURCE
+#include <crypt.h> // gcc -lcrypt
+#include <errno.h>
+#include <stdio.h>
+
+int main(int argc, char **argv){
+    if(argc != 3){
+        printf("USAGE:\n\t%s <password> <salt>\n", argv[0]);
+        printf("Example:\n\t%s password \\$6\\$salt  -> SHA512 with given salt (see man 3 crypt)\n", argv[0]);
+        return 1;
+    }
+    struct crypt_data x;
+    x.initialized = 0;
+
+    char *c = crypt_r(argv[1], argv[2], &x);
+    int e = errno;
+    if(!c){
+        switch(e){
+            case EINVAL:
+                fprintf(stderr, "Wrong SALT format!\n");
+            break;
+            case ENOSYS:
+                fprintf(stderr, "crypt() not realized!\n");
+            break;
+            case EPERM:
+                fprintf(stderr, "/proc/sys/crypto/fips_enabled prohibits this encryption method!\n");
+            default:
+                fprintf(stderr, "Unknown error!\n");
+        }
+        return e;
+    }
+    printf("Result of crypt_r():\n%s\n", c);
+    return 0;
+}

_snippets/factor.c

@@ -0,0 +1,144 @@
+#include <stdio.h>
+#include <stdint.h>
+#include <omp.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <usefull_macros.h>
+
+#define THREADNO    8
+
+typedef struct _num{
+    uint64_t n;
+    uint8_t isprime;
+    struct _num *next;
+    struct _num *prev;
+} num;
+
+static num *listhead = NULL, *listtail = NULL;
+
+static void ins(uint64_t n, uint8_t isprime){//  printf("ins %zd %d\n", n, isprime);
+    num *N = calloc(1, sizeof(num));
+    N->n = n; N->isprime = isprime;
+    if(listhead){
+        N->prev = listtail;
+        listtail->next = N;
+        listtail = N;
+    }else{
+        listhead = N;
+        listtail = N;
+    }
+}
+
+static num *del(num *el){// printf("del %zd %d\n", el->n, el->isprime);
+    num *nxt = el->next, *prev = el->prev;
+    if(el == listhead) listhead = nxt;
+    if(el == listtail) listtail = prev;
+    if(el->prev) el->prev->next = nxt;
+    if(el->next) el->next->prev = prev;
+    free(el);
+    return nxt;
+}
+
+static uint64_t __1st[2] = {2,3};
+static void factorize(uint64_t n){
+    if(n < 2) return;
+    for(int i = 0; i < 2; ++i){
+        if(n % __1st[i] == 0){
+            ins(__1st[i], 1);
+            n /= __1st[i];
+        }
+    }
+    if(n < 2) return;
+    uint8_t prime = 1;
+    register uint64_t last = (2 + sqrt(n))/6;
+    #pragma omp parallel for shared(n, prime)
+    for(uint64_t i = 1; i < last; ++i){
+        uint64_t nmb = i*6 - 1;
+        for(uint64_t x = 0; x < 3; x += 2){
+        nmb += x;
+        if(n % nmb == 0){
+            #pragma omp critical
+            {
+                ins(nmb, 0);
+                prime = 0;
+            }
+        }}
+    }
+    if(prime){ // prime number
+        ins(n, 1);
+    }
+}
+
+static void calc_fact(uint64_t numb){
+    num *n = NULL;
+    if(listhead){ // free previous list
+        n = listhead;
+        do{
+            num *nxt = n->next;
+            free(n); n = nxt;
+        }while(n);
+        listhead = listtail = NULL;
+    }
+    factorize(numb);
+    n = listhead;
+    do{
+        if(n->n == numb){
+            n = del(n);
+            continue;
+        }
+        if(!n->isprime){
+            uint64_t p = n->n;
+            n = del(n);
+            factorize(p);
+            continue;
+        }
+        n = n->next;
+    }while(n);
+}
+
+int main(int argc, char **argv){
+    if(argc != 2){
+        printf("Usage: %s number or -t for test\n", argv[0]);
+        return 1;
+    }
+    omp_set_num_threads(THREADNO);
+    if(strcmp(argv[1], "-t") == 0){ // run test
+        uint64_t nmax = 0; double tmax = 0.;
+        for(uint64_t x = UINT64_MAX; x > 1ULL<<58; --x){
+            double t0 = dtime();
+            calc_fact(x);
+            double t = dtime() - t0;
+            if(t > tmax){
+                tmax = t;
+                nmax = x;
+                printf("n = %lu, t=%g\n", nmax, tmax);
+            }
+        }
+        return 0;
+    }
+    uint64_t numb = atoll(argv[1]), norig = numb;
+    calc_fact(numb);
+    printf("Number: %zd\n", numb);
+    if(!listhead || (listhead == listtail && listhead->n == numb)){
+        printf("Prime number %zd\n", numb);
+        return 0;
+    }
+    num *n = listhead;
+    printf("Fact: ");
+    uint64_t test = 1;
+    do{
+        while(numb % n->n == 0){
+            printf("%zd ", n->n);
+            numb /= n->n;
+            test *= n->n;
+        }
+        n = n->next;
+    }while(n);
+    if(numb != 1){
+        printf("%zd", numb);
+        test *= numb;
+    }
+    printf("\n\nTest: %zd == %zd ?\n", norig, test);
+    return 0;
+}

_snippets/fixed_sort_example.c

@@ -0,0 +1,98 @@
+/*
+ * fixed_sort_example.c
+ *
+ * Copyright 2015 Edward V. Emelianov <eddy@sao.ru, edward.emelianoff@gmail.com>
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ */
+
+#include <stdio.h>
+
+static __inline__ unsigned long long rdtsc(void){
+	unsigned hi, lo;
+	__asm__ __volatile__ ("rdtsc" : "=a"(lo), "=d"(hi));
+	return ( (unsigned long long)lo)|( ((unsigned long long)hi)<<32 );
+}
+int arra[5][3] = {{3,2,1}, {3,6,4}, {7,5,8}, {1024,5543,9875}, {1001,-1001,1002}};
+int arrb[5][3] = {{3,2,1}, {3,6,4}, {7,5,8}, {1024,5543,9875}, {1001,-1001,1002}};
+int arrc[5][3] = {{3,2,1}, {3,6,4}, {7,5,8}, {1024,5543,9875}, {1001,-1001,1002}};
+
+/*
+ * It seems that this manner of sorting would be less productive than sort3b or sort3c
+ * but even on -O1 it gives better results (median by 1000 seq):
+ *
+ * -Ox  timing a    timing b   timing c
+ *   0   453        402        366
+ *   1   111        174        147
+ *   2   126        165        144
+ *   3   117        171        141
+ */
+static inline void sort3a(int *d){
+#define min(x, y) (x<y?x:y)
+#define max(x, y) (x<y?y:x)
+#define SWAP(x,y) { const int a = min(d[x], d[y]); const int b = max(d[x], d[y]); d[x] = a; d[y] = b;}
+	SWAP(0, 1);
+	SWAP(1, 2);
+	SWAP(0, 1);
+#undef SWAP
+#undef min
+#undef max
+}
+
+static inline void sort3b(int *d){
+#define CMPSWAP(x,y) {if(d[x] > d[y]){const int a = d[x]; d[x] = d[y]; d[y] = a;}}
+	CMPSWAP(0, 1);
+	CMPSWAP(1, 2);
+	CMPSWAP(0, 1);
+#undef CMPSWAP
+}
+
+static inline void sort3c(int *d){
+	register int a, b;
+#define CMPSWAP(x,y) {a = d[x]; b = d[y]; if(a > b){d[x] = b; d[y] = a;}}
+	CMPSWAP(0, 1);
+	CMPSWAP(1, 2);
+	CMPSWAP(0, 1);
+#undef CMPSWAP
+}
+
+int main(){
+	unsigned long long time1, time2, time3;
+	int i;
+	time1 = rdtsc();
+	for(i = 0; i < 5 ; i++){
+		sort3a(arra[i]);
+	}
+	time1 = rdtsc() - time1;
+	time2 = rdtsc();
+	for(i = 0; i < 5 ; i++){
+		sort3b(arrb[i]);
+	}
+	time2 = rdtsc() - time2;
+	time3 = rdtsc();
+	for(i = 0; i < 5 ; i++){
+		sort3c(arrc[i]);
+	}
+	time3 = rdtsc() - time3;
+	printf("%llu, %llu, %llu; ", time1, time2, time3);
+/*	printf("Time is %llu (A) & %llu (B)\nData:\n\n", time1, time2);
+	for(i = 0; i < 5; i++)
+		printf("A[%d]: %d, %d, %d\n", i, arra[i][0], arra[i][1], arra[i][2]);
+	printf("\n");
+	for(i = 0; i < 5; i++)
+		printf("B[%d]: %d, %d, %d\n", i, arrb[i][0], arrb[i][1], arrb[i][2]);*/
+}
+

_snippets/get_avail_mem.c

@@ -0,0 +1,14 @@
+#include <stdio.h>
+#include <unistd.h>
+
+static unsigned long long get_available_mem(){
+    return sysconf(_SC_AVPHYS_PAGES) * (unsigned long long) sysconf(_SC_PAGE_SIZE);
+}
+
+int main(){
+	unsigned long long m = get_available_mem();
+	printf("MEM: %llu == %lluGB\n", m, m/1024/1024/1024);
+	return 0;
+}
+
+// Never allocate memory by big pieces with malloc! Only mmap with MAP_POPULATE!!!!!!!!!11111111111

_snippets/inotify.c

@@ -0,0 +1,118 @@
+#include <errno.h>
+#include <poll.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/inotify.h>
+#include <unistd.h>
+
+/* Read all available inotify events from the file descriptor 'fd'.
+  wd is the table of watch descriptors for the directories in argv.
+  argc is the length of wd and argv.
+  argv is the list of watched directories.
+  Entry 0 of wd and argv is unused. */
+
+static void
+handle_events(int fd, int wd, int argc, char* argv[])
+{
+   /* Some systems cannot read integer variables if they are not
+      properly aligned. On other systems, incorrect alignment may
+      decrease performance. Hence, the buffer used for reading from
+      the inotify file descriptor should have the same alignment as
+      struct inotify_event. */
+    struct inotify_event buf;
+   int i;
+   ssize_t len;
+   char *ptr;
+
+   /* Loop while events can be read from inotify file descriptor. */
+
+   for (;;) {
+
+       /* Read some events. */
+
+       len = read(fd, &buf, sizeof buf);
+       if (len == -1 && errno != EAGAIN) {
+           perror("read");
+           exit(EXIT_FAILURE);
+       }
+
+       /* If the nonblocking read() found no events to read, then
+          it returns -1 with errno set to EAGAIN. In that case,
+          we exit the loop. */
+
+       if (len <= 0)
+           break;
+       if (buf.mask & IN_CLOSE_WRITE)
+           printf("IN_CLOSE_WRITE: ");
+       /* Print the name of the watched directory */
+       for (i = 1; i < argc; ++i) {
+           if (wd == buf.wd) {
+               printf("%s/", argv[i]);
+               break;
+           }
+       }
+
+       /* Print the name of the file */
+
+       if (buf.len)
+           printf("%s", buf.name);
+
+       /* Print type of filesystem object */
+
+       if (buf.mask & IN_ISDIR)
+           printf(" [directory]\n");
+       else
+           printf(" [file]\n");
+   }
+
+}
+
+int
+main(int argc, char* argv[])
+{
+    int fd, i, poll_num;
+    int wd;
+    struct pollfd fds;
+
+    if (argc != 2) {
+       printf("Usage: %s file\n", argv[0]);
+       exit(EXIT_FAILURE);
+    }
+    fd = inotify_init1(IN_NONBLOCK);
+    if(fd == -1) {
+       perror("inotify_init1");
+       exit(EXIT_FAILURE);
+    }
+    // check if there's a file
+    FILE* f = fopen(argv[1], "r");
+    if(!f){perror("open"); return 1;}
+    char buf[256];
+    if(fread(buf, 1, 256, f) < 1){ perror("read"); return 1;}
+    fclose(f);
+
+    if(-1 == (wd = inotify_add_watch(fd, argv[1], IN_CLOSE_WRITE))){
+        perror("inotify_add_watch");
+        exit(EXIT_FAILURE);
+    }
+    /* Inotify input */
+    fds.fd = fd;
+    fds.events = POLLIN;
+    /* Wait for events */
+    printf("Listening for events.\n");
+    while(1){
+       poll_num = poll(&fds, 1, -1);
+       if(poll_num == -1){
+           if (errno == EINTR)
+               continue;
+           perror("poll");
+           exit(EXIT_FAILURE);
+       }
+       if (poll_num > 0) {
+           if(fds.revents & POLLIN){
+               handle_events(fd, wd, argc, argv);
+           }
+       }
+    }
+    close(fd);
+    exit(EXIT_SUCCESS);
+}

_snippets/med.c

@@ -0,0 +1,245 @@
+//Copyright (c) 2011 ashelly.myopenid.com under <http://www.opensource.org/licenses/mit-license>
+#include <stdlib.h>
+#include <sys/time.h>
+#include <string.h>
+
+#define OMP_NUM_THREADS 4
+#define Stringify(x) #x
+#define OMP_FOR(x) _Pragma(Stringify(omp parallel for x))
+
+//Customize for your data Item type
+typedef int Item;
+#define ItemLess(a,b) ((a)<(b))
+#define ItemMean(a,b) (((a)+(b))/2)
+
+typedef struct Mediator_t{
+	Item* data; // circular queue of values
+	int* pos;   // index into `heap` for each value
+	int* heap;  // max/median/min heap holding indexes into `data`.
+	int N;      // allocated size.
+	int idx;    // position in circular queue
+	int ct;     // count of items in queue
+} Mediator;
+
+/*--- Helper Functions ---*/
+
+#define minCt(m) (((m)->ct-1)/2) //count of items in minheap
+#define maxCt(m) (((m)->ct)/2) //count of items in maxheap
+
+//returns 1 if heap[i] < heap[j]
+inline int mmless(Mediator* m, int i, int j){
+	return ItemLess(m->data[m->heap[i]],m->data[m->heap[j]]);
+}
+
+//swaps items i&j in heap, maintains indexes
+inline int mmexchange(Mediator* m, int i, int j){
+	int t = m->heap[i];
+	m->heap[i] = m->heap[j];
+	m->heap[j] = t;
+	m->pos[m->heap[i]] = i;
+	m->pos[m->heap[j]] = j;
+	return 1;
+}
+
+//swaps items i&j if i<j; returns true if swapped
+inline int mmCmpExch(Mediator* m, int i, int j){
+	return (mmless(m,i,j) && mmexchange(m,i,j));
+}
+
+//maintains minheap property for all items below i/2.
+void minSortDown(Mediator* m, int i){
+	for(; i <= minCt(m); i*=2){
+		if(i>1 && i < minCt(m) && mmless(m, i+1, i)) ++i;
+		if(!mmCmpExch(m,i,i/2)) break;
+	}
+}
+
+//maintains maxheap property for all items below i/2. (negative indexes)
+void maxSortDown(Mediator* m, int i){
+	for(; i >= -maxCt(m); i*=2){
+		if(i<-1 && i > -maxCt(m) && mmless(m, i, i-1)) --i;
+	if(!mmCmpExch(m,i/2,i)) break;
+	}
+}
+
+//maintains minheap property for all items above i, including median
+//returns true if median changed
+int minSortUp(Mediator* m, int i){
+	while (i > 0 && mmCmpExch(m, i, i/2)) i /= 2;
+	return (i == 0);
+}
+
+//maintains maxheap property for all items above i, including median
+//returns true if median changed
+int maxSortUp(Mediator* m, int i){
+	while (i < 0 && mmCmpExch(m, i/2, i)) i /= 2;
+	return (i == 0);
+}
+
+/*--- Public Interface ---*/
+
+
+//creates new Mediator: to calculate `nItems` running median.
+//mallocs single block of memory, caller must free.
+Mediator* MediatorNew(int nItems){
+	int size = sizeof(Mediator) + nItems*(sizeof(Item)+sizeof(int)*2);
+	Mediator* m = malloc(size);
+	m->data = (Item*)(m + 1);
+	m->pos = (int*) (m->data + nItems);
+	m->heap = m->pos + nItems + (nItems / 2); //points to middle of storage.
+	m->N = nItems;
+	m->ct = m->idx = 0;
+	while (nItems--){ //set up initial heap fill pattern: median,max,min,max,...
+		m->pos[nItems] = ((nItems+1)/2) * ((nItems&1)? -1 : 1);
+		m->heap[m->pos[nItems]] = nItems;
+	}
+	return m;
+}
+
+
+//Inserts item, maintains median in O(lg nItems)
+void MediatorInsert(Mediator* m, Item v){
+	int isNew=(m->ct<m->N);
+	int p = m->pos[m->idx];
+	Item old = m->data[m->idx];
+	m->data[m->idx]=v;
+	m->idx = (m->idx+1) % m->N;
+	m->ct+=isNew;
+	if(p>0){ //new item is in minHeap
+		if (!isNew && ItemLess(old,v)) minSortDown(m,p*2);
+		else if (minSortUp(m,p)) maxSortDown(m,-1);
+	}else if (p<0){ //new item is in maxheap
+		if (!isNew && ItemLess(v,old)) maxSortDown(m,p*2);
+		else if (maxSortUp(m,p)) minSortDown(m, 1);
+	}else{ //new item is at median
+		if (maxCt(m)) maxSortDown(m,-1);
+		if (minCt(m)) minSortDown(m, 1);
+	}
+}
+
+//returns median item (or average of 2 when item count is even)
+Item MediatorMedian(Mediator* m, Item *minval, Item *maxval){
+	Item v= m->data[m->heap[0]];
+	if ((m->ct&1) == 0) v = ItemMean(v,m->data[m->heap[-1]]);
+	Item min, max;
+	if(minval){
+		min = v;
+		int i;
+		for(i = -maxCt(m); i < 0; ++i){
+			int v = m->data[m->heap[i]];
+			if(v < min) min = v;
+		}
+		*minval = min;
+	}
+	if(maxval){
+		max = v;
+		int i;
+		for(i = 0; i <= minCt(m); ++i){
+			int v = m->data[m->heap[i]];
+			if(v > max) max = v;
+		}
+		*maxval = v;
+	}
+	return v;
+}
+
+
+/*--- Test Code ---*/
+#include <stdio.h>
+void PrintMaxHeap(Mediator* m){
+	int i;
+	if(maxCt(m))
+	printf("Max: %3d",m->data[m->heap[-1]]);
+	for (i=2;i<=maxCt(m);++i){
+		printf("|%3d ",m->data[m->heap[-i]]);
+		if(++i<=maxCt(m)) printf("%3d",m->data[m->heap[-i]]);
+	}
+	printf("\n");
+}
+void PrintMinHeap(Mediator* m){
+	int i;
+	if(minCt(m))
+	printf("Min: %3d",m->data[m->heap[1]]);
+	for (i=2;i<=minCt(m);++i){
+		printf("|%3d ",m->data[m->heap[i]]);
+		if(++i<=minCt(m)) printf("%3d",m->data[m->heap[i]]);
+	}
+	printf("\n");
+}
+
+void ShowTree(Mediator* m){
+	PrintMaxHeap(m);
+	printf("Mid: %3d\n",m->data[m->heap[0]]);
+	PrintMinHeap(m);
+	printf("\n");
+}
+
+double dtime(){
+	double t;
+	struct timeval tv;
+	gettimeofday(&tv, NULL);
+	t = tv.tv_sec + ((double)tv.tv_usec)/1e6;
+	return t;
+}
+
+#define SZ (4000)
+int main(int argc, char* argv[])
+{
+	int siz = SZ * SZ;
+	int i, *v = malloc(siz * sizeof(int)), *o = malloc(siz*sizeof(int)),
+		*small = malloc(siz * sizeof(int)), *large = malloc(siz * sizeof(int));
+
+for (i=0; i < siz; ++i){
+	int s = (rand() & 0xf) + 0xf00, d = rand() & 0xffff;
+	if(d > 0xefff) s+= 0xf000;
+	else if(d < 0xfff) s = 0;
+	v[i] = s;
+}
+double t0;
+void printmas(int *arr){
+	int  x, y, fst = SZ/2-5, lst = SZ/2+4;
+	for(y = fst; y < lst; ++y){
+		for(x = fst; x < lst; ++x)
+			printf("%5d ", arr[x + y*SZ]);
+		printf("\n");
+	}
+}
+void printels(int hsz){
+	int fulls = hsz*2+1;
+	printf("median %dx%d for %dx%d image: %g seconds\n", fulls, fulls, SZ, SZ, dtime()-t0);
+	printf("mid 10 values\n");
+	printf("ori:\n");
+	printmas(v);
+	printf("\nmed:\n");
+	printmas(o);
+	printf("\nmin:\n");
+	printmas(small);
+	printf("\nmax:\n");
+	printmas(large);
+	printf("\n\n");
+}
+
+int hs;
+for(hs = 1; hs < 5; ++hs){
+	t0 = dtime();
+	int blksz = hs*2+1, fullsz = blksz * blksz;
+	OMP_FOR(shared(o, v, small, large))
+	for(int x = hs; x < SZ - hs; ++x){
+		int xx, yy, xm = x+hs+1, y;
+		Mediator* m = MediatorNew(fullsz);
+		// initial fill
+		for(yy = 0; yy < blksz - 1; ++yy) for(xx = x-hs; xx < xm; ++xx)  MediatorInsert(m, v[xx+yy*SZ]);
+		for(y = hs; y < SZ - hs; ++y){
+			for(xx = x-hs; xx < xm; ++xx) MediatorInsert(m, v[xx+(y+hs)*SZ]);
+			int curpos = x+y*SZ;
+			o[curpos] = MediatorMedian(m, &small[curpos], &large[curpos]);
+			//ShowTree(m);
+			//if(y == 2) return 0;
+		}
+		free(m);
+	}
+	printels(hs);
+}
+free(o); free(v);
+	free(small); free(large);
+}

_snippets/powerof2.c

@@ -0,0 +1,3 @@
+int IsPowerOfTwo(uint64_t x){
+    return (x != 0) && ((x & (x - 1)) == 0);
+}

_snippets/pwd.c

@@ -0,0 +1,132 @@
+#include <sys/stat.h>
+#include <dirent.h>
+#include <limits.h>
+#include <stdio.h>
+#include <string.h>
+#include <unistd.h>
+#include <stdlib.h>
+#include <mntent.h>
+
+#define ROOTINO     (2)
+#define BUFSZ       (128)
+
+typedef struct _mntdef{
+	char *mntpoint;
+	struct _mntdef *next;
+} mntdef;
+
+mntdef *mntlist = NULL;
+
+void addmntdef(char *name){
+	mntdef *m, *last = NULL;
+	last = mntlist;
+	if(last){
+		do{
+			if(strcmp(last->mntpoint, name) == 0) return;
+			if(last->next) last = last->next;
+			else break;
+		}while(1);
+	}
+	m = calloc(1, sizeof(mntdef));
+	m->mntpoint = strdup(name);
+	if(!last) mntlist = m;
+	else last->next = m;
+}
+
+
+void push_dir(char **old, char* new, size_t *len, size_t *L){
+	size_t dlen = strlen(new) + 2;
+	size_t plen = strlen(*old);
+	if(plen + dlen < *L){
+		*old = realloc(*old, plen + dlen);
+		*L = plen + dlen;
+	}
+	if(len)
+		memmove(*old + dlen - 1, *old, plen+1);
+	memcpy(*old + 1, new, dlen-2);
+	if(plen == 0) *(*old + dlen - 1) = 0;
+	**old = '/';
+}
+
+struct stat original_st;
+char *chkmounted(char *path){
+	mntdef *m = mntlist;
+	char p[1024];
+	struct stat st;
+	do{
+		snprintf(p, 1023, "%s%s", m->mntpoint, path);
+		stat(p, &st);
+		if(st.st_ino == original_st.st_ino){
+			if(strcmp("/", m->mntpoint))
+				return m->mntpoint + 1;
+			else
+				return NULL;
+		}
+		if(m->next) m = m->next;
+		else break;
+	}while(1);
+	return NULL;
+}
+
+char *pwd(char **path, size_t *len){
+	size_t L = 0;
+	struct dirent *dp;
+	struct stat st;
+	stat(".", &st);
+	if(!path) return NULL;
+	if(!*path){
+		*path = malloc(BUFSZ + 1);
+		**path = 0;
+		L = BUFSZ;
+		if(len) *len = BUFSZ;
+	}else{
+		if(len) L = *len;
+		else return NULL;
+	}
+	if(st.st_ino == ROOTINO){
+		char *m = NULL;
+		if((m = chkmounted(*path))) push_dir(path, m, len, &L);
+		if(**path != '/'){
+			**path = '/';
+			*(*path+1) = 0;
+		}
+		return *path;
+	}
+	DIR *dirp = opendir("..");
+	while((dp = readdir(dirp))){
+			if(strcmp(dp->d_name, ".") == 0 || strcmp(dp->d_name, "..") == 0)
+				continue;
+			if(st.st_ino == dp->d_ino){
+					push_dir(path, dp->d_name, len, &L);
+					chdir("..");
+					return pwd(path, &L);
+			}
+	}
+	chdir("..");
+	return pwd(path, &L);
+	closedir(dirp);
+	return *path;
+}
+
+
+int main(int argc, char **argv){
+	char *path = NULL;
+	struct stat st;
+	struct mntent *m;
+	FILE *f;
+	f = setmntent(_PATH_MOUNTED, "r");
+	while((m = getmntent(f))){
+		stat(m->mnt_dir, &st);
+		addmntdef(m->mnt_dir);
+	}
+	endmntent(f);
+	if(argc == 2){
+		if(chdir(argv[1])){
+			perror("chdir()");
+			return -1;
+		}
+	}
+	stat(".", &original_st);
+	puts(pwd(&path, NULL));
+	return 0;
+}

_snippets/sendfile.c

@@ -0,0 +1,122 @@
+#include <unistd.h>
+#include <sys/sendfile.h>
+#include <fcntl.h>
+#include <sys/stat.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <sys/time.h>
+#include <sys/mman.h>
+
+/*
+./sendfile Titanik.mp4 
+Copied by sendfile, time: 30.1055s
+Copied by mmap, time: 35.6469s
+
+du Titanik.mp4 
+2.6GTitanik.mp4
+*/
+
+
+double dtime(){
+	double t;
+	struct timeval tv;
+	gettimeofday(&tv, NULL);
+	t = tv.tv_sec + ((double)tv.tv_usec)/1e6;
+	return t;
+}
+
+typedef struct{
+	char *data;
+	size_t len;
+} mmapbuf;
+
+void My_munmap(mmapbuf **b){
+	if(munmap((*b)->data, (*b)->len))
+		perror("Can't munmap");
+	free(*b);
+	*b = NULL;
+}
+
+mmapbuf *My_mmap(char *filename){
+	int fd;
+	char *ptr = NULL;
+	size_t Mlen;
+	struct stat statbuf;
+	if((fd = open(filename, O_RDONLY)) < 0){
+		perror("Can't open file for reading");
+		goto ret;
+	}
+	if(fstat (fd, &statbuf) < 0){
+		perror("Can't stat file");
+		goto ret;
+	}
+	Mlen = statbuf.st_size;
+	if((ptr = mmap (0, Mlen, PROT_READ, MAP_PRIVATE, fd, 0)) == MAP_FAILED){
+		perror("Mmap error for input");
+		goto ret;
+	}
+	mmapbuf *ret = calloc(sizeof(mmapbuf), 1);
+	if(ret){
+		ret->data = ptr;
+		ret->len = Mlen;
+	}else munmap(ptr, Mlen);
+ret:
+	if(close(fd)) perror("Can't close mmap'ed file");
+	return  ret;
+}
+
+
+int main(int argc, char **argv){
+	struct stat st;
+	double T0;
+	if(argc != 2){
+		printf("Usage: %s file\n", argv[0]);
+		return 2;
+	}
+	if (stat(argv[1], &st) == 0 && S_ISREG(st.st_mode)){
+		off_t off = 0;
+		int fd = open(argv[1], O_RDONLY);
+		int fdo = open("tmpoutput", O_WRONLY | O_CREAT,  S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH);
+		if (fd < 0 || fdo < 0){
+			perror("can't open");
+			return 1;
+		}
+		T0 = dtime();
+		posix_fadvise(fd, 0, 0, POSIX_FADV_SEQUENTIAL);
+		while (off < st.st_size) {
+			ssize_t bytes = sendfile(fdo, fd, &off, st.st_size - off);
+			if (bytes <= 0){
+				perror("can't sendfile");
+			}
+		}
+		close(fd);
+		close(fdo);
+		printf("Copied by sendfile, time: %gs\n", dtime()-T0);
+		T0 = dtime();
+		mmapbuf *map = My_mmap(argv[1]);
+		if(map){
+			fdo = open("tmpoutput", O_WRONLY | O_CREAT,  S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH);
+			if (fd < 0 || fdo < 0){
+				perror("can't open");
+				return 1;
+			}
+			size_t written = 0, towrite = map->len;
+			char *ptr = map->data;
+			do{
+				ssize_t wr = write(fdo, ptr, towrite);
+				if(wr <= 0) break;
+				written += wr;
+				towrite -= wr;
+				ptr += wr;
+			}while(towrite);
+			if(written != map->len){
+				printf("err: writed only %zd byted of %zd\n", written, map->len);
+			}
+			close(fdo);
+			My_munmap(&map);
+			printf("Copied by mmap, time: %gs\n", dtime()-T0);
+		}
+	}else
+		perror("Can't stat file");
+	return 0;
+}

_snippets/simple_Eratosfen_sieve.c

@@ -0,0 +1,57 @@
+#include <stdio.h>
+#include <stdint.h>
+#include <string.h>
+#include <stdlib.h>
+#include <omp.h>
+#include <usefull_macros.h>
+
+// 1GB of RAM
+#define MASKSZ      (1073741824L)
+#define MNUMBERS    (MASKSZ*8)
+static uint8_t *masks;
+
+static inline __attribute__((always_inline)) uint8_t get(uint64_t idx){
+    register uint64_t i = idx >> 3;
+    register uint8_t j = idx - (i<<3);
+    return masks[i] & (1<<j);
+}
+
+static inline __attribute__((always_inline)) void reset(uint64_t idx){
+    register uint64_t i = idx >> 3;
+    register uint8_t j = idx - (i<<3);
+    masks[i] &= ~(1<<j);
+}
+
+int main(){
+    masks = malloc(MASKSZ);
+    double t0 = dtime();
+    memset(masks, 0b10101010, MASKSZ); // without even numbers
+    masks[0] = 0b10101100; // with 2 and withoun 0 and 1
+    for(uint64_t i = 3; i < 64; i += 2){
+        if(get(i)){
+            for(uint64_t j = i*2; j < MNUMBERS; j += i){
+                reset(j);
+            }
+        }
+    }
+    printf("\tfirst 64: %g\n", dtime()-t0);
+#define THREADNO	8
+    omp_set_num_threads(THREADNO);
+#pragma omp parallel for shared(masks)
+    for(uint64_t i = 0; i < THREADNO; ++i){
+        for(uint64_t idx = 65 + 2*omp_get_thread_num(); idx < MNUMBERS; idx += THREADNO*2){
+        if(get(idx)){
+            for(uint64_t j = idx*2; j < MNUMBERS; j += idx){
+                reset(j);
+            }
+        }}
+    }
+    printf("\tfull table: %g\n", dtime()-t0);
+    for(uint64_t i = 0; i < 1000; ++i){
+        if(get(i)) printf("%zd\n", i);
+    }
+	for(uint64_t last = MNUMBERS-1; last > 1000; --last){
+		if(get(last)){ printf("last in list: %zd\n", last); break;}
+	}
+    return 0;
+}

_snippets/sorting.c

@@ -0,0 +1,80 @@
+#include <stdio.h>
+#include <time.h>
+#include <stdlib.h>
+
+#define NMEMB 10000000
+
+/*
+ * Timings (in seconds):
+ * -Ox    sort3a    sort3c     qsort
+ * 0     0.22787    0.24655    0.54751
+ * 1     0.034301   0.120855   0.543839
+ * 2     0.032332   0.102875   0.515567
+ * 3     0.032472   0.102684   0.514301
+ */
+
+static inline void sort3a(int *d){
+#define min(x, y) (x<y?x:y)
+#define max(x, y) (x<y?y:x)
+#define SWAP(x,y) { const int a = min(d[x], d[y]); const int b = max(d[x], d[y]); d[x] = a; d[y] = b;}
+	SWAP(0, 1);
+	SWAP(1, 2);
+	SWAP(0, 1);
+#undef SWAP
+#undef min
+#undef max
+}
+
+static inline void sort3c(int *d){
+	register int a, b;
+#define CMPSWAP(x,y) {a = d[x]; b = d[y]; if(a > b){d[x] = b; d[y] = a;}}
+	CMPSWAP(0, 1);
+	CMPSWAP(1, 2);
+	CMPSWAP(0, 1);
+#undef CMPSWAP
+}
+
+static int cmpints(const void *p1, const void *p2){
+	return (*((int*)p1) - *((int*)p2));
+}
+
+double timediff(struct timespec *tstart, struct timespec *tend){
+	return (((double)tend->tv_sec + 1.0e-9*tend->tv_nsec) -
+		   ((double)tstart->tv_sec + 1.0e-9*tstart->tv_nsec));
+}
+
+int main(){
+	double t1, t2, t3;
+	int *arra, *arrb, *arrc;
+	struct timespec tstart, tend;
+	int i, nmemb = 3*NMEMB;
+	// generate
+	i = nmemb * sizeof(int);
+	arra = malloc(i);
+	arrb = malloc(i);
+	arrc = malloc(i);
+	srand(time(NULL));
+	for(i = 0; i < nmemb; i++){
+		arra[i] = arrb[i] = arrc[i] = rand();
+	}
+	clock_gettime(CLOCK_MONOTONIC, &tstart);
+	for(i = 0; i < nmemb ; i+=3){
+		sort3a(&arra[i]);
+	}
+	clock_gettime(CLOCK_MONOTONIC, &tend);
+	t1 = timediff(&tstart, &tend);
+	clock_gettime(CLOCK_MONOTONIC, &tstart);
+	for(i = 0; i < nmemb ; i+=3){
+		sort3c(&arrb[i]);
+	}
+	clock_gettime(CLOCK_MONOTONIC, &tend);
+	t2 = timediff(&tstart, &tend);
+	clock_gettime(CLOCK_MONOTONIC, &tstart);
+	for(i = 0; i < nmemb ; i+=3){
+		qsort(&arrc[i], 3, sizeof(int), cmpints);
+	}
+	clock_gettime(CLOCK_MONOTONIC, &tend);
+	t3 = timediff(&tstart, &tend);
+	printf("%g, %g, %g; ", t1, t2, t3);
+}
+

_snippets/tmout.c

@@ -0,0 +1,114 @@
+/*
+ *                                                                                                  geany_encoding=koi8-r
+ * tmout.c
+ *
+ * Copyright 2018 Edward V. Emelianov <eddy@sao.ru, edward.emelianoff@gmail.com>
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+ * MA 02110-1301, USA.
+ *
+ */
+
+// define something else
+#define WARN(...)    do{fprintf(stderr, __VA_ARGS__);}while(0)
+#define _(arg)      arg
+
+// content of header file (tmout.h):
+#include <unistd.h>
+
+extern int verbose;
+extern char indi[];
+
+#define PRINT(...) do{if(verbose) printf(__VA_ARGS__);}while(0)
+
+#define WAIT_EVENT(evt, max_delay)  do{int __ = 0; set_timeout(max_delay); \
+        char *iptr = indi; PRINT(" "); while(!tmout && !(evt)){ \
+        usleep(100000); if(!*(++iptr)) iptr = indi; if(++__%10==0) PRINT("\b. "); \
+        PRINT("\b%c", *iptr);}; PRINT("\n");}while(0)
+
+void set_timeout(double delay);
+extern volatile int tmout;
+
+// content of c file (tmout.c):
+
+#include <stdio.h>
+#include <pthread.h>
+#include <errno.h>
+#include <sys/select.h>
+#include <signal.h>
+
+char indi[] = "|/-\\";
+volatile int tmout = 0;
+
+void *tmout_thread(void *buf){
+    int selfd = -1;
+    struct timeval *tv = (struct timeval *) buf;
+    errno = 0;
+    while(selfd < 0){
+        selfd = select(0, NULL, NULL, NULL, tv);
+        if(selfd < 0 && errno != EINTR){
+            WARN(_("Error while select()"));
+            tmout = 1;
+            return NULL;
+        }
+    }
+    tmout = 1;
+    return NULL;
+}
+
+/**
+ * run thread with pause [delay] (in seconds), at its end set variable tmout
+ */
+void set_timeout(double delay){
+    static int run = 0;
+    static pthread_t athread;
+    static struct timeval tv; // should be static to send this as argument of tmout_thread
+    if(delay < 0.){
+        tmout = 1;
+        return;
+    }
+    if(run && (pthread_kill(athread, 0) != ESRCH)){ // another timeout process detected - kill it
+        pthread_cancel(athread);
+        pthread_join(athread, NULL);
+    }
+    tmout = 0;
+    run = 1;
+    tv.tv_sec  = (time_t) delay;
+    tv.tv_usec = (suseconds_t)((delay-(double)tv.tv_sec)*1e6);
+    if(pthread_create(&athread, NULL, tmout_thread, (void*)&tv)){
+        WARN(_("Can't create timeout thread!"));
+        tmout = 1;
+        return;
+    }
+}
+
+// main.c
+#include <string.h>
+
+int verbose = 0;
+
+int main(int argc, char **argv){
+    // check `verbose`
+    if(argc == 2 && strcmp(argv[1], "-v") == 0){
+        printf("Will show indicator\n");
+        verbose = 1;
+    }
+    setbuf(stdout, NULL); // without this string there will be no indicator!
+    printf("wait for 2.574 seconds\n");
+    WAIT_EVENT(0, 2.574);
+    if(tmout) printf("Timeout\n");
+    else printf("Never reached\n");
+    return 0;
+}

_snippets/translateCintoASM

@@ -0,0 +1 @@
+translate c into asm: gcc -S -fverbose-asm file.c

_snippets/xor.c

@@ -0,0 +1,33 @@
+#include <stdio.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
+#include <stdlib.h>
+
+void xor(char *buf, char *pattern, size_t l){
+	size_t i;
+	for(i = 0; i < l; ++i)
+		printf("%c", buf[i] ^ pattern[i]);
+}
+
+int main(int argc, char **argv){
+	if(argc != 3){
+		printf("Usage: %s <infile> <password>\n", argv[0]);
+		return -1;
+	}
+	char *key = strdup(argv[2]);
+	size_t n = 0, keylen = strlen(key);
+	char *buff = malloc(keylen);
+	int f = open(argv[1], O_RDONLY);
+	if(f < 0){
+		perror("Can't open file");
+		return f;
+	}
+	do{
+		n = read(f, buff, keylen);
+		xor(buff, key, n);
+	}while(n);
+	return 0;
+}