fork(1) download
copy 
  1. // Generate nth Hamming number, for n < 1.5e12
  2. // Joe Knapp  knapp.167@osu.edu
  3. //
  4. // To generate "band" of triples in the vicinity of the nth number,
  5. // uses basic scheme of Louis Klauder described in:
  6. // http://w...content-available-to-author-only...s.com/architecture-and-design/hamming-problem/228700538
  7. // 
  8. // Once the lower triple in the band is identified, and the offset
  9. // from that number to the nth number determined, the table in
  10. // hamtab3.h is used to generate iterates up to the nth number.
  11. //
  12. #define TRUE 1
  13. #define FALSE 0
  14.  
  15. #define LOG2 log(2.0)
  16. #define LOG3 log(3.0)
  17. #define LOG5 log(5.0)
  18. #define THIRD (1.0/3.0)
  19.  
  20. #define LOWOFFSET 1.703
  21. #define HIOFFSET 1.693
  22.  
  23. #define NMIN 50000   // below NMIN, just get the Nth number by chugging
  24.                      // through the whole sequence
  25.  
  26. #include <unistd.h>
  27. #include <stdlib.h>
  28. #include <stdio.h>
  29. #include <math.h>
  30.  
  31. // table of intervals to compute the next Hamming number (2^p2*3^p3*5^p5)
  32. // The interval to use is the first one in the table twhere multiplication
  33. // results in positive powers for all three prime factors 2, 3, and 5.
  34. // 
  35. // These factors cover the first 1.59 trillion (1.59e12) Hamming numbers,
  36. // yielding values up to about 10^10000.
  37. #define NTAB 119  // number of table entries
  38. struct {
  39.     short p2 ; // power of 2
  40.     short p3 ; // power of 3
  41.     short p5 ; // power of 5
  42. } hamtab[NTAB] = {
  43. {16875,7633,-12478}, // 1.000000000053896
  44. {-22109,20086,-4189}, // 1.000000001373946
  45. {14639,-17295,5501}, // 1.000000002068432
  46. {31514,-9662,-6977}, // 1.000000002122329
  47. {-24345,-4842,13790}, // 1.000000003388482
  48. {-7470,2791,1312}, // 1.000000003442378
  49. {9405,10424,-11166}, // 1.000000003496275
  50. {7169,-14504,6813}, // 1.000000005510811
  51. {24044,-6871,-5665}, // 1.000000005564707
  52. {1935,13215,-9854}, // 1.000000006938654
  53. {-301,-11713,8125}, // 1.000000008953190
  54. {16574,-4080,-4353}, // 1.000000009007086
  55. {-5535,16006,-8542}, // 1.000000010381033
  56. {9104,-1289,-3041}, // 1.000000012449465
  57. {1634,1502,-1729}, // 1.000000015891844
  58. {-5836,4293,-417}, // 1.000000019334223
  59. {8803,-13002,5084}, // 1.000000021402655
  60. {1333,-10211,6396}, // 1.000000024845034
  61. {-6137,-7420,7708}, // 1.000000028287413
  62. {2967,-8709,4667}, // 1.000000040736879
  63. {-4503,-5918,5979}, // 1.000000044179258
  64. {4601,-7207,2938}, // 1.000000056628724
  65. {-2869,-4416,4250}, // 1.000000060071103
  66. {6235,-5705,1209}, // 1.000000072520570
  67. {-1235,-2914,2521}, // 1.000000075962949
  68. {7869,-4203,-520}, // 1.000000088412415
  69. {399,-1412,792}, // 1.000000091854795
  70. {-7071,1379,2104}, // 1.000000095297174
  71. {2033,90,-937}, // 1.000000107746641
  72. {-5437,2881,375}, // 1.000000111189020
  73. {2432,-1322,-145}, // 1.000000199601446
  74. {-5038,1469,1167}, // 1.000000203043825
  75. {-3404,2971,-562}, // 1.000000218935673
  76. {-4639,57,1959}, // 1.000000294898639
  77. {-3005,1559,230}, // 1.000000310790488
  78. {-2606,147,1022}, // 1.000000402645312
  79. {-972,1649,-707}, // 1.000000418537163
  80. {-2207,-1265,1814}, // 1.000000494500144
  81. {-573,237,85}, // 1.000000510391996
  82. {-174,-1175,877}, // 1.000000602246838
  83. {1460,327,-852}, // 1.000000618138692
  84. {1859,-1085,-60}, // 1.000000709993544
  85. {887,564,-767}, // 1.000001128531004
  86. {1286,-848,25}, // 1.000001220385903
  87. {314,801,-682}, // 1.000001638923577
  88. {713,-611,110}, // 1.000001730778523
  89. {-259,1038,-597}, // 1.000002149316410
  90. {140,-374,195}, // 1.000002241171403
  91. {-433,-137,280}, // 1.000002751564543
  92. {1600,-47,-657}, // 1.000002859311481
  93. {1027,190,-572}, // 1.000003369704937
  94. {454,427,-487}, // 1.000003880098653
  95. {-119,664,-402}, // 1.000004390492630
  96. {1167,-184,-377}, // 1.000005610883892
  97. {594,53,-292}, // 1.000006121278753
  98. {21,290,-207}, // 1.000006631673873
  99. {-552,527,-122}, // 1.000007142069255
  100. {734,-321,-97}, // 1.000008362463875
  101. {161,-84,-12}, // 1.000008872860139
  102. {-412,153,73}, // 1.000009383256665
  103. {-272,-221,268}, // 1.000011624449097
  104. {-391,443,-134}, // 1.000016014992765
  105. {-251,69,61}, // 1.000018256200061
  106. {-111,-305,256}, // 1.000020497412379
  107. {-230,359,-146}, // 1.000024887995004
  108. {-90,-15,49}, // 1.000027129222185
  109. {-69,275,-158}, // 1.000033761075971
  110. {71,-99,37}, // 1.000036002323039
  111. {92,191,-170}, // 1.000042634235669
  112. {-159,260,-109}, // 1.000060891214069
  113. {2,176,-121}, // 1.000069764614488
  114. {-249,245,-60}, // 1.000088022088186
  115. {-88,161,-72}, // 1.000096895729334
  116. {73,77,-84}, // 1.000105769449216
  117. {-178,146,-23}, // 1.000124027580225
  118. {-17,62,-35}, // 1.000132901540844
  119. {144,-22,-47}, // 1.000141775580201
  120. {-107,47,14}, // 1.000160034368546
  121. {54,-37,2}, // 1.000168908648648
  122. {127,40,-82}, // 1.000274695963239
  123. {37,25,-33}, // 1.000301832637713
  124. {-53,10,16}, // 1.000328970048383
  125. {91,-12,-31}, // 1.000470792268504
  126. {1,-27,18}, // 1.000497934262918
  127. {-16,35,-17}, // 1.000630901979994
  128. {38,-2,-15}, // 1.000799917193443
  129. {-52,-17,34}, // 1.000827068116760
  130. {-15,8,1}, // 1.001129150390625
  131. {22,33,-32}, // 1.001431323842778
  132. {-14,-19,19}, // 1.001627646896210
  133. {23,6,-14}, // 1.001929970810880
  134. {24,-21,4}, // 1.002428866072391
  135. {8,14,-13}, // 1.003061300428800
  136. {9,-13,5}, // 1.003560759018091
  137. {-7,22,-12}, // 1.004193907488000
  138. {-6,-5,6}, // 1.004693930041152
  139. {32,-7,-9}, // 1.005497601989940
  140. {17,1,-8}, // 1.006632960000000
  141. {2,9,-7}, // 1.007769600000000
  142. {-13,17,-6}, // 1.008907523437500
  143. {26,-12,-3}, // 1.010217337390227
  144. {11,-4,-2}, // 1.011358024691358
  145. {-4,4,-1}, // 1.012500000000000
  146. {5,-9,4}, // 1.016105268505817
  147. {-10,-1,5}, // 1.017252604166666
  148. {7,0,-3}, // 1.024000000000000
  149. {1,-5,3}, // 1.028806584362139
  150. {-14,3,4}, // 1.029968261718750
  151. {-3,-1,2}, // 1.041666666666666
  152. {-7,3,1}, // 1.054687500000000
  153. {4,-1,-1}, // 1.066666666666666
  154. {0,3,-2}, // 1.080000000000000
  155. {1,-2,1}, // 1.111111111111111
  156. {-3,2,0}, // 1.125000000000000
  157. {1,1,-1}, // 1.200000000000000
  158. {-2,0,1}, // 1.250000000000000
  159. {2,-1,0}, // 1.333333333333333
  160. {-1,1,0}, // 1.500000000000000
  161. {1,0,0} // 2.000000000000000
  162. } ;
  163.  
  164. void hammahead(unsigned long n, int *i0, int *j0, int *k0)  ;
  165.  
  166. int main(int argc, char **argv) {
  167.     double l2 = LOG2, l3 = LOG3, l5 = LOG5 ; // precomputed logarithms
  168.     double logM, logMlo, logMhi ; // per the Dr. Dobbs article, 
  169.                                   // logMlo < i*l2 + j*l3 + k*l5 < logMhi
  170.     int i, j, k ;  //loop counters 
  171.     int imin, imax, jmax, kmax ;  // limits for i,j,k loops
  172.     double logMmin = 1e9 ;  // minimum logM seen, defaults to big number
  173.  
  174.     unsigned long n ;  // want nth Hamming number
  175.     unsigned long totn = 0 ;  // total i,j,k combinations below logMlo
  176.  
  177.     int i0, j0, k0 ; // i,j,k values for minimum logM in band
  178.     unsigned long nahead ;  // number of values to skip ahead from i0, j0, k0
  179.  
  180.     unsigned char found = FALSE ;
  181.  
  182.     n = 1000000000000 ;
  183.  
  184.     if (n >= NMIN) {
  185.         // get range of logM
  186.         logM = pow((6.0*l2*l3*l5)*n,THIRD) ;
  187.         logMlo = logM - LOWOFFSET ;
  188.         logMhi = logM - HIOFFSET ;
  189.  
  190.         kmax = floor(logMhi/l5) ;
  191.         for (k = 0 ; k <= kmax ; k++) {
  192.             jmax = floor((logMhi - k*l5)/l3) ;
  193.             for (j = 0 ; j <= jmax ; j++) {
  194.                 imin = ceil((logMlo - k*l5 - j*l3)/l2) ;
  195.                 imax = floor((logMhi - k*l5 - j*l3)/l2) ;
  196.                 totn += imin ;
  197.                 for (i = imin ; i <= imax ; i++) {
  198.                     logM = i*l2 + j*l3 + k*l5 ;
  199.                     if (logM < logMmin) {
  200.                         logMmin = logM ;
  201.                         i0 = i ;
  202.                         j0 = j ;
  203.                         k0 = k ;
  204.                     }
  205.                 }
  206.             }
  207.         }
  208.         nahead = n - totn - 1 ;
  209.         hammahead(nahead,&i0,&j0,&k0) ;
  210.     }
  211.     else {
  212. 		int m, in ;
  213. 		i0 = 0 ;
  214. 		j0 = 0 ;
  215. 		k0 = 0 ;
  216. 		for (in = 0 ; in < n-1 ; in++) {
  217. 		    m = 0 ;
  218. 		    found = FALSE ;
  219. 		    while (!found && m < NTAB) {
  220. 			    if (i0 + hamtab[m].p2 >= 0)
  221. 			        if (j0 + hamtab[m].p3 >= 0)
  222. 			            if (k0 + hamtab[m].p5 >= 0)
  223. 						    found = TRUE ;
  224. 			    m++ ;
  225. 		    }
  226. 			m-- ;
  227. 			i0 += hamtab[m].p2 ;
  228. 			j0 += hamtab[m].p3 ;
  229. 			k0 += hamtab[m].p5 ;
  230.         }
  231.     }
  232.     printf("%d %d %d\n",i0,j0,k0) ;  // output nth i,j,k values
  233. }
  234.  
  235. // From the Hamming number 2^i0*3^j0*5^k0, skip ahead
  236. // nhead numbers and update i0, j0 and k0 accordingly.
  237. void hammahead(unsigned long nahead, int *i0, int *j0, int *k0) {
  238.     unsigned char k ;
  239.     unsigned long j ;
  240.     unsigned char found ;
  241.  
  242.     for (j = 0 ; j < nahead ; j++) {
  243.         // for each Hamming number, find the minimum interval in hamtab
  244.         // i.e., the minimum entry where multiplication would result
  245.         // in positive powers for i, j & k
  246.         k = 0 ;
  247.         found = FALSE ;
  248.         while (!found && k < NTAB) {
  249.             if (*i0 + hamtab[k].p2 >= 0)
  250.                 if (*j0 + hamtab[k].p3 >= 0)
  251.                     if (*k0 + hamtab[k].p5 >= 0)
  252.                             found = TRUE ;
  253.             if (!found)
  254.                 k++ ;
  255.         }
  256.  
  257.         if (!found) {
  258.             fprintf(stderr,"error: suitable table entry not found, j=%ld\n",j) ;
  259.             exit(1) ;
  260.         }
  261.         else {
  262.             // got the next interval, update counts
  263.             *i0 += hamtab[k].p2 ;
  264.             *j0 += hamtab[k].p3 ;
  265.             *k0 += hamtab[k].p5 ;
  266.         }
  267.     }
  268. }
  269.  
Success #stdin #stdout 0.51s 9432KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
1126 16930 40
https://ideone.com/r5lCLS
language:
C (gcc 8.3)
created:
8 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!