Sunday, April 16, 2017

Digest for comp.lang.c++@googlegroups.com - 25 updates in 4 topics

Marcel Mueller <news.5.maazl@spamgourmet.org>: Apr 13 04:31PM +0200

AFAIK the address of an external symbol with static linkage can be used
as constexpr (e.g. a method reference as template parameter).
 
Is it also possible to import /absolute/ linker symbols this way?
 
Example
 
shader.o:
00000000 R shader
00000410 A shader_size
 
shader.h:
// import the shader symbol from shader.o, OK
extern unsigned shader[];
 
// Try to import shader_size - does not work
extern const int shader_size;
constexpr int shader_size_proxy = (int)&shader_size;
#define shader_size shader_size_proxy
 
// without constexpr
const int shader_size_proxy2 = (int)&shader_size;
#undef shader_size
#define shader_size shader_size_proxy2
 
 
The basic problem is that C++ only imports addresses but not absolute
values like the global symbol shader_size. And the following required
conversion from int* to int is no longer a valid constexpr.
 
Is there anything to come around this restriction?
 
 
Marcel
"Alf P. Steinbach" <alf.p.steinbach+usenet@gmail.com>: Apr 13 06:24PM +0200

On 13-Apr-17 4:31 PM, Marcel Mueller wrote:
> values like the global symbol shader_size. And the following required
> conversion from int* to int is no longer a valid constexpr.
 
> Is there anything to come around this restriction?
 
Why can't you just provide the definition in a header file?
 
 
Cheers!,
 
- Alf
"Chris M. Thomasson" <invalid@invalid.invalid>: Apr 12 09:51PM -0700

Here is a quick and crude example of using the ct_roots function
described in the following thread:
 
"n-ary roots from complex numbers..."
https://groups.google.com/d/topic/comp.lang.c++/05XwgswUnDg/discussion
 
to actually store data. The following code stores my name "CHRIS" from
the following symbol table:
 
"0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"
 
in a single complex number via root map. For storing data, it maps each
symbol index to a root in the complex number z during iteration. For
loading data, it raises z back up through the root chain. Loading needs
the complex number, the origin number (for termination), and the correct
power, and it can recreate the stored data. It locates each root by
doing a little search in the n-ary roots returned by ct_roots. The
experimental function is ct_try_find. Take note of it because it is key
wrt being able to load the stored data... There has to be a more
efficient way of doing the root search. Humm... For termination, it
compares z iterates to the original number used for storing the information.
 
 
the function ct_store stores n-ary tokens in a single complex number.
the function ct_load loads n-ary tokens from a complex number.
 
The termination condition is iffy and really sensitive to floating point
errors. Its better to store data in strict blocks, so a load knows to
iterate a block length, and then quite. Now it has to compare the
iterates of z to the damn origin. Humm... The can be improved upon.
 
We cannot store too much data here, or one with get a data-incoherent
error, via assertion. I can store the symbols DEADBEEF for sure. ;^)
 
This is where floating point errors can really break things down.
 
 
Btw, I am wondering if you can run it and perhaps post some of the
output? It would be interesting to compare the results I get on to other
systems.
 
 
Before I go on, can anybody run this code? Also, I think this deserved
its own thread? Will get back to the group tomorrow.
 
Thanks everybody.
 
__________________________________________
// 4/12/2017: n-ary complex storage example by Chris M. Thomasson
 
#include <complex>
#include <iostream>
#include <vector>
#include <limits>
#include <algorithm> // reverse
#include <cstdint> // want to work with 64-bit numbers
#include <cassert> // want to sanity check run time...
#include <cstring> // faster than iostream's?
 
 
// Well, I need some consistent typenames... ;^)
typedef std::int64_t ct_int;
typedef std::uint64_t ct_uint;
typedef double ct_float;
typedef std::numeric_limits<ct_float> ct_float_nlim;
typedef std::complex<ct_float> ct_complex;
typedef std::complex<ct_uint> ct_complex_uint;
typedef std::vector<ct_complex> ct_complex_vec;
 
#define CT_PI 3.14159265358979323846
 
 
// Round up and convert the real and imaginary
// parts of z to unsigned integers of type ct_uint
// return a complex number with unsigned integer parts
ct_complex_uint
ct_round_uint(
ct_complex const& z
) {
ct_uint re = (ct_uint)std::floor(std::abs(z.real()) + .5);
ct_uint im = (ct_uint)std::floor(std::abs(z.imag()) + .5);
return ct_complex_uint(re, im);
}
 
 
// the integer p shall not be zero
// create abs(p) roots of z wrt z^(1/p);
// store them in out, and return the average error.
ct_float
ct_roots(
ct_complex const& z,
ct_int p,
ct_complex_vec& out
) {
assert(p != 0);
 
// Gain the basics
ct_float radius = std::pow(std::abs(z), 1.0 / p);
ct_float angle_base = std::arg(z) / p;
ct_float angle_step = (CT_PI * 2.0) / p;
 
// Setup the iteration
ct_uint n = std::abs(p);
ct_float avg_err = 0.0;
 
// Calculate the n roots...
for (ct_uint i = 0; i < n; ++i)
{
// our angle
ct_float angle = angle_step * i;
 
// our point
ct_complex c = {
std::cos(angle_base + angle) * radius,
std::sin(angle_base + angle) * radius
};
 
// output data
out.push_back(c);
 
// Raise our root the the power...
ct_complex raised = std::pow(c, p);
 
// Sum up the Go% damn floating point errors!
avg_err = avg_err + std::abs(raised - z);
}
 
// gain the average error sum... ;^o
return avg_err / n;
}
 
 
// Try's to find the target root z out of roots using
// eps, return the index of the root, or -1 for failure.
int
ct_try_find(
ct_complex const& z,
ct_complex_vec const& roots,
ct_float eps
) {
std::size_t n = roots.size();
 
for (std::size_t i = 0; i < n; ++i)
{
ct_complex const& root = roots[i];
 
ct_float adif = std::abs(root - z);
 
if (adif < eps)
{
return i;
}
}
 
return -1;
}
 
 
 
 
 
// The Token Table
// Will deal with scrambled token vectors in further posts.
// This is global for now, easy to convert to per store/load
// pairs
static std::string const g_tokens_str =
"0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
 
// Gains the power from the largest token position
// from tokens in g_tokens_str
ct_int
ct_gain_power(
std::string const& tokens
) {
ct_uint n = tokens.length();
 
std::size_t pmax = 0;
 
for (ct_uint i = 0; i < n; ++i)
{
std::size_t fridx = g_tokens_str.find_first_of(tokens[i]);
assert(fridx != std::string::npos);
pmax = std::max(pmax, fridx);
}
 
return (ct_int)(pmax + 1);
}
 
 
// Store tokens using a power of p starting in z-origin
// Return the complex number holding said tokens.
ct_complex
ct_store(
ct_complex const& z_origin,
ct_int p,
std::string const& tokens
) {
ct_uint n = tokens.length();
 
ct_complex z = z_origin;
ct_float store_avg_err = 0.0;
 
std::cout << "Storing Data..." << "\n";
std::cout << "stored:z_origin:" << z_origin << "\n";
 
for (ct_uint i = 0; i < n; ++i)
{
// Gain all of the roots, and running store error
ct_complex_vec roots;
ct_float avg_err = ct_roots(z, p, roots);
store_avg_err = store_avg_err + avg_err;
 
// reference our root
std::size_t fridx = g_tokens_str.find_first_of(tokens[i]);
assert(fridx != std::string::npos);
 
z = roots[fridx];
 
std::cout << "stored[" << i << "]:" << z << "\n";
}
 
store_avg_err = store_avg_err / n;
std::cout << "store_avg_err:" << store_avg_err << "\n";
 
return z;
}
 
 
// Load our tokens from z_store, power of p,
// stopping at z_target using eps, storing tokens
// in out_tokens, and the resulting z in out_z
ct_float
ct_load(
ct_complex const& z_store,
ct_complex const& z_target,
ct_int p,
ct_float eps, // epsilon
std::string& out_tokens,
ct_complex& out_z
) {
ct_complex z = z_store;
ct_uint n = 128; // max iter
ct_float load_err_sum = 0.0;
 
std::cout << "Loading Data..." << "\n";
for (ct_uint i = 0; i < n; ++i)
{
// raise...
ct_complex z_next = std::pow(z, p);
 
// Gain all of the roots, and running load error
ct_complex_vec roots;
ct_float avg_err = ct_roots(z_next, p, roots);
load_err_sum += avg_err;
 
// try to find our root...
int root_idx = ct_try_find(z, roots, eps);
if (root_idx < 0 ||
(ct_uint)root_idx >= g_tokens_str.length()) break;
 
std::cout << "loaded[" << i << "]:" << z << "\n";
 
out_tokens += g_tokens_str[root_idx];
 
// advance
z = z_next;
 
// check for termination condition...
if (std::abs(z - z_target) < eps)
{
std::cout << "fin detected!:[" << i << "]:" << z << "\n";
break;
}
}
 
// reverse our tokens
std::reverse(out_tokens.begin(), out_tokens.end());
 
out_z = z;
return load_err_sum;
}
 
 
int main()
{
std::cout.precision(ct_float_nlim::max_digits10);
 
std::cout << "g_tokens_str:" << g_tokens_str << "\n\n";
 
{
ct_complex z_origin = { -.75, .06 };
 
// The original data to be stored
std::string stored = "CHRIS";
ct_int power = ct_gain_power(stored);
 
std::cout << "stored:" << stored << "\n";
std::cout << "power:" << power << "\n\n";
std::cout << "________________________________________\n";
 
// STORE
ct_complex z_stored = ct_store(z_origin, power, stored);
 
std::cout << "________________________________________\n";
 
 
std::cout << "\nSTORED POINT:" << z_stored << "\n";
 
 
std::cout << "________________________________________\n";
 
// The data loaded from the stored.
std::string loaded;
ct_complex z_loaded;
ct_float eps = .001; // epsilon
 
// LOAD
ct_float load_err_sum =
ct_load(z_stored, z_origin, power, eps, loaded, z_loaded);
 
std::cout << "________________________________________\n";
 
std::cout << "\nORIGIN POINT:" << z_origin << "\n";
std::cout << "LOADED POINT:" << z_loaded << "\n";
 
std::cout << "\nloaded:" << loaded << "\n"
"load_err_sum:" << load_err_sum << "\n";
 
// make sure everything is okay...
if (stored == loaded)
{
std::cout << "\n\nDATA COHERENT! :^D" << "\n";
}
 
else
{
std::cout << "\n\n***** DATA CORRUPTED!!! Shi%! *****" << "\n";
assert(stored == loaded);
}
}
 
// Fin...
std::cout << "\n\nFin, hit <ENTER> to exit...\n";
std::fflush(stdout);
std::cin.get();
 
return 0;
}
__________________________________________
"Alf P. Steinbach" <alf.p.steinbach+usenet@gmail.com>: Apr 13 07:52AM +0200

On 13-Apr-17 6:51 AM, Chris M. Thomasson wrote:
 
> in a single complex number via root map. For storing data, it maps each
> symbol index to a root in the complex number z during iteration. For
> loading data, it raises z back up through the root chain.
 
Uhm, this sounds exquisitely without practical application. :)
 
 
Cheers!,
 
- Alf
"Chris M. Thomasson" <invalid@invalid.invalid>: Apr 13 09:24AM -0700

On 4/12/2017 10:52 PM, Alf P. Steinbach wrote:
>> symbol index to a root in the complex number z during iteration. For
>> loading data, it raises z back up through the root chain.
 
> Uhm, this sounds exquisitely without practical application. :)
 
Maybe... ;^)
 
Fwiw, I am currently experimenting with using the technique in fractal
encryption. The form used is z^n+c to load data, and (z-c)^(1/n) to
store data in its roots. I personally think its interesting how we can
store and load data using inverse iteration of an n-ary Julia set.
 
https://en.wikipedia.org/wiki/Julia_set#Using_backwards_.28inverse.29_iteration_.28IIM.29
red floyd <dont.bother@its.invalid>: Apr 11 09:28AM -0700

On 4/11/2017 8:58 AM, bitrex wrote:
> if any of the following is true:
 
> T has a non-static data member of non-class type (or array thereof) that
> is const;
 
[redacted]
 
I'm no language guru, but it seems to me that _ref2 is a non-static data
member of non-class type that is const.
bitrex <bitrex@de.lete.earthlink.net>: Apr 11 11:58AM -0400

Suppose I have a class such as:
 
class Foo
{
public:
Foo(std::stack<Bar>& ref1, Baz *const ref2) :
_ref1(ref1), _ref2(ref2)
{}
 
// etc.
 
private:
std::stack<Bar>& _ref1;
Baz *const _ref2;
};
 
The standard says (for C++11):
 
"A defaulted copy assignment operator for class T is defined as deleted
if any of the following is true:
 
T has a non-static data member of non-class type (or array thereof) that
is const;
 
T has a non-static data member of a reference type;
 
T has a non-static data member or a direct or virtual base class that
cannot be copy-assigned (overload resolution for the copy assignment
fails, or selects a deleted or inaccessible function);
 
T is a union-like class, and has a variant member whose corresponding
assignment operator is non-trivial."
 
It also says:
 
"The copy assignment operator for class T is trivial if all of the
following is true:
 
it is not user-provided (meaning, it is implicitly-defined or
defaulted), , and if it is defaulted, its signature is the same as
implicitly-defined (until C++14);
 
T has no virtual member functions;
 
T has no virtual base classes;
 
the copy assignment operator selected for every direct base of T is trivial;
 
the copy assignment operator selected for every non-static class type
(or array of class type) member of T is trivial.
 
A trivial copy assignment operator makes a copy of the object
representation as if by std::memmove. All data types compatible with the
C language (POD types) are trivially copy-assignable."
 
 
 
Which of these situations takes precedence for a class like my example?
It does contain a reference, yet given that it has no other fields other
than what are essentially pointers which aren't allowed to change, a
std::memmove should work fine for copying it.
Bonita Montero <Bonita.Montero@gmail.com>: Apr 11 08:08PM +0200

> std::memmove should work fine for copying it.
 
LOL!
bitrex <bitrex@de.lete.earthlink.net>: Apr 12 11:42AM -0400

On 04/11/2017 02:27 PM, Alf P. Steinbach wrote:
 
 
> Lols. :)
 
> Cheers & hth.,
 
> - Alf
 
This seems to work fine (at least on my architecture and compiler.)
Though is probably undefined behavior. At least from a "C++ outsider's"
perspective, why you couldn't at least from a theoretical implementation
standpoint memmove a class whose fields consisted solely of immutable
pointers or aliases for such and otherwise met the requirements for
being trivially copyable escapes me.
 
#include <stack>
#include <cstring>
#include <iostream>
 
struct Foo
{
Foo(std::stack<int>& thing1,
int *const thing2) :
_thing1(thing1),
_thing2(thing2)
{}
 
std::stack<int>& _thing1;
int *const _thing2;
};
 
int main()
{
std::stack<int> a;
std::stack<int> b;
int c = 5;
int *const d = &c;
int e = 6;
int *const f = &e;
 
a.push(c);
b.push(e);
 
auto bar = Foo(a, d);
auto baz = Foo(b, f);
 
std::cout << *bar._thing2 << std::endl;
std::cout << *baz._thing2 << std::endl;
std::cout << bar._thing1.top() << std::endl;
std::cout << baz._thing1.top() << std::endl;
 
std::memmove(&bar, &baz, sizeof(Foo));
 
std::cout << *bar._thing2 << std::endl;
std::cout << *baz._thing2 << std::endl;
std::cout << bar._thing1.top() << std::endl;
std::cout << baz._thing1.top() << std::endl;
}
"Alf P. Steinbach" <alf.p.steinbach+usenet@gmail.com>: Apr 11 08:27PM +0200

On 11-Apr-17 5:58 PM, bitrex wrote:
 
> Which of these situations takes precedence for a class like my example?
> It does contain a reference, yet given that it has no other fields other
> than what are essentially pointers which aren't allowed to change, a
 
[c:\Temp]
> type foo.cpp
#include <stdio.h>
 
struct X
{
int& ref;
};
 
auto main()
-> int
{
int gosh = 42;
X a{ gosh }, b{ gosh };
a = b;
}
 
[c:\Temp]
> g++ foo.cpp
foo.cpp: In function 'int main()':
foo.cpp:13:9: error: use of deleted function 'X& X::operator=(const X&)'
a = b;
^
foo.cpp:3:8: note: 'X& X::operator=(const X&)' is implicitly deleted
because the default definition would be ill-formed:
struct X
^
foo.cpp:3:8: error: non-static reference member 'int& X::ref', can't use
default assignment operator
 
[c:\Temp]
> _
 
 
> std::memmove should work fine for copying it.
 
Lols. :)
 
 
Cheers & hth.,
 
- Alf
"Chris M. Thomasson" <invalid@invalid.invalid>: Apr 10 04:09PM -0700

The goal of this crude sample code is to take the n-ary roots of a
complex number. Iterating a complex number z to all of its z^(1/n)
roots. This code derives the polar form from z, and iterates all of the
roots. This is an essential part of providing the ability of storing
n-ary tokens in a set of iterated roots. We can map the n roots, to a
token table, and store data. I will explain more later, however, please
take a look at the code, run it, tear it apart, and give me pointers on
how to make it better, please? Thank you all for your time. The function
at the heart of this is the function ct_roots that breaks apart a
complex number into its n-ary roots. Here is the code:
 
______________________________
#include <complex>
#include <iostream>
#include <vector>
#include <limits>
#include <cstdint> // want to work with 64-bit numbers
#include <cassert> // want to sanity check run time...
 
 
// Well, I need some consistent typenames... ;^)
typedef std::int64_t ct_int;
typedef std::uint64_t ct_uint;
typedef double ct_float;
typedef std::numeric_limits<ct_float> ct_float_nlim;
typedef std::complex<ct_float> ct_complex;
typedef std::complex<ct_uint> ct_complex_uint;
typedef std::vector<ct_complex> ct_complex_vec;
 
#define CT_PI 3.14159265358979323846
 
 
// Round up and convert the real and imaginary
// parts of z to unsigned integers of type ct_uint
// return a complex number with unsigned integer parts
ct_complex_uint
ct_round_uint(
ct_complex const& z
) {
ct_uint re = (ct_uint)std::floor(std::abs(z.real()) + .5);
ct_uint im = (ct_uint)std::floor(std::abs(z.imag()) + .5);
return ct_complex_uint(re, im);
}
 
 
// the integer p shall not be zero
// create abs(p) roots of z wrt z^(1/p);
// store them in out, and return the average error.
ct_float
ct_roots(
ct_complex const& z,
ct_int p,
ct_complex_vec& out
) {
assert(p != 0);
 
// Gain the basics
ct_float radius = std::pow(std::abs(z), 1.0 / p);
ct_float angle_base = std::arg(z) / p;
ct_float angle_step = (CT_PI * 2.0) / p;
 
// Setup the iteration
ct_uint n = std::abs(p);
ct_float avg_err = 0.0;
 
// Calculate the n roots...
for (ct_uint i = 0; i < n; ++i)
{
// our angle
ct_float angle = angle_step * i;
 
// our point
ct_complex c = {
std::cos(angle_base + angle) * radius,
std::sin(angle_base + angle) * radius
};
 
// output data
out.push_back(c);
 
// Raise our root the the power...
ct_complex raised = std::pow(c, p);
ct_complex_uint raised_rounded = ct_round_uint(raised);
 
// Output...
std::cout << "ct_roots::raised[" << i << "]:" <<
raised << "\n";
std::cout << "ct_roots::raised_rounded[" << i << "]:" <<
raised_rounded << "\n";
std::cout << "___________________________________\n";
 
 
// Sum up the Go% damn floating point errors!
avg_err = avg_err + std::abs(raised - z);
 
}
 
// gain the average error sum... ;^o
return avg_err / n;
}
 
 
 
 
int main()
{
std::cout.precision(ct_float_nlim::max_digits10);
 
{
// Our origin point z
ct_complex z = { 94143178827, 0 };
 
std::cout << "z:" << z << " - " << std::abs(z) << "\n";
std::cout << "**********************************************\n\n";
 
// The call to ct_roots...
ct_complex_vec roots;
ct_float avg_err = ct_roots(z, 23, roots);
std::cout << "\n\navg_err:" << avg_err << "\n";
 
// Output
std::cout << "___________________________________\n";
for (unsigned int i = 0; i < roots.size(); ++i)
{
ct_complex& root = roots[i];
 
std::cout << "root[" << i << "]:" << root << " - " <<
std::abs(root) << "\n";
}
std::cout << "___________________________________\n";
}
 
// Fin...
std::cout << "\n\nFin, hit <ENTER> to exit...\n";
std::fflush(stdout);
std::cin.get();
 
return 0;
}
______________________________
 
 
Can you run this? BTW, what about possibly reducing the average error
ratio? Any ideas?
 
Will post more context later on tonight or tommorw.
 
:^)
"Chris M. Thomasson" <invalid@invalid.invalid>: Apr 10 07:35PM -0700


> I'll give it a go tomorrow! Right now I'm working on a project to check out
> nearly a thousand proposed fractal formulas *grin* I'm writing a parser for
> them to rewrite as Ulta Fractal formulas which allows a quick worthiness check…
 
Nice! I need a good parser for fractal formulas. :^)
"Alf P. Steinbach" <alf.p.steinbach+usenet@gmail.com>: Apr 11 08:23AM +0200

On 11-Apr-17 7:41 AM, Christian Gollwitzer wrote:
> functions. Concerning accuracy, your algorithm is more accurate if you
> compute roots of very high order (say 1e10) - probably not for useful
> cases.
 
I still don't understand how rotation, and thereby the structure
(metric) of our physical environment, emerges in complex arithmetic.
 
It's amazing.
 
It's simple to show that `i` multiplied by itself goes round through
four positions (i, -1, -i, 1, ...), and then that with e.g. the point 45
degrees up right, the powers give a rotation through eight positions,
including `i` so that it can be identified with the square root of `i`,
and so on. But what /causes/ this behavior eludes me. And even though
I'm old and pretty stupid now, I have been good at math, and I believe
I've been thinking about this also in my younger days.
 
 
Cheers!,
 
- Alf (baffled, again :) )
"Chris M. Thomasson" <invalid@invalid.invalid>: Apr 11 01:39PM -0700

On 4/10/2017 10:41 PM, Christian Gollwitzer wrote:
>> ct_float radius = std::pow(std::abs(z), 1.0 / p);
>> ct_float angle_base = std::arg(z) / p;
>> ct_float angle_step = (CT_PI * 2.0) / p;
[..]
> myroots.push_back(baseroot);
 
> ct_complex root = baseroot;
> ct_complex phasor = baseroot / std::abs(baseroot);
 
Afaict, phasor will always be 1+0i.
 
 
Christian Gollwitzer <auriocus@gmx.de>: Apr 11 10:45PM +0200

Am 11.04.17 um 21:13 schrieb Chris M. Thomasson:
> [...]
 
> The shows all of the roots. Yours only shows the primary root.
 
> Please correct me if I am totally wrong.
 
puh. this endless listing of code and results makes it hard to find the
bug. The problem is that the phasor is wrong; it is not z/|z| as I
posted (a "unit" complex pointing to z) but a unit complex with the
angle 2*pi/n.
Try:
 
const ct_complex i = ct_complex(0,1);
ct_complex phasor = std::exp(2 * i * M_PI / p);
 
(untested, again).
 
Sorry for the mistake.
 
Christian
Christian Gollwitzer <auriocus@gmx.de>: Apr 11 10:58PM +0200

Am 11.04.17 um 22:39 schrieb Chris M. Thomasson:
>> ct_complex root = baseroot;
>> ct_complex phasor = baseroot / std::abs(baseroot);
 
> Afaict, phasor will always be 1+0i.
 
For positive real input, yes - try a complex number like 1+2i to see
what it gives. Of course, the correct solution is a different number.
 
Christian
"Chris M. Thomasson" <invalid@invalid.invalid>: Apr 11 02:00PM -0700

On 4/11/2017 1:45 PM, Christian Gollwitzer wrote:
 
>> Please correct me if I am totally wrong.
 
> puh. this endless listing of code and results makes it hard to find the
> bug.
 
sorry about that.
 
 
> const ct_complex i = ct_complex(0,1);
> ct_complex phasor = std::exp(2 * i * M_PI / p);
 
> (untested, again).
 
Works like a charm. Thank you Christian. :^)
 
As you mentioned, the average error is greater at:
 
cg_roots = avg_err:0.0010976249994551542
 
vs. the slow trig original:
 
ct_roots = avg_err:0.00056815732700208645
 
but it works fine. Thanks again.
 
 
> Sorry for the mistake.
 
No problem at all! :^)
"Chris M. Thomasson" <invalid@invalid.invalid>: Apr 11 02:15PM -0700

On 4/11/2017 1:58 PM, Christian Gollwitzer wrote:
 
>> Afaict, phasor will always be 1+0i.
 
> For positive real input, yes - try a complex number like 1+2i to see
> what it gives. Of course, the correct solution is a different number.
 
Doh! Of course you are correct here. Sorry.
"Chris M. Thomasson" <invalid@invalid.invalid>: Apr 11 02:20PM -0700

On 4/11/2017 11:54 AM, Stefan Ram wrote:
> generator of rotations around the z axis. (The following
> pages then explain how finite rotations can be obtained
> from such "infinitesimal rotations".)
 
Nice.
Christian Gollwitzer <auriocus@gmx.de>: Apr 11 07:41AM +0200

Am 11.04.17 um 01:09 schrieb Chris M. Thomasson:
> std::cos(angle_base + angle) * radius,
> std::sin(angle_base + angle) * radius
> };
 
 
While this is correct, I think that the following is easier and should
run faster:
 
std::vector<ct_complex> myroots;
ct_complex baseroot = std::pow(z, 1.0/p);
myroots.push_back(baseroot);
 
ct_complex root = baseroot;
ct_complex phasor = baseroot / std::abs(baseroot);
for (int i=1; i<n; i++) {
root *= phasor;
myroots.push_back(root);
}
 
(untested)
 
This algorithm saves the computation of 2*(n-1) transcendental
functions. Concerning accuracy, your algorithm is more accurate if you
compute roots of very high order (say 1e10) - probably not for useful
cases.
 
Christian
"Chris M. Thomasson" <invalid@invalid.invalid>: Apr 11 02:25PM -0700

On 4/10/2017 11:23 PM, Alf P. Steinbach wrote:
> and so on. But what /causes/ this behavior eludes me. And even though
> I'm old and pretty stupid now, I have been good at math, and I believe
> I've been thinking about this also in my younger days.
[...]
 
Once one converts a complex number z to polar form, raising to a power
becomes clear to me. Christian's use of raising e to the power of
 
2.0 * i * CT_PI / (ct_float)p
 
_________________
const ct_complex i = ct_complex(0, 1);
ct_complex phasor = std::exp(2.0 * i * CT_PI / (ct_float)p);
_________________
 
is right on par! ;^)
 
A little less accurate, but much much faster!
Marcel Mueller <news.5.maazl@spamgourmet.org>: Apr 12 08:42AM +0200

On 11.04.17 07.41, Christian Gollwitzer wrote:
> root *= phasor;
> myroots.push_back(root);
> }
 
I did not ckeck whether your code gives correct results, but this kind
of calculation results in amazingly large floating point errors. A
similar method is used in the GPU FFT algorithm for the Raspberry Pi and
has errors in the order of 10^-3 for n of a few hundred. This even bad
for single precision floats.
The iterative multiplication tends to create cumulative errors for small
angles. Chris's method is by far more accurate.
 
There is a workaround possible by dealing with residuals of the complex
multiplication, more precisely the additions and subtractions in the
complex multiplication are the root of the errors. This workaround
increases accuracy but not to the level of the transcendental functions.
 
 
Marcel
Manfred <noname@invalid.add>: Apr 12 02:57PM +0200

On 4/11/2017 1:09 AM, Chris M. Thomasson wrote:
> std::cos(angle_base + angle) * radius,
> std::sin(angle_base + angle) * radius
> };
 
The pair of cos+sin calculations can be efficiently combined in a single
call on platforms that allow that. This makes use of the power of
numerical algorithms to calculate both values simultaneously in
significantly less time than the two separate calls.
I think it is a pity that sincos() was not included (was it considered?)
in the standard, anyway these are the options I know of:
- GNU glibc exports a sincos() function
- on x87 a FSINCOS instruction is built in the FPU which can easily be
used on 32-bit x86
- on x86_64 it requires more work, but still the core numerical
algorithm (e.g. a Newton expansion) can be easily adapted to this.
Marcel Mueller <news.5.maazl@spamgourmet.org>: Apr 12 09:09AM +0200

On 11.04.17 01.09, Chris M. Thomasson wrote:
> std::cos(angle_base + angle) * radius,
> std::sin(angle_base + angle) * radius
> };
 
Depending on your portability and speed requirements you might want to
calculate sin and cos at once. Unfortunately sincos() never made it in
the standard. However,
ct_complex c = std::polar(radius, angle_base + angle);
might do the job for you if it is appropriately optimized at your target
platform. At least it is more readable.
 
 
> }
 
> // gain the average error sum... ;^o
> return avg_err / n;
 
Is there a reason why you calculated the /absolute/ floating point error
rather than the RMS value as usual?
 
 
> Can you run this? BTW, what about possibly reducing the average error
> ratio? Any ideas?
 
Hmm, there is probably no much you can do to reduce floating point
errors except if you do your own calculations at higher precision. E.g.
the complex angle calculations could be done with integer arithmetic.
This gives you about 63 significant bits rather than 52 bits of IEEE 754
double. Since you are operating at the unit circle fixed point is no big
deal. But reasonably implementing the transcendental functions is no fun
either.
 
 
Marcel
"Chris M. Thomasson" <invalid@invalid.invalid>: Apr 13 10:25AM -0700

On 4/12/2017 12:09 AM, Marcel Mueller wrote:
> ct_complex c = std::polar(radius, angle_base + angle);
> might do the job for you if it is appropriately optimized at your target
> platform. At least it is more readable.
 
Agreed. std::polar is easier to read, and can be faster.
 
 
>> return avg_err / n;
 
> Is there a reason why you calculated the /absolute/ floating point error
> rather than the RMS value as usual?
 
Not really. I just wanted to get a quick idea about the errors that are
going on.
 
> double. Since you are operating at the unit circle fixed point is no big
> deal. But reasonably implementing the transcendental functions is no fun
> either.
 
Right. I am going to use a high-precision floating point library in
later versions. 64-bit doubles are pretty good for this experimental phase.
You received this digest because you're subscribed to updates for this group. You can change your settings on the group membership page.
To unsubscribe from this group and stop receiving emails from it send an email to comp.lang.c+++unsubscribe@googlegroups.com.
Posted by at

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)