122

Community / Exhibition / Gallery / Re: Arduino generating prime numbers

on: November 26, 2013, 01:13:55 am

Eleven bad. That should be a one. If that statement ever executed, it should have caused an error; nevertheless, a test with five million primes ran without a hitch last night. Here's what I did to test it:  Added a "print something and stop" inside those brackets. The Due went to 120K primes found, and that didn't execute.
 Generated a whole bunch of uint64_t random pairs, multiplied them, and checked about 40K results in Perl with Math::GMPz. Again, the print and stop didn't execute. All the products were correct.
 Ran 100K random pairs, this time with the most significant nybble set to 0xF. Again, no print and stop, and products checked out.
 Rewrote the routine to operate on bytes, work on nybbles, and generate integer products; ran it against all 65K combinations of inputs, and checked. Again, no print and stop, and all the products checked.
 Multiplied 0xFFFFFFFFFFFFFFFF by itself. Nothing went wrong.
Apparently, that statement never executes  there's never a carry after the second addition to the intermediate product. That product captures the products of the high half of one input multiplied by the low half of the other. That happens twice, and those products are added, and that may generate a carry  that's the first test. Then the high half of the product of the low halves of the inputs is added to the result. I'd expect that addition to sometimes generate a carry, but it appears that it never does. I can see why it won't when the inputs are restricted to 54 bits, but I can't see how it never gets a carry with arbitrary inputs. If I can figure out why that test is unnecessary, I can remove it. [Edit: Changed "11ULL" to "1ULL" in file attached to reply #87; fixed spelling.]



123

Community / Exhibition / Gallery / Re: Arduino generating prime numbers

on: November 25, 2013, 04:35:19 pm

Indeed, the MillerRabin test calculates a whole lot of (x*y) % n's, where n is the number under test. There's benefit to be had from doing something first with n if it will simplify the mod operations. Here's "Barrett reduction," and it does just that, though without making a table: http://en.wikipedia.org/wiki/Barrett_reduction. It lets us calculate an integer function of the number under test, about the same size, and use it to get modulo results by multiplication, shifting, and subtraction. It adds complexity, but reduces the effort of multiple modulo calculations with the same modulus. The problem is that for this project, the multiplication yields a 3uint64_t product. Routines for multiple uint64_t's are inefficient compared to calculating with native data types, as they can't see the carry bit. Barrett provides a speed improvement, but not the whopping acceleration that I hoped for. It takes about 60% as long to calculate a 16digit prime, average, as the version from reply #81. It beats one second per prime, average, on the Uno. Barrett numbers k and m are added the cheesy way  as global variables, just before the MillerRabin tests. There are a bunch of tweaks that would improve the speed a little bit, but the ones that I can see are in routines that run once per candidate, and they're unlikely to make a big difference. The code spends twothirds of its time in function multll(), which multiplies two uint64_t's; I can't see way to improve that routine much. There may be ways that I can't see, or ways to eliminate calls to it. Barrett reduction was something of a chore to implement. I suspect that it would be a lot easier with BigNumber, but it doing it would require us to get into the raisemod() function in number.c. I'm not sure I want to try it, and I suspect that it wouldn't beat this speed, anyway. The code is so long that I can't fit it in with this post. It's attached, [Edit: attached file, changed this] along with the included file. Here's output for the representative primes that have become the staple of this thread, from the Uno: 1000000000054141 is prime. Took 849 mS. Found 4 primes. 1000000000054147 is prime. Took 432 mS. Found 5 primes. 1000000000054259 is prime. Took 995 mS. Found 6 primes.
The Due, of course, ran very fast. Times were about 10 ms for integers in the high end of the 16digit range. Accuracy looks great: five million primes, from the Due, exactly matched the 46th through 50th million primes from this page: http://primes.utm.edu/lists/small/millions/. [Edit: Fixed error in attached file.]



125

Community / Exhibition / Gallery / Re: Arduino generating prime numbers

on: November 19, 2013, 07:50:44 pm

One more tweak. This one doesn't use BigNumber, and it runs nearly twice as fast. We used BigNumber to get around the fact that intermediate results can overflow an unsigned long long in the MillerRabin test. That test requires calculating a ^{b} % n for large numbers. Using numbers that are less than 10 ^{16}  54 bits  the results always fit into a long long, but the intermediate products that generate the exponentials can be bigger than 2 ^{64}. BigNumber gets around the restriction, but at some cost in speed and storage. This code uses function "mulmodll()," that produces the product of two unsigned long long's, modulo a third unsigned long long, without overflowing on intermediate products. It generates the product internally as two unsigned long long's, and calculates the modulo by repeated shifts and mods. It has its own overhead  it has to take <something> % <something else> as many as seven times. But, it's working in native number formats, rather than in BCD, and it goes faster. This code depends heavily on the fact that numbers bigger than 10 ^{16} aren't investigated. That fact simplifies the 128bit multiplication, in that certain intermediate results can't overflow, and it simplifies the mod operations, too. The code isn't general. It's assumptions are that the lower half of any argument times the upper half of any argument will be less than 2 ^{63}, and that the modulus will always be less than 2 ^{56}. With 10 ^{16} = 2 ^{53.15}, the problem fits, with a bit or two to spare. Results: Running on the Uno, 100000 primes starting at 961,748,929, and comparing to the fiftieth million from this list  http://primes.utm.edu/lists/small/millions/  found an exact match, with calculation time averaging about 170 ms per prime. Using the tests of reply #80, calculation time for 16digit numbers added up to 134,578 ms, for an improvement of a factor of a little better than 5 over using BigNumber, and about 1.9 times as fast as using Numbler4. For 12digit numbers, calculation time was 52,472 ms, for a similar speed improvement. For numbers near the top of the range, starting at 9,999,999,999,990,001, the average calculation time for 100 primes was around 1.5 seconds  about three times what Nick was looking for. The Due was ferocious. It averaged 10.6 ms for 1,000 12digit primes starting at 100000000001; 25.8 ms for 16digit starting at 1000000000000001, and 33.4 starting at 9900000000036791. My Due ran these tests, and then stopped working  I'm hoping that it's just bored. There's room for improvement here. I think there's a sweet spot for the number of trial divisions, I think it depends to some extent on the size of the numbers under test, and I don't think I've found it here. There may be a way to reduce the number of modulo calls in function powmodll(), but I can't think of any yet. Here's the code. The included file is attached at reply #70: PROGMEM prog_uchar primes[] = { 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,101,103,107,109,113,127,131,137, 139,149,151,157,163,167,173,179,181,191,193, 197,199,211,223,227,229,233,239,241,251 }; const byte nPrimes = sizeof(primes)/sizeof(primes[0]);
#include "Pseudoprimes.h"; const uint16_t nPSP = (sizeof(psps)/sizeof(psps[0])) >> 3; uint64_t candidate;
unsigned long start; long found = 0; // Number of primes found
void setup() { Serial.begin(115200); Serial.println ("Starting.");
candidate = 961748927ULL; } // end setup()
void loop() { uint64_t nextPSP; int nextPSPindex; for (int i=0;i<nPSP;i++) { nextPSP = fetchPSP(i); if ((nextPSP == candidate)  (nextPSP > candidate)) { nextPSPindex = i; break; } nextPSP = 0; nextPSPindex = 0; } start = millis(); for ( ; candidate < 9999999999999999ULL; candidate += 2) { if (nextPSP != candidate) { if (isPrime(candidate)) { found++; showPrime (candidate); } } else { nextPSPindex++; if (nextPSPindex < nPSP) { nextPSP = fetchPSP(nextPSPindex); } } } found = 1; showPrime(1); found++; showPrime(2); candidate = 3; // back to start } // end loop()
bool Miller(uint64_t source) { if(source <= 3) { return true; } if((source & 1) == 0) { return false; }
uint64_t d = source  1; int s = 0;
while((d & 1) == 0) { d >>= 1; s ++; }
uint8_t a; const byte bases[] = {2, 3, 5, 7,};// 11, 13, 17, 19, 23 }; const byte nbases = sizeof(bases)/sizeof(bases[0]); for(int i = 0; i < nbases; i++) { a = bases[i]; uint64_t x = powmodll(a, d, source); if((x == 1)  (x == (source  1))) { continue; }
bool stillOK = false; for(int r = 1; r < s; r++) { x = mulmodll(x, x, source); if(x == 1) return false; if(x == (source  1)) { stillOK = true; break; } }
if(!stillOK) return false; }
return true; } // end Miller()
bool isPrime(uint64_t n) { if (!check(n)) { return false; } if (n > 0x10000) { if (!Miller (n)) { return false; } } return true; }
void showPrime (uint64_t n) { printll(n); Serial.print(" is prime. Took "); Serial.print(millis ()  start); Serial.print (" mS. Found "); Serial.print (found); Serial.println (" primes."); start = millis ();
} // end showIfPrime()
bool check(uint64_t n) { if (n <= 3) { return true; } if (!(n & 1)) { return false; } uint8_t i = 1; uint8_t f = pgm_read_byte_near(primes); while ((i < nPrimes) && (f < n)) { if ((n % f) == 0) { return false; } f = pgm_read_byte_near(primes+i); i++; } return true; } // end check()
uint64_t fetchPSP(int i) { i <<= 3; uint64_t n = 0; for (int8_t k=0;k<sizeof(n);k++) { n = n <<= 8; n = pgm_read_byte_near(psps+i+k); } return n; } // end fetchPSP()
void printll(uint64_t n) { // print an unsigned long long if (n == 0) { Serial.print('0'); return; } const byte DIGITS = 21; char s[DIGITS+1]; for (byte i=0;i<DIGITS;i++) { s[i] = 'X'; } s[DIGITS] = 0; byte i = DIGITS1; while (n > 0) { s[i] = (n % 10)  '0'; n = n / 10; i; } i++; char *z = &s[i]; Serial.print(z); } // end printll()
uint64_t powmodll(uint8_t a, uint64_t b, uint64_t c) { // Calculate a**b % c, uint64_t x = 1; uint64_t y = a;
while(b > 0) { if(b & 1) { x = mulmodll(x, y, c); } y = mulmodll(y, y, c); b >>= 1; } return x; } // end powmodll()
uint64_t mulmodll(uint64_t a, uint64_t b, uint64_t c) { // Calculate a*b %c uint64_t imh; // Intermediate results for multiplication uint64_t iml; uint64_t imm;
// Multiply a and b // Routine is not general: // relies on fact that MSB of inputs is zero // Product stored as two unsigned long long's imm = (a & 0xffffffffUL) * (b >> 32); iml = (a >> 32) * (b & 0xffffffffUL); imm += iml; // overflow for general inputs iml = (a & 0xffffffffUL) * (b & 0xffffffffUL); imm += iml >> 32; iml = (iml & 0xffffffffUL) + (imm << 32); imh = (a >> 32) * (b >> 32); imh += imm >> 32;
if (imh == 0) { return (iml % c); } for (uint8_t i=0;i<8;i++) { imh = imh << 8; imh = iml >> 56; iml = iml << 8; if ((imh >> 56) != 0) { if (imh > c) { imh %= c; } else { if (imh == c) { imh = 0; } } } } if (imh > c) { imh %= c; } return imh; } // end mulmodll()
Here's Uno output for comparison to reply #72, with clock speed the standard 16 MHz: 1000000000054141 is prime. Took 1433 mS. Found 1 primes. 1000000000054147 is prime. Took 778 mS. Found 2 primes. 1000000000054259 is prime. Took 1594 mS. Found 3 primes.



126

Community / Exhibition / Gallery / Re: Arduino generating prime numbers

on: November 13, 2013, 01:41:41 am

I don't understand why you seem to think that BigNumber is so quick. No offense intended. BigNumber is quick compared to my expectations. I'm surprised to see how well it does compared to operating with native integers. I had expected it to have significant overhead. OK, here are results with BigNumber, and with Numbler4 as posted at http://forum.arduino.cc/index.php?topic=85692.msg1416890#msg1416890. I transplanted powMod() from BigNumber into Numbler4. Tests ran on an Uno and a Due, and, just for fun, with a TI Stellaris Launchpad, compiled under the opensource Arduinolike IDE, "Energia. " Calculated primes agreed for all three devices. Results are total time in milliseconds for calculation of the first 100 primes after 10 ^{N1}, where N is 12 and 16: 12 digit 16 digit Uno: 102002 263207 254004 688223 Due: 14431 13790 34571 32051 Stellaris: 10349 10281 24643 24058
It ran significantly faster with Numbler4 on the Uno, with execution time less than 40% of the BigNumber version. It ran a little slower than BigNumber on the ARMs. On the Due, Numbler4 ran about 58% slower. There was less difference with the Stellaris, but it was still bit slower. I can qualitatively understand how eliminating a lot of modulo operations can speed things up on the Uno; my understanding of the internals of the ARMs is too thin to help me explain the results with the Due and Stellaris. I'm fool enough to speculate, though: TI says that the Stellaris' CortexM3 has hardware divide, and Atmel says the same thing about the SAM3X. My guess is that division on those processors is fast enough to be competitive with the alternative. Uno and Due programs were identical, generally the code from reply 66, maybe with a minor change or two. The unused function rng(), for generating random bases for the MillerRabin test, was deleted. The IDE v1.5.4 wouldn't compile it for the Due  it didn't like the random() function without arguments, preferring rand() instead. v1.5.2 was OK with random(). I don't even want to think about about the machinations I had to go through to get the Stellaris to run the program; just as well, since that kind of talk doesn't belong on an Arduino forum.



127

Community / Exhibition / Gallery / Re: Arduino generating prime numbers

on: November 11, 2013, 11:26:17 pm

One last tweak. This one includes a list of numbers that fool the MillerRabin test for bases 2, 3, 5, and 7. It compares each odd number to the next item in the list, and skips it if it finds them equal. The previous version, in reply #66, uses nine bases, as required to eliminate pseudoprimes below 10 ^{16}. This one uses four and a table lookup, so it's a lot faster. It kind of feels like cheating, though  Nick's versions knew that 2 an 3 are prime, while this program has a lot of arcane knowledge about special numbers. I was surprised to see BigNumber work so quickly compared to native number formats, as I'd expect it to carry substantial overhead. I suspect that its progenitors were developed, at least in part, for dealing with calculations of primes, and that may account for its fitness for this use. My thanks to Nick for posting this project, and for coming up with the BigNumber library. Hunting for primes is easy; speeding up the process is really tricky. And, I learned how to spell, "pseudo," I think, so my ability to eruditely describe things that aren't quite what they seem is much improved. A list of strong psuedo fake primes to bases 2 and 3 is here  http://www.belllabs.com/user/bleichen/diss/thesis.html. The CPAN module Math::Primality provided the engine for extracting the pseudoprimes that penetrate tests for bases 5 and 7. The included file containing the bytewide bits of the list is attached  "Pseudoprimes.h". The code has its oddities. Because there seems to be no way to store or access uint64_t values with PROGMEM, and because I couldn't get BigNumber to play nice with variables bigger than one byte, pseudoprimes are stored as 8 bytes, bigendian, and decoded into BigNumbers. Here's the sketch: #include <BigNumber.h> PROGMEM prog_uchar primes[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,101,103,107,109,113,127,131,137, 139,149,151,157,163,167,173,179,181,191,193, 197,199,211,223,227,229,233,239,241,251 }; const byte nPrimes = sizeof(primes)/sizeof(primes[0]);
#include "Pseudoprimes.h"; const uint16_t nPSP = (sizeof(psps)/sizeof(psps[0])) >> 3; BigNumber candidate; BigNumber one (1); BigNumber two (2);
unsigned long start; long found = 0; // Number of primes found
void setup() { Serial.begin(115200); Serial.println ("Starting.");
BigNumber::begin ();
candidate = BigNumber ("100000000003"); } // end of setup
void loop() { BigNumber nextPSP; int nextPSPindex; for (int i=0;i<nPSP;i++) { nextPSP = fetchPSP(i); if ((nextPSP == candidate)  (nextPSP > candidate)) { nextPSPindex = i; break; } nextPSP = 0; nextPSPindex = 0; }
start = millis(); for ( ; candidate < BigNumber("9999999999999999"); candidate += 2) { if (nextPSP != candidate) { showIfPrime (candidate); } else { nextPSPindex++; if (nextPSPindex < nPSP) { nextPSP = fetchPSP(nextPSPindex); } } } candidate = 3; // back to start! } // end of loop
bool Miller(BigNumber source) {
if(source == 2  source == 3) return true; if(source < two  source % two == 0) return false;
BigNumber d = source  one; int s = 0;
while(d % two == 0) { d /= two; s += one; }
BigNumber a; const byte bases[] = {2, 3, 5, 7,}; const byte nbases = sizeof(bases)/sizeof(bases[0]); for(int i = 0; i < nbases; i++) { a = BigNumber (bases[i]); BigNumber x = a.powMod (d, source); if(x == one  x == source  one) continue;
bool stillOK = false; for(int r = 1; r < s; r++) { x = x.powMod(two, source); if(x == one) return false; if(x == source  one) { stillOK = true; break; } }
if(!stillOK) return false; }
return true; } // end of Miller
void showIfPrime (BigNumber n) { if (!check(n)) return; if (!Miller (n)) return; found++; Serial.print(n); Serial.print(" is prime. Took "); Serial.print(millis ()  start); Serial.print (" mS. Found "); Serial.print (found); Serial.println (" primes."); start = millis ();
} // end of showIfPrime
bool check(BigNumber n) { BigNumber f; for (byte i = 0;i < nPrimes; i++) { f = BigNumber(pgm_read_byte_near(primes+i)); if ((n % f) == 0) { return false; } } return true; }
BigNumber fetchPSP(int i) { // Progmem array psps contains strong pseudoprimes // 8byte numbers, arranged as bytes, bigendian i <<= 3; BigNumber x256 ("256"); BigNumber n = 0; for (byte k=0;k<8;k++) { n = n * x256; uint8_t b = pgm_read_byte_near(psps+i+k); n = n + BigNumber(b); } return n; } Here's some output: 100000000003 is prime. Took 1676 mS. Found 1 primes. 100000000019 is prime. Took 2798 mS. Found 2 primes. 100000000057 is prime. Took 3539 mS. Found 3 primes. 100000000063 is prime. Took 1770 mS. Found 4 primes. Here's the output for a Due running the same code: 100000000003 is prime. Took 86 mS. Found 1 primes. 100000000019 is prime. Took 145 mS. Found 2 primes. 100000000057 is prime. Took 189 mS. Found 3 primes. 100000000063 is prime. Took 89 mS. Found 4 primes. [Edit: Add this] And, the Due in the 16digit zone: 1000000000000037 is prime. Took 416 mS. Found 1 primes. 1000000000000091 is prime. Took 236 mS. Found 2 primes. 1000000000000159 is prime. Took 529 mS. Found 3 primes. 1000000000000187 is prime. Took 261 mS. Found 4 primes.



128

Using Arduino / Audio / Re: Amplified Signal to Arduino Mega

on: November 07, 2013, 05:22:42 pm

... automatic gain control amplifiers, I assume this IC's advertised automatic level control is the same thing? Sameish. Automatic level control, as used in the datasheet, implies that the output will be within a vary narrow range. "...voltage input range of 1000." Do you mean over a range of 1000dB, 1000V, 1000x, or what? A range of 1000x. If the datasheet's claims are accurate, you should be able to design a circuit with that will deliver a constant AC voltage at the output while the input signal varies by a factor of 1000. I have drawn a schematic based on the diagram in the data sheet and attached. Does it look correct? I don't know. I'm not an expert in using the SA575. I have an application for this IC myself, and I was hoping that you'd design and debug something, and then tell me how to use it. I think that you're in for some experimentation. My suggestion is to implement the circuit shown in the datasheet, and see how it performs. Armed with that information, you'll be able to calculate how to manage the IC's output to use it as input to something else. I'd recommend considering the output voltage of your amp to be more like sqrt(2) * sqrt(50*6), and then add a safety factor, and then protect your inputs to avoid damage from overvoltage. The output power quoted for the amplifier is probably imprecise, and it's based on calculations that might be different from the ones that you and I would make. After you know more, then you can hook it up to the next delicate circuit.



129

Community / Exhibition / Gallery / Re: Arduino generating prime numbers

on: November 05, 2013, 01:14:09 pm

Got a speed improvement of an average factor of about 2.3 over the Nick's BigNumberbased code shown in reply #43, starting at 1000000000000001, 16 digits, by making these changes:  Added a trial division test using bytesized primes, before other testing.
 Modified the MillerRabin test to use 2, 3, 5, 7, 11, 13, 17, 19, and 23, rather than random numbers. This Wikipedia link  http://en.wikipedia.org/wiki/Strong_pseudoprime#Examples  says that there's no number less than 3825123056546413051, 19 digits, about 1.6 x 2^{61}, that will fool every singlepass MillerRabin tests with these bases.
 [Edit: Add this] Modified loop parameters to wield only integers = 6k + 1 and 6k  1, for integer k.
Most of the speed improvement appears to come from eliminating MillerRabin tests for numbers that divide by small primes. That suggests that there may be a sweet spot for the number of prime divisors in the trial division test that depends on the speed of the MillerRabin test, and it may be different from the one I chose. Lesser improvement comes from eliminating selection of random BigNumbers, and reducing the number of MillerRabin tests from ten to nine. [Edit: Add this] The 6K + 1 trick doesn't do much, since it only eliminates multiples of 3, and those fail the trial division test almost immediately. Average displayed time was about 11.5 mS for the primes shown in replies #44 and #45. It also ran about 3.5% slower than the values shown in reply #64  calculation times were sometimes wildly different, but on average, they were almost as good. Here's the code, based on the sketch in reply #43, minimally changed, autoformatted: #include <BigNumber.h> PROGMEM prog_uchar primes[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,101,103,107,109,113,127,131,137, 139,149,151,157,163,167,173,179,181,191,193, 197,199,211,223,227,229,233,239,241,251 }; const byte nPrimes = sizeof(primes)/sizeof(primes[0]);
const char FIRST_PRIME_TO_USE [] = "1000000000000001";
BigNumber candidate; BigNumber one (1); BigNumber two (2);
byte bump = 2;
void setup() { BigNumber five (5); BigNumber six (6);
Serial.begin(115200); while (!Serial) { }
BigNumber::begin ();
Serial.println (); Serial.println (); Serial.println ("Starting.");
candidate = BigNumber (FIRST_PRIME_TO_USE);
BigNumber rem = candidate % six; while ((rem != one) && (rem != five)) { rem = rem + one; rem = rem % six; candidate = candidate + one; } if (rem == one) { bump = 4; } } // end of setup
void loop() {
for ( ; candidate < BigNumber ("9999999999999999"); candidate += bump, bump = 6  bump) showIfPrime (candidate);
candidate = 3; // back to start! } // end of loop
void rng (BigNumber & result) { result = BigNumber (random ()); } // end of rng
bool Miller(BigNumber source) {
if(source == 2  source == 3) return true; if(source < two  source % two == 0) return false;
BigNumber d = source  one; int s = 0;
while(d % two == 0) { d /= two; s += one; }
const byte bases[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23 }; const byte nbases = sizeof(bases)/sizeof(bases[0]); for(int i = 0; i < nbases; i++) { // do // { // rng (a); // } // while(a < two  a >= source  two);
BigNumber a = BigNumber (bases[i]); BigNumber x = a.powMod (d, source); if(x == one  x == source  one) continue;
for(int r = 1; r < s; r++) { x = x.powMod(two, source); if(x == one) return false; if(x == source  one) break; }
if(x != source  one) return false; }
return true; } // end of Miller
const byte DIGITS = 16;
unsigned long start; unsigned long lastShown; long found = 1; // Number we found
void showIfPrime (BigNumber num) { if (!check(num)) return; if (!Miller (num)) return;
char buf [20] = { 0 }; long long temp = num; byte pos = 0;
// make printable for (pos = 0; pos < DIGITS; pos++) { if (temp == 0) break; char digit = temp % 10; buf [DIGITS  1  pos] = digit  '0'; temp /= 10; }
// insert leading spaces for (; pos < DIGITS; pos++) buf [DIGITS  1  pos] = ' '; // = 15;
Serial.print(num); Serial.print(" is prime. Took "); Serial.print(millis ()  start); Serial.print (" mS. Found "); Serial.print (found); Serial.println (" primes."); start = millis ();
found++;
lastShown = millis ();
} // end of showIfPrime
bool check(BigNumber n) { for (byte i = 0;i < nPrimes; i++) { BigNumber f = BigNumber(pgm_read_byte_near(primes + i)); if ((n % f) == 0) { return false; } } return true; }



130

Using Arduino / Audio / Re: Amplified Signal to Arduino Mega

on: October 31, 2013, 09:01:37 pm

This gizmo claims that it can operate as a limiter, and accommodate 60 dB on its inputs with only about 0.5 dB change on the outputs. http://www.onsemi.com//pub/Collateral/SA575D.PDFIt's available in a throughhole mount for something around $5.50 from Mouser, though the throughhole version is marked "endoflife." If that's true, you should be able to get the a flat output response over a voltage input range of 1000. Using the "+10 dB ~ twice as loud" rule of thumb, that's an apparent loudness range factor of 64. That might do it.



131

Using Arduino / Programming Questions / Re: kurtosis calculations on an analog signal

on: October 29, 2013, 07:32:50 pm

Every time loop() executes, the counter variable  i  is incremented, whether it's time to take a sample or not. loop() executes many times between samples, so the counter reaches its limit quite quickly, and the output data prints. Since the counter is supposed to be counting samples, it needs to be incremented only when a sample is taken. It needs to be inside the brackets of the "if (currentMillis ... " test. [Edit: Add this] After a closer look:  Proper indentation will make your code easier to read and understand, for others and for you. You can manage that yourself, or you can use "Auto Format" under the "Tools" menu in the Arduino IDE.
 I recommend that you refrain from adding complexity to this code until you get the kurtosis calculation working. Leave xbee stuff out of your program until you get kurtosis working.
 Patience, Grasshopper.



133

Using Arduino / Programming Questions / Re: kurtosis calculations on an analog signal

on: October 29, 2013, 08:15:10 am

Ingenious. That's a lot of reasonablelooking code to come from a guy who didn't know where to start a few days ago. Here are some suggestions:  [Edit: Add this]
Something's wrong with the timing and program flow. The code decides when to take a sample, but it processes a sample regardless of whether or not it takes one. Examine it: loops starts; the test passes, it takes a sample, and processes it; loop terminates, and starts over; the test fails, and it doesn't take a sample, but it processes the last sample again anyway. Most of that code needs to move into the "if(currentMillis ..." section, and the rest should only be executed when it's time to do final calculations and print.  You're trying to calculate a function of a series of values, and you don't yet know whether the calculation is correct. You also don't know what the values are, since they're analog readings, and they aren't displayed. Troubleshooting the calculation is difficult enough without having to wonder what the input data was. While you're troubleshooting the calculation, use some convenient calculated data whose kurtosis is wellknown, instead of analog readings. Shutting off sample timing during troubleshooting will make the results show up faster. Here are two suggestions:
xi = i % 10 // uniform distribution, expected excess kurtosis = 1.2 xi = sin(i*0.001*2*3.1415927*100) + 1 // sine function, 100 Hz, sampled at 1/millisecond, expected excess kurtosis = 1.5  Decide whether you want to calculate kurtosis the oldfashioned way, or whether you want to calculate excess kurtosis.
 Print some debugging information, like the value of the running average, fourth moment, and variance. Be prepared to get deeper into the calculations if those values don't help.
 Your "running average" technique, where newValue = alpha*sampleValue + (1alpha)*oldValue, is a singlepole digital lowpass filter. Its output value tends toward the mean of the input signal. When you filter variance, you get something that looks like the statistical variance  the mean of the squares of the differences between the samples and the mean of the signal. When you filter the fourth moment, you get something that looks like the fourth moment  the mean of the fourth power of those same differences. Division by the number of samples is essentially accomplished in filtering those values. Consequently, I'll recommend that you reconsider whether this line needs the multiplication by the number of samples:
denominator = sigma4*numberSamples;  When you start getting output data, experiment with larger values of alpha. I don't think it's necessary that alpha be the same as 1.0/numberSamples. You might want to define beta = 1.0  alpha.
 This looks to me to be an ingenious solution. My experience with ingenious solutions to statistical calculations is that they often yield unpextected results, usually because there's some nugget of probability theory that I've overlooked that makes them invalid. I can't see the flaw here, though  the level of filtering is pretty intense, and the signal mean will get pretty close to the actual value; once it stabilizes, I think that the variance and fourth moment calculations will be valid. Nevertheless, I'd recommend testing it with several kinds of distributions to see that it gives the expected results, and be prepared for disappointment.
 After you're confident of the kurtosis calculation, think about ways to mange the calculation of the signal mean. When I run this code with known input data and print the value of runningaverage, it takes a while to stabilize. Until it does, the output data will not be valid. If the DC offset of the signal won't change during the measurement, consider instead taking some time to determine the signal mean while outputting no values, and then using that mean to calculate the second and fourth moments. There's no loss of performance  there's a delay between startup and valid data either way. With the mean known, the calculations are pretty simple. Once the signal mean is established, I think that there's no necessity to filter the second and fourth moments  the calculation can be straightforward from the formula.
 Consider ACcoupling the input signal. The signal mean will then be the value of the analog offset voltage. That'll give control of the signal mean to you, rather than to some setting on a sensor amplifier. There are a lot of options for how to manage getting that voltage into the calculation, some involving hardware, some not. A simple one would be to use a resistor divider to set it at half scale, and use that value as the initial value of runningaverage, rather than starting with zero. That would reduce the time it takes for data to become valid, though possibly not by all that much  the second and fourth moments are still filtered.
 As PaulS suggests, don't make sensorPin a float. It seems to work, but it's just too weird. sensorPin doesn't change, so it could easily be const int.



134

Using Arduino / Programming Questions / Re: kurtosis calculations on an analog signal

on: October 28, 2013, 04:25:12 pm

It's good to see the original poster put up some code. It suggests that he's actually interested in a solution.
The posted code doesn't compile. There are repeated references to kurtvoltages, as if it were a variable. But kurtvoltages is an array, and, if it's going to be treated like a single value, it needs an index.
There are logic errors, too, but they can wait. If it doesn't compile, you'll never get any results at all.
[Edit: Add this] And, until you know what you're doing, please don't try to hook up an analog source  like a signal generator or a biological monitor  to your Arduino. It can be quite easy to damage an analog input when the input voltage goes out of bounds. I suspect that your signal source is centered at zero volts, and that it's negative for at least some of the time. A negative voltage on an analog pin can damage it.
This sounds like an academic assignment. Is it?



