1 #ifdef HAVE_CONFIG
2 #include "fsudace_config.h"
3 #endif
4 # ifdef HAVE_STD
5 # include <cstdlib>
6 # include <cmath>
7 # include <cstring>
8 # include <ctime>
9 # else
10 # include <stdlib.h>
11 # include <math.h>
12 # include <string.h>
13 # include <time.h>
14 # endif
15
16 # include <iostream>
17 # include <fstream>
18 # include <iomanip>
19
20 using namespace std;
21
22 # include "fsu.H"
23
24 //*******************************************************************************
25
ch_cap(char c)26 char ch_cap ( char c )
27
28 //*******************************************************************************
29 //
30 // Purpose:
31 //
32 // CH_CAP capitalizes a single character.
33 //
34 // License:
35 //
36 // Copyright (C) 2004 John Burkardt and Max Gunzburger
37 //
38 // This library is free software; you can redistribute it and/or
39 // modify it under the terms of the GNU Lesser General Public
40 // License as published by the Free Software Foundation; either
41 // version 2.1 of the License, or (at your option) any later version.
42 //
43 // This library is distributed in the hope that it will be useful,
44 // but WITHOUT ANY WARRANTY; without even the implied warranty of
45 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
46 // Lesser General Public License for more details.
47 //
48 // You should have received a copy of the GNU Lesser General Public
49 // License along with this library; if not, write to the Free Software
50 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
51 //
52 // Discussion:
53 //
54 // This routine should be equivalent to the library "toupper" function.
55 //
56 // Modified:
57 //
58 // 30 September 2004
59 //
60 // Author:
61 //
62 // John Burkardt
63 //
64 // Parameters:
65 //
66 // Input, char C, the character to capitalize.
67 //
68 // Output, char CH_CAP, the capitalized character.
69 //
70 {
71 if ( 97 <= c && c <= 122 )
72 {
73 c = c - 32;
74 }
75
76 return c;
77 }
78 //*******************************************************************************
79
ch_eqi(char c1,char c2)80 bool ch_eqi ( char c1, char c2 )
81
82 //*******************************************************************************
83 //
84 // Purpose:
85 //
86 // CH_EQI is true if two characters are equal, disregarding case.
87 //
88 // License:
89 //
90 // Copyright (C) 2004 John Burkardt and Max Gunzburger
91 //
92 // This library is free software; you can redistribute it and/or
93 // modify it under the terms of the GNU Lesser General Public
94 // License as published by the Free Software Foundation; either
95 // version 2.1 of the License, or (at your option) any later version.
96 //
97 // This library is distributed in the hope that it will be useful,
98 // but WITHOUT ANY WARRANTY; without even the implied warranty of
99 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
100 // Lesser General Public License for more details.
101 //
102 // You should have received a copy of the GNU Lesser General Public
103 // License along with this library; if not, write to the Free Software
104 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
105 //
106 // Modified:
107 //
108 // 13 June 2003
109 //
110 // Author:
111 //
112 // John Burkardt
113 //
114 // Parameters:
115 //
116 // Input, char C1, C2, the characters to compare.
117 //
118 // Output, bool CH_EQI, is true if the two characters are equal,
119 // disregarding case.
120 //
121 {
122 if ( 97 <= c1 && c1 <= 122 )
123 {
124 c1 = c1 - 32;
125 }
126 if ( 97 <= c2 && c2 <= 122 )
127 {
128 c2 = c2 - 32;
129 }
130
131 return ( c1 == c2 );
132 }
133 //******************************************************************************
134
ch_to_digit(char c)135 int ch_to_digit ( char c )
136
137 //******************************************************************************
138 //
139 // Purpose:
140 //
141 // CH_TO_DIGIT returns the integer value of a base 10 digit.
142 //
143 // License:
144 //
145 // Copyright (C) 2004 John Burkardt and Max Gunzburger
146 //
147 // This library is free software; you can redistribute it and/or
148 // modify it under the terms of the GNU Lesser General Public
149 // License as published by the Free Software Foundation; either
150 // version 2.1 of the License, or (at your option) any later version.
151 //
152 // This library is distributed in the hope that it will be useful,
153 // but WITHOUT ANY WARRANTY; without even the implied warranty of
154 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
155 // Lesser General Public License for more details.
156 //
157 // You should have received a copy of the GNU Lesser General Public
158 // License along with this library; if not, write to the Free Software
159 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
160 //
161 // Example:
162 //
163 // C DIGIT
164 // --- -----
165 // '0' 0
166 // '1' 1
167 // ... ...
168 // '9' 9
169 // ' ' 0
170 // 'X' -1
171 //
172 // Modified:
173 //
174 // 13 June 2003
175 //
176 // Author:
177 //
178 // John Burkardt
179 //
180 // Parameters:
181 //
182 // Input, char C, the decimal digit, '0' through '9' or blank are legal.
183 //
184 // Output, int CH_TO_DIGIT, the corresponding integer value. If C was
185 // 'illegal', then DIGIT is -1.
186 //
187 {
188 int digit;
189
190 if ( '0' <= c && c <= '9' )
191 {
192 digit = c - '0';
193 }
194 else if ( c == ' ' )
195 {
196 digit = 0;
197 }
198 else
199 {
200 digit = -1;
201 }
202
203 return digit;
204 }
205 //******************************************************************************
206
d_epsilon(void)207 double d_epsilon ( void )
208
209 //******************************************************************************
210 //
211 // Purpose:
212 //
213 // D_EPSILON returns the round off unit for double precision arithmetic.
214 //
215 // License:
216 //
217 // Copyright (C) 2004 John Burkardt and Max Gunzburger
218 //
219 // This library is free software; you can redistribute it and/or
220 // modify it under the terms of the GNU Lesser General Public
221 // License as published by the Free Software Foundation; either
222 // version 2.1 of the License, or (at your option) any later version.
223 //
224 // This library is distributed in the hope that it will be useful,
225 // but WITHOUT ANY WARRANTY; without even the implied warranty of
226 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
227 // Lesser General Public License for more details.
228 //
229 // You should have received a copy of the GNU Lesser General Public
230 // License along with this library; if not, write to the Free Software
231 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
232 //
233 // Discussion:
234 //
235 // D_EPSILON is a number R which is a power of 2 with the property that,
236 // to the precision of the computer's arithmetic,
237 // 1 < 1 + R
238 // but
239 // 1 = ( 1 + R / 2 )
240 //
241 // Modified:
242 //
243 // 17 September 2004
244 //
245 // Author:
246 //
247 // John Burkardt
248 //
249 // Parameters:
250 //
251 // Output, double D_EPSILON, the double precision round-off unit.
252 //
253 {
254 double r;
255
256 r = 1.0E+00;
257
258 while ( 1.0E+00 < ( double ) ( 1.0E+00 + r ) )
259 {
260 r = r / 2.0E+00;
261 }
262
263 return ( 2.0E+00 * r );
264 }
265 //******************************************************************************
266
d_huge(void)267 double d_huge ( void )
268
269 //******************************************************************************
270 //
271 // Purpose:
272 //
273 // D_HUGE returns a "huge" double precision value.
274 //
275 // License:
276 //
277 // Copyright (C) 2004 John Burkardt and Max Gunzburger
278 //
279 // This library is free software; you can redistribute it and/or
280 // modify it under the terms of the GNU Lesser General Public
281 // License as published by the Free Software Foundation; either
282 // version 2.1 of the License, or (at your option) any later version.
283 //
284 // This library is distributed in the hope that it will be useful,
285 // but WITHOUT ANY WARRANTY; without even the implied warranty of
286 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
287 // Lesser General Public License for more details.
288 //
289 // You should have received a copy of the GNU Lesser General Public
290 // License along with this library; if not, write to the Free Software
291 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
292 //
293 // Discussion:
294 //
295 // HUGE_VAL is the largest representable legal double precision number, and is usually
296 // defined in math.h, or sometimes in stdlib.h.
297 //
298 // Modified:
299 //
300 // 31 August 2004
301 //
302 // Author:
303 //
304 // John Burkardt
305 //
306 // Parameters:
307 //
308 // Output, double D_HUGE, a "huge" double precision value.
309 //
310 {
311 return HUGE_VAL;
312 }
313 //*********************************************************************
314
d_max(double x,double y)315 double d_max ( double x, double y )
316
317 //*********************************************************************
318 //
319 // Purpose:
320 //
321 // D_MAX returns the maximum of two double precision values.
322 //
323 // License:
324 //
325 // Copyright (C) 2004 John Burkardt and Max Gunzburger
326 //
327 // This library is free software; you can redistribute it and/or
328 // modify it under the terms of the GNU Lesser General Public
329 // License as published by the Free Software Foundation; either
330 // version 2.1 of the License, or (at your option) any later version.
331 //
332 // This library is distributed in the hope that it will be useful,
333 // but WITHOUT ANY WARRANTY; without even the implied warranty of
334 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
335 // Lesser General Public License for more details.
336 //
337 // You should have received a copy of the GNU Lesser General Public
338 // License along with this library; if not, write to the Free Software
339 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
340 //
341 // Modified:
342 //
343 // 18 August 2004
344 //
345 // Author:
346 //
347 // John Burkardt
348 //
349 // Parameters:
350 //
351 // Input, double X, Y, the quantities to compare.
352 //
353 // Output, double D_MAX, the maximum of X and Y.
354 //
355 {
356 if ( y < x )
357 {
358 return x;
359 }
360 else
361 {
362 return y;
363 }
364 }
365 //*********************************************************************
366
d_min(double x,double y)367 double d_min ( double x, double y )
368
369 //*********************************************************************
370 //
371 // Purpose:
372 //
373 // D_MIN returns the minimum of two double precision values.
374 //
375 // License:
376 //
377 // Copyright (C) 2004 John Burkardt and Max Gunzburger
378 //
379 // This library is free software; you can redistribute it and/or
380 // modify it under the terms of the GNU Lesser General Public
381 // License as published by the Free Software Foundation; either
382 // version 2.1 of the License, or (at your option) any later version.
383 //
384 // This library is distributed in the hope that it will be useful,
385 // but WITHOUT ANY WARRANTY; without even the implied warranty of
386 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
387 // Lesser General Public License for more details.
388 //
389 // You should have received a copy of the GNU Lesser General Public
390 // License along with this library; if not, write to the Free Software
391 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
392 //
393 // Modified:
394 //
395 // 31 August 2004
396 //
397 // Author:
398 //
399 // John Burkardt
400 //
401 // Parameters:
402 //
403 // Input, double X, Y, the quantities to compare.
404 //
405 // Output, double D_MIN, the minimum of X and Y.
406 //
407 {
408 if ( y < x )
409 {
410 return y;
411 }
412 else
413 {
414 return x;
415 }
416 }
417 //******************************************************************************
418
dge_det(int n,double a[])419 double dge_det ( int n, double a[] )
420
421 //******************************************************************************
422 //
423 // Purpose:
424 //
425 // DGE_DET computes the determinant of a square matrix in DGE storage.
426 //
427 // License:
428 //
429 // Copyright (C) 2004 John Burkardt and Max Gunzburger
430 //
431 // This library is free software; you can redistribute it and/or
432 // modify it under the terms of the GNU Lesser General Public
433 // License as published by the Free Software Foundation; either
434 // version 2.1 of the License, or (at your option) any later version.
435 //
436 // This library is distributed in the hope that it will be useful,
437 // but WITHOUT ANY WARRANTY; without even the implied warranty of
438 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
439 // Lesser General Public License for more details.
440 //
441 // You should have received a copy of the GNU Lesser General Public
442 // License along with this library; if not, write to the Free Software
443 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
444 //
445 // Discussion:
446 //
447 // The DGE storage format is used for a general M by N matrix. A storage
448 // space is made for each logical entry. The two dimensional logical
449 // array is mapped to a vector, in which storage is by columns.
450 //
451 // Modified:
452 //
453 // 05 October 2004
454 //
455 // Author:
456 //
457 // John Burkardt
458 //
459 // Reference:
460 //
461 // Dongarra, Bunch, Moler, Stewart,
462 // LINPACK User's Guide,
463 // SIAM, 1979
464 //
465 // Parameters:
466 //
467 // Input, int N, the order of the matrix.
468 // N must be positive.
469 //
470 // Input/output, double A[N*N], the matrix to be analyzed.
471 // On output, the matrix has been overwritten by factorization information.
472 //
473 // Output, double DGE_DET, the determinant of the matrix.
474 //
475 {
476 double det;
477 int i;
478 int j;
479 int k;
480 int l;
481 double t;
482
483 det = 1.0;
484
485 for ( k = 1; k <= n-1; k++ )
486 {
487 //
488 // Find L, the index of the pivot row.
489 //
490 l = k;
491 for ( i = k+1; i <= n; i++ )
492 {
493 if ( fabs ( a[(l-1)+(k-1)*n] ) < fabs ( a[(i-1)+(k-1)*n] ) )
494 {
495 l = i;
496 }
497 }
498
499 det = det * a[(l-1)+(k-1)*n];
500
501 if ( a[(l-1)+(k-1)*n] == 0.0 )
502 {
503 return det;
504 }
505 //
506 // Interchange rows L and K if necessary.
507 //
508 if ( l != k )
509 {
510 t = a[(l-1)+(k-1)*n];
511 a[(l-1)+(k-1)*n] = a[(k-1)+(k-1)*n];
512 a[(k-1)+(k-1)*n] = t;
513 }
514 //
515 // Normalize the values that lie below the pivot entry A(K,K).
516 //
517 for ( i = k+1; i <= n; i++ )
518 {
519 a[(i-1)+(k-1)*n] = -a[(i-1)+(k-1)*n] / a[(k-1)+(k-1)*n];
520 }
521 //
522 // Row elimination with column indexing.
523 //
524 for ( j = k+1; j <= n; j++ )
525 {
526 if ( l != k )
527 {
528 t = a[(l-1)+(j-1)*n];
529 a[(l-1)+(j-1)*n] = a[(k-1)+(j-1)*n];
530 a[(k-1)+(j-1)*n] = t;
531 }
532
533 for ( i = k+1; i <= n; i++ )
534 {
535 a[(i-1)+(j-1)*n] = a[(i-1)+(j-1)*n]
536 + a[(i-1)+(k-1)*n] * a[(k-1)+(j-1)*n];
537 }
538 }
539 }
540
541 det = det * a[(n-1)+(n-1)*n];
542
543 return det;
544 }
545 //******************************************************************************
546
digit_to_ch(int i)547 char digit_to_ch ( int i )
548
549 //******************************************************************************
550 //
551 // Purpose:
552 //
553 // DIGIT_TO_CH returns the base 10 digit character corresponding to a digit.
554 //
555 // License:
556 //
557 // Copyright (C) 2004 John Burkardt and Max Gunzburger
558 //
559 // This library is free software; you can redistribute it and/or
560 // modify it under the terms of the GNU Lesser General Public
561 // License as published by the Free Software Foundation; either
562 // version 2.1 of the License, or (at your option) any later version.
563 //
564 // This library is distributed in the hope that it will be useful,
565 // but WITHOUT ANY WARRANTY; without even the implied warranty of
566 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
567 // Lesser General Public License for more details.
568 //
569 // You should have received a copy of the GNU Lesser General Public
570 // License along with this library; if not, write to the Free Software
571 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
572 //
573 // Example:
574 //
575 // I C
576 // ----- ---
577 // 0 '0'
578 // 1 '1'
579 // ... ...
580 // 9 '9'
581 // 10 '*'
582 // -83 '*'
583 //
584 // Modified:
585 //
586 // 17 September 2004
587 //
588 // Author:
589 //
590 // John Burkardt
591 //
592 // Parameters:
593 //
594 // Input, int I, the digit, which should be between 0 and 9.
595 //
596 // Output, char DIGIT_TO_CH, the appropriate character '0' through '9' or '*'.
597 //
598 {
599 char c;
600
601 if ( 0 <= i && i <= 9 )
602 {
603 c = '0' + i;
604 }
605 else
606 {
607 c = '*';
608 }
609
610 return c;
611 }
612 //******************************************************************************
613
dmat_in_01(int m,int n,double a[])614 bool dmat_in_01 ( int m, int n, double a[] )
615
616 //******************************************************************************
617 //
618 // Purpose:
619 //
620 // DMAT_IN_01 is TRUE if the entries of an array are in the range [0,1].
621 //
622 // Modified:
623 //
624 // 06 October 2004
625 //
626 // Author:
627 //
628 // John Burkardt
629 //
630 // Parameters:
631 //
632 // Input, int M, the number of rows in A.
633 //
634 // Input, int N, the number of columns in A.
635 //
636 // Input, double A[M*N], the matrix.
637 //
638 // Output, bool DMAT_IN_01, is TRUE if every entry of A is
639 // between 0 and 1.
640 //
641 {
642 int i;
643 int j;
644
645 for ( j = 0; j < n; j++ )
646 {
647 for ( i = 0; i < m; i++ )
648 {
649 if ( a[i+j*m] < 0.0 || 1.0 < a[i+j*m] )
650 {
651 return false;
652 }
653 }
654 }
655
656 return true;
657 }
658 //******************************************************************************
659
dmat_transpose_print(int m,int n,double a[],const char * title)660 void dmat_transpose_print ( int m, int n, double a[], const char *title )
661
662 //******************************************************************************
663 //
664 // Purpose:
665 //
666 // DMAT_TRANSPOSE_PRINT prints a real matrix, transposed.
667 //
668 // License:
669 //
670 // Copyright (C) 2004 John Burkardt
671 //
672 // This library is free software; you can redistribute it and/or
673 // modify it under the terms of the GNU Lesser General Public
674 // License as published by the Free Software Foundation; either
675 // version 2.1 of the License, or (at your option) any later version.
676 //
677 // This library is distributed in the hope that it will be useful,
678 // but WITHOUT ANY WARRANTY; without even the implied warranty of
679 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
680 // Lesser General Public License for more details.
681 //
682 // You should have received a copy of the GNU Lesser General Public
683 // License along with this library; if not, write to the Free Software
684 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
685 //
686 // Modified:
687 //
688 // 11 August 2004
689 //
690 // Author:
691 //
692 // John Burkardt
693 //
694 // Parameters:
695 //
696 // Input, int M, N, the number of rows and columns.
697 //
698 // Input, double A[M*N], an M by N matrix to be printed.
699 //
700 // Input, char *TITLE, an optional title.
701 //
702 {
703 dmat_transpose_print_some ( m, n, a, 1, 1, m, n, title );
704
705 return;
706 }
707 //******************************************************************************
708
dmat_transpose_print_some(int m,int n,double a[],int ilo,int jlo,int ihi,int jhi,const char * title)709 void dmat_transpose_print_some ( int m, int n, double a[], int ilo, int jlo,
710 int ihi, int jhi, const char *title )
711
712 //******************************************************************************
713 //
714 // Purpose:
715 //
716 // DMAT_TRANSPOSE_PRINT_SOME prints some of a real matrix, transposed.
717 //
718 // License:
719 //
720 // Copyright (C) 2004 John Burkardt
721 //
722 // This library is free software; you can redistribute it and/or
723 // modify it under the terms of the GNU Lesser General Public
724 // License as published by the Free Software Foundation; either
725 // version 2.1 of the License, or (at your option) any later version.
726 //
727 // This library is distributed in the hope that it will be useful,
728 // but WITHOUT ANY WARRANTY; without even the implied warranty of
729 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
730 // Lesser General Public License for more details.
731 //
732 // You should have received a copy of the GNU Lesser General Public
733 // License along with this library; if not, write to the Free Software
734 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
735 //
736 // Modified:
737 //
738 // 11 August 2004
739 //
740 // Author:
741 //
742 // John Burkardt
743 //
744 // Parameters:
745 //
746 // Input, int M, N, the number of rows and columns.
747 //
748 // Input, double A[M*N], an M by N matrix to be printed.
749 //
750 // Input, int ILO, JLO, the first row and column to print.
751 //
752 // Input, int IHI, JHI, the last row and column to print.
753 //
754 // Input, char *TITLE, an optional title.
755 //
756 {
757 # define INCX 5
758
759 int i;
760 int i2;
761 int i2hi;
762 int i2lo;
763 int inc;
764 int j;
765 int j2hi;
766 int j2lo;
767
768 if ( 0 < s_len_trim ( title ) )
769 {
770 cout << "\n";
771 cout << title << "\n";
772 }
773
774 for ( i2lo = i_max ( ilo, 1 ); i2lo <= i_min ( ihi, m ); i2lo = i2lo + INCX )
775 {
776 i2hi = i2lo + INCX - 1;
777 i2hi = i_min ( i2hi, m );
778 i2hi = i_min ( i2hi, ihi );
779
780 inc = i2hi + 1 - i2lo;
781
782 cout << "\n";
783 cout << " Row: ";
784 for ( i = i2lo; i <= i2hi; i++ )
785 {
786 cout << setw(7) << i << " ";
787 }
788 cout << "\n";
789 cout << " Col\n";
790
791 j2lo = i_max ( jlo, 1 );
792 j2hi = i_min ( jhi, n );
793
794 for ( j = j2lo; j <= j2hi; j++ )
795 {
796 cout << setw(5) << j << " ";
797 for ( i2 = 1; i2 <= inc; i2++ )
798 {
799 i = i2lo - 1 + i2;
800 cout << setw(14) << a[(i-1)+(j-1)*m];
801 }
802 cout << "\n";
803 }
804 }
805 cout << "\n";
806
807 return;
808 # undef INCX
809 }
810 //******************************************************************************
811
dmat_uniform_01(int m,int n,int * seed,double r[])812 void dmat_uniform_01 ( int m, int n, int *seed, double r[] )
813
814 //******************************************************************************
815 //
816 // Purpose:
817 //
818 // DMAT_UNIFORM_01 fills a double precision array with pseudorandom values.
819 //
820 // License:
821 //
822 // Copyright (C) 2004 John Burkardt and Max Gunzburger
823 //
824 // This library is free software; you can redistribute it and/or
825 // modify it under the terms of the GNU Lesser General Public
826 // License as published by the Free Software Foundation; either
827 // version 2.1 of the License, or (at your option) any later version.
828 //
829 // This library is distributed in the hope that it will be useful,
830 // but WITHOUT ANY WARRANTY; without even the implied warranty of
831 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
832 // Lesser General Public License for more details.
833 //
834 // You should have received a copy of the GNU Lesser General Public
835 // License along with this library; if not, write to the Free Software
836 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
837 //
838 // Discussion:
839 //
840 // This routine implements the recursion
841 //
842 // seed = 16807 * seed mod ( 2**31 - 1 )
843 // unif = seed / ( 2**31 - 1 )
844 //
845 // The integer arithmetic never requires more than 32 bits,
846 // including a sign bit.
847 //
848 // Modified:
849 //
850 // 11 August 2004
851 //
852 // Reference:
853 //
854 // Paul Bratley, Bennett Fox, L E Schrage,
855 // A Guide to Simulation,
856 // Springer Verlag, pages 201-202, 1983.
857 //
858 // Bennett Fox,
859 // Algorithm 647:
860 // Implementation and Relative Efficiency of Quasirandom
861 // Sequence Generators,
862 // ACM Transactions on Mathematical Software,
863 // Volume 12, Number 4, pages 362-376, 1986.
864 //
865 // Parameters:
866 //
867 // Input, int M, N, the number of rows and columns.
868 //
869 // Input/output, int *SEED, the "seed" value. Normally, this
870 // value should not be 0. On output, SEED has
871 // been updated.
872 //
873 // Output, double DMAT_UNIFORM_01[M*N], a matrix of pseudorandom values.
874 //
875 {
876 int i;
877 int j;
878 int k;
879
880 for ( j = 0; j < n; j++ )
881 {
882 for ( i = 0; i < m; i++ )
883 {
884 k = *seed / 127773;
885
886 *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;
887
888 if ( *seed < 0 )
889 {
890 *seed = *seed + 2147483647;
891 }
892 r[i+j*m] = ( double ) ( *seed ) * 4.656612875E-10;
893 }
894 }
895
896 return;
897 }
898 //******************************************************************************
899
dtable_data_read(const char * input_filename,int m,int n)900 double *dtable_data_read ( const char *input_filename, int m, int n )
901
902 //******************************************************************************
903 //
904 // Purpose:
905 //
906 // DTABLE_DATA_READ reads the data from a real TABLE file.
907 //
908 // License:
909 //
910 // Copyright (C) 2004 John Burkardt and Max Gunzburger
911 //
912 // This library is free software; you can redistribute it and/or
913 // modify it under the terms of the GNU Lesser General Public
914 // License as published by the Free Software Foundation; either
915 // version 2.1 of the License, or (at your option) any later version.
916 //
917 // This library is distributed in the hope that it will be useful,
918 // but WITHOUT ANY WARRANTY; without even the implied warranty of
919 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
920 // Lesser General Public License for more details.
921 //
922 // You should have received a copy of the GNU Lesser General Public
923 // License along with this library; if not, write to the Free Software
924 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
925 //
926 // Discussion:
927 //
928 // The file is assumed to contain one record per line.
929 //
930 // Records beginning with the '#' character are comments, and are ignored.
931 // Blank lines are also ignored.
932 //
933 // Each line that is not ignored is assumed to contain exactly (or at least)
934 // M real numbers, representing the coordinates of a point.
935 //
936 // There are assumed to be exactly (or at least) N such records.
937 //
938 // Modified:
939 //
940 // 11 October 2003
941 //
942 // Author:
943 //
944 // John Burkardt
945 //
946 // Parameters:
947 //
948 // Input, char *INPUT_FILENAME, the name of the input file.
949 //
950 // Input, int M, the number of spatial dimensions.
951 //
952 // Input, int N, the number of points. The program
953 // will stop reading data once N values have been read.
954 //
955 // Output, double DTABLE_DATA_READ[M*N], the table data.
956 //
957 {
958 bool error;
959 ifstream input;
960 int i;
961 int j;
962 char line[255];
963 double *table;
964 double *x;
965
966 input.open ( input_filename );
967
968 if ( !input )
969 {
970 cout << "\n";
971 cout << "DTABLE_DATA_READ - Fatal error!\n";
972 cout << " Could not open the input file: \"" << input_filename << "\"\n";
973 exit ( 1 );
974 }
975
976 table = new double[m*n];
977
978 x = new double[m];
979
980 j = 0;
981
982 while ( j < n )
983 {
984 input.getline ( line, sizeof ( line ) );
985
986 if ( input.eof ( ) )
987 {
988 break;
989 }
990
991 if ( line[0] == '#' || s_len_trim ( line ) == 0 )
992 {
993 continue;
994 }
995
996 error = s_to_dvec ( line, m, x );
997
998 if ( error )
999 {
1000 continue;
1001 }
1002
1003 for ( i = 0; i < m; i++ )
1004 {
1005 table[i+j*m] = x[i];
1006 }
1007 j = j + 1;
1008
1009 }
1010
1011 input.close ( );
1012
1013 delete [] x;
1014
1015 return table;
1016 }
1017 //******************************************************************************
1018
dtable_data_write(int m,int n,double table[],ofstream & output)1019 void dtable_data_write ( int m, int n, double table[], ofstream &output )
1020
1021 //******************************************************************************
1022 //
1023 // Purpose:
1024 //
1025 // DTABLE_DATA_WRITE writes real table data to a TABLE file.
1026 //
1027 // License:
1028 //
1029 // Copyright (C) 2004 John Burkardt and Max Gunzburger
1030 //
1031 // This library is free software; you can redistribute it and/or
1032 // modify it under the terms of the GNU Lesser General Public
1033 // License as published by the Free Software Foundation; either
1034 // version 2.1 of the License, or (at your option) any later version.
1035 //
1036 // This library is distributed in the hope that it will be useful,
1037 // but WITHOUT ANY WARRANTY; without even the implied warranty of
1038 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1039 // Lesser General Public License for more details.
1040 //
1041 // You should have received a copy of the GNU Lesser General Public
1042 // License along with this library; if not, write to the Free Software
1043 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
1044 //
1045 // Modified:
1046 //
1047 // 10 October 2004
1048 //
1049 // Author:
1050 //
1051 // John Burkardt
1052 //
1053 // Parameters:
1054 //
1055 // Input, int M, the spatial dimension.
1056 //
1057 // Input, int N, the number of points.
1058 //
1059 // Input, double TABLE[M*N], the table data.
1060 //
1061 // Input, ofstream &OUTPUT, a pointer to the output stream.
1062 {
1063 int i;
1064 int j;
1065 char *s;
1066
1067 for ( j = 0; j < n; j++ )
1068 {
1069 for ( i = 0; i < m; i++ )
1070 {
1071 output << setw(10) << table[i+j*m] << " ";
1072 }
1073 output << "\n";
1074 }
1075
1076 output.close ( );
1077
1078 return;
1079 }
1080 //******************************************************************************
1081
dtable_header_read(const char * input_filename,int * m,int * n)1082 void dtable_header_read ( const char *input_filename, int *m, int *n )
1083
1084 //******************************************************************************
1085 //
1086 // Purpose:
1087 //
1088 // DTABLE_HEADER_READ reads the header from a real TABLE file.
1089 //
1090 // License:
1091 //
1092 // Copyright (C) 2004 John Burkardt and Max Gunzburger
1093 //
1094 // This library is free software; you can redistribute it and/or
1095 // modify it under the terms of the GNU Lesser General Public
1096 // License as published by the Free Software Foundation; either
1097 // version 2.1 of the License, or (at your option) any later version.
1098 //
1099 // This library is distributed in the hope that it will be useful,
1100 // but WITHOUT ANY WARRANTY; without even the implied warranty of
1101 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1102 // Lesser General Public License for more details.
1103 //
1104 // You should have received a copy of the GNU Lesser General Public
1105 // License along with this library; if not, write to the Free Software
1106 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
1107 //
1108 // Modified:
1109 //
1110 // 04 June 2004
1111 //
1112 // Author:
1113 //
1114 // John Burkardt
1115 //
1116 // Parameters:
1117 //
1118 // Input, char *INPUT_FILENAME, the name of the input file.
1119 //
1120 // Output, int *M, the number of spatial dimensions.
1121 //
1122 // Output, int *N, the number of points
1123 //
1124 {
1125 *m = file_column_count ( input_filename );
1126
1127 if ( *m <= 0 )
1128 {
1129 cout << "\n";
1130 cout << "DTABLE_HEADER_READ - Fatal error!\n";
1131 cout << " FILE_COLUMN_COUNT failed.\n";
1132 *n = -1;
1133 return;
1134 }
1135
1136 *n = file_row_count ( input_filename );
1137
1138 if ( *n <= 0 )
1139 {
1140 cout << "\n";
1141 cout << "DTABLE_HEADER_READ - Fatal error!\n";
1142 cout << " FILE_ROW_COUNT failed.\n";
1143 return;
1144 }
1145
1146 return;
1147 }
1148 //******************************************************************************
1149
dtable_header_write(int m,int n,char * output_filename,ofstream & output)1150 void dtable_header_write ( int m, int n, char *output_filename, ofstream &output )
1151
1152 //******************************************************************************
1153 //
1154 // Purpose:
1155 //
1156 // DTABLE_HEADER_WRITE writes the header of a TABLE file.
1157 //
1158 // License:
1159 //
1160 // Copyright (C) 2004 John Burkardt and Max Gunzburger
1161 //
1162 // This library is free software; you can redistribute it and/or
1163 // modify it under the terms of the GNU Lesser General Public
1164 // License as published by the Free Software Foundation; either
1165 // version 2.1 of the License, or (at your option) any later version.
1166 //
1167 // This library is distributed in the hope that it will be useful,
1168 // but WITHOUT ANY WARRANTY; without even the implied warranty of
1169 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1170 // Lesser General Public License for more details.
1171 //
1172 // You should have received a copy of the GNU Lesser General Public
1173 // License along with this library; if not, write to the Free Software
1174 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
1175 //
1176 // Modified:
1177 //
1178 // 10 October 2004
1179 //
1180 // Author:
1181 //
1182 // John Burkardt
1183 //
1184 // Parameters:
1185 //
1186 // Input, int M, the spatial dimension.
1187 //
1188 // Input, int N, the number of points.
1189 //
1190 // Input, char *OUTPUT_FILENAME, the output filename.
1191 //
1192 // Input, ofstream &OUTPUT, the output stream.
1193 //
1194 {
1195 double roundoff;
1196 char *s;
1197
1198 s = timestring ( );
1199 roundoff = d_epsilon ( );
1200
1201 output << "# " << output_filename << "\n";
1202 output << "# created at " << s << "\n";
1203 output << "#\n";
1204 output << "# Spatial dimension M = " << m << "\n";
1205 output << "# Number of points N = " << n << "\n";
1206 output << "# EPSILON (unit roundoff) = " << roundoff << "\n";
1207 output << "#\n";
1208
1209 return;
1210 }
1211 //******************************************************************************
1212
dvec_sort_heap_index_a(int n,double a[])1213 int *dvec_sort_heap_index_a ( int n, double a[] )
1214
1215 //******************************************************************************
1216 //
1217 // Purpose:
1218 //
1219 // DVEC_SORT_HEAP_INDEX_A does an indexed heap ascending sort of a real vector.
1220 //
1221 // License:
1222 //
1223 // Copyright (C) 2004 John Burkardt and Max Gunzburger
1224 //
1225 // This library is free software; you can redistribute it and/or
1226 // modify it under the terms of the GNU Lesser General Public
1227 // License as published by the Free Software Foundation; either
1228 // version 2.1 of the License, or (at your option) any later version.
1229 //
1230 // This library is distributed in the hope that it will be useful,
1231 // but WITHOUT ANY WARRANTY; without even the implied warranty of
1232 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1233 // Lesser General Public License for more details.
1234 //
1235 // You should have received a copy of the GNU Lesser General Public
1236 // License along with this library; if not, write to the Free Software
1237 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
1238 //
1239 // Discussion:
1240 //
1241 // The sorting is not actually carried out. Rather an index array is
1242 // created which defines the sorting. This array may be used to sort
1243 // or index the array, or to sort or index related arrays keyed on the
1244 // original array.
1245 //
1246 // Once the index array is computed, the sorting can be carried out
1247 // "implicitly:
1248 //
1249 // A(INDX(I)), I = 1 to N is sorted,
1250 //
1251 // after which A(I), I = 1 to N is sorted.
1252 //
1253 // Modified:
1254 //
1255 // 10 October 2004
1256 //
1257 // Author:
1258 //
1259 // John Burkardt
1260 //
1261 // Parameters:
1262 //
1263 // Input, int N, the number of entries in the array.
1264 //
1265 // Input, double A[N], an array to be index-sorted.
1266 //
1267 // Output, int DVEC_SORT_HEAP_INDEX_A[N], contains the sort index. The
1268 // I-th element of the sorted array is A(INDX(I)).
1269 //
1270 {
1271 double aval;
1272 int i;
1273 int *indx;
1274 int indxt;
1275 int ir;
1276 int j;
1277 int l;
1278
1279 indx = new int[n];
1280
1281 for ( i = 1; i <= n; i++ )
1282 {
1283 indx[i-1] = i;
1284 }
1285
1286 l = n / 2 + 1;
1287 ir = n;
1288
1289 for ( ; ; )
1290 {
1291 if ( 1 < l )
1292 {
1293 l = l - 1;
1294 indxt = indx[l-1];
1295 aval = a[indxt-1];
1296 }
1297 else
1298 {
1299 indxt = indx[ir-1];
1300 aval = a[indxt-1];
1301 indx[ir-1] = indx[0];
1302 ir = ir - 1;
1303
1304 if ( ir == 1 )
1305 {
1306 indx[0] = indxt;
1307 for ( i = 0; i < n; i++ )
1308 {
1309 indx[i] = indx[i] - 1;
1310 }
1311 return indx;
1312 }
1313
1314 }
1315
1316 i = l;
1317 j = l + l;
1318
1319 while ( j <= ir )
1320 {
1321 if ( j < ir )
1322 {
1323 if ( a[indx[j-1]-1] < a[indx[j]-1] )
1324 {
1325 j = j + 1;
1326 }
1327 }
1328
1329 if ( aval < a[indx[j-1]-1] )
1330 {
1331 indx[i-1] = indx[j-1];
1332 i = j;
1333 j = j + j;
1334 }
1335 else
1336 {
1337 j = ir + 1;
1338 }
1339 }
1340 indx[i-1] = indxt;
1341 }
1342 }
1343 //******************************************************************************
1344
dvec_uniform_01(int n,int * seed,double r[])1345 void dvec_uniform_01 ( int n, int *seed, double r[] )
1346
1347 //******************************************************************************
1348 //
1349 // Purpose:
1350 //
1351 // DVEC_UNIFORM_01 fills a double precision vector with pseudorandom values.
1352 //
1353 // License:
1354 //
1355 // Copyright (C) 2004 John Burkardt and Max Gunzburger
1356 //
1357 // This library is free software; you can redistribute it and/or
1358 // modify it under the terms of the GNU Lesser General Public
1359 // License as published by the Free Software Foundation; either
1360 // version 2.1 of the License, or (at your option) any later version.
1361 //
1362 // This library is distributed in the hope that it will be useful,
1363 // but WITHOUT ANY WARRANTY; without even the implied warranty of
1364 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1365 // Lesser General Public License for more details.
1366 //
1367 // You should have received a copy of the GNU Lesser General Public
1368 // License along with this library; if not, write to the Free Software
1369 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
1370 //
1371 // Discussion:
1372 //
1373 // This routine implements the recursion
1374 //
1375 // seed = 16807 * seed mod ( 2**31 - 1 )
1376 // unif = seed / ( 2**31 - 1 )
1377 //
1378 // The integer arithmetic never requires more than 32 bits,
1379 // including a sign bit.
1380 //
1381 // Modified:
1382 //
1383 // 19 August 2004
1384 //
1385 // Reference:
1386 //
1387 // Paul Bratley, Bennett Fox, L E Schrage,
1388 // A Guide to Simulation,
1389 // Springer Verlag, pages 201-202, 1983.
1390 //
1391 // Bennett Fox,
1392 // Algorithm 647:
1393 // Implementation and Relative Efficiency of Quasirandom
1394 // Sequence Generators,
1395 // ACM Transactions on Mathematical Software,
1396 // Volume 12, Number 4, pages 362-376, 1986.
1397 //
1398 // Parameters:
1399 //
1400 // Input, int N, the number of entries in the vector.
1401 //
1402 // Input/output, int *SEED, the "seed" value. Normally, this
1403 // value should not be 0. On output, SEED has been updated.
1404 //
1405 // Output, double R[N], the vector of pseudorandom values.
1406 //
1407 {
1408 int i;
1409 int k;
1410
1411 for ( i = 0; i < n; i++ )
1412 {
1413 k = *seed / 127773;
1414
1415 *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;
1416
1417 if ( *seed < 0 )
1418 {
1419 *seed = *seed + 2147483647;
1420 }
1421
1422 r[i] = ( double ) ( *seed ) * 4.656612875E-10;
1423 }
1424
1425 return;
1426 }
1427 //******************************************************************************
1428
file_column_count(const char * input_filename)1429 int file_column_count ( const char *input_filename )
1430
1431 //******************************************************************************
1432 //
1433 // Purpose:
1434 //
1435 // FILE_COLUMN_COUNT counts the number of columns in the first line of a file.
1436 //
1437 // License:
1438 //
1439 // Copyright (C) 2004 John Burkardt and Max Gunzburger
1440 //
1441 // This library is free software; you can redistribute it and/or
1442 // modify it under the terms of the GNU Lesser General Public
1443 // License as published by the Free Software Foundation; either
1444 // version 2.1 of the License, or (at your option) any later version.
1445 //
1446 // This library is distributed in the hope that it will be useful,
1447 // but WITHOUT ANY WARRANTY; without even the implied warranty of
1448 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1449 // Lesser General Public License for more details.
1450 //
1451 // You should have received a copy of the GNU Lesser General Public
1452 // License along with this library; if not, write to the Free Software
1453 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
1454 //
1455 // Discussion:
1456 //
1457 // The file is assumed to be a simple text file.
1458 //
1459 // Most lines of the file is presumed to consist of COLUMN_NUM words, separated
1460 // by spaces. There may also be some blank lines, and some comment lines,
1461 // which have a "#" in column 1.
1462 //
1463 // The routine tries to find the first non-comment non-blank line and
1464 // counts the number of words in that line.
1465 //
1466 // If all lines are blanks or comments, it goes back and tries to analyze
1467 // a comment line.
1468 //
1469 // Modified:
1470 //
1471 // 13 June 2003
1472 //
1473 // Author:
1474 //
1475 // John Burkardt
1476 //
1477 // Parameters:
1478 //
1479 // Input, char *INPUT_FILENAME, the name of the file.
1480 //
1481 // Output, int FILE_COLUMN_COUNT, the number of columns assumed
1482 // to be in the file.
1483 //
1484 {
1485 int column_num;
1486 ifstream input;
1487 bool got_one;
1488 char line[256];
1489 //
1490 // Open the file.
1491 //
1492 input.open ( input_filename );
1493
1494 if ( !input )
1495 {
1496 column_num = -1;
1497 cout << "\n";
1498 cout << "FILE_COLUMN_COUNT - Fatal error!\n";
1499 cout << " Could not open the file:\n";
1500 cout << " \"" << input_filename << "\"\n";
1501 return column_num;
1502 }
1503 //
1504 // Read one line, but skip blank lines and comment lines.
1505 //
1506 got_one = false;
1507
1508 for ( ; ; )
1509 {
1510 input.getline ( line, sizeof ( line ) );
1511
1512 if ( input.eof ( ) )
1513 {
1514 break;
1515 }
1516
1517 if ( s_len_trim ( line ) == 0 )
1518 {
1519 continue;
1520 }
1521
1522 if ( line[0] == '#' )
1523 {
1524 continue;
1525 }
1526
1527 got_one = true;
1528 break;
1529
1530 }
1531
1532 if ( !got_one )
1533 {
1534 input.close ( );
1535
1536 input.open ( input_filename );
1537
1538 for ( ; ; )
1539 {
1540 input.getline ( line, sizeof ( line ) );
1541
1542 if ( input.eof ( ) )
1543 {
1544 break;
1545 }
1546
1547 if ( s_len_trim ( line ) == 0 )
1548 {
1549 continue;
1550 }
1551
1552 got_one = true;
1553 break;
1554
1555 }
1556
1557 }
1558
1559 input.close ( );
1560
1561 if ( !got_one )
1562 {
1563 cout << "\n";
1564 cout << "FILE_COLUMN_COUNT - Warning!\n";
1565 cout << " The file does not seem to contain any data.\n";
1566 return -1;
1567 }
1568
1569 column_num = word_count ( line );
1570
1571 return column_num;
1572 }
1573 //******************************************************************************
1574
file_name_ext_get(const char * file_name,int * i,int * j)1575 void file_name_ext_get ( const char *file_name, int *i, int *j )
1576
1577 //******************************************************************************
1578 //
1579 // Purpose:
1580 //
1581 // FILE_NAME_EXT_GET determines the "extension" of a file name.
1582 //
1583 // License:
1584 //
1585 // Copyright (C) 2004 John Burkardt and Max Gunzburger
1586 //
1587 // This library is free software; you can redistribute it and/or
1588 // modify it under the terms of the GNU Lesser General Public
1589 // License as published by the Free Software Foundation; either
1590 // version 2.1 of the License, or (at your option) any later version.
1591 //
1592 // This library is distributed in the hope that it will be useful,
1593 // but WITHOUT ANY WARRANTY; without even the implied warranty of
1594 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1595 // Lesser General Public License for more details.
1596 //
1597 // You should have received a copy of the GNU Lesser General Public
1598 // License along with this library; if not, write to the Free Software
1599 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
1600 //
1601 // Definition:
1602 //
1603 // The "extension" of a filename is the string of characters
1604 // that appears after the LAST period in the name. A file
1605 // with no period, or with a period as the last character
1606 // in the name, has a "null" extension.
1607 //
1608 // Note:
1609 //
1610 // Blanks are unusual in filenames. This routine ignores all
1611 // trailing blanks, but will treat initial or internal blanks
1612 // as regular characters acceptable in a file name.
1613 //
1614 // Examples:
1615 //
1616 // FILE_NAME I J
1617 //
1618 // bob.for 3 6
1619 // N.B.C.D 5 6
1620 // Naomi. 5 5
1621 // Arthur -1 -1
1622 // .com 0 0
1623 //
1624 // Modified:
1625 //
1626 // 10 October 2004
1627 //
1628 // Author:
1629 //
1630 // John Burkardt
1631 //
1632 // Parameters:
1633 //
1634 // Input, char *FILE_NAME, a file name to be examined.
1635 //
1636 // Output, int *I, *J, the indices of the first and last characters
1637 // in the file extension.
1638 //
1639 // If no period occurs in FILE_NAME, then
1640 // I = J = -1;
1641 // Otherwise,
1642 // I is the position of the LAST period in FILE_NAME, and J is the
1643 // position of the last nonblank character following the period.
1644 //
1645 {
1646 *i = s_index_last_c ( file_name, '.' );
1647
1648 if ( *i == -1 )
1649 {
1650 *j = -1;
1651 }
1652 else
1653 {
1654 *j = s_len_trim ( file_name ) - 1;
1655 }
1656
1657 return;
1658 }
1659 //******************************************************************************
1660
file_name_ext_swap(const char * file_name,const char * ext)1661 char *file_name_ext_swap ( const char *file_name, const char *ext )
1662
1663 //******************************************************************************
1664 //
1665 // Purpose:
1666 //
1667 // FILE_NAME_EXT_SWAP replaces the current "extension" of a file name.
1668 //
1669 // License:
1670 //
1671 // Copyright (C) 2004 John Burkardt and Max Gunzburger
1672 //
1673 // This library is free software; you can redistribute it and/or
1674 // modify it under the terms of the GNU Lesser General Public
1675 // License as published by the Free Software Foundation; either
1676 // version 2.1 of the License, or (at your option) any later version.
1677 //
1678 // This library is distributed in the hope that it will be useful,
1679 // but WITHOUT ANY WARRANTY; without even the implied warranty of
1680 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1681 // Lesser General Public License for more details.
1682 //
1683 // You should have received a copy of the GNU Lesser General Public
1684 // License along with this library; if not, write to the Free Software
1685 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
1686 //
1687 // Definition:
1688 //
1689 // The "extension" of a filename is the string of characters
1690 // that appears after the LAST period in the name. A file
1691 // with no period, or with a period as the last character
1692 // in the name, has a "null" extension.
1693 //
1694 // Examples:
1695 //
1696 // Input Output
1697 // ================ ==================
1698 // FILE_NAME EXT FILE_NAME_EXT_SWAP
1699 //
1700 // bob.for obj bob.obj
1701 // bob.bob.bob txt bob.bob.txt
1702 // bob yak bob.yak
1703 //
1704 // Modified:
1705 //
1706 // 10 October 2004
1707 //
1708 // Author:
1709 //
1710 // John Burkardt
1711 //
1712 // Parameters:
1713 //
1714 // Input, char *FILE_NAME, a file name.
1715 //
1716 // Input, char *EXT, the extension to be used on the output
1717 // copy of FILE_NAME, replacing the current extension if any.
1718 //
1719 // Output, char *FILE_NAME_EXT_SWAP, the file name with the new extension.
1720 //
1721 {
1722 int ext_len;
1723 int file_name_len;
1724 char *file_name2;
1725 int i;
1726 int j;
1727 int period;
1728 char *s;
1729 const char *t;
1730 //
1731 ext_len = s_len_trim ( ext );
1732 file_name_len = s_len_trim ( file_name );
1733 //
1734 // Call FILE_NAME_EXT_GET to get I, the index of the period character
1735 // beginning the file extension.
1736 //
1737 file_name_ext_get ( file_name, &period, &j );
1738 //
1739 // If there is no extension, point to where the "." should be.
1740 //
1741 if ( period == -1 )
1742 {
1743 period = file_name_len;
1744 }
1745
1746 file_name2 = new char[ period + 1 + ext_len + 1];
1747
1748 s = file_name2;
1749 t = file_name;
1750 for ( i = 0; i < period; i++ )
1751 {
1752 *s = *t;
1753 s++;
1754 t++;
1755 }
1756 *s = '.';
1757 s++;
1758
1759 t = ext;
1760 for ( i = 0; i < ext_len; i++ )
1761 {
1762 *s = *t;
1763 s++;
1764 t++;
1765 }
1766
1767 *s = '\0';
1768
1769 return file_name2;
1770 }
1771 //******************************************************************************
1772
file_row_count(const char * input_filename)1773 int file_row_count ( const char *input_filename )
1774
1775 //******************************************************************************
1776 //
1777 // Purpose:
1778 //
1779 // FILE_ROW_COUNT counts the number of row records in a file.
1780 //
1781 // License:
1782 //
1783 // Copyright (C) 2004 John Burkardt and Max Gunzburger
1784 //
1785 // This library is free software; you can redistribute it and/or
1786 // modify it under the terms of the GNU Lesser General Public
1787 // License as published by the Free Software Foundation; either
1788 // version 2.1 of the License, or (at your option) any later version.
1789 //
1790 // This library is distributed in the hope that it will be useful,
1791 // but WITHOUT ANY WARRANTY; without even the implied warranty of
1792 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1793 // Lesser General Public License for more details.
1794 //
1795 // You should have received a copy of the GNU Lesser General Public
1796 // License along with this library; if not, write to the Free Software
1797 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
1798 //
1799 // Discussion:
1800 //
1801 // It does not count lines that are blank, or that begin with a
1802 // comment symbol '#'.
1803 //
1804 // Modified:
1805 //
1806 // 13 June 2003
1807 //
1808 // Author:
1809 //
1810 // John Burkardt
1811 //
1812 // Parameters:
1813 //
1814 // Input, char *INPUT_FILENAME, the name of the input file.
1815 //
1816 // Output, int FILE_ROW_COUNT, the number of rows found.
1817 //
1818 {
1819 int bad_num;
1820 int comment_num;
1821 ifstream input;
1822 int i;
1823 char line[100];
1824 int record_num;
1825 int row_num;
1826
1827 row_num = 0;
1828 comment_num = 0;
1829 record_num = 0;
1830 bad_num = 0;
1831
1832 input.open ( input_filename );
1833
1834 if ( !input )
1835 {
1836 cout << "\n";
1837 cout << "FILE_ROW_COUNT - Fatal error!\n";
1838 cout << " Could not open the input file: \"" << input_filename << "\"\n";
1839 exit ( 1 );
1840 }
1841
1842 for ( ; ; )
1843 {
1844 input.getline ( line, sizeof ( line ) );
1845
1846 if ( input.eof ( ) )
1847 {
1848 break;
1849 }
1850
1851 record_num = record_num + 1;
1852
1853 if ( line[0] == '#' )
1854 {
1855 comment_num = comment_num + 1;
1856 continue;
1857 }
1858
1859 if ( s_len_trim ( line ) == 0 )
1860 {
1861 comment_num = comment_num + 1;
1862 continue;
1863 }
1864
1865 row_num = row_num + 1;
1866
1867 }
1868
1869 input.close ( );
1870
1871 return row_num;
1872 }
1873 //******************************************************************************
1874
find_closest(int ndim,int n,int sample_num,double s[],double r[],int nearest[])1875 void find_closest ( int ndim, int n, int sample_num, double s[], double r[],
1876 int nearest[] )
1877
1878 //******************************************************************************
1879 //
1880 // Purpose:
1881 //
1882 // FIND_CLOSEST finds the nearest R point to each S point.
1883 //
1884 // License:
1885 //
1886 // Copyright (C) 2004 John Burkardt and Max Gunzburger
1887 //
1888 // This library is free software; you can redistribute it and/or
1889 // modify it under the terms of the GNU Lesser General Public
1890 // License as published by the Free Software Foundation; either
1891 // version 2.1 of the License, or (at your option) any later version.
1892 //
1893 // This library is distributed in the hope that it will be useful,
1894 // but WITHOUT ANY WARRANTY; without even the implied warranty of
1895 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1896 // Lesser General Public License for more details.
1897 //
1898 // You should have received a copy of the GNU Lesser General Public
1899 // License along with this library; if not, write to the Free Software
1900 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
1901 //
1902 // Discussion:
1903 //
1904 // This routine finds the closest Voronoi cell generator by checking every
1905 // one. For problems with many cells, this process can take the bulk
1906 // of the CPU time. Other approaches, which group the cell generators into
1907 // bins, can run faster by a large factor.
1908 //
1909 // HUGE_VAL is the largest representable legal real number, and is usually
1910 // defined in math.h, or sometimes in stdlib.h.
1911 //
1912 // Modified:
1913 //
1914 // 02 August 2004
1915 //
1916 // Author:
1917 //
1918 // John Burkardt
1919 //
1920 // Parameters:
1921 //
1922 // Input, int NDIM, the spatial dimension.
1923 //
1924 // Input, int N, the number of cell generators.
1925 //
1926 // Input, int SAMPLE_NUM, the number of sample points.
1927 //
1928 // Input, double S[NDIM*SAMPLE_NUM], the points to be checked.
1929 //
1930 // Input, double R[NDIM*N], the cell generators.
1931 //
1932 // Output, int NEAREST[SAMPLE_NUM], the index of the nearest
1933 // cell generators.
1934 //
1935 {
1936 double dist_sq_min;
1937 double dist_sq;
1938 int i;
1939 int jr;
1940 int js;
1941 //
1942 for ( js = 0; js < sample_num; js++ )
1943 {
1944 dist_sq_min = HUGE_VAL;
1945 nearest[js] = -1;
1946
1947 for ( jr = 0; jr < n; jr++ )
1948 {
1949 dist_sq = 0.0;
1950 for ( i = 0; i < ndim; i++ )
1951 {
1952 dist_sq = dist_sq + ( s[i+js*ndim] - r[i+jr*ndim] )
1953 * ( s[i+js*ndim] - r[i+jr*ndim] );
1954 }
1955
1956 if ( jr == 0 || dist_sq < dist_sq_min )
1957 {
1958 dist_sq_min = dist_sq;
1959 nearest[js] = jr;
1960 }
1961
1962 }
1963
1964 }
1965
1966 return;
1967 }
1968 //******************************************************************************
1969
get_seed(void)1970 int get_seed ( void )
1971
1972 //******************************************************************************
1973 //
1974 // Purpose:
1975 //
1976 // GET_SEED returns a random seed for the random number generator.
1977 //
1978 // License:
1979 //
1980 // Copyright (C) 2004 John Burkardt and Max Gunzburger
1981 //
1982 // This library is free software; you can redistribute it and/or
1983 // modify it under the terms of the GNU Lesser General Public
1984 // License as published by the Free Software Foundation; either
1985 // version 2.1 of the License, or (at your option) any later version.
1986 //
1987 // This library is distributed in the hope that it will be useful,
1988 // but WITHOUT ANY WARRANTY; without even the implied warranty of
1989 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1990 // Lesser General Public License for more details.
1991 //
1992 // You should have received a copy of the GNU Lesser General Public
1993 // License along with this library; if not, write to the Free Software
1994 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
1995 //
1996 // Modified:
1997 //
1998 // 15 September 2003
1999 //
2000 // Author:
2001 //
2002 // John Burkardt
2003 //
2004 // Parameters:
2005 //
2006 // Output, int GET_SEED, a random seed value.
2007 //
2008 {
2009 # define I_MAX 2147483647
2010 time_t clock;
2011 int i;
2012 int ihour;
2013 int imin;
2014 int isec;
2015 int seed;
2016 struct tm *lt;
2017 time_t tloc;
2018 //
2019 // If the internal seed is 0, generate a value based on the time.
2020 //
2021 clock = time ( &tloc );
2022 lt = localtime ( &clock );
2023 //
2024 // Hours is 1, 2, ..., 12.
2025 //
2026 ihour = lt->tm_hour;
2027
2028 if ( 12 < ihour )
2029 {
2030 ihour = ihour - 12;
2031 }
2032 //
2033 // Move Hours to 0, 1, ..., 11
2034 //
2035 ihour = ihour - 1;
2036
2037 imin = lt->tm_min;
2038
2039 isec = lt->tm_sec;
2040
2041 seed = isec + 60 * ( imin + 60 * ihour );
2042 //
2043 // We want values in [1,43200], not [0,43199].
2044 //
2045 seed = seed + 1;
2046 //
2047 // Remap ISEED from [1,43200] to [1,IMAX].
2048 //
2049 seed = ( int )
2050 ( ( ( double ) seed )
2051 * ( ( double ) I_MAX ) / ( 60.0E+00 * 60.0E+00 * 12.0E+00 ) );
2052 //
2053 // Never use a seed of 0.
2054 //
2055 if ( seed == 0 )
2056 {
2057 seed = 1;
2058 }
2059
2060 return seed;
2061 #undef I_MAX
2062 }
2063 //**********************************************************************
2064
halham_leap_check(int ndim,int leap[])2065 bool halham_leap_check ( int ndim, int leap[] )
2066
2067 //**********************************************************************
2068 //
2069 // Purpose:
2070 //
2071 // HALHAM_LEAP_CHECK checks LEAP for a Halton or Hammersley sequence.
2072 //
2073 // License:
2074 //
2075 // Copyright (C) 2004 John Burkardt and Max Gunzburger
2076 //
2077 // This library is free software; you can redistribute it and/or
2078 // modify it under the terms of the GNU Lesser General Public
2079 // License as published by the Free Software Foundation; either
2080 // version 2.1 of the License, or (at your option) any later version.
2081 //
2082 // This library is distributed in the hope that it will be useful,
2083 // but WITHOUT ANY WARRANTY; without even the implied warranty of
2084 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2085 // Lesser General Public License for more details.
2086 //
2087 // You should have received a copy of the GNU Lesser General Public
2088 // License along with this library; if not, write to the Free Software
2089 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
2090 //
2091 // Modified:
2092 //
2093 // 17 September 2004
2094 //
2095 // Author:
2096 //
2097 // John Burkardt
2098 //
2099 // Parameters:
2100 //
2101 // Input, int LEAP[HALTON_NDIM], the successive jumps in the sequence.
2102 // Each entry must be greater than 0.
2103 //
2104 // Output, bool HALHAM_LEAP_CHECK, is true if LEAP is legal.
2105 //
2106 {
2107 int i;
2108 bool value;
2109
2110 value = true;
2111
2112 for ( i = 0; i < ndim; i++ )
2113 {
2114 if ( leap[i] < 1 )
2115 {
2116 cout << "\n";
2117 cout << "HALHAM_LEAP_CHECK - Fatal error!\n";
2118 cout << " Leap entries must be greater than 0.\n";
2119 cout << " leap[" << i << "] = " << leap[i] << "\n";
2120 value = false;
2121 break;
2122 }
2123 }
2124
2125 return value;
2126 }
2127 //**********************************************************************
2128
halham_n_check(int n)2129 bool halham_n_check ( int n )
2130
2131 //**********************************************************************
2132 //
2133 // Purpose:
2134 //
2135 // HALHAM_N_CHECK checks N for a Halton or Hammersley sequence.
2136 //
2137 // License:
2138 //
2139 // Copyright (C) 2004 John Burkardt and Max Gunzburger
2140 //
2141 // This library is free software; you can redistribute it and/or
2142 // modify it under the terms of the GNU Lesser General Public
2143 // License as published by the Free Software Foundation; either
2144 // version 2.1 of the License, or (at your option) any later version.
2145 //
2146 // This library is distributed in the hope that it will be useful,
2147 // but WITHOUT ANY WARRANTY; without even the implied warranty of
2148 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2149 // Lesser General Public License for more details.
2150 //
2151 // You should have received a copy of the GNU Lesser General Public
2152 // License along with this library; if not, write to the Free Software
2153 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
2154 //
2155 // Modified:
2156 //
2157 // 17 September 2004
2158 //
2159 // Author:
2160 //
2161 // John Burkardt
2162 //
2163 // Parameters:
2164 //
2165 // Input, int N, the number of elements of the subsequence.
2166 // N must be positive.
2167 //
2168 // Output, bool HALHAM_N_CHECK, is true if N is legal.
2169 //
2170 {
2171 bool value;
2172
2173 if ( n < 1 )
2174 {
2175 cout << "\n";
2176 cout << "HALHAM_N_CHECK - Fatal error!\n";
2177 cout << " N < 0.";
2178 cout << " N = " << n << "\n";
2179 value = false;
2180 }
2181 else
2182 {
2183 value = true;
2184 }
2185
2186 return value;
2187 }
2188 //**********************************************************************
2189
halham_ndim_check(int ndim)2190 bool halham_ndim_check ( int ndim )
2191
2192 //**********************************************************************
2193 //
2194 // Purpose:
2195 //
2196 // HALHAM_NDIM_CHECK checks NDIM for a Halton or Hammersley sequence.
2197 //
2198 // License:
2199 //
2200 // Copyright (C) 2004 John Burkardt and Max Gunzburger
2201 //
2202 // This library is free software; you can redistribute it and/or
2203 // modify it under the terms of the GNU Lesser General Public
2204 // License as published by the Free Software Foundation; either
2205 // version 2.1 of the License, or (at your option) any later version.
2206 //
2207 // This library is distributed in the hope that it will be useful,
2208 // but WITHOUT ANY WARRANTY; without even the implied warranty of
2209 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2210 // Lesser General Public License for more details.
2211 //
2212 // You should have received a copy of the GNU Lesser General Public
2213 // License along with this library; if not, write to the Free Software
2214 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
2215 //
2216 // Modified:
2217 //
2218 // 17 September 2004
2219 //
2220 // Author:
2221 //
2222 // John Burkardt
2223 //
2224 // Parameters:
2225 //
2226 // Input, int NDIM, the spatial dimension.
2227 // NDIM must be positive.
2228 //
2229 // Output, bool HALHAM_NDIM_CHECK, is true if NDIM is legal.
2230 //
2231 {
2232 bool value;
2233
2234 if ( ndim < 1 )
2235 {
2236 cout << "\n";
2237 cout << "HALHAM_NDIM_CHECK - Fatal error!\n";
2238 cout << " NDIM < 0.";
2239 cout << " NDIM = " << ndim << "\n";
2240 value = false;
2241 }
2242 else
2243 {
2244 value = true;
2245 }
2246
2247 return value;
2248 }
2249 //**********************************************************************
2250
halham_seed_check(int ndim,int seed[])2251 bool halham_seed_check ( int ndim, int seed[] )
2252
2253 //**********************************************************************
2254 //
2255 // Purpose:
2256 //
2257 // HALHAM_SEED_CHECK checks SEED for a Halton or Hammersley sequence.
2258 //
2259 // License:
2260 //
2261 // Copyright (C) 2004 John Burkardt and Max Gunzburger
2262 //
2263 // This library is free software; you can redistribute it and/or
2264 // modify it under the terms of the GNU Lesser General Public
2265 // License as published by the Free Software Foundation; either
2266 // version 2.1 of the License, or (at your option) any later version.
2267 //
2268 // This library is distributed in the hope that it will be useful,
2269 // but WITHOUT ANY WARRANTY; without even the implied warranty of
2270 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2271 // Lesser General Public License for more details.
2272 //
2273 // You should have received a copy of the GNU Lesser General Public
2274 // License along with this library; if not, write to the Free Software
2275 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
2276 //
2277 // Modified:
2278 //
2279 // 17 September 2004
2280 //
2281 // Author:
2282 //
2283 // John Burkardt
2284 //
2285 // Parameters:
2286 //
2287 // Input, int SEED[NDIM], the sequence index
2288 // corresponding to STEP = 0. Each entry must be 0 or greater.
2289 //
2290 // Output, bool HALHAM_SEED_CHECK is true if SEED is legal.
2291 //
2292 {
2293 int i;
2294 bool value;
2295
2296 value = true;
2297
2298 for ( i = 0; i < ndim; i++ )
2299 {
2300 if ( seed[i] < 0 )
2301 {
2302 cout << "\n";
2303 cout << "HALHAM_SEED_CHECK - Fatal error!\n";
2304 cout << " SEED entries must be nonnegative.\n";
2305 cout << " seed[" << i << "] = " << seed[i] << "\n";
2306 value = false;
2307 break;
2308 }
2309 }
2310
2311 return value;
2312 }
2313 //**********************************************************************
2314
halham_step_check(int step)2315 bool halham_step_check ( int step )
2316
2317 //**********************************************************************
2318 //
2319 // Purpose:
2320 //
2321 // HALHAM_STEP_CHECK checks STEP for a Halton or Hammersley sequence.
2322 //
2323 // License:
2324 //
2325 // Copyright (C) 2004 John Burkardt and Max Gunzburger
2326 //
2327 // This library is free software; you can redistribute it and/or
2328 // modify it under the terms of the GNU Lesser General Public
2329 // License as published by the Free Software Foundation; either
2330 // version 2.1 of the License, or (at your option) any later version.
2331 //
2332 // This library is distributed in the hope that it will be useful,
2333 // but WITHOUT ANY WARRANTY; without even the implied warranty of
2334 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2335 // Lesser General Public License for more details.
2336 //
2337 // You should have received a copy of the GNU Lesser General Public
2338 // License along with this library; if not, write to the Free Software
2339 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
2340 //
2341 // Modified:
2342 //
2343 // 17 September 2004
2344 //
2345 // Author:
2346 //
2347 // John Burkardt
2348 //
2349 // Parameters:
2350 //
2351 // Input, int STEP, the index of the subsequence element.
2352 // STEP must be 1 or greater.
2353 //
2354 // Output, bool HALHAM_STEP_CHECK is true if STEP is legal.
2355 //
2356 {
2357 int i;
2358 bool value;
2359
2360 if ( step < 0 )
2361 {
2362 cout << "\n";
2363 cout << "HALHAM_STEP_CHECK - Fatal error!\n";
2364 cout << " STEP < 0.";
2365 cout << " STEP = " << step << "\n";
2366 value = false;
2367 }
2368 else
2369 {
2370 value = true;
2371 }
2372
2373 return value;
2374 }
2375 //******************************************************************************
2376
halham_write(int ndim,int n,int step,int seed[],int leap[],int base[],double r[],char * file_out_name)2377 void halham_write ( int ndim, int n, int step, int seed[], int leap[],
2378 int base[], double r[], char *file_out_name )
2379
2380 //******************************************************************************
2381 //
2382 // Purpose:
2383 //
2384 // HALHAM_WRITE writes a Halton or Hammersley dataset to a file.
2385 //
2386 // License:
2387 //
2388 // Copyright (C) 2004 John Burkardt and Max Gunzburger
2389 //
2390 // This library is free software; you can redistribute it and/or
2391 // modify it under the terms of the GNU Lesser General Public
2392 // License as published by the Free Software Foundation; either
2393 // version 2.1 of the License, or (at your option) any later version.
2394 //
2395 // This library is distributed in the hope that it will be useful,
2396 // but WITHOUT ANY WARRANTY; without even the implied warranty of
2397 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2398 // Lesser General Public License for more details.
2399 //
2400 // You should have received a copy of the GNU Lesser General Public
2401 // License along with this library; if not, write to the Free Software
2402 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
2403 //
2404 // Discussion:
2405 //
2406 // The initial lines of the file are comments, which begin with a
2407 // "#" character.
2408 //
2409 // Thereafter, each line of the file contains the NDIM-dimensional
2410 // components of the next entry of the dataset.
2411 //
2412 // Modified:
2413 //
2414 // 17 September 2004
2415 //
2416 // Author:
2417 //
2418 // John Burkardt
2419 //
2420 // Parameters:
2421 //
2422 // Input, int NDIM, the spatial dimension.
2423 //
2424 // Input, int N, the number of elements in the subsequence.
2425 //
2426 // Input, int STEP, the index of the subsequence element.
2427 // 0 <= STEP is required.
2428 //
2429 // Input, int SEED[NDIM], the sequence index for STEP = 0.
2430 //
2431 // Input, int LEAP[NDIM], the successive jumps in the sequence.
2432 //
2433 // Input, int BASE[NDIM], the bases.
2434 //
2435 // Input, double R[NDIM*N], the points.
2436 //
2437 // Input, char *FILE_OUT_NAME, the name of the output file.
2438 //
2439 {
2440 ofstream file_out;
2441 int i;
2442 int j;
2443 int mhi;
2444 int mlo;
2445 char *s;
2446
2447 file_out.open ( file_out_name );
2448
2449 if ( !file_out )
2450 {
2451 cout << "\n";
2452 cout << "HALHAM_WRITE - Fatal error!\n";
2453 cout << " Could not open the output file.\n";
2454 exit ( 1 );
2455 }
2456
2457 s = timestring ( );
2458
2459 file_out << "# " << file_out_name << "\n";
2460 file_out << "# created by routine HALHAM_WRITE.CC" << "\n";
2461 file_out << "# at " << s << "\n";
2462 file_out << "#\n";
2463 file_out << "# NDIM = " << setw(12) << ndim << "\n";
2464 file_out << "# N = " << setw(12) << n << "\n";
2465 file_out << "# STEP = " << setw(12) << step << "\n";
2466 for ( mlo = 1; mlo <= ndim; mlo = mlo + 5 )
2467 {
2468 mhi = i_min ( mlo + 5 - 1, ndim );
2469 if ( mlo == 1 )
2470 {
2471 file_out << "# SEED = ";
2472 }
2473 else
2474 {
2475 file_out << "# ";
2476 }
2477 for ( i = mlo; i <= mhi; i++ )
2478 {
2479 file_out << setw(12) << seed[i-1];
2480 }
2481 file_out << "\n";
2482 }
2483 for ( mlo = 1; mlo <= ndim; mlo = mlo + 5 )
2484 {
2485 mhi = i_min ( mlo + 5 - 1, ndim );
2486 if ( mlo == 1 )
2487 {
2488 file_out << "# LEAP = ";
2489 }
2490 else
2491 {
2492 file_out << "# ";
2493 }
2494 for ( i = mlo; i <= mhi; i++ )
2495 {
2496 file_out << setw(12) << leap[i-1];
2497 }
2498 file_out << "\n";
2499 }
2500 for ( mlo = 1; mlo <= ndim; mlo = mlo + 5 )
2501 {
2502 mhi = i_min ( mlo + 5 - 1, ndim );
2503 if ( mlo == 1 )
2504 {
2505 file_out << "# BASE = ";
2506 }
2507 else
2508 {
2509 file_out << "# ";
2510 }
2511 for ( i = mlo; i <= mhi; i++ )
2512 {
2513 file_out << setw(12) << base[i-1];
2514 }
2515 file_out << "\n";
2516 }
2517 file_out << "# EPSILON (unit roundoff) = " << d_epsilon ( ) << "\n";
2518 file_out << "#\n";
2519
2520 for ( j = 0; j < n; j++ )
2521 {
2522 for ( i = 0; i < ndim; i++ )
2523 {
2524 file_out << setw(10) << r[i+j*ndim] << " ";
2525 }
2526 file_out << "\n";
2527 }
2528
2529 file_out.close ( );
2530
2531 return;
2532 }
2533 //******************************************************************************
2534
i_log_10(int i)2535 int i_log_10 ( int i )
2536
2537 //******************************************************************************
2538 //
2539 // Purpose:
2540 //
2541 // I_LOG_10 returns the whole part of the logarithm base 10 of an integer.
2542 //
2543 // License:
2544 //
2545 // Copyright (C) 2004 John Burkardt and Max Gunzburger
2546 //
2547 // This library is free software; you can redistribute it and/or
2548 // modify it under the terms of the GNU Lesser General Public
2549 // License as published by the Free Software Foundation; either
2550 // version 2.1 of the License, or (at your option) any later version.
2551 //
2552 // This library is distributed in the hope that it will be useful,
2553 // but WITHOUT ANY WARRANTY; without even the implied warranty of
2554 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2555 // Lesser General Public License for more details.
2556 //
2557 // You should have received a copy of the GNU Lesser General Public
2558 // License along with this library; if not, write to the Free Software
2559 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
2560 //
2561 // Discussion:
2562 //
2563 // It should be the case that 10^I_LOG_10(I) <= |I| < 10^(I_LOG_10(I)+1).
2564 // (except for I = 0).
2565 //
2566 // The number of decimal digits in I is I_LOG_10(I) + 1.
2567 //
2568 // Examples:
2569 //
2570 // I I_LOG_10(I)
2571 //
2572 // 0 0
2573 // 1 0
2574 // 2 0
2575 //
2576 // 9 0
2577 // 10 1
2578 // 11 1
2579 //
2580 // 99 1
2581 // 100 2
2582 // 101 2
2583 //
2584 // 999 2
2585 // 1000 3
2586 // 1001 3
2587 //
2588 // 9999 3
2589 // 10000 4
2590 // 10001 4
2591 //
2592 // Modified:
2593 //
2594 // 17 September 2004
2595 //
2596 // Author:
2597 //
2598 // John Burkardt
2599 //
2600 // Parameters:
2601 //
2602 // Input, int I, the integer.
2603 //
2604 // Output, int I_LOG_10, the whole part of the logarithm of abs ( I ).
2605 //
2606 {
2607 int ten_pow;
2608 int value;
2609
2610 i = abs ( i );
2611
2612 ten_pow = 10;
2613 value = 0;
2614
2615 while ( ten_pow <= i )
2616 {
2617 ten_pow = ten_pow * 10;
2618 value = value + 1;
2619 }
2620
2621 return value;
2622 }
2623 //****************************************************************************
2624
i_max(int i1,int i2)2625 int i_max ( int i1, int i2 )
2626
2627 //****************************************************************************
2628 //
2629 // Purpose:
2630 //
2631 // I_MAX returns the maximum of two integers.
2632 //
2633 // License:
2634 //
2635 // Copyright (C) 2004 John Burkardt and Max Gunzburger
2636 //
2637 // This library is free software; you can redistribute it and/or
2638 // modify it under the terms of the GNU Lesser General Public
2639 // License as published by the Free Software Foundation; either
2640 // version 2.1 of the License, or (at your option) any later version.
2641 //
2642 // This library is distributed in the hope that it will be useful,
2643 // but WITHOUT ANY WARRANTY; without even the implied warranty of
2644 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2645 // Lesser General Public License for more details.
2646 //
2647 // You should have received a copy of the GNU Lesser General Public
2648 // License along with this library; if not, write to the Free Software
2649 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
2650 //
2651 // Modified:
2652 //
2653 // 13 October 1998
2654 //
2655 // Author:
2656 //
2657 // John Burkardt
2658 //
2659 // Parameters:
2660 //
2661 // Input, int I1, I2, are two integers to be compared.
2662 //
2663 // Output, int I_MAX, the larger of I1 and I2.
2664 //
2665 //
2666 {
2667 if ( i2 < i1 )
2668 {
2669 return i1;
2670 }
2671 else
2672 {
2673 return i2;
2674 }
2675
2676 }
2677 //****************************************************************************
2678
i_min(int i1,int i2)2679 int i_min ( int i1, int i2 )
2680
2681 //****************************************************************************
2682 //
2683 // Purpose:
2684 //
2685 // I_MIN returns the smaller of two integers.
2686 //
2687 // License:
2688 //
2689 // Copyright (C) 2004 John Burkardt and Max Gunzburger
2690 //
2691 // This library is free software; you can redistribute it and/or
2692 // modify it under the terms of the GNU Lesser General Public
2693 // License as published by the Free Software Foundation; either
2694 // version 2.1 of the License, or (at your option) any later version.
2695 //
2696 // This library is distributed in the hope that it will be useful,
2697 // but WITHOUT ANY WARRANTY; without even the implied warranty of
2698 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2699 // Lesser General Public License for more details.
2700 //
2701 // You should have received a copy of the GNU Lesser General Public
2702 // License along with this library; if not, write to the Free Software
2703 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
2704 //
2705 // Modified:
2706 //
2707 // 17 September 2004
2708 //
2709 // Author:
2710 //
2711 // John Burkardt
2712 //
2713 // Parameters:
2714 //
2715 // Input, int I1, I2, two integers to be compared.
2716 //
2717 // Output, int I_MIN, the smaller of I1 and I2.
2718 //
2719 //
2720 {
2721 if ( i1 < i2 )
2722 {
2723 return i1;
2724 }
2725 else
2726 {
2727 return i2;
2728 }
2729
2730 }
2731 //******************************************************************************
2732
i_to_s(int i)2733 char *i_to_s ( int i )
2734
2735 //******************************************************************************
2736 //
2737 // Purpose:
2738 //
2739 // I_TO_S converts an integer to a string.
2740 //
2741 // License:
2742 //
2743 // Copyright (C) 2004 John Burkardt and Max Gunzburger
2744 //
2745 // This library is free software; you can redistribute it and/or
2746 // modify it under the terms of the GNU Lesser General Public
2747 // License as published by the Free Software Foundation; either
2748 // version 2.1 of the License, or (at your option) any later version.
2749 //
2750 // This library is distributed in the hope that it will be useful,
2751 // but WITHOUT ANY WARRANTY; without even the implied warranty of
2752 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2753 // Lesser General Public License for more details.
2754 //
2755 // You should have received a copy of the GNU Lesser General Public
2756 // License along with this library; if not, write to the Free Software
2757 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
2758 //
2759 // Examples:
2760 //
2761 // INTVAL S
2762 //
2763 // 1 1
2764 // -1 -1
2765 // 0 0
2766 // 1952 1952
2767 // 123456 123456
2768 // 1234567 1234567
2769 //
2770 // Modified:
2771 //
2772 // 17 September 2004
2773 //
2774 // Author:
2775 //
2776 // John Burkardt
2777 //
2778 // Parameters:
2779 //
2780 // Input, int I, an integer to be converted.
2781 //
2782 // Output, char *I_TO_S, the representation of the integer.
2783 //
2784 {
2785 int digit;
2786 int j;
2787 int length;
2788 int ten_power;
2789 char *s;
2790
2791 length = i_log_10 ( i );
2792
2793 ten_power = ( int ) pow ( ( double ) 10, ( double ) length );
2794
2795 if ( i < 0 )
2796 {
2797 length = length + 1;
2798 }
2799 //
2800 // Add one position for the trailing null.
2801 //
2802 length = length + 1;
2803
2804 s = new char[length];
2805
2806 if ( i == 0 )
2807 {
2808 s[0] = '0';
2809 s[1] = '\0';
2810 return s;
2811 }
2812 //
2813 // Now take care of the sign.
2814 //
2815 j = 0;
2816 if ( i < 0 )
2817 {
2818 s[j] = '-';
2819 j = j + 1;
2820 i = abs ( i );
2821 }
2822 //
2823 // Find the leading digit of I, strip it off, and stick it into the string.
2824 //
2825 while ( 0 < ten_power )
2826 {
2827 digit = i / ten_power;
2828 s[j] = digit_to_ch ( digit );
2829 j = j + 1;
2830 i = i - digit * ten_power;
2831 ten_power = ten_power / 10;
2832 }
2833 //
2834 // Tack on the trailing NULL.
2835 //
2836 s[j] = '\0';
2837 j = j + 1;
2838
2839 return s;
2840 }
2841 //******************************************************************************
2842
ivec_transpose_print(int n,int a[],const char * title)2843 void ivec_transpose_print ( int n, int a[], const char *title )
2844
2845 //******************************************************************************
2846 //
2847 // Purpose:
2848 //
2849 // IVEC_TRANSPOSE_PRINT prints an integer vector "transposed".
2850 //
2851 // License:
2852 //
2853 // Copyright (C) 2004 John Burkardt and Max Gunzburger
2854 //
2855 // This library is free software; you can redistribute it and/or
2856 // modify it under the terms of the GNU Lesser General Public
2857 // License as published by the Free Software Foundation; either
2858 // version 2.1 of the License, or (at your option) any later version.
2859 //
2860 // This library is distributed in the hope that it will be useful,
2861 // but WITHOUT ANY WARRANTY; without even the implied warranty of
2862 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2863 // Lesser General Public License for more details.
2864 //
2865 // You should have received a copy of the GNU Lesser General Public
2866 // License along with this library; if not, write to the Free Software
2867 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
2868 //
2869 // Example:
2870 //
2871 // A = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 }
2872 // TITLE = "My vector: "
2873 //
2874 // My vector: 1 2 3 4 5
2875 // 6 7 8 9 10
2876 // 11
2877 //
2878 // Modified:
2879 //
2880 // 17 September 2004
2881 //
2882 // Author:
2883 //
2884 // John Burkardt
2885 //
2886 // Parameters:
2887 //
2888 // Input, int N, the number of components of the vector.
2889 //
2890 // Input, int A[N], the vector to be printed.
2891 //
2892 // Input, char *TITLE, a title to be printed first.
2893 // TITLE may be blank or NULL.
2894 //
2895 {
2896 int i;
2897 int ihi;
2898 int ilo;
2899 int title_len;
2900
2901 if ( 0 < s_len_trim ( title ) )
2902 {
2903 title_len = strlen ( title );
2904
2905 for ( ilo = 1; ilo <= n; ilo = ilo + 5 )
2906 {
2907 ihi = i_min ( ilo + 5 - 1, n );
2908 if ( ilo == 1 )
2909 {
2910 cout << title;
2911 }
2912 else
2913 {
2914 for ( i = 1; i <= title_len; i++ )
2915 {
2916 cout << " ";
2917 }
2918 }
2919 for ( i = ilo; i <= ihi; i++ )
2920 {
2921 cout << setw(12) << a[i-1];
2922 }
2923 cout << "\n";
2924 }
2925 }
2926 else
2927 {
2928 for ( ilo = 1; ilo <= n; ilo = ilo + 5 )
2929 {
2930 ihi = i_min ( ilo + 5 - 1, n );
2931 for ( i = ilo; i <= ihi; i++ )
2932 {
2933 cout << setw(12) << a[i-1];
2934 }
2935 cout << "\n";
2936 }
2937 }
2938
2939 return;
2940 }
2941 //******************************************************************************
2942
prime(int n)2943 int prime ( int n )
2944
2945 //******************************************************************************
2946 //
2947 // Purpose:
2948 //
2949 // PRIME returns any of the first PRIME_MAX prime numbers.
2950 //
2951 // License:
2952 //
2953 // Copyright (C) 2004 John Burkardt and Max Gunzburger
2954 //
2955 // This library is free software; you can redistribute it and/or
2956 // modify it under the terms of the GNU Lesser General Public
2957 // License as published by the Free Software Foundation; either
2958 // version 2.1 of the License, or (at your option) any later version.
2959 //
2960 // This library is distributed in the hope that it will be useful,
2961 // but WITHOUT ANY WARRANTY; without even the implied warranty of
2962 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2963 // Lesser General Public License for more details.
2964 //
2965 // You should have received a copy of the GNU Lesser General Public
2966 // License along with this library; if not, write to the Free Software
2967 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
2968 //
2969 // Note:
2970 //
2971 // PRIME_MAX is 1500, and the largest prime stored is 12553.
2972 //
2973 // Modified:
2974 //
2975 // 17 September 2004
2976 //
2977 // Author:
2978 //
2979 // John Burkardt
2980 //
2981 // Reference:
2982 //
2983 // Milton Abramowitz and Irene Stegun,
2984 // Handbook of Mathematical Functions,
2985 // US Department of Commerce, 1964, pages 870-873.
2986 //
2987 // Daniel Zwillinger,
2988 // CRC Standard Mathematical Tables and Formulae,
2989 // 30th Edition,
2990 // CRC Press, 1996, pages 95-98.
2991 //
2992 // Parameters:
2993 //
2994 // Input, int N, the index of the desired prime number.
2995 // N = -1 returns PRIME_MAX, the index of the largest prime available.
2996 // N = 0 is legal, returning PRIME = 1.
2997 // It should generally be true that 0 <= N <= PRIME_MAX.
2998 //
2999 // Output, int PRIME, the N-th prime. If N is out of range, PRIME
3000 // is returned as -1.
3001 //
3002 {
3003 #define PRIME_MAX 1500
3004
3005 int npvec[PRIME_MAX] = {
3006 2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
3007 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
3008 73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
3009 127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
3010 179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
3011 233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
3012 283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
3013 353, 359, 367, 373, 379, 383, 389, 397, 401, 409,
3014 419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
3015 467, 479, 487, 491, 499, 503, 509, 521, 523, 541,
3016 547, 557, 563, 569, 571, 577, 587, 593, 599, 601,
3017 607, 613, 617, 619, 631, 641, 643, 647, 653, 659,
3018 661, 673, 677, 683, 691, 701, 709, 719, 727, 733,
3019 739, 743, 751, 757, 761, 769, 773, 787, 797, 809,
3020 811, 821, 823, 827, 829, 839, 853, 857, 859, 863,
3021 877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
3022 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013,
3023 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069,
3024 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151,
3025 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223,
3026 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291,
3027 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373,
3028 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451,
3029 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511,
3030 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583,
3031 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657,
3032 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733,
3033 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811,
3034 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889,
3035 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987,
3036 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053,
3037 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129,
3038 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213,
3039 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287,
3040 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357,
3041 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423,
3042 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531,
3043 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617,
3044 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687,
3045 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741,
3046 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819,
3047 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903,
3048 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999,
3049 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079,
3050 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181,
3051 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257,
3052 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331,
3053 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413,
3054 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511,
3055 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571,
3056 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643,
3057 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727,
3058 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821,
3059 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907,
3060 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989,
3061 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057,
3062 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139,
3063 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231,
3064 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297,
3065 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409,
3066 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493,
3067 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583,
3068 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657,
3069 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751,
3070 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831,
3071 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937,
3072 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003,
3073 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087,
3074 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179,
3075 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279,
3076 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387,
3077 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443,
3078 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521,
3079 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639,
3080 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693,
3081 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791,
3082 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857,
3083 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939,
3084 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053,
3085 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133,
3086 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221,
3087 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301,
3088 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367,
3089 6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473,
3090 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571,
3091 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673,
3092 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761,
3093 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833,
3094 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917,
3095 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997,
3096 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103,
3097 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207,
3098 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297,
3099 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411,
3100 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499,
3101 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561,
3102 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643,
3103 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723,
3104 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829,
3105 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919,
3106 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017,
3107 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111,
3108 8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219,
3109 8221, 8231, 8233, 8237, 8243, 8263, 8269, 8273, 8287, 8291,
3110 8293, 8297, 8311, 8317, 8329, 8353, 8363, 8369, 8377, 8387,
3111 8389, 8419, 8423, 8429, 8431, 8443, 8447, 8461, 8467, 8501,
3112 8513, 8521, 8527, 8537, 8539, 8543, 8563, 8573, 8581, 8597,
3113 8599, 8609, 8623, 8627, 8629, 8641, 8647, 8663, 8669, 8677,
3114 8681, 8689, 8693, 8699, 8707, 8713, 8719, 8731, 8737, 8741,
3115 8747, 8753, 8761, 8779, 8783, 8803, 8807, 8819, 8821, 8831,
3116 8837, 8839, 8849, 8861, 8863, 8867, 8887, 8893, 8923, 8929,
3117 8933, 8941, 8951, 8963, 8969, 8971, 8999, 9001, 9007, 9011,
3118 9013, 9029, 9041, 9043, 9049, 9059, 9067, 9091, 9103, 9109,
3119 9127, 9133, 9137, 9151, 9157, 9161, 9173, 9181, 9187, 9199,
3120 9203, 9209, 9221, 9227, 9239, 9241, 9257, 9277, 9281, 9283,
3121 9293, 9311, 9319, 9323, 9337, 9341, 9343, 9349, 9371, 9377,
3122 9391, 9397, 9403, 9413, 9419, 9421, 9431, 9433, 9437, 9439,
3123 9461, 9463, 9467, 9473, 9479, 9491, 9497, 9511, 9521, 9533,
3124 9539, 9547, 9551, 9587, 9601, 9613, 9619, 9623, 9629, 9631,
3125 9643, 9649, 9661, 9677, 9679, 9689, 9697, 9719, 9721, 9733,
3126 9739, 9743, 9749, 9767, 9769, 9781, 9787, 9791, 9803, 9811,
3127 9817, 9829, 9833, 9839, 9851, 9857, 9859, 9871, 9883, 9887,
3128 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973,10007,
3129 10009,10037,10039,10061,10067,10069,10079,10091,10093,10099,
3130 10103,10111,10133,10139,10141,10151,10159,10163,10169,10177,
3131 10181,10193,10211,10223,10243,10247,10253,10259,10267,10271,
3132 10273,10289,10301,10303,10313,10321,10331,10333,10337,10343,
3133 10357,10369,10391,10399,10427,10429,10433,10453,10457,10459,
3134 10463,10477,10487,10499,10501,10513,10529,10531,10559,10567,
3135 10589,10597,10601,10607,10613,10627,10631,10639,10651,10657,
3136 10663,10667,10687,10691,10709,10711,10723,10729,10733,10739,
3137 10753,10771,10781,10789,10799,10831,10837,10847,10853,10859,
3138 10861,10867,10883,10889,10891,10903,10909,19037,10939,10949,
3139 10957,10973,10979,10987,10993,11003,11027,11047,11057,11059,
3140 11069,11071,11083,11087,11093,11113,11117,11119,11131,11149,
3141 11159,11161,11171,11173,11177,11197,11213,11239,11243,11251,
3142 11257,11261,11273,11279,11287,11299,11311,11317,11321,11329,
3143 11351,11353,11369,11383,11393,11399,11411,11423,11437,11443,
3144 11447,11467,11471,11483,11489,11491,11497,11503,11519,11527,
3145 11549,11551,11579,11587,11593,11597,11617,11621,11633,11657,
3146 11677,11681,11689,11699,11701,11717,11719,11731,11743,11777,
3147 11779,11783,11789,11801,11807,11813,11821,11827,11831,11833,
3148 11839,11863,11867,11887,11897,11903,11909,11923,11927,11933,
3149 11939,11941,11953,11959,11969,11971,11981,11987,12007,12011,
3150 12037,12041,12043,12049,12071,12073,12097,12101,12107,12109,
3151 12113,12119,12143,12149,12157,12161,12163,12197,12203,12211,
3152 12227,12239,12241,12251,12253,12263,12269,12277,12281,12289,
3153 12301,12323,12329,12343,12347,12373,12377,12379,12391,12401,
3154 12409,12413,12421,12433,12437,12451,12457,12473,12479,12487,
3155 12491,12497,12503,12511,12517,12527,12539,12541,12547,12553 };
3156
3157 if ( n == -1 )
3158 {
3159 return PRIME_MAX;
3160 }
3161 else if ( n == 0 )
3162 {
3163 return 1;
3164 }
3165 else if ( n < PRIME_MAX )
3166 {
3167 return npvec[n-1];
3168 }
3169 else
3170 {
3171 return -1;
3172 }
3173
3174 #undef PRIME_MAX
3175 }
3176 //******************************************************************************
3177
random_initialize(int seed)3178 unsigned long random_initialize ( int seed )
3179
3180 //******************************************************************************
3181 //
3182 // Purpose:
3183 //
3184 // RANDOM_INITIALIZE initializes the RANDOM random number generator.
3185 //
3186 // License:
3187 //
3188 // Copyright (C) 2004 John Burkardt and Max Gunzburger
3189 //
3190 // This library is free software; you can redistribute it and/or
3191 // modify it under the terms of the GNU Lesser General Public
3192 // License as published by the Free Software Foundation; either
3193 // version 2.1 of the License, or (at your option) any later version.
3194 //
3195 // This library is distributed in the hope that it will be useful,
3196 // but WITHOUT ANY WARRANTY; without even the implied warranty of
3197 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3198 // Lesser General Public License for more details.
3199 //
3200 // You should have received a copy of the GNU Lesser General Public
3201 // License along with this library; if not, write to the Free Software
3202 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
3203 //
3204 // Discussion:
3205 //
3206 // If you don't initialize RANDOM, the random number generator,
3207 // it will behave as though it were seeded with value 1.
3208 // This routine will either take a user-specified seed, or
3209 // (if the user passes a 0) make up a "random" one. In either
3210 // case, the seed is passed to SRAND (the appropriate routine
3211 // to call when setting the seed for RANDOM). The seed is also
3212 // returned to the user as the value of the function.
3213 //
3214 // Modified:
3215 //
3216 // 07 December 2004
3217 //
3218 // Author:
3219 //
3220 // John Burkardt
3221 //
3222 // Parameters:
3223 //
3224 // Input, int SEED, is either 0, which means that the user
3225 // wants this routine to come up with a seed, or nonzero, in which
3226 // case the user has supplied the seed.
3227 //
3228 // Output, unsigned long RANDOM_INITIALIZE, is the value of the seed
3229 // passed to SRAND, which is either the user's input value, or if
3230 // that was zero, the value selected by this routine.
3231 //
3232 {
3233 unsigned long ul_seed;
3234
3235 if ( seed == 0 )
3236 {
3237 seed = get_seed ( );
3238 }
3239 //
3240 // Now set the seed.
3241 //
3242 ul_seed = ( unsigned long ) seed;
3243
3244 srand ( ul_seed );
3245
3246 return ul_seed;
3247 }
3248 //******************************************************************************
3249
s_eqi(const char * s1,const char * s2)3250 bool s_eqi ( const char *s1, const char *s2 )
3251
3252 //******************************************************************************
3253 //
3254 // Purpose:
3255 //
3256 // S_EQI reports whether two strings are equal, ignoring case.
3257 //
3258 // License:
3259 //
3260 // Copyright (C) 2004 John Burkardt and Max Gunzburger
3261 //
3262 // This library is free software; you can redistribute it and/or
3263 // modify it under the terms of the GNU Lesser General Public
3264 // License as published by the Free Software Foundation; either
3265 // version 2.1 of the License, or (at your option) any later version.
3266 //
3267 // This library is distributed in the hope that it will be useful,
3268 // but WITHOUT ANY WARRANTY; without even the implied warranty of
3269 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3270 // Lesser General Public License for more details.
3271 //
3272 // You should have received a copy of the GNU Lesser General Public
3273 // License along with this library; if not, write to the Free Software
3274 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
3275 //
3276 // Modified:
3277 //
3278 // 05 May 2003
3279 //
3280 // Author:
3281 //
3282 // John Burkardt
3283 //
3284 // Parameters:
3285 //
3286 // Input, char *S1, char *S2, pointers to two strings.
3287 //
3288 // Output, bool S_EQI, is true if the strings are equal.
3289 //
3290 {
3291 int i;
3292 int nchar;
3293 int nchar1;
3294 int nchar2;
3295
3296 nchar1 = strlen ( s1 );
3297 nchar2 = strlen ( s2 );
3298 nchar = i_min ( nchar1, nchar2 );
3299
3300 //
3301 // The strings are not equal if they differ over their common length.
3302 //
3303 for ( i = 0; i < nchar; i++ )
3304 {
3305
3306 if ( ch_cap ( s1[i] ) != ch_cap ( s2[i] ) )
3307 {
3308 return false;
3309 }
3310 }
3311 //
3312 // The strings are not equal if the longer one includes nonblanks
3313 // in the tail.
3314 //
3315 if ( nchar < nchar1 )
3316 {
3317 for ( i = nchar; i < nchar1; i++ )
3318 {
3319 if ( s1[i] != ' ' )
3320 {
3321 return false;
3322 }
3323 }
3324 }
3325 else if ( nchar < nchar2 )
3326 {
3327 for ( i = nchar; i < nchar2; i++ )
3328 {
3329 if ( s2[i] != ' ' )
3330 {
3331 return false;
3332 }
3333 }
3334 }
3335
3336 return true;
3337
3338 }
3339 //******************************************************************************
3340
s_index_last_c(const char * s,char c)3341 int s_index_last_c ( const char *s, char c )
3342
3343 //******************************************************************************
3344 //
3345 // Purpose:
3346 //
3347 // S_INDEX_LAST_C returns a pointer to the last occurrence of a given character.
3348 //
3349 // License:
3350 //
3351 // Copyright (C) 2004 John Burkardt and Max Gunzburger
3352 //
3353 // This library is free software; you can redistribute it and/or
3354 // modify it under the terms of the GNU Lesser General Public
3355 // License as published by the Free Software Foundation; either
3356 // version 2.1 of the License, or (at your option) any later version.
3357 //
3358 // This library is distributed in the hope that it will be useful,
3359 // but WITHOUT ANY WARRANTY; without even the implied warranty of
3360 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3361 // Lesser General Public License for more details.
3362 //
3363 // You should have received a copy of the GNU Lesser General Public
3364 // License along with this library; if not, write to the Free Software
3365 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
3366 //
3367 // Modified:
3368 //
3369 // 10 October 2004
3370 //
3371 // Author:
3372 //
3373 // John Burkardt
3374 //
3375 // Parameters:
3376 //
3377 // Input, char *S, a pointer to a string.
3378 //
3379 // Input, char C, the character to search for.
3380 //
3381 // Output, int S_INDEX_LAST_C, the index in S of the last occurrence of the character,
3382 // or -1 if it does not occur.
3383 //
3384 {
3385 int n;
3386 const char *t;
3387
3388 n = strlen ( s ) - 1;
3389 t = s + strlen ( s ) - 1;
3390
3391 while ( 0 <= n )
3392 {
3393 if ( *t == c )
3394 {
3395 return n;
3396 }
3397 t--;
3398 n--;
3399 }
3400
3401 return (-1);
3402 }
3403 //******************************************************************************
3404
s_len_trim(const char * s)3405 int s_len_trim ( const char *s )
3406
3407 //******************************************************************************
3408 //
3409 // Purpose:
3410 //
3411 // S_LEN_TRIM returns the length of a string to the last nonblank.
3412 //
3413 // License:
3414 //
3415 // Copyright (C) 2004 John Burkardt and Max Gunzburger
3416 //
3417 // This library is free software; you can redistribute it and/or
3418 // modify it under the terms of the GNU Lesser General Public
3419 // License as published by the Free Software Foundation; either
3420 // version 2.1 of the License, or (at your option) any later version.
3421 //
3422 // This library is distributed in the hope that it will be useful,
3423 // but WITHOUT ANY WARRANTY; without even the implied warranty of
3424 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3425 // Lesser General Public License for more details.
3426 //
3427 // You should have received a copy of the GNU Lesser General Public
3428 // License along with this library; if not, write to the Free Software
3429 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
3430 //
3431 // Modified:
3432 //
3433 // 17 September 2004
3434 //
3435 // Author:
3436 //
3437 // John Burkardt
3438 //
3439 // Parameters:
3440 //
3441 // Input, char *S, a pointer to a string.
3442 //
3443 // Output, int S_LEN_TRIM, the length of the string to the last nonblank.
3444 // If S_LEN_TRIM is 0, then the string is entirely blank.
3445 //
3446 {
3447 int n;
3448 const char* t;
3449
3450 n = strlen ( s );
3451 t = s + strlen ( s ) - 1;
3452
3453 while ( 0 < n )
3454 {
3455 if ( *t != ' ' )
3456 {
3457 return n;
3458 }
3459 t--;
3460 n--;
3461 }
3462
3463 return n;
3464 }
3465 //******************************************************************************
3466
s_to_d(char * s,int * lchar,bool * error)3467 double s_to_d ( char *s, int *lchar, bool *error )
3468
3469 //******************************************************************************
3470 //
3471 // Purpose:
3472 //
3473 // S_TO_D reads a real number from a string.
3474 //
3475 // License:
3476 //
3477 // Copyright (C) 2004 John Burkardt and Max Gunzburger
3478 //
3479 // This library is free software; you can redistribute it and/or
3480 // modify it under the terms of the GNU Lesser General Public
3481 // License as published by the Free Software Foundation; either
3482 // version 2.1 of the License, or (at your option) any later version.
3483 //
3484 // This library is distributed in the hope that it will be useful,
3485 // but WITHOUT ANY WARRANTY; without even the implied warranty of
3486 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3487 // Lesser General Public License for more details.
3488 //
3489 // You should have received a copy of the GNU Lesser General Public
3490 // License along with this library; if not, write to the Free Software
3491 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
3492 //
3493 // Discussion:
3494 //
3495 // This routine will read as many characters as possible until it reaches
3496 // the end of the string, or encounters a character which cannot be
3497 // part of the real number.
3498 //
3499 // Legal input is:
3500 //
3501 // 1 blanks,
3502 // 2 '+' or '-' sign,
3503 // 2.5 spaces
3504 // 3 integer part,
3505 // 4 decimal point,
3506 // 5 fraction part,
3507 // 6 'E' or 'e' or 'D' or 'd', exponent marker,
3508 // 7 exponent sign,
3509 // 8 exponent integer part,
3510 // 9 exponent decimal point,
3511 // 10 exponent fraction part,
3512 // 11 blanks,
3513 // 12 final comma or semicolon.
3514 //
3515 // with most quantities optional.
3516 //
3517 // Examples:
3518 //
3519 // S R
3520 //
3521 // '1' 1.0
3522 // ' 1 ' 1.0
3523 // '1A' 1.0
3524 // '12,34,56' 12.0
3525 // ' 34 7' 34.0
3526 // '-1E2ABCD' -100.0
3527 // '-1X2ABCD' -1.0
3528 // ' 2E-1' 0.2
3529 // '23.45' 23.45
3530 // '-4.2E+2' -420.0
3531 // '17d2' 1700.0
3532 // '-14e-2' -0.14
3533 // 'e2' 100.0
3534 // '-12.73e-9.23' -12.73 * 10.0**(-9.23)
3535 //
3536 // Modified:
3537 //
3538 // 07 August 2003
3539 //
3540 // Author:
3541 //
3542 // John Burkardt
3543 //
3544 // Parameters:
3545 //
3546 // Input, char *S, the string containing the
3547 // data to be read. Reading will begin at position 1 and
3548 // terminate at the end of the string, or when no more
3549 // characters can be read to form a legal real. Blanks,
3550 // commas, or other nonnumeric data will, in particular,
3551 // cause the conversion to halt.
3552 //
3553 // Output, int *LCHAR, the number of characters read from
3554 // the string to form the number, including any terminating
3555 // characters such as a trailing comma or blanks.
3556 //
3557 // Output, bool *ERROR, is true if an error occurred.
3558 //
3559 // Output, double S_TO_D, the real value that was read from the string.
3560 //
3561 {
3562 char c;
3563 int ihave;
3564 int isgn;
3565 int iterm;
3566 int jbot;
3567 int jsgn;
3568 int jtop;
3569 int nchar;
3570 int ndig;
3571 double r;
3572 double rbot;
3573 double rexp;
3574 double rtop;
3575 char TAB = 9;
3576
3577 nchar = s_len_trim ( s );
3578 *error = false;
3579 r = 0.0;
3580 *lchar = -1;
3581 isgn = 1;
3582 rtop = 0.0;
3583 rbot = 1.0;
3584 jsgn = 1;
3585 jtop = 0;
3586 jbot = 1;
3587 ihave = 1;
3588 iterm = 0;
3589
3590 for ( ; ; )
3591 {
3592 c = s[*lchar+1];
3593 *lchar = *lchar + 1;
3594 //
3595 // Blank or TAB character.
3596 //
3597 if ( c == ' ' || c == TAB )
3598 {
3599 if ( ihave == 2 )
3600 {
3601 }
3602 else if ( ihave == 6 || ihave == 7 )
3603 {
3604 iterm = 1;
3605 }
3606 else if ( 1 < ihave )
3607 {
3608 ihave = 11;
3609 }
3610 }
3611 //
3612 // Comma.
3613 //
3614 else if ( c == ',' || c == ';' )
3615 {
3616 if ( ihave != 1 )
3617 {
3618 iterm = 1;
3619 ihave = 12;
3620 *lchar = *lchar + 1;
3621 }
3622 }
3623 //
3624 // Minus sign.
3625 //
3626 else if ( c == '-' )
3627 {
3628 if ( ihave == 1 )
3629 {
3630 ihave = 2;
3631 isgn = -1;
3632 }
3633 else if ( ihave == 6 )
3634 {
3635 ihave = 7;
3636 jsgn = -1;
3637 }
3638 else
3639 {
3640 iterm = 1;
3641 }
3642 }
3643 //
3644 // Plus sign.
3645 //
3646 else if ( c == '+' )
3647 {
3648 if ( ihave == 1 )
3649 {
3650 ihave = 2;
3651 }
3652 else if ( ihave == 6 )
3653 {
3654 ihave = 7;
3655 }
3656 else
3657 {
3658 iterm = 1;
3659 }
3660 }
3661 //
3662 // Decimal point.
3663 //
3664 else if ( c == '.' )
3665 {
3666 if ( ihave < 4 )
3667 {
3668 ihave = 4;
3669 }
3670 else if ( 6 <= ihave && ihave <= 8 )
3671 {
3672 ihave = 9;
3673 }
3674 else
3675 {
3676 iterm = 1;
3677 }
3678 }
3679 //
3680 // Exponent marker.
3681 //
3682 else if ( ch_eqi ( c, 'E' ) || ch_eqi ( c, 'D' ) )
3683 {
3684 if ( ihave < 6 )
3685 {
3686 ihave = 6;
3687 }
3688 else
3689 {
3690 iterm = 1;
3691 }
3692 }
3693 //
3694 // Digit.
3695 //
3696 else if ( ihave < 11 && '0' <= c && c <= '9' )
3697 {
3698 if ( ihave <= 2 )
3699 {
3700 ihave = 3;
3701 }
3702 else if ( ihave == 4 )
3703 {
3704 ihave = 5;
3705 }
3706 else if ( ihave == 6 || ihave == 7 )
3707 {
3708 ihave = 8;
3709 }
3710 else if ( ihave == 9 )
3711 {
3712 ihave = 10;
3713 }
3714
3715 ndig = ch_to_digit ( c );
3716
3717 if ( ihave == 3 )
3718 {
3719 rtop = 10.0E+00 * rtop + ( double ) ndig;
3720 }
3721 else if ( ihave == 5 )
3722 {
3723 rtop = 10.0E+00 * rtop + ( double ) ndig;
3724 rbot = 10.0E+00 * rbot;
3725 }
3726 else if ( ihave == 8 )
3727 {
3728 jtop = 10 * jtop + ndig;
3729 }
3730 else if ( ihave == 10 )
3731 {
3732 jtop = 10 * jtop + ndig;
3733 jbot = 10 * jbot;
3734 }
3735
3736 }
3737 //
3738 // Anything else is regarded as a terminator.
3739 //
3740 else
3741 {
3742 iterm = 1;
3743 }
3744 //
3745 // If we haven't seen a terminator, and we haven't examined the
3746 // entire string, go get the next character.
3747 //
3748 if ( iterm == 1 || nchar <= *lchar + 1 )
3749 {
3750 break;
3751 }
3752
3753 }
3754 //
3755 // If we haven't seen a terminator, and we have examined the
3756 // entire string, then we're done, and LCHAR is equal to NCHAR.
3757 //
3758 if ( iterm != 1 && (*lchar) + 1 == nchar )
3759 {
3760 *lchar = nchar;
3761 }
3762 //
3763 // Number seems to have terminated. Have we got a legal number?
3764 // Not if we terminated in states 1, 2, 6 or 7!
3765 //
3766 if ( ihave == 1 || ihave == 2 || ihave == 6 || ihave == 7 )
3767 {
3768 *error = true;
3769 return r;
3770 }
3771 //
3772 // Number seems OK. Form it.
3773 //
3774 if ( jtop == 0 )
3775 {
3776 rexp = 1.0E+00;
3777 }
3778 else
3779 {
3780 if ( jbot == 1 )
3781 {
3782 rexp = pow ( 10.0E+00, jsgn * jtop );
3783 }
3784 else
3785 {
3786 rexp = jsgn * jtop;
3787 rexp = rexp / jbot;
3788 rexp = pow ( 10.0E+00, rexp );
3789 }
3790
3791 }
3792
3793 r = isgn * rexp * rtop / rbot;
3794
3795 return r;
3796 }
3797 //******************************************************************************
3798
s_to_dvec(char * s,int n,double rvec[])3799 bool s_to_dvec ( char *s, int n, double rvec[] )
3800
3801 //******************************************************************************
3802 //
3803 // Purpose:
3804 //
3805 // S_TO_DVEC reads a real vector from a string.
3806 //
3807 // License:
3808 //
3809 // Copyright (C) 2004 John Burkardt and Max Gunzburger
3810 //
3811 // This library is free software; you can redistribute it and/or
3812 // modify it under the terms of the GNU Lesser General Public
3813 // License as published by the Free Software Foundation; either
3814 // version 2.1 of the License, or (at your option) any later version.
3815 //
3816 // This library is distributed in the hope that it will be useful,
3817 // but WITHOUT ANY WARRANTY; without even the implied warranty of
3818 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3819 // Lesser General Public License for more details.
3820 //
3821 // You should have received a copy of the GNU Lesser General Public
3822 // License along with this library; if not, write to the Free Software
3823 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
3824 //
3825 // Modified:
3826 //
3827 // 19 February 2001
3828 //
3829 // Author:
3830 //
3831 // John Burkardt
3832 //
3833 // Parameters:
3834 //
3835 // Input, char *S, the string to be read.
3836 //
3837 // Input, int N, the number of values expected.
3838 //
3839 // Output, double RVEC[N], the values read from the string.
3840 //
3841 // Output, bool S_TO_DVEC, is true if an error occurred.
3842 //
3843 {
3844 bool error;
3845 int i;
3846 int lchar;
3847 double x;
3848
3849 for ( i = 0; i < n; i++ )
3850 {
3851 rvec[i] = s_to_d ( s, &lchar, &error );
3852
3853 if ( error )
3854 {
3855 return error;
3856 }
3857
3858 s = s + lchar;
3859
3860 }
3861
3862 return error;
3863 }
3864 //**********************************************************************
3865
timestamp(void)3866 void timestamp ( void )
3867
3868 //**********************************************************************
3869 //
3870 // Purpose:
3871 //
3872 // TIMESTAMP prints the current YMDHMS date as a time stamp.
3873 //
3874 // License:
3875 //
3876 // Copyright (C) 2004 John Burkardt and Max Gunzburger
3877 //
3878 // This library is free software; you can redistribute it and/or
3879 // modify it under the terms of the GNU Lesser General Public
3880 // License as published by the Free Software Foundation; either
3881 // version 2.1 of the License, or (at your option) any later version.
3882 //
3883 // This library is distributed in the hope that it will be useful,
3884 // but WITHOUT ANY WARRANTY; without even the implied warranty of
3885 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3886 // Lesser General Public License for more details.
3887 //
3888 // You should have received a copy of the GNU Lesser General Public
3889 // License along with this library; if not, write to the Free Software
3890 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
3891 //
3892 // Example:
3893 //
3894 // May 31 2001 09:45:54 AM
3895 //
3896 // Modified:
3897 //
3898 // 17 September 2004
3899 //
3900 // Author:
3901 //
3902 // John Burkardt
3903 //
3904 // Parameters:
3905 //
3906 // None
3907 //
3908 {
3909 #define TIME_SIZE 40
3910
3911 static char time_buffer[TIME_SIZE];
3912 const struct tm *tm;
3913 size_t len;
3914 time_t now;
3915
3916 now = time ( NULL );
3917 tm = localtime ( &now );
3918
3919 len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
3920
3921 cout << time_buffer << "\n";
3922
3923 return;
3924 #undef TIME_SIZE
3925 }
3926 //**********************************************************************
3927
timestring(void)3928 char *timestring ( void )
3929
3930 //**********************************************************************
3931 //
3932 // Purpose:
3933 //
3934 // TIMESTRING returns the current YMDHMS date as a string.
3935 //
3936 // License:
3937 //
3938 // Copyright (C) 2004 John Burkardt and Max Gunzburger
3939 //
3940 // This library is free software; you can redistribute it and/or
3941 // modify it under the terms of the GNU Lesser General Public
3942 // License as published by the Free Software Foundation; either
3943 // version 2.1 of the License, or (at your option) any later version.
3944 //
3945 // This library is distributed in the hope that it will be useful,
3946 // but WITHOUT ANY WARRANTY; without even the implied warranty of
3947 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3948 // Lesser General Public License for more details.
3949 //
3950 // You should have received a copy of the GNU Lesser General Public
3951 // License along with this library; if not, write to the Free Software
3952 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
3953 //
3954 // Example:
3955 //
3956 // May 31 2001 09:45:54 AM
3957 //
3958 // Modified:
3959 //
3960 // 17 September 2004
3961 //
3962 // Author:
3963 //
3964 // John Burkardt
3965 //
3966 // Parameters:
3967 //
3968 // Output, char *TIMESTRING, a string containing the current YMDHMS date.
3969 //
3970 {
3971 #define TIME_SIZE 40
3972
3973 const struct tm *tm;
3974 size_t len;
3975 time_t now;
3976 char *s;
3977
3978 now = time ( NULL );
3979 tm = localtime ( &now );
3980
3981 s = new char[TIME_SIZE];
3982
3983 len = strftime ( s, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
3984
3985 return s;
3986 #undef TIME_SIZE
3987 }
3988 //******************************************************************************
3989
tuple_next_fast(int m,int n,int rank,int x[])3990 void tuple_next_fast ( int m, int n, int rank, int x[] )
3991
3992 //******************************************************************************
3993 //
3994 // Purpose:
3995 //
3996 // TUPLE_NEXT_FAST computes the next element of a tuple space, "fast".
3997 //
3998 // License:
3999 //
4000 // Copyright (C) 2004 John Burkardt and Max Gunzburger
4001 //
4002 // This library is free software; you can redistribute it and/or
4003 // modify it under the terms of the GNU Lesser General Public
4004 // License as published by the Free Software Foundation; either
4005 // version 2.1 of the License, or (at your option) any later version.
4006 //
4007 // This library is distributed in the hope that it will be useful,
4008 // but WITHOUT ANY WARRANTY; without even the implied warranty of
4009 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4010 // Lesser General Public License for more details.
4011 //
4012 // You should have received a copy of the GNU Lesser General Public
4013 // License along with this library; if not, write to the Free Software
4014 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
4015 //
4016 // Discussion:
4017 //
4018 // The elements are N vectors. Each entry is constrained to lie
4019 // between 1 and M. The elements are produced one at a time.
4020 // The first element is
4021 // (1,1,...,1)
4022 // and the last element is
4023 // (M,M,...,M)
4024 // Intermediate elements are produced in lexicographic order.
4025 //
4026 // Example:
4027 //
4028 // N = 2,
4029 // M = 3
4030 //
4031 // INPUT OUTPUT
4032 // ------- -------
4033 // Rank X
4034 // ---- ----
4035 // -1 -1 -1
4036 //
4037 // 0 1 1
4038 // 1 1 2
4039 // 2 1 3
4040 // 3 2 1
4041 // 4 2 2
4042 // 5 2 3
4043 // 6 3 1
4044 // 7 3 2
4045 // 8 3 3
4046 // 9 1 1
4047 //
4048 // Modified:
4049 //
4050 // 11 August 2004
4051 //
4052 // Author:
4053 //
4054 // John Burkardt
4055 //
4056 // Parameters:
4057 //
4058 // Input, int M, the maximum entry in each component.
4059 // M must be greater than 0.
4060 //
4061 // Input, int N, the number of components.
4062 // N must be greater than 0.
4063 //
4064 // Input, integer RANK, indicates the rank of the tuples.
4065 // Typically, 0 <= RANK < N**M; values larger than this are legal
4066 // and meaningful, and are equivalent to the corresponding value
4067 // MOD N**M. If RANK < 0, this indicates that this is the first
4068 // call for the given values of (M,N). Initialization is done,
4069 // and X is set to a dummy value.
4070 //
4071 // Output, int X[N], the next tuple, or a dummy value if initialization
4072 // is being done.
4073 //
4074 {
4075 static int *base = NULL;
4076 int i;
4077 //
4078 if ( rank < 0 )
4079 {
4080 if ( m <= 0 )
4081 {
4082 cout << "\n";
4083 cout << "TUPLE_NEXT_FAST - Fatal error!\n";
4084 cout << " The value M <= 0 is not legal.\n";
4085 cout << " M = " << m << "\n";
4086 exit ( 1 );
4087 }
4088 if ( n <= 0 )
4089 {
4090 cout << "\n";
4091 cout << "TUPLE_NEXT_FAST - Fatal error!\n";
4092 cout << " The value N <= 0 is not legal.\n";
4093 cout << " N = " << n << "\n";
4094 exit ( 1 );
4095 }
4096
4097 if ( base )
4098 {
4099 delete [] base;
4100 }
4101 base = new int[n];
4102
4103 base[n-1] = 1;
4104 for ( i = n-2; 0 <= i; i-- )
4105 {
4106 base[i] = base[i+1] * m;
4107 }
4108 for ( i = 0; i < n; i++ )
4109 {
4110 x[i] = -1;
4111 }
4112 }
4113 else
4114 {
4115 for ( i = 0; i < n; i++ )
4116 {
4117 x[i] = ( ( rank / base[i] ) % m ) + 1;
4118 }
4119 }
4120 return;
4121 }
4122
4123 //******************************************************************************
4124
word_count(char * s)4125 int word_count ( char *s )
4126
4127 //******************************************************************************
4128 //
4129 // Purpose:
4130 //
4131 // WORD_COUNT counts the number of "words" in a string.
4132 //
4133 // License:
4134 //
4135 // Copyright (C) 2004 John Burkardt and Max Gunzburger
4136 //
4137 // This library is free software; you can redistribute it and/or
4138 // modify it under the terms of the GNU Lesser General Public
4139 // License as published by the Free Software Foundation; either
4140 // version 2.1 of the License, or (at your option) any later version.
4141 //
4142 // This library is distributed in the hope that it will be useful,
4143 // but WITHOUT ANY WARRANTY; without even the implied warranty of
4144 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4145 // Lesser General Public License for more details.
4146 //
4147 // You should have received a copy of the GNU Lesser General Public
4148 // License along with this library; if not, write to the Free Software
4149 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
4150 //
4151 // Modified:
4152 //
4153 // 13 June 2003
4154 //
4155 // Author:
4156 //
4157 // John Burkardt
4158 //
4159 // Parameters:
4160 //
4161 // Input, char *S, the string to be examined.
4162 //
4163 // Output, int WORD_COUNT, the number of "words" in the string.
4164 // Words are presumed to be separated by one or more blanks.
4165 //
4166 {
4167 bool blank;
4168 int i;
4169 int nword;
4170
4171 nword = 0;
4172 blank = true;
4173
4174 /* 20140120, BMA: changed to avoid compiler warning:
4175 while ( *s != NULL )
4176 */
4177 while ( *s != '\0' )
4178 {
4179 if ( *s == ' ' )
4180 {
4181 blank = true;
4182 }
4183 else if ( blank )
4184 {
4185 nword = nword + 1;
4186 blank = false;
4187 }
4188 /* 20140120, BMA: right to left precedence for ++ and * --> * not needed
4189 but can't overrun string, so leaving on least-change principle */
4190 *s++;
4191 }
4192
4193 return nword;
4194 }
4195