[JT] Re: [Comm] Re: MathCAD, MathLab

[Michael Shigorin]

On Thu, Nov 03, 2005 at 09:43:30AM -0000, Aleksander N. Gorohovski wrote:
> >>>>У него свой язык, типа программирования, но, по моему, лучше
> >>>>уж писать на С.
> >>>Математику?  Ну-ну.
> >>Сурьезную математику, а не беллетристику.
> >Сурьёзную математику на Це писать сильно дорого выходит.
> >Разберёте такую мелочь, как аттачик (в смысле "что оно делает")
> Миша, я до конца не понял ваш язык изопа?

Вообще-то это был язык C. :]

> What is "аттачик"?

an attach.  Цепляю ещё раз (хотя и тогда вроде не забыл?).

> :-)

-- 
 ---- WBR, Michael Shigorin <mike на altlinux.ru>
  ------ Linux.Kiev http://www.linux.kiev.ua/
----------- следующая часть -----------
/*
 * ALCANOL -- the alcanol isomeres number calculations
 * Version 0.20 of May 10 2000 by Michael Shigorin <michael.shigorin на usa.net>
 *
 * THIS IS FREE SOFTWARE - USE FREELY BUT ON YOUR OWN RISK (CISK;)
 * MAY BE REDISTRIBUTED/MODIFIED UNDER THE TERMS OF GNU GPL ][
 *
 */
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <stdio.h>
#include <string.h>
#include "gmp.h"

#define counter_t unsigned int		/* used for counters -- 32bit for me */

//#define TIMER				/* timing support */
//#define DEBUG				/* more verbosity */
#define MAX_DEGREE	4096		/* CPU/RAM/time's the limit :) */
#define MAX_LEN		255		/* max length of iofile name */
#define FORMAT		"%d:\t"		/* ^data file format excl. HUGE numbers$ */
#ifdef TIMER
  #define LOGFORMAT	"%d\t%f"	/* number <tab> seconds */
  #define TIMEXT	".log"		/* extension to add for log */
  #include <sys/time.h>
  #include <unistd.h>
#endif

mpz_t source [MAX_DEGREE], factor, buffer;
mpz_t temp;

FILE *iofile;
char filename [MAX_LEN];
counter_t max;

#ifdef TIMER
  char logname [MAX_LEN];
  FILE *logfile;
  float timer;
  struct timeval raw_timer;
  struct timezone tz;		/* stub */

  float alk_time (void);
#endif

void init (void);
void parse (void);
void calculate (void);
void print (void);
void usage (void);
void shutdown (void);

int main (int argc, char *argv[]) {
	if( argc != 2 ) usage ();
	strncpy (filename, argv [1], strlen (argv [1])+1);

#ifdef TIMER
	snprintf (logname, strlen (filename)+strlen (TIMEXT)+1, "%s%s", filename, TIMEXT);
#endif
	init ();
	parse ();
	calculate ();
	print ();
	shutdown ();
	return (0);
}

#ifdef TIMER
float alk_time (void) {
	gettimeofday (&raw_timer, &tz);
	return (float)raw_timer.tv_sec + (float)raw_timer.tv_usec/1000000;
}
#endif

void usage (void) {
	fprintf (stderr, "I need one argument (a file to proceed)!\n");
	exit (1);
}

void init (void) {
#ifdef DEBUG
	fprintf (stderr, "\nReading data from <%s>...\n", filename);
#endif

	if( !(iofile = fopen (filename, "r")) ) {
		fprintf (stderr, "Error: cannot Open Source file <%s>\n", filename);
		exit (1);
	}
}

void parse (void) {
	mpz_init (temp);

#ifdef DEBUG
	fprintf (stderr, "Parsing...");
#endif
	while( !feof(iofile) ) {
		if( fscanf(iofile, FORMAT, &max) ) {
			if( mpz_inp_str (temp, iofile, 10) ) mpz_set (source [max], temp);
		} else {
			fprintf (stderr, "Error: malformed line #%d", max);
			exit (2);
		}
	}
#ifdef DEBUG
	fprintf (stderr, " ok, maximal degree is %d.\n", max);
#endif
}

void calculate (void) {
	counter_t i, j;

#ifdef TIMER
	timer = alk_time ();
#endif
#ifdef DEBUG
	fprintf (stderr, "Calculating factor by degree %d.\n", max+1);
#endif

	mpz_init_set_ui (factor, 0);
	/* first part -- C^3(x). Take only those parts of degree == max */
	for( i=0; i<=max; i++ )
		for( j=0; j<=max-i; j++ ) {
			mpz_mul (temp, source [i], source [j]);
			mpz_mul (temp, source [max-i-j], temp);
			mpz_add (factor, temp, factor);
		}
	/* next part -- 3*C(x)*C(x^2). A bit tricky... */
	for( i=max%2; i<=max; i+=2 ) {			/* evens and odds */
		mpz_mul_ui (temp, source [i], 3);
		mpz_mul (temp, source [(max-i)>>1], temp);
		mpz_add (factor, temp, factor);
	}
	/* last part -- 2*C(x^3). */
	if( !(max % 3) ) {
		mpz_mul_ui (temp, source [max/3], 2);
		mpz_add (factor, temp, factor);
	}
	/* and finally, divide factor by 6. If it doesn't........ */
	mpz_tdiv_qr_ui (factor, temp, factor, 6);
	if( mpz_cmp_ui (temp, 0) ) {
		fprintf (stderr, "INTERNAL ERROR OR BAD DATA!\n");
		fprintf (stderr, " is not divisible by 6 (#%d)\n", max+1);
		exit (3);				/* ...... */
	}
#ifdef TIMER
	timer = alk_time () - timer;
#endif
}

void print (void) {
	if( !(iofile = fopen (filename, "a")) ) {
		fprintf (stderr, "Error: cannot open output file <%s>\n", filename);
		fclose (iofile);
		exit (1);
	}

#ifdef TIMER
	if( !(logfile = fopen (logname, "a")) ) {
		fprintf (stderr, "Warning: cannot open log file <%s>\n", logname);
	} else {
		fprintf (logfile, LOGFORMAT, max+1, timer);
		putc ('\n', logfile);
	}
#endif
#ifdef DEBUG
	fprintf (stderr, "Writing to <%s>.\n", filename);
#endif

	fprintf (iofile, FORMAT, max+1);
	mpz_out_str (iofile, 10, factor);
	putc ('\n', iofile);
}

void shutdown (void) {
	fclose (iofile);
#ifdef TIMER
	fclose (logfile);
#endif
}