--- /dev/null
+/* The Computer Language Shootout Benchmarks
+ http://shootout.alioth.debian.org/
+
+ contributed by Kevin Carson
+ compilation:
+ gcc -O3 -fomit-frame-pointer -funroll-loops -static binary-trees.c -lm
+ icc -O3 -ip -unroll -static binary-trees.c -lm
+*/
+
+#include <malloc.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+
+typedef struct tn {
+ struct tn* left;
+ struct tn* right;
+ long item;
+} treeNode;
+
+
+treeNode* NewTreeNode(treeNode* left, treeNode* right, long item)
+{
+ treeNode* new;
+
+ new = (treeNode*)malloc(sizeof(treeNode));
+
+ new->left = left;
+ new->right = right;
+ new->item = item;
+
+ return new;
+} /* NewTreeNode() */
+
+
+long ItemCheck(treeNode* tree)
+{
+ if (tree->left == NULL)
+ return tree->item;
+ else
+ return tree->item + ItemCheck(tree->left) - ItemCheck(tree->right);
+} /* ItemCheck() */
+
+
+treeNode* BottomUpTree(long item, unsigned depth)
+{
+ if (depth > 0)
+ return NewTreeNode
+ (
+ BottomUpTree(2 * item - 1, depth - 1),
+ BottomUpTree(2 * item, depth - 1),
+ item
+ );
+ else
+ return NewTreeNode(NULL, NULL, item);
+} /* BottomUpTree() */
+
+
+void DeleteTree(treeNode* tree)
+{
+ if (tree->left != NULL)
+ {
+ DeleteTree(tree->left);
+ DeleteTree(tree->right);
+ }
+
+ free(tree);
+} /* DeleteTree() */
+
+
+int main(int argc, char* argv[])
+{
+ unsigned N, depth, minDepth, maxDepth, stretchDepth;
+ treeNode *stretchTree, *longLivedTree, *tempTree;
+
+ if(argc < 2) {
+ N = 10;
+ printf("Using default input %d\n", N);
+ } else {
+ N = atol(argv[1]);
+ }
+
+ minDepth = 4;
+
+ if ((minDepth + 2) > N)
+ maxDepth = minDepth + 2;
+ else
+ maxDepth = N;
+
+ stretchDepth = maxDepth + 1;
+
+ stretchTree = BottomUpTree(0, stretchDepth);
+ printf
+ (
+ "stretch tree of depth %u\t check: %li\n",
+ stretchDepth,
+ ItemCheck(stretchTree)
+ );
+
+ DeleteTree(stretchTree);
+
+ longLivedTree = BottomUpTree(0, maxDepth);
+
+ for (depth = minDepth; depth <= maxDepth; depth += 2)
+ {
+ long i, iterations, check;
+
+ iterations = pow(2, maxDepth - depth + minDepth);
+
+ check = 0;
+
+ for (i = 1; i <= iterations; i++)
+ {
+ tempTree = BottomUpTree(i, depth);
+ check += ItemCheck(tempTree);
+ DeleteTree(tempTree);
+
+ tempTree = BottomUpTree(-i, depth);
+ check += ItemCheck(tempTree);
+ DeleteTree(tempTree);
+ } /* for(i = 1...) */
+
+ printf
+ (
+ "%li\t trees of depth %u\t check: %li\n",
+ iterations * 2,
+ depth,
+ check
+ );
+ } /* for(depth = minDepth...) */
+
+ printf
+ (
+ "long lived tree of depth %u\t check: %li\n",
+ maxDepth,
+ ItemCheck(longLivedTree)
+ );
+
+ return 0;
+} /* main() */
--- /dev/null
+/*
+ * The Computer Lannguage Shootout
+ * http://shootout.alioth.debian.org/
+ * Contributed by Heiner Marxen
+ *
+ * "fannkuch" for C gcc
+ *
+ * $Id$
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#define Int int
+#define Aint int
+
+ static long
+fannkuch( int n )
+{
+ Aint* perm;
+ Aint* perm1;
+ Aint* count;
+ long flips;
+ long flipsMax;
+ Int r;
+ Int i;
+ Int k;
+ Int didpr;
+ const Int n1 = n - 1;
+
+ if( n < 1 ) return 0;
+
+ perm = calloc(n, sizeof(*perm ));
+ perm1 = calloc(n, sizeof(*perm1));
+ count = calloc(n, sizeof(*count));
+
+ for( i=0 ; i<n ; ++i ) perm1[i] = i; /* initial (trivial) permu */
+
+ r = n; didpr = 0; flipsMax = 0;
+ for(;;) {
+ if( didpr < 30 ) {
+ for( i=0 ; i<n ; ++i ) printf("%d", (int)(1+perm1[i]));
+ printf("\n");
+ ++didpr;
+ }
+ for( ; r!=1 ; --r ) {
+ count[r-1] = r;
+ }
+
+#define XCH(x,y) { Aint t_mp; t_mp=(x); (x)=(y); (y)=t_mp; }
+
+ if( ! (perm1[0]==0 || perm1[n1]==n1) ) {
+ flips = 0;
+ for( i=1 ; i<n ; ++i ) { /* perm = perm1 */
+ perm[i] = perm1[i];
+ }
+ k = perm1[0]; /* cache perm[0] in k */
+ do { /* k!=0 ==> k>0 */
+ Int j;
+ for( i=1, j=k-1 ; i<j ; ++i, --j ) {
+ XCH(perm[i], perm[j])
+ }
+ ++flips;
+ /*
+ * Now exchange k (caching perm[0]) and perm[k]... with care!
+ * XCH(k, perm[k]) does NOT work!
+ */
+ j=perm[k]; perm[k]=k ; k=j;
+ }while( k );
+ if( flipsMax < flips ) {
+ flipsMax = flips;
+ }
+ }
+
+ for(;;) {
+ if( r == n ) {
+ return flipsMax;
+ }
+ /* rotate down perm[0..r] by one */
+ {
+ Int perm0 = perm1[0];
+ i = 0;
+ while( i < r ) {
+ k = i+1;
+ perm1[i] = perm1[k];
+ i = k;
+ }
+ perm1[r] = perm0;
+ }
+ if( (count[r] -= 1) > 0 ) {
+ break;
+ }
+ ++r;
+ }
+ }
+}
+
+int
+main( int argc, char* argv[] )
+{
+ int n = (argc>1) ? atoi(argv[1]) : 10;
+
+ printf("Pfannkuchen(%d) = %ld\n", n, fannkuch(n));
+ return 0;
+}
--- /dev/null
+/* -*- mode: c -*-
+ * $Id$
+ * http://shootout.alioth.debian.org/
+ *
+ * by Paul Hsieh
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#define IM 139968
+#define IA 3877
+#define IC 29573
+
+double gen_random (double max) {
+ static long last = 42;
+ return max * (last = (last * IA + IC) % IM) / IM;
+}
+
+struct aminoacids {
+ char c;
+ double p;
+};
+
+/* Weighted selection from alphabet */
+
+void makeCumulative (struct aminoacids * genelist, int count) {
+ double cp = 0.0;
+ int i;
+
+ for (i=0; i < count; i++) {
+ cp += genelist[i].p;
+ genelist[i].p = cp;
+ }
+}
+
+char selectRandom (const struct aminoacids * genelist, int count) {
+ double r = gen_random (1);
+ int i, lo, hi;
+
+ if (r < genelist[0].p) return genelist[0].c;
+
+ lo = 0;
+ hi = count-1;
+
+ while (hi > lo+1) {
+ i = (hi + lo) / 2;
+ if (r < genelist[i].p) hi = i; else lo = i;
+ }
+ return genelist[hi].c;
+}
+
+/* Generate and write FASTA format */
+
+#define LINE_LENGTH (60)
+
+void makeRandomFasta (const char * id, const char * desc, const struct
+aminoacids * genelist, int count, int n) {
+ int todo = n;
+ int i, m;
+
+ printf (">%s %s\n", id, desc);
+
+ for (; todo > 0; todo -= LINE_LENGTH) {
+ char pick[LINE_LENGTH+1];
+ if (todo < LINE_LENGTH) m = todo; else m = LINE_LENGTH;
+ for (i=0; i < m; i++) pick[i] = selectRandom (genelist, count);
+ pick[m] = '\0';
+ puts (pick);
+ }
+}
+
+void makeRepeatFasta (const char * id, const char * desc, const char *
+s, int n) {
+ char * ss;
+ int todo = n, k = 0, kn = strlen (s);
+ int m;
+
+ ss = (char *) malloc (kn + 1);
+ memcpy (ss, s, kn+1);
+
+ printf (">%s %s\n", id, desc);
+
+ for (; todo > 0; todo -= LINE_LENGTH) {
+ if (todo < LINE_LENGTH) m = todo; else m = LINE_LENGTH;
+
+ while (m >= kn - k) {
+ printf ("%s", s+k);
+ m -= kn - k;
+ k = 0;
+ }
+
+ ss[k + m] = '\0';
+ puts (ss+k);
+ ss[k + m] = s[m+k];
+ k += m;
+ }
+
+ free (ss);
+}
+
+/* Main -- define alphabets, make 3 fragments */
+
+struct aminoacids iub[] = {
+ { 'a', 0.27 },
+ { 'c', 0.12 },
+ { 'g', 0.12 },
+ { 't', 0.27 },
+
+ { 'B', 0.02 },
+ { 'D', 0.02 },
+ { 'H', 0.02 },
+ { 'K', 0.02 },
+ { 'M', 0.02 },
+ { 'N', 0.02 },
+ { 'R', 0.02 },
+ { 'S', 0.02 },
+ { 'V', 0.02 },
+ { 'W', 0.02 },
+ { 'Y', 0.02 }
+};
+
+#define IUB_LEN (sizeof (iub) / sizeof (struct aminoacids))
+
+struct aminoacids homosapiens[] = {
+ { 'a', 0.3029549426680 },
+ { 'c', 0.1979883004921 },
+ { 'g', 0.1975473066391 },
+ { 't', 0.3015094502008 },
+};
+
+#define HOMOSAPIENS_LEN (sizeof (homosapiens) / sizeof (struct aminoacids))
+
+char * alu =
+ "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG"
+ "GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA"
+ "CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT"
+ "ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA"
+ "GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG"
+ "AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC"
+ "AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA";
+
+int main (int argc, char * argv[]) {
+ int n = 1000;
+
+ if (argc > 1) sscanf (argv[1], "%d", &n);
+
+ makeCumulative (iub, IUB_LEN);
+ makeCumulative (homosapiens, HOMOSAPIENS_LEN);
+
+ makeRepeatFasta ("ONE", "Homo sapiens alu", alu, n*2);
+ makeRandomFasta ("TWO", "IUB ambiguity codes", iub, IUB_LEN, n*3);
+ makeRandomFasta ("THREE", "Homo sapiens frequency", homosapiens,
+HOMOSAPIENS_LEN, n*5);
+
+ return 0;
+}
--- /dev/null
+/*
+ * The Great Computer Language Shootout
+ * http://shootout.alioth.debian.org/
+ *
+ * contributed by Christoph Bauer
+ *
+ */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#define pi 3.141592653589793
+#define solar_mass (4 * pi * pi)
+#define days_per_year 365.24
+
+struct planet {
+ double x, y, z;
+ double vx, vy, vz;
+ double mass;
+};
+
+void advance(int nbodies, struct planet * bodies, double dt)
+{
+ int i, j;
+
+ for (i = 0; i < nbodies; i++) {
+ struct planet * b = &(bodies[i]);
+ for (j = i + 1; j < nbodies; j++) {
+ struct planet * b2 = &(bodies[j]);
+ double dx = b->x - b2->x;
+ double dy = b->y - b2->y;
+ double dz = b->z - b2->z;
+ double distance = sqrt(dx * dx + dy * dy + dz * dz);
+ double mag = dt / (distance * distance * distance);
+ b->vx -= dx * b2->mass * mag;
+ b->vy -= dy * b2->mass * mag;
+ b->vz -= dz * b2->mass * mag;
+ b2->vx += dx * b->mass * mag;
+ b2->vy += dy * b->mass * mag;
+ b2->vz += dz * b->mass * mag;
+ }
+ }
+ for (i = 0; i < nbodies; i++) {
+ struct planet * b = &(bodies[i]);
+ b->x += dt * b->vx;
+ b->y += dt * b->vy;
+ b->z += dt * b->vz;
+ }
+}
+
+double energy(int nbodies, struct planet * bodies)
+{
+ double e;
+ int i, j;
+
+ e = 0.0;
+ for (i = 0; i < nbodies; i++) {
+ struct planet * b = &(bodies[i]);
+ e += 0.5 * b->mass * (b->vx * b->vx + b->vy * b->vy + b->vz * b->vz);
+ for (j = i + 1; j < nbodies; j++) {
+ struct planet * b2 = &(bodies[j]);
+ double dx = b->x - b2->x;
+ double dy = b->y - b2->y;
+ double dz = b->z - b2->z;
+ double distance = sqrt(dx * dx + dy * dy + dz * dz);
+ e -= (b->mass * b2->mass) / distance;
+ }
+ }
+ return e;
+}
+
+void offset_momentum(int nbodies, struct planet * bodies)
+{
+ double px = 0.0, py = 0.0, pz = 0.0;
+ int i;
+ for (i = 0; i < nbodies; i++) {
+ px += bodies[i].vx * bodies[i].mass;
+ py += bodies[i].vy * bodies[i].mass;
+ pz += bodies[i].vz * bodies[i].mass;
+ }
+ bodies[0].vx = - px / solar_mass;
+ bodies[0].vy = - py / solar_mass;
+ bodies[0].vz = - pz / solar_mass;
+}
+
+#define NBODIES 5
+struct planet bodies[NBODIES] = {
+ { /* sun */
+ 0, 0, 0, 0, 0, 0, solar_mass
+ },
+ { /* jupiter */
+ 4.84143144246472090e+00,
+ -1.16032004402742839e+00,
+ -1.03622044471123109e-01,
+ 1.66007664274403694e-03 * days_per_year,
+ 7.69901118419740425e-03 * days_per_year,
+ -6.90460016972063023e-05 * days_per_year,
+ 9.54791938424326609e-04 * solar_mass
+ },
+ { /* saturn */
+ 8.34336671824457987e+00,
+ 4.12479856412430479e+00,
+ -4.03523417114321381e-01,
+ -2.76742510726862411e-03 * days_per_year,
+ 4.99852801234917238e-03 * days_per_year,
+ 2.30417297573763929e-05 * days_per_year,
+ 2.85885980666130812e-04 * solar_mass
+ },
+ { /* uranus */
+ 1.28943695621391310e+01,
+ -1.51111514016986312e+01,
+ -2.23307578892655734e-01,
+ 2.96460137564761618e-03 * days_per_year,
+ 2.37847173959480950e-03 * days_per_year,
+ -2.96589568540237556e-05 * days_per_year,
+ 4.36624404335156298e-05 * solar_mass
+ },
+ { /* neptune */
+ 1.53796971148509165e+01,
+ -2.59193146099879641e+01,
+ 1.79258772950371181e-01,
+ 2.68067772490389322e-03 * days_per_year,
+ 1.62824170038242295e-03 * days_per_year,
+ -9.51592254519715870e-05 * days_per_year,
+ 5.15138902046611451e-05 * solar_mass
+ }
+};
+
+int main(int argc, char ** argv)
+{
+ int n = atoi(argv[1]);
+ int i;
+
+ offset_momentum(NBODIES, bodies);
+ printf ("%.9f\n", energy(NBODIES, bodies));
+ for (i = 1; i <= n; i++)
+ advance(NBODIES, bodies, 0.01);
+ printf ("%.9f\n", energy(NBODIES, bodies));
+ return 0;
+}
--- /dev/null
+/*
+ * The Great Computer Language Shootout
+ * http://shootout.alioth.debian.org/
+ *
+ * Written by Dima Dorfman, 2004
+ * Compile: gcc -std=c99 -O2 -o nsieve_bits_gcc nsieve_bits.c
+ */
+
+#include <limits.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <string.h>
+
+typedef uint_fast8_t bits;
+#define NBITS (CHAR_BIT * sizeof(bits))
+
+static uintmax_t
+nsieve(uintmax_t m)
+{
+ uintmax_t count, i, j;
+ bits a[m / NBITS];
+
+ memset(a, (1 << CHAR_BIT) - 1, sizeof(a));
+ count = 0;
+ for (i = 2; i < m; ++i)
+ if (a[i / NBITS] & (1 << i % NBITS)) {
+ for (j = i + i; j < m; j += i)
+ a[j / NBITS] &= ~(1 << j % NBITS);
+ ++count;
+ }
+ return (count);
+}
+
+static void
+test(unsigned long n)
+{
+ uintmax_t count, m;
+
+ m = (1 << n) * 10000;
+ count = nsieve(m);
+ printf("Primes up to %8ju %8ju\n", m, count);
+}
+
+int
+main(int ac, char **av)
+{
+ unsigned long n;
+ char *cp;
+
+ if (ac < 2) {
+usage: fprintf(stderr, "usage: nsieve N\n");
+ exit(2);
+ }
+ n = strtoul(av[1], &cp, 10);
+ if (*av[1] == '\0' || *cp != '\0' || n == ULONG_MAX)
+ goto usage;
+ test(n);
+ if (n >= 1)
+ test(n - 1);
+ if (n >= 2)
+ test(n - 2);
+ exit(0);
+}
--- /dev/null
+// The Computer Language Shootout
+// http://shootout.alioth.debian.org/
+// Precedent C entry modified by bearophile for speed and size, 31 Jan 2006
+// Compile with: -O3 -s -std=c99 -fomit-frame-pointer
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+typedef unsigned char boolean;
+
+
+static void nsieve(int m) {
+ unsigned int count = 0, i, j;
+ boolean * flags = (boolean *) malloc(m * sizeof(boolean));
+ memset(flags, 1, m);
+
+ for (i = 2; i < m; ++i)
+ if (flags[i]) {
+ ++count;
+ for (j = i << 1; j < m; j += i)
+ if (flags[j]) flags[j] = 0;
+ }
+
+ free(flags);
+ printf("Primes up to %8u %8u\n", m, count);
+}
+
+int main(int argc, char * argv[]) {
+ int m = atoi(argv[1]);
+ for (int i = 0; i < 3; i++)
+ nsieve(10000 << (m-i));
+ return 0;
+}
--- /dev/null
+/*
+** The Computer Language Shootout
+** http://shootout.alioth.debian.org/
+** contributed by Mike Pall
+** de-optimized by Isaac Gouy
+**
+** compile with:
+** gcc -O3 -fomit-frame-pointer -ffast-math -o partialsums partialsums.c -lm
+** Adding -march=<yourcpu> may help, too.
+** On a P4/K8 or later try adding: --march=<yourcpu> -mfpmath=sse -msse2
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+int main(int argc, char **argv)
+{
+ int k;
+ int n = 100000;
+ double sum, a;
+
+ if(argc > 1)
+ n = atoi(argv[1]);
+
+/*
+** Yes, I (Mike Pall) tried using a double as a primary or secondary loop variable.
+** But the x86 ABI requires a cleared x87 FPU stack before every call
+** (e.g. to sin()) which nullifies any performance gains.
+**
+** Combining all loops does not pay off because the x87 FPU has to shuffle
+** stack slots and/or runs out of registers. This may not be entirely true
+** for SSE2 with fully inlined FPU code (-ffast-math required). Dito for
+** other CPUs with a register-based FPU and a sane FP ABI.
+**
+** Auto vectorization may be a bit easier with separate loops, too.
+*/
+
+#define kd ((double)k)
+
+ sum = 0.0;
+ for (k = 0; k <= n; k++) sum += pow(2.0/3.0, kd);
+ printf("%.9f\t(2/3)^k\n", sum);
+
+ sum = 0.0;
+ for (k = 1 ; k <= n; k++) sum += 1/sqrt(kd); /* aka pow(kd, -0.5) */
+ printf("%.9f\tk^-0.5\n", sum);
+
+ sum = 0.0;
+ for (k = 1; k <= n; k++) sum += 1.0/(kd*(kd+1.0));
+ printf("%.9f\t1/k(k+1)\n", sum);
+
+ sum = 0.0;
+ for (k = 1; k <= n; k++) {
+ double sk = sin(kd);
+ sum += 1.0/(kd*kd*kd*sk*sk);
+ }
+ printf("%.9f\tFlint Hills\n", sum);
+
+ sum = 0.0;
+ for (k = 1; k <= n; k++) {
+ double ck = cos(kd);
+ sum += 1.0/(kd*kd*kd*ck*ck);
+ }
+ printf("%.9f\tCookson Hills\n", sum);
+
+ sum = 0.0;
+ for (k = 1; k <= n; k++) sum += 1.0/kd;
+ printf("%.9f\tHarmonic\n", sum);
+
+ sum = 0.0;
+ for (k = 1; k <= n; k++) sum += 1.0/(kd*kd);
+ printf("%.9f\tRiemann Zeta\n", sum);
+
+ sum = 0.0; a = -1.0;
+ for (k = 1; k <= n; k++) sum += (a = -a)/kd;
+ printf("%.9f\tAlternating Harmonic\n", sum);
+
+ sum = 0.0; a = -1.0;
+ for (k = 1; k <= n; k++) sum += (a = -a)/(2.0*kd - 1.0);
+ printf("%.9f\tGregory\n", sum);
+
+ return 0;
+}
--- /dev/null
+// The Computer Language Shootout
+// http://shootout.alioth.debian.org/
+// recursive test, by bearophile, Jan 24 2006
+// Compile with: -O3 -s -fomit-frame-pointer -funroll-loops
+
+#include <stdio.h>
+
+int Ack(int x, int y) {
+ if (x == 0)
+ return y+1;
+ if (y == 0)
+ return Ack(x-1, 1);
+ return Ack(x-1, Ack(x, y-1));
+}
+
+int Fib(int n) {
+ if (n < 2)
+ return 1;
+ return Fib(n-2) + Fib(n-1);
+}
+
+double FibFP(double n) {
+ if (n < 2.0)
+ return 1.0;
+ return FibFP(n-2.0) + FibFP(n-1.0);
+}
+
+int Tak(int x, int y, int z) {
+ if (y < x)
+ return Tak( Tak(x-1, y, z), Tak(y-1, z, x), Tak(z-1, x, y) );
+ return z;
+}
+
+double TakFP(double x, double y, double z) {
+ if (y < x)
+ return TakFP( TakFP(x-1.0, y, z), TakFP(y-1.0, z, x), TakFP(z-1.0, x, y) );
+ return z;
+}
+
+int main(int argc, char **argv) {
+ int n = 7;
+
+ if(argc > 1)
+ n = atoi(argv[1]) - 1;
+ printf("Ack(3,%d): %d\n", n+1, Ack(3, n+1));
+ printf("Fib(%.1f): %.1f\n", 28.0+n, FibFP(28.0+n));
+ printf("Tak(%d,%d,%d): %d\n", 3*n, 2*n, n, Tak(3*n, 2*n, n));
+ printf("Fib(3): %d\n", Fib(3));
+ printf("Tak(3.0,2.0,1.0): %.1f\n", TakFP(3.0, 2.0, 1.0));
+ return 0;
+}
--- /dev/null
+/* -*- mode: c -*-
+ *
+ * The Great Computer Language Shootout
+ * http://shootout.alioth.debian.org/
+ *
+ * Contributed by Sebastien Loisel
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+double eval_A(int i, int j) { return 1.0/((i+j)*(i+j+1)/2+i+1); }
+
+void eval_A_times_u(int N, const double u[], double Au[])
+{
+ int i,j;
+ for(i=0;i<N;i++)
+ {
+ Au[i]=0;
+ for(j=0;j<N;j++) Au[i]+=eval_A(i,j)*u[j];
+ }
+}
+
+void eval_At_times_u(int N, const double u[], double Au[])
+{
+ int i,j;
+ for(i=0;i<N;i++) {
+ Au[i]=0;
+ for(j=0;j<N;j++)
+ Au[i] += eval_A(j,i)*u[j];
+ }
+}
+
+void eval_AtA_times_u(int N, const double u[], double AtAu[])
+{
+ double *v[N];
+ eval_A_times_u(N,u,v);
+ eval_At_times_u(N,v,AtAu);
+}
+
+int main(int argc, char *argv[])
+{
+ int i;
+ int N = ((argc == 2) ? atoi(argv[1]) : 2000);
+
+ double *u = alloca(sizeof(u[0]) * N);
+ double *v = alloca(sizeof(v[0]) * N);
+ double vBv, vv;
+
+ for(i=0;i<N;i++)
+ u[i]=1;
+
+ for(i=0;i<10;i++) {
+ eval_AtA_times_u(N,u,v);
+ eval_AtA_times_u(N,v,u);
+ }
+
+ vBv = vv = 0;
+
+ for(i=0;i<N;i++) {
+ vBv+=u[i]*v[i];
+ vv+=v[i]*v[i];
+ }
+
+ printf("%0.9f\n", sqrt(vBv/vv));
+ return 0;
+}
projects / libfirm / commitdiff