/*                                                    B Bunk 4/2001
  Tabelle von Quadrat- und Kubikwurzeln berechnen     rev    3/2012
      und in File schreiben
*/
#include <stdlib.h>                 /* fuer exit() */
#include <stdio.h>                  /* fuer printf, scanf, fopen etc */
#include <math.h>                   /* nicht noetig */

#define NMAX 100
static float tol = 1.e-6;

float w2(float z);                  /* Prototypen der weiter unten */
void sub3(float z, float *w3);      /* definierten Unterprogramme */

int main(void){

  float tabelle[NMAX][3];       /* Index-Bereiche 0:(NMAX-1), 0:2 */
  float x0, dx, x, w3;
  int   n, i;
  FILE  *tabfile;

  printf("x0 dx n? ");                /* Parameter abfragen */
  scanf("%f%f%i", &x0, &dx, &n);      /* und einlesen */

  if (n > NMAX){
    printf("Fehler: Anzahl zu gross fuer Tabelle\n");
    exit(1);
  }

  for (i = 0; i < n; i++){       /* Schleife ueber Tabellenzeilen */
    x = x0 + i*dx;

    tabelle[i][0] = x;           /* Tabellenzeile berechnen */
    tabelle[i][1] = w2(x);       /* Funktionsaufruf */

    sub3(x, &w3);                /* Aufruf von Unterprogramm */
                                 /* Adresse von w3 uebergeben */
                                 /* fuer Rueckgabe des Ergebnisses */
    tabelle[i][2] = w3;

    /* Kontrollausdruck */
    printf("%10f %10f %10f\n",
           tabelle[i][0], tabelle[i][1], tabelle[i][2]);
  }

  tabfile = fopen("wurzeln.tab","w"); /* File zum Schreiben oeffnen */
  for (i = 0; i < n; i++){
    /* Zeile schreiben */
    fprintf(tabfile,"%10f %10f %10f\n",
            tabelle[i][0], tabelle[i][1], tabelle[i][2]);
  }
  fclose(tabfile);                    /* File schliessen */
  exit(0);
}

/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */

float w2(float z){

  /* w2 <- sqrt(z) berechnen */

  float x, xalt;

  if (z < 0.){                      /* Argument pruefen */
    printf("Fehler in w2: negatives Argument\n");
    exit(2);
  }

  xalt = 0.;
  x = 1.;                           /* Startwert */
  while (fabs(x-xalt) > tol){       /* Konvergenztest */
    xalt = x;
    x = .5 * (x + z/x);             /* Iteration */
  }
  return x;
}

/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */

void sub3(float z, float *w3){

  /* w3 <- z**(1/3) berechnen */

  float w, walt;

  if (z < 0.){                      /* Argument pruefen */
    printf("Fehler in sub3: negatives Argument\n");
    exit(3);
  }

  w = 1.;                           /* Startwert */
  do    {
    walt = w;
    w = (2*w + z/(w*w))/3;          /* Iteration */
  }
  while (fabs(w-walt) > tol);       /* Konvergenztest */
  *w3 = w;
  return;
}