1 #ifndef NUMSUP_H
2 
3 /* Numerical routine general support declarations */
4 
5 /*
6  * Copyright 2000-2010 Graeme W. Gill
7  * All rights reserved.
8  *
9  * This material is licenced under the GNU GENERAL PUBLIC LICENSE Version 2 or later :-
10  * see the License2.txt file for licencing details.
11  */
12 
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <stdarg.h>
16 #include <string.h>
17 #include <float.h>
18 #include <math.h>
19 
20 #ifdef NT
21 #include <basetsd.h>		/* So jpg header doesn't define INT32 */
22 # if !defined(_WIN32_WINNT) || _WIN32_WINNT < 0x0501
23 #  undef _WIN32_WINNT
24 #  define _WIN32_WINNT 0x0501
25 # endif
26 # define WIN32_LEAN_AND_MEAN
27 # include <windows.h>
28 #endif
29 #ifdef UNIX
30 # include <pthread.h>
31 #endif
32 
33 #ifdef __cplusplus
34 	extern "C" {
35 #endif
36 
37 /* =========================================================== */
38 /* Platform specific primitive defines. */
39 /* This really needs checking for each different platform. */
40 /* Using C99 and MSC covers a lot of cases, */
41 /* and the fallback default is pretty reliable with modern compilers and machines. */
42 /* Note that MSWin is LLP64 == 32 bit long, while OS X/Linux is LP64 == 64 bit long. */
43 /* so long shouldn't really be used in any code.... */
44 /* (duplicated in icc.h) */
45 
46 /* Use __P64__ as cross platform 64 bit pointer #define */
47 #if defined(__LP64__) || defined(__ILP64__) || defined(__LLP64__) || defined(_WIN64)
48 # define __P64__ 1
49 #endif
50 
51 #ifndef ORD32
52 
53 #if (__STDC_VERSION__ >= 199901L)	/* C99 */		\
54  || defined(_STDINT_H_) || defined(_STDINT_H)		\
55  || defined(_SYS_TYPES_H)
56 
57 #include <stdint.h>
58 
59 #define INR8   int8_t		/* 8 bit signed */
60 #define INR16  int16_t		/* 16 bit signed */
61 #define INR32  int32_t		/* 32 bit signed */
62 #define INR64  int64_t		/* 64 bit signed */
63 #define ORD8   uint8_t		/* 8 bit unsigned */
64 #define ORD16  uint16_t		/* 16 bit unsigned */
65 #define ORD32  uint32_t		/* 32 bit unsigned */
66 #define ORD64  uint64_t		/* 64 bit unsigned */
67 
68 #define PNTR intptr_t
69 
70 #define PF64PREC "ll"		/* printf format precision specifier */
71 #define CF64PREC "LL"		/* Constant precision specifier */
72 
73 #ifndef ATTRIBUTE_NORETURN
74 # ifdef _MSC_VER
75 #  define ATTRIBUTE_NORETURN __declspec(noreturn)
76 # else
77 #  define ATTRIBUTE_NORETURN __attribute__((noreturn))
78 # endif
79 #endif
80 
81 #else  /* !__STDC_VERSION__ */
82 
83 #ifdef _MSC_VER
84 
85 #define INR8   __int8				/* 8 bit signed */
86 #define INR16  __int16				/* 16 bit signed */
87 #define INR32  __int32				/* 32 bit signed */
88 #define INR64  __int64				/* 64 bit signed */
89 #define ORD8   unsigned __int8		/* 8 bit unsigned */
90 #define ORD16  unsigned __int16		/* 16 bit unsigned */
91 #define ORD32  unsigned __int32		/* 32 bit unsigned */
92 #define ORD64  unsigned __int64		/* 64 bit unsigned */
93 
94 #define PNTR UINT_PTR
95 
96 #define PF64PREC "I64"				/* printf format precision specifier */
97 #define CF64PREC "LL"				/* Constant precision specifier */
98 
99 #ifndef ATTRIBUTE_NORETURN
100 #define ATTRIBUTE_NORETURN __attribute__((noreturn))
101 #endif
102 
103 #else  /* !_MSC_VER */
104 
105 /* The following works on a lot of modern systems, including */
106 /* LLP64 and LP64 models, but won't work with ILP64 which needs int32 */
107 
108 #define INR8   signed char		/* 8 bit signed */
109 #define INR16  signed short		/* 16 bit signed */
110 #define INR32  signed int		/* 32 bit signed */
111 #define ORD8   unsigned char	/* 8 bit unsigned */
112 #define ORD16  unsigned short	/* 16 bit unsigned */
113 #define ORD32  unsigned int		/* 32 bit unsigned */
114 
115 #ifdef __GNUC__
116 # ifdef __LP64__	/* long long could be 128 bit ? */
117 #  define INR64  long				/* 64 bit signed */
118 #  define ORD64  unsigned long		/* 64 bit unsigned */
119 #  define PF64PREC "l"			/* printf format precision specifier */
120 #  define CF64PREC "L"			/* Constant precision specifier */
121 # else
122 #  define INR64  long long			/* 64 bit signed */
123 #  define ORD64  unsigned long long	/* 64 bit unsigned */
124 #  define PF64PREC "ll"			/* printf format precision specifier */
125 #  define CF64PREC "LL"			/* Constant precision specifier */
126 # endif /* !__LP64__ */
127 #endif /* __GNUC__ */
128 
129 #define PNTR unsigned long
130 
131 #ifndef ATTRIBUTE_NORETURN
132 # define ATTRIBUTE_NORETURN __attribute__((noreturn))
133 #endif
134 
135 #endif /* !_MSC_VER */
136 #endif /* !__STDC_VERSION__ */
137 #endif /* !ORD32 */
138 
139 /* =========================================================== */
140 /* System compatibility #defines */
141 #if defined (NT)
142 
143 #ifndef sys_stat
144 # define sys_stat _stat
145 #endif
146 #ifndef sys_mkdir
147 # define sys_mkdir _mkdir
148 #endif
149 #ifndef sys_read
150 # define sys_read _read
151 #endif
152 #ifndef sys_utime
153 # define sys_utime _utime
154 # define sys_utimbuf _utimbuf
155 #endif
156 #ifndef sys_access
157 # define sys_access _access
158 #endif
159 
160 #ifndef snprintf
161 # define snprintf _snprintf
162 # define vsnprintf _vsnprintf
163 #endif
164 #ifndef stricmp
165 # define stricmp _stricmp
166 #endif
167 
168 #endif	/* NT */
169 
170 #if defined (UNIX)
171 
172 #ifndef sys_stat
173 # define sys_stat stat
174 #endif
175 #ifndef sys_mkdir
176 # define sys_mkdir mkdir
177 #endif
178 #ifndef sys_read
179 # define sys_read read
180 #endif
181 #ifndef sys_utime
182 # define sys_utime utime
183 # define sys_utimbuf utimbuf
184 #endif
185 #ifndef sys_access
186 # define sys_access access
187 #endif
188 
189 #ifndef stricmp
190 # define stricmp strcasecmp
191 #endif
192 
193 #endif	/* UNIX */
194 
195 /* =========================================================== */
196 /* Some default math limits and constants */
197 #ifndef DBL_EPSILON
198 #define DBL_EPSILON     2.2204460492503131e-016		/* 1.0+DBL_EPSILON != 1.0 */
199 #endif
200 #ifndef DBL_MIN
201 #define DBL_MIN         2.2250738585072014e-308		/* IEEE 64 bit min value */
202 #endif
203 #ifndef DBL_MAX
204 #define DBL_MAX         1.7976931348623158e+308		/* IEEE 64 bit max value */
205 #endif
206 #ifndef DBL_PI
207 #define DBL_PI         3.1415926535897932384626433832795
208 #endif
209 
210 /* =========================================================== */
211 /* General verbose, debug, warning and error logging object and functions */
212 
213 /*
214 
215 	Verbose is simply informational for the end user during normal
216 	operation, typically on stdout or in a text information log.
217 	This will be logged to the verbose log stream.  Level:
218 
219 	 0   = Off
220 	 1   = Brief progress and informational messages
221 	 2   = Much more verbose information
222 	 3+  = Very verbose detail
223 
224 	Debug output is for the developer, to trace what has happened so
225 	as to help diagnose a fault. This will be logged to the debug log stream.
226 	Level:
227 
228 	 0   = Off
229 	 1   = Brief debug info and any internal errors that may be significant
230 	       if something subsequently fails.
231 	 1-5 = Application internals at increasing level of detail
232 	 2-6 = Driver level.(overlaps app & coms)
233 	 6-7 = high level communications
234 	 8-9 = low level communications.
235 
236 	Warning is a serious internal fault that is going to be ignored at the
237 	point it is noticed, but may explain any unexpected behaviour.
238 	It will be reported to the vebose, debug and error log streams.
239 
240 	An error is something that is regarded as fatal at the point it
241 	is noticed. It will be reported to the vebose, debug and error log
242 	streams. The error code and error message will be latched within the
243 	log object if it is the first error logged. It can (theoretically) be
244 	treated as a warning at a higher calling level by calling
245 	a1logue (unlatch error) to reset the error code and message.
246 
247  */
248 
249 #define A1_LOG_BUFSIZE 500
250 
251 
252 struct _a1log {
253 	int refc;					/* Reference count */
254 	char *tag;					/* Optional tag name */
255 	int verb;					/* Current verbosity level (public) */
256 	int debug;					/* Current debug level (public) */
257 
258 	void *cntx;					/* Context to provide to log functions */
259 	void (*logv)(void *cntx, struct _a1log *p, char *fmt, va_list args);
260 								/* Implementation of log verbose */
261 	void (*logd)(void *cntx, struct _a1log *p, char *fmt, va_list args);
262 								/* Implementation of log debug */
263 	void (*loge)(void *cntx, struct _a1log *p, char *fmt, va_list args);
264 								/* Implementation of log warning/error */
265 
266 	int errc; 					/* error code */
267 	char errm[A1_LOG_BUFSIZE];	/* error message (public) */
268 
269 #ifdef NT
270 	CRITICAL_SECTION lock;
271 #endif
272 #ifdef UNIX
273 	pthread_mutex_t lock;
274 #endif
275 }; typedef struct _a1log a1log;
276 
277 
278 
279 /* If log NULL, allocate a new log and return it, */
280 /* otherwise increment reference count and return existing log, */
281 /* exit() if malloc fails. */
282 a1log *new_a1log(
283 	a1log *log,						/* Existing log to reference, NULL if none */
284 	int verb,						/* Verbose level to set */
285 	int debug,						/* Debug level to set */
286 	void *cntx,						/* Function context value */
287 	/* Debug log function to call - stdout if NULL */
288 	void (*logv)(void *cntx, a1log *p, char *fmt, va_list args),
289 	/* Debug log function to call - stderr if NULL */
290 	void (*logd)(void *cntx, a1log *p, char *fmt, va_list args),
291 	/* Warning/error Log function to call - stderr if NULL */
292 	void (*loge)(void *cntx, a1log *p, char *fmt, va_list args)
293 );
294 
295 /* Same as above but set default functions */
296 a1log *new_a1log_d(a1log *log);
297 
298 /* Decrement reference count and free log. */
299 /* Returns NULL */
300 a1log *del_a1log(a1log *log);
301 
302 /* Set the tag. Note that the tag string is NOT copied, just referenced */
303 void a1log_tag(a1log *log, char *tag);
304 
305 /* Log a verbose message if level >= verb */
306 void a1logv(a1log *log, int level, char *fmt, ...);
307 
308 /* Log a debug message if level >= debug */
309 void a1logd(a1log *log, int level, char *fmt, ...);
310 
311 /* log a warning message to the error output  */
312 void a1logw(a1log *log, char *fmt, ...);
313 
314 /* log an error message,  */
315 /* ecode = system, icoms or instrument error */
316 void a1loge(a1log *log, int ecode, char *fmt, ...);
317 
318 /* Unlatch an error message. */
319 /* This resets errc and errm */
320 void a1logue(a1log *log);
321 
322 /* Print bytes as hex to debug log */
323 /* base is the base of the displayed offset */
324 void adump_bytes(a1log *log, char *pfx, unsigned char *buf, int base, int len);
325 
326 /* =========================================================== */
327 /* Globals used to hold certain information */
328 extern char *exe_path;		/* Path leading to executable, not including exe name itself. */
329 							/* Always uses '/' separator */
330 							/* Malloce'd - won't be freed on exit (intended leak) */
331 
332 extern a1log *g_log;		/* Default log */
333 
334 /* These legacy functions that now call through to the default log */
335 #define error_program g_log->tag
336 extern void set_exe_path(char *arg0);
337 
338 extern void ATTRIBUTE_NORETURN error(char *fmt, ...);
339 extern void warning(char *fmt, ...);
340 extern void verbose(int level, char *fmt, ...);
341 
342 extern int ret_null_on_malloc_fail;
343 
344 extern void check_if_not_interactive();
345 extern int not_interactive;
346 extern char cr_char;
347 
348 /* =========================================================== */
349 
350 
351 /* =========================================================== */
352 
353 /* reallocate and clear new allocation */
354 void *recalloc(		/* Return new address */
355 void *ptr,			/* Current address */
356 size_t cnum,		/* Current number and unit size */
357 size_t csize,
358 size_t nnum,		/* New number and unit size */
359 size_t nsize
360 );
361 
362 /* =========================================================== */
363 
364 #if defined(__APPLE__)
365 
366 /* Tell App Nap that this is user initiated */
367 void osx_userinitiated_start();
368 
369 /* Done with user initiated */
370 void osx_userinitiated_end();
371 
372 /* Tell App Nap that this is latency critical */
373 void osx_latencycritical_start();
374 
375 /* Done with latency critical */
376 void osx_latencycritical_end();
377 
378 #endif	/* __APPLE__ */
379 
380 /* =========================================================== */
381 
382 /* Numerical recipes vector/matrix support functions */
383 /* Note that the index arguments are the inclusive low and high values */
384 
385 /* Double */
386 double *dvector(int nl,int nh);
387 double *dvectorz(int nl,int nh);
388 void free_dvector(double *v,int nl,int nh);
389 
390 double **dmatrix(int nrl, int nrh, int ncl, int nch);
391 double **dmatrixz(int nrl, int nrh, int ncl, int nch);
392 void free_dmatrix(double **m, int nrl, int nrh, int ncl, int nch);
393 
394 double **dhmatrix(int nrl, int nrh, int ncl, int nch);
395 double **dhmatrixz(int nrl, int nrh, int ncl, int nch);
396 void free_dhmatrix(double **m, int nrl, int nrh, int ncl, int nch);
397 
398 void copy_dmatrix(double **dst, double **src, int nrl, int nrh, int ncl, int nch);
399 void copy_dmatrix_to3x3(double dst[3][3], double **src, int nrl, int nrh, int ncl, int nch);
400 
401 double **convert_dmatrix(double *a,int nrl,int nrh,int ncl,int nch);
402 void free_convert_dmatrix(double **m,int nrl,int nrh,int ncl,int nch);
403 
404 
405 /* Float */
406 float *fvector(int nl,int nh);
407 float *fvectorz(int nl,int nh);
408 void free_fvector(float *v,int nl,int nh);
409 
410 float **fmatrix(int nrl, int nrh, int ncl, int nch);
411 float **fmatrixz(int nrl, int nrh, int ncl, int nch);
412 void free_fmatrix(float **m, int nrl, int nrh, int ncl, int nch);
413 
414 
415 /* Int */
416 int *ivector(int nl,int nh);
417 int *ivectorz(int nl,int nh);
418 void free_ivector(int *v,int nl,int nh);
419 
420 int **imatrix(int nrl,int nrh,int ncl,int nch);
421 int **imatrixz(int nrl,int nrh,int ncl,int nch);
422 void free_imatrix(int **m,int nrl,int nrh,int ncl,int nch);
423 
424 
425 /* Short */
426 short *svector(int nl, int nh);
427 short *svectorz(int nl, int nh);
428 void free_svector(short *v, int nl, int nh);
429 
430 short **smatrix(int nrl,int nrh,int ncl,int nch);
431 short **smatrixz(int nrl,int nrh,int ncl,int nch);
432 void free_smatrix(short **m,int nrl,int nrh,int ncl,int nch);
433 
434 /* ----------------------------------------------------------- */
435 /* Basic matrix operations */
436 
437 /* Transpose a 0 base matrix */
438 void matrix_trans(double **d, double **s, int nr,  int nc);
439 
440 /* Matrix multiply 0 based matricies */
441 int matrix_mult(
442 	double **d,  int nr,  int nc,
443 	double **s1, int nr1, int nc1,
444 	double **s2, int nr2, int nc2
445 );
446 
447 /* Diagnostic */
448 void matrix_print(char *c, double **a, int nr,  int nc);
449 
450 /* =========================================================== */
451 
452 /* Cast a native double to an IEEE754 encoded single precision value, */
453 /* in a platform independent fashion. */
454 ORD32 doubletoIEEE754(double d);
455 
456 /* Cast an IEEE754 encoded single precision value to a native double, */
457 /* in a platform independent fashion. */
458 double IEEE754todouble(ORD32 ip);
459 
460 
461 /* Cast a native double to an IEEE754 encoded double precision value, */
462 /* in a platform independent fashion. */
463 ORD64 doubletoIEEE754_64(double d);
464 
465 /* Cast an IEEE754 encoded double precision value to a native double, */
466 /* in a platform independent fashion. */
467 double IEEE754_64todouble(ORD64 ip);
468 
469 /* Return a string representation of a 32 bit ctime. */
470 /* A static buffer is used. There is no \n at the end */
471 char *ctime_32(const INR32 *timer);
472 
473 /* Return a string representation of a 64 bit ctime. */
474 /* A static buffer is used. There is no \n at the end */
475 char *ctime_64(const INR64 *timer);
476 
477 /*******************************************/
478 /* Native to/from byte buffer functions    */
479 /*******************************************/
480 
481 /* No overflow detection is done - numbers are clipped */
482 
483 /* be = Big Endian */
484 /* le = Little Endian */
485 
486 /* Unsigned 8 bit */
487 unsigned int read_ORD8(ORD8 *p);
488 void write_ORD8(ORD8 *p, unsigned int d);
489 
490 /* Signed 8 bit */
491 int read_INR8(ORD8 *p);
492 void write_INR8(ORD8 *p, int d);
493 
494 /* Unsigned 16 bit */
495 unsigned int read_ORD16_be(ORD8 *p);
496 unsigned int read_ORD16_le(ORD8 *p);
497 void write_ORD16_be(ORD8 *p, unsigned int d);
498 void write_ORD16_le(ORD8 *p, unsigned int d);
499 
500 /* Signed 16 bit */
501 int read_INR16_be(ORD8 *p);
502 int read_INR16_le(ORD8 *p);
503 void write_INR16_be(ORD8 *p, int d);
504 void write_INR16_le(ORD8 *p, int d);
505 
506 /* Unsigned 32 bit */
507 unsigned int read_ORD32_be(ORD8 *p);
508 unsigned int read_ORD32_le(ORD8 *p);
509 void write_ORD32_be(ORD8 *p, unsigned int d);
510 void write_ORD32_le(ORD8 *p, unsigned int d);
511 
512 /* Signed 32 bit */
513 int read_INR32_be(ORD8 *p);
514 int read_INR32_le(ORD8 *p);
515 void write_INR32_be(ORD8 *p, int d);
516 void write_INR32_le(ORD8 *p, int d);
517 
518 /* Unsigned 64 bit */
519 ORD64 read_ORD64_be(ORD8 *p);
520 ORD64 read_ORD64_le(ORD8 *p);
521 void write_ORD64_be(ORD8 *p, ORD64 d);
522 void write_ORD64_le(ORD8 *p, ORD64 d);
523 
524 /* Signed 64 bit */
525 INR64 read_INR64_be(ORD8 *p);
526 INR64 read_INR64_le(ORD8 *p);
527 void write_INR64_be(ORD8 *p, INR64 d);
528 void write_INR64_le(ORD8 *p, INR64 d);
529 
530 /*******************************************/
531 
532 /* Sleep for the given number of msec */
533 void msec_sleep(unsigned int msec);
534 
535 /* Return the current time in msec since */
536 /* the first invokation of msec_time() */
537 unsigned int msec_time();
538 
539 /* Return the current time in usec */
540 /* (The first invokation of usec_time() returns zero) */
541 double usec_time();
542 
543 /*******************************************/
544 /* Debug convenience functions (duplicated in icc) */
545 
546 /* Print an int vector to a string. */
547 /* Returned static buffer is re-used every 5 calls. */
548 char *debPiv(int di, int *p);
549 
550 /* Print a double color vector to a string. */
551 /* Returned static buffer is re-used every 5 calls. */
552 char *debPdv(int di, double *p);
553 
554 /* Print a float color vector to a string. */
555 /* Returned static buffer is re-used every 5 calls. */
556 char *debPfv(int di, float *p);
557 
558 /*******************************************/
559 /* Numerical diagnostics */
560 
561 #ifndef isNan
562 #define isNan(x) ((x) != (x))
563 #define isFinite(x) ((x) == 0.0 || (x) * 1.0000001 != (x))
564 #define isNFinite(x) ((x) != 0.0 && (x) * 1.0000001 == (x))
565 #endif
566 
567 
568 #ifdef __cplusplus
569 	}
570 #endif
571 
572 #define NUMSUP_H
573 #endif /* NUMSUP_H */
574