1 /*<html><pre> -<a href="index.htm#TOC"
2 >-------------------------------</a><a name="TOP">-</a>
3
4 rboxlib.c
5 Generate input points
6
7 notes:
8 For documentation, see prompt[] of rbox.c
9 50 points generated for 'rbox D4'
10
11 WARNING:
12 incorrect range if qh_RANDOMmax is defined wrong (user.h)
13 */
14
15 #include "random.h"
16 #include "libqhull.h"
17
18 #include <ctype.h>
19 #include <limits.h>
20 #include <math.h>
21 #include <setjmp.h>
22 #include <string.h>
23 #include <time.h>
24 #include <stdio.h>
25 #include <stdlib.h>
26
27 #ifdef _MSC_VER /* Microsoft Visual C++ */
28 #pragma warning( disable : 4706) /* assignment within conditional expression. */
29 #pragma warning( disable : 4996) /* this function (strncat,sprintf,strcpy) or variable may be unsafe. */
30 #endif
31
32 #define MAXdim 200
33
34 /* ------------------------------ prototypes ----------------*/
35 int roundi( double a);
36 void out1( double a);
37 void out2n( double a, double b);
38 void out3n( double a, double b, double c);
39
40 void qh_fprintf_rbox(FILE *fp, int msgcode, const char *fmt, ... );
41 void qh_free(void *mem);
42 void *qh_malloc(size_t size);
43 int qh_rand( void);
44 void qh_srand( int seed);
45
46
47 /* ------------------------------ globals -------------------*/
48
49 /* No state is carried between rbox requests */
50 typedef struct rboxT rboxT;
51 struct rboxT {
52 FILE *fout;
53 FILE *ferr;
54 int isinteger;
55 double out_offset;
56 jmp_buf errexit; /* exit label for rboxpoints, defined by setjmp(), called by qh_errexit_rbox() */
57 };
58
59
60 int rbox_inuse= 0;
61 rboxT rbox;
62
63 /*-<a href="qh-qhull.htm#TOC"
64 >-------------------------------</a><a name="rboxpoints">-</a>
65
66 qh_rboxpoints( fout, ferr, rbox_command )
67 Generate points to fout according to rbox options
68 Report errors on ferr
69
70 returns:
71 0 (qh_ERRnone) on success
72 1 (qh_ERRinput) on input error
73 4 (qh_ERRmem) on memory error
74 5 (qh_ERRqhull) on internal error
75
76 notes:
77 To avoid stdio, redefine qh_malloc, qh_free, and qh_fprintf_rbox (user.c)
78 Rbox is not multithreaded.
79
80 design:
81 Straight line code (consider defining a struct and functions):
82
83 Parse arguments into variables
84 Determine the number of points
85 Generate the points
86 */
qh_rboxpoints(FILE * fout,FILE * ferr,char * rbox_command)87 int qh_rboxpoints(FILE* fout, FILE* ferr, char* rbox_command) {
88 int i,j,k;
89 int gendim;
90 int cubesize, diamondsize, seed=0, count, apex;
91 int dim=3 , numpoints= 0, totpoints, addpoints=0;
92 int issphere=0, isaxis=0, iscdd= 0, islens= 0, isregular=0, iswidth=0, addcube=0;
93 int isgap=0, isspiral=0, NOcommand= 0, adddiamond=0;
94 int israndom=0, istime=0;
95 int isbox=0, issimplex=0, issimplex2=0, ismesh=0;
96 double width=0.0, gap=0.0, radius= 0.0;
97 double coord[MAXdim], offset, meshm=3.0, meshn=4.0, meshr=5.0;
98 double *simplex= NULL, *simplexp;
99 int nthroot, mult[MAXdim];
100 double norm, factor, randr, rangap, lensangle= 0, lensbase= 1;
101 double anglediff, angle, x, y, cube= 0.0, diamond= 0.0;
102 double box= qh_DEFAULTbox; /* scale all numbers before output */
103 double randmax= qh_RANDOMmax;
104 char command[200], seedbuf[200];
105 char *s= command, *t, *first_point= NULL;
106 time_t timedata = 0;
107 int exitcode;
108
109 if (rbox_inuse) {
110 qh_fprintf_rbox(rbox.ferr, 6188, "rbox error: rbox in use by another process. Please lock calls to rbox.\n");
111 return qh_ERRqhull;
112 }
113 rbox_inuse= True;
114 rbox.ferr= ferr;
115 rbox.fout= fout;
116
117 exitcode= setjmp(rbox.errexit);
118 if (exitcode) {
119 /* same code for error exit and normal return */
120 if (simplex)
121 qh_free(simplex);
122 rbox_inuse= False;
123 return exitcode;
124 }
125
126 *command= '\0';
127 strncat(command, rbox_command, sizeof(command)-strlen(command)-1);
128
129 while (*s && !isspace(*s)) /* skip program name */
130 s++;
131 while (*s) {
132 while (*s && isspace(*s))
133 s++;
134 if (*s == '-')
135 s++;
136 if (!*s)
137 break;
138 if (isdigit(*s)) {
139 numpoints= qh_strtol(s, &s);
140 continue;
141 }
142 /* ============= read flags =============== */
143 switch (*s++) {
144 case 'c':
145 addcube= 1;
146 t= s;
147 while (isspace(*t))
148 t++;
149 if (*t == 'G')
150 cube= qh_strtod(++t, &s);
151 break;
152 case 'd':
153 adddiamond= 1;
154 t= s;
155 while (isspace(*t))
156 t++;
157 if (*t == 'G')
158 diamond= qh_strtod(++t, &s);
159 break;
160 case 'h':
161 iscdd= 1;
162 break;
163 case 'l':
164 isspiral= 1;
165 break;
166 case 'n':
167 NOcommand= 1;
168 break;
169 case 'r':
170 isregular= 1;
171 break;
172 case 's':
173 issphere= 1;
174 break;
175 case 't':
176 istime= 1;
177 if (isdigit(*s)) {
178 seed= qh_strtol(s, &s);
179 israndom= 0;
180 }else
181 israndom= 1;
182 break;
183 case 'x':
184 issimplex= 1;
185 break;
186 case 'y':
187 issimplex2= 1;
188 break;
189 case 'z':
190 rbox.isinteger= 1;
191 break;
192 case 'B':
193 box= qh_strtod(s, &s);
194 isbox= 1;
195 break;
196 case 'D':
197 dim= qh_strtol(s, &s);
198 if (dim < 1
199 || dim > MAXdim) {
200 qh_fprintf_rbox(rbox.ferr, 6189, "rbox error: dimension, D%d, out of bounds (>=%d or <=0)", dim, MAXdim);
201 qh_errexit_rbox(qh_ERRinput);
202 }
203 break;
204 case 'G':
205 if (isdigit(*s))
206 gap= qh_strtod(s, &s);
207 else
208 gap= 0.5;
209 isgap= 1;
210 break;
211 case 'L':
212 if (isdigit(*s))
213 radius= qh_strtod(s, &s);
214 else
215 radius= 10;
216 islens= 1;
217 break;
218 case 'M':
219 ismesh= 1;
220 if (*s)
221 meshn= qh_strtod(s, &s);
222 if (*s == ',') {
223 ++s;
224 meshm= qh_strtod(s, &s);
225 }else
226 meshm= 0.0;
227 if (*s == ',') {
228 ++s;
229 meshr= qh_strtod(s, &s);
230 }else
231 meshr= sqrt(meshn*meshn + meshm*meshm);
232 if (*s && !isspace(*s)) {
233 qh_fprintf_rbox(rbox.ferr, 7069, "rbox warning: assuming 'M3,4,5' since mesh args are not integers or reals\n");
234 meshn= 3.0, meshm=4.0, meshr=5.0;
235 }
236 break;
237 case 'O':
238 rbox.out_offset= qh_strtod(s, &s);
239 break;
240 case 'P':
241 if (!first_point)
242 first_point= s-1;
243 addpoints++;
244 while (*s && !isspace(*s)) /* read points later */
245 s++;
246 break;
247 case 'W':
248 width= qh_strtod(s, &s);
249 iswidth= 1;
250 break;
251 case 'Z':
252 if (isdigit(*s))
253 radius= qh_strtod(s, &s);
254 else
255 radius= 1.0;
256 isaxis= 1;
257 break;
258 default:
259 qh_fprintf_rbox(rbox.ferr, 7070, "rbox error: unknown flag at %s.\nExecute 'rbox' without arguments for documentation.\n", s);
260 qh_errexit_rbox(qh_ERRinput);
261 }
262 if (*s && !isspace(*s)) {
263 qh_fprintf_rbox(rbox.ferr, 7071, "rbox error: missing space between flags at %s.\n", s);
264 qh_errexit_rbox(qh_ERRinput);
265 }
266 }
267
268 /* ============= defaults, constants, and sizes =============== */
269 if (rbox.isinteger && !isbox)
270 box= qh_DEFAULTzbox;
271 if (addcube) {
272 cubesize= (int)floor(ldexp(1.0,dim)+0.5);
273 if (cube == 0.0)
274 cube= box;
275 }else
276 cubesize= 0;
277 if (adddiamond) {
278 diamondsize= 2*dim;
279 if (diamond == 0.0)
280 diamond= box;
281 }else
282 diamondsize= 0;
283 if (islens) {
284 if (isaxis) {
285 qh_fprintf_rbox(rbox.ferr, 6190, "rbox error: can not combine 'Ln' with 'Zn'\n");
286 qh_errexit_rbox(qh_ERRinput);
287 }
288 if (radius <= 1.0) {
289 qh_fprintf_rbox(rbox.ferr, 6191, "rbox error: lens radius %.2g should be greater than 1.0\n",
290 radius);
291 qh_errexit_rbox(qh_ERRinput);
292 }
293 lensangle= asin(1.0/radius);
294 lensbase= radius * cos(lensangle);
295 }
296
297 if (!numpoints) {
298 if (issimplex2)
299 ; /* ok */
300 else if (isregular + issimplex + islens + issphere + isaxis + isspiral + iswidth + ismesh) {
301 qh_fprintf_rbox(rbox.ferr, 6192, "rbox error: missing count\n");
302 qh_errexit_rbox(qh_ERRinput);
303 }else if (adddiamond + addcube + addpoints)
304 ; /* ok */
305 else {
306 numpoints= 50; /* ./rbox D4 is the test case */
307 issphere= 1;
308 }
309 }
310 if ((issimplex + islens + isspiral + ismesh > 1)
311 || (issimplex + issphere + isspiral + ismesh > 1)) {
312 qh_fprintf_rbox(rbox.ferr, 6193, "rbox error: can only specify one of 'l', 's', 'x', 'Ln', or 'Mn,m,r' ('Ln s' is ok).\n");
313 qh_errexit_rbox(qh_ERRinput);
314 }
315
316 /* ============= print header with total points =============== */
317 if (issimplex || ismesh)
318 totpoints= numpoints;
319 else if (issimplex2)
320 totpoints= numpoints+dim+1;
321 else if (isregular) {
322 totpoints= numpoints;
323 if (dim == 2) {
324 if (islens)
325 totpoints += numpoints - 2;
326 }else if (dim == 3) {
327 if (islens)
328 totpoints += 2 * numpoints;
329 else if (isgap)
330 totpoints += 1 + numpoints;
331 else
332 totpoints += 2;
333 }
334 }else
335 totpoints= numpoints + isaxis;
336 totpoints += cubesize + diamondsize + addpoints;
337
338 /* ============= seed randoms =============== */
339 if (istime == 0) {
340 for (s=command; *s; s++) {
341 if (issimplex2 && *s == 'y') /* make 'y' same seed as 'x' */
342 i= 'x';
343 else
344 i= *s;
345 seed= 11*seed + i;
346 }
347 }else if (israndom) {
348 seed= (int)time(&timedata);
349 sprintf(seedbuf, " t%d", seed); /* appends an extra t, not worth removing */
350 strncat(command, seedbuf, sizeof(command)-strlen(command)-1);
351 t= strstr(command, " t ");
352 if (t)
353 strcpy(t+1, t+3); /* remove " t " */
354 } /* else, seed explicitly set to n */
355 qh_RANDOMseed_(seed);
356
357 /* ============= print header =============== */
358
359 if (iscdd)
360 qh_fprintf_rbox(rbox.fout, 9391, "%s\nbegin\n %d %d %s\n",
361 NOcommand ? "" : command,
362 totpoints, dim+1,
363 rbox.isinteger ? "integer" : "real");
364 else if (NOcommand)
365 qh_fprintf_rbox(rbox.fout, 9392, "%d\n%d\n", dim, totpoints);
366 else
367 qh_fprintf_rbox(rbox.fout, 9393, "%d %s\n%d\n", dim, command, totpoints);
368
369 /* ============= explicit points =============== */
370 if ((s= first_point)) {
371 while (s && *s) { /* 'P' */
372 count= 0;
373 if (iscdd)
374 out1( 1.0);
375 while (*++s) {
376 out1( qh_strtod(s, &s));
377 count++;
378 if (isspace(*s) || !*s)
379 break;
380 if (*s != ',') {
381 qh_fprintf_rbox(rbox.ferr, 6194, "rbox error: missing comma after coordinate in %s\n\n", s);
382 qh_errexit_rbox(qh_ERRinput);
383 }
384 }
385 if (count < dim) {
386 for (k=dim-count; k--; )
387 out1( 0.0);
388 }else if (count > dim) {
389 qh_fprintf_rbox(rbox.ferr, 6195, "rbox error: %d coordinates instead of %d coordinates in %s\n\n",
390 count, dim, s);
391 qh_errexit_rbox(qh_ERRinput);
392 }
393 qh_fprintf_rbox(rbox.fout, 9394, "\n");
394 while ((s= strchr(s, 'P'))) {
395 if (isspace(s[-1]))
396 break;
397 }
398 }
399 }
400
401 /* ============= simplex distribution =============== */
402 if (issimplex+issimplex2) {
403 if (!(simplex= (double*)qh_malloc( dim * (dim+1) * sizeof(double)))) {
404 qh_fprintf_rbox(rbox.ferr, 6196, "rbox error: insufficient memory for simplex\n");
405 qh_errexit_rbox(qh_ERRmem); /* qh_ERRmem */
406 }
407 simplexp= simplex;
408 if (isregular) {
409 for (i=0; i<dim; i++) {
410 for (k=0; k<dim; k++)
411 *(simplexp++)= i==k ? 1.0 : 0.0;
412 }
413 for (k=0; k<dim; k++)
414 *(simplexp++)= -1.0;
415 }else {
416 for (i=0; i<dim+1; i++) {
417 for (k=0; k<dim; k++) {
418 randr= qh_RANDOMint;
419 *(simplexp++)= 2.0 * randr/randmax - 1.0;
420 }
421 }
422 }
423 if (issimplex2) {
424 simplexp= simplex;
425 for (i=0; i<dim+1; i++) {
426 if (iscdd)
427 out1( 1.0);
428 for (k=0; k<dim; k++)
429 out1( *(simplexp++) * box);
430 qh_fprintf_rbox(rbox.fout, 9395, "\n");
431 }
432 }
433 for (j=0; j<numpoints; j++) {
434 if (iswidth)
435 apex= qh_RANDOMint % (dim+1);
436 else
437 apex= -1;
438 for (k=0; k<dim; k++)
439 coord[k]= 0.0;
440 norm= 0.0;
441 for (i=0; i<dim+1; i++) {
442 randr= qh_RANDOMint;
443 factor= randr/randmax;
444 if (i == apex)
445 factor *= width;
446 norm += factor;
447 for (k=0; k<dim; k++) {
448 simplexp= simplex + i*dim + k;
449 coord[k] += factor * (*simplexp);
450 }
451 }
452 for (k=0; k<dim; k++)
453 coord[k] /= norm;
454 if (iscdd)
455 out1( 1.0);
456 for (k=0; k < dim; k++)
457 out1( coord[k] * box);
458 qh_fprintf_rbox(rbox.fout, 9396, "\n");
459 }
460 isregular= 0; /* continue with isbox */
461 numpoints= 0;
462 }
463
464 /* ============= mesh distribution =============== */
465 if (ismesh) {
466 nthroot= (int)(pow((double)numpoints, 1.0/dim) + 0.99999);
467 for (k=dim; k--; )
468 mult[k]= 0;
469 for (i=0; i < numpoints; i++) {
470 for (k=0; k < dim; k++) {
471 if (k == 0)
472 out1( mult[0] * meshn + mult[1] * (-meshm));
473 else if (k == 1)
474 out1( mult[0] * meshm + mult[1] * meshn);
475 else
476 out1( mult[k] * meshr );
477 }
478 qh_fprintf_rbox(rbox.fout, 9397, "\n");
479 for (k=0; k < dim; k++) {
480 if (++mult[k] < nthroot)
481 break;
482 mult[k]= 0;
483 }
484 }
485 }
486 /* ============= regular points for 's' =============== */
487 else if (isregular && !islens) {
488 if (dim != 2 && dim != 3) {
489 qh_fprintf_rbox(rbox.ferr, 6197, "rbox error: regular points can be used only in 2-d and 3-d\n\n");
490 qh_errexit_rbox(qh_ERRinput);
491 }
492 if (!isaxis || radius == 0.0) {
493 isaxis= 1;
494 radius= 1.0;
495 }
496 if (dim == 3) {
497 if (iscdd)
498 out1( 1.0);
499 out3n( 0.0, 0.0, -box);
500 if (!isgap) {
501 if (iscdd)
502 out1( 1.0);
503 out3n( 0.0, 0.0, box);
504 }
505 }
506 angle= 0.0;
507 anglediff= 2.0 * M_PI/numpoints;
508 for (i=0; i < numpoints; i++) {
509 angle += anglediff;
510 x= radius * cos(angle);
511 y= radius * sin(angle);
512 if (dim == 2) {
513 if (iscdd)
514 out1( 1.0);
515 out2n( x*box, y*box);
516 }else {
517 norm= sqrt(1.0 + x*x + y*y);
518 if (iscdd)
519 out1( 1.0);
520 out3n( box*x/norm, box*y/norm, box/norm);
521 if (isgap) {
522 x *= 1-gap;
523 y *= 1-gap;
524 norm= sqrt(1.0 + x*x + y*y);
525 if (iscdd)
526 out1( 1.0);
527 out3n( box*x/norm, box*y/norm, box/norm);
528 }
529 }
530 }
531 }
532 /* ============= regular points for 'r Ln D2' =============== */
533 else if (isregular && islens && dim == 2) {
534 double cos_0;
535
536 angle= lensangle;
537 anglediff= 2 * lensangle/(numpoints - 1);
538 cos_0= cos(lensangle);
539 for (i=0; i < numpoints; i++, angle -= anglediff) {
540 x= radius * sin(angle);
541 y= radius * (cos(angle) - cos_0);
542 if (iscdd)
543 out1( 1.0);
544 out2n( x*box, y*box);
545 if (i != 0 && i != numpoints - 1) {
546 if (iscdd)
547 out1( 1.0);
548 out2n( x*box, -y*box);
549 }
550 }
551 }
552 /* ============= regular points for 'r Ln D3' =============== */
553 else if (isregular && islens && dim != 2) {
554 if (dim != 3) {
555 qh_fprintf_rbox(rbox.ferr, 6198, "rbox error: regular points can be used only in 2-d and 3-d\n\n");
556 qh_errexit_rbox(qh_ERRinput);
557 }
558 angle= 0.0;
559 anglediff= 2* M_PI/numpoints;
560 if (!isgap) {
561 isgap= 1;
562 gap= 0.5;
563 }
564 offset= sqrt(radius * radius - (1-gap)*(1-gap)) - lensbase;
565 for (i=0; i < numpoints; i++, angle += anglediff) {
566 x= cos(angle);
567 y= sin(angle);
568 if (iscdd)
569 out1( 1.0);
570 out3n( box*x, box*y, 0.0);
571 x *= 1-gap;
572 y *= 1-gap;
573 if (iscdd)
574 out1( 1.0);
575 out3n( box*x, box*y, box * offset);
576 if (iscdd)
577 out1( 1.0);
578 out3n( box*x, box*y, -box * offset);
579 }
580 }
581 /* ============= apex of 'Zn' distribution + gendim =============== */
582 else {
583 if (isaxis) {
584 gendim= dim-1;
585 if (iscdd)
586 out1( 1.0);
587 for (j=0; j < gendim; j++)
588 out1( 0.0);
589 out1( -box);
590 qh_fprintf_rbox(rbox.fout, 9398, "\n");
591 }else if (islens)
592 gendim= dim-1;
593 else
594 gendim= dim;
595 /* ============= generate random point in unit cube =============== */
596 for (i=0; i < numpoints; i++) {
597 norm= 0.0;
598 for (j=0; j < gendim; j++) {
599 randr= qh_RANDOMint;
600 coord[j]= 2.0 * randr/randmax - 1.0;
601 norm += coord[j] * coord[j];
602 }
603 norm= sqrt(norm);
604 /* ============= dim-1 point of 'Zn' distribution ========== */
605 if (isaxis) {
606 if (!isgap) {
607 isgap= 1;
608 gap= 1.0;
609 }
610 randr= qh_RANDOMint;
611 rangap= 1.0 - gap * randr/randmax;
612 factor= radius * rangap / norm;
613 for (j=0; j<gendim; j++)
614 coord[j]= factor * coord[j];
615 /* ============= dim-1 point of 'Ln s' distribution =========== */
616 }else if (islens && issphere) {
617 if (!isgap) {
618 isgap= 1;
619 gap= 1.0;
620 }
621 randr= qh_RANDOMint;
622 rangap= 1.0 - gap * randr/randmax;
623 factor= rangap / norm;
624 for (j=0; j<gendim; j++)
625 coord[j]= factor * coord[j];
626 /* ============= dim-1 point of 'Ln' distribution ========== */
627 }else if (islens && !issphere) {
628 if (!isgap) {
629 isgap= 1;
630 gap= 1.0;
631 }
632 j= qh_RANDOMint % gendim;
633 if (coord[j] < 0)
634 coord[j]= -1.0 - coord[j] * gap;
635 else
636 coord[j]= 1.0 - coord[j] * gap;
637 /* ============= point of 'l' distribution =============== */
638 }else if (isspiral) {
639 if (dim != 3) {
640 qh_fprintf_rbox(rbox.ferr, 6199, "rbox error: spiral distribution is available only in 3d\n\n");
641 longjmp(rbox.errexit,qh_ERRinput);
642 }
643 coord[0]= cos(2*M_PI*i/(numpoints - 1));
644 coord[1]= sin(2*M_PI*i/(numpoints - 1));
645 coord[2]= 2.0*(double)i/(double)(numpoints-1) - 1.0;
646 /* ============= point of 's' distribution =============== */
647 }else if (issphere) {
648 factor= 1.0/norm;
649 if (iswidth) {
650 randr= qh_RANDOMint;
651 factor *= 1.0 - width * randr/randmax;
652 }
653 for (j=0; j<dim; j++)
654 coord[j]= factor * coord[j];
655 }
656 /* ============= project 'Zn s' point in to sphere =============== */
657 if (isaxis && issphere) {
658 coord[dim-1]= 1.0;
659 norm= 1.0;
660 for (j=0; j<gendim; j++)
661 norm += coord[j] * coord[j];
662 norm= sqrt(norm);
663 for (j=0; j<dim; j++)
664 coord[j]= coord[j] / norm;
665 if (iswidth) {
666 randr= qh_RANDOMint;
667 coord[dim-1] *= 1 - width * randr/randmax;
668 }
669 /* ============= project 'Zn' point onto cube =============== */
670 }else if (isaxis && !issphere) { /* not very interesting */
671 randr= qh_RANDOMint;
672 coord[dim-1]= 2.0 * randr/randmax - 1.0;
673 /* ============= project 'Ln' point out to sphere =============== */
674 }else if (islens) {
675 coord[dim-1]= lensbase;
676 for (j=0, norm= 0; j<dim; j++)
677 norm += coord[j] * coord[j];
678 norm= sqrt(norm);
679 for (j=0; j<dim; j++)
680 coord[j]= coord[j] * radius/ norm;
681 coord[dim-1] -= lensbase;
682 if (iswidth) {
683 randr= qh_RANDOMint;
684 coord[dim-1] *= 1 - width * randr/randmax;
685 }
686 if (qh_RANDOMint > randmax/2)
687 coord[dim-1]= -coord[dim-1];
688 /* ============= project 'Wn' point toward boundary =============== */
689 }else if (iswidth && !issphere) {
690 j= qh_RANDOMint % gendim;
691 if (coord[j] < 0)
692 coord[j]= -1.0 - coord[j] * width;
693 else
694 coord[j]= 1.0 - coord[j] * width;
695 }
696 /* ============= write point =============== */
697 if (iscdd)
698 out1( 1.0);
699 for (k=0; k < dim; k++)
700 out1( coord[k] * box);
701 qh_fprintf_rbox(rbox.fout, 9399, "\n");
702 }
703 }
704
705 /* ============= write cube vertices =============== */
706 if (addcube) {
707 for (j=0; j<cubesize; j++) {
708 if (iscdd)
709 out1( 1.0);
710 for (k=dim-1; k>=0; k--) {
711 if (j & ( 1 << k))
712 out1( cube);
713 else
714 out1( -cube);
715 }
716 qh_fprintf_rbox(rbox.fout, 9400, "\n");
717 }
718 }
719
720 /* ============= write diamond vertices =============== */
721 if (adddiamond) {
722 for (j=0; j<diamondsize; j++) {
723 if (iscdd)
724 out1( 1.0);
725 for (k=dim-1; k>=0; k--) {
726 if (j/2 != k)
727 out1( 0.0);
728 else if (j & 0x1)
729 out1( diamond);
730 else
731 out1( -diamond);
732 }
733 qh_fprintf_rbox(rbox.fout, 9401, "\n");
734 }
735 }
736
737 if (iscdd)
738 qh_fprintf_rbox(rbox.fout, 9402, "end\nhull\n");
739
740 /* same code for error exit and normal return */
741 if (simplex)
742 qh_free(simplex);
743 rbox_inuse= False;
744 return qh_ERRnone;
745 } /* rboxpoints */
746
747 /*------------------------------------------------
748 outxxx - output functions
749 */
roundi(double a)750 int roundi( double a) {
751 if (a < 0.0) {
752 if (a - 0.5 < INT_MIN) {
753 qh_fprintf_rbox(rbox.ferr, 6200, "rbox input error: negative coordinate %2.2g is too large. Reduce 'Bn'\n", a);
754 qh_errexit_rbox(qh_ERRinput);
755 }
756 return (int)(a - 0.5);
757 }else {
758 if (a + 0.5 > INT_MAX) {
759 qh_fprintf_rbox(rbox.ferr, 6201, "rbox input error: coordinate %2.2g is too large. Reduce 'Bn'\n", a);
760 qh_errexit_rbox(qh_ERRinput);
761 }
762 return (int)(a + 0.5);
763 }
764 } /* roundi */
765
out1(double a)766 void out1(double a) {
767
768 if (rbox.isinteger)
769 qh_fprintf_rbox(rbox.fout, 9403, "%d ", roundi( a+rbox.out_offset));
770 else
771 qh_fprintf_rbox(rbox.fout, 9404, qh_REAL_1, a+rbox.out_offset);
772 } /* out1 */
773
out2n(double a,double b)774 void out2n( double a, double b) {
775
776 if (rbox.isinteger)
777 qh_fprintf_rbox(rbox.fout, 9405, "%d %d\n", roundi(a+rbox.out_offset), roundi(b+rbox.out_offset));
778 else
779 qh_fprintf_rbox(rbox.fout, 9406, qh_REAL_2n, a+rbox.out_offset, b+rbox.out_offset);
780 } /* out2n */
781
out3n(double a,double b,double c)782 void out3n( double a, double b, double c) {
783
784 if (rbox.isinteger)
785 qh_fprintf_rbox(rbox.fout, 9407, "%d %d %d\n", roundi(a+rbox.out_offset), roundi(b+rbox.out_offset), roundi(c+rbox.out_offset));
786 else
787 qh_fprintf_rbox(rbox.fout, 9408, qh_REAL_3n, a+rbox.out_offset, b+rbox.out_offset, c+rbox.out_offset);
788 } /* out3n */
789
qh_errexit_rbox(int exitcode)790 void qh_errexit_rbox(int exitcode)
791 {
792 longjmp(rbox.errexit, exitcode);
793 } /* rbox_errexit */
794