Cela ne répond pas à l'origine de poser la question, mais si vous êtes un programmeur qui travaille avec des problèmes similaires, cette réponse pourrait vous aider.
Je ne vois vraiment pas où la perception de la difficulté. Fournir strictes de la norme IEEE-754 binary64 sémantique, tout en étant limité à 80387 mathématiques à virgule flottante, et de la retenue de 80 bits de long double calcul, semble suivre bien spécifié C99 les règles de casting avec GCC-4.6.3 et clang-3.0 (basé sur LLVM 3.0).
Edité pour ajouter: Pourtant, Pascal Cuoq est correcte: ni gcc 4.6.3 ou clang-llvm-3.0 que l'application effective de ces règles correctement pour '387 mathématiques à virgule flottante. Si les options du compilateur, les règles sont correctement appliquées à des expressions évaluées au moment de la compilation, mais pas pour le moment de l'exécution des expressions. Il existe des solutions de contournement, énumérés après le saut ci-dessous.
Je ne simulation de dynamique moléculaire de code, et je suis très familier avec la répétabilité et la prévisibilité des exigences et aussi avec le désir de conserver un maximum de précision disponibles lorsque cela est possible, donc, je ne demande que je sais de quoi je parle ici. Cette réponse devrait montrer que les outils existent et sont simples à utiliser; les problèmes se posent de ne pas être au courant ou non à l'aide de ces outils.
(Un exemple préféré que j'aime, est le Kahan algorithme de sommation. Avec le C99 et bon casting (ajout de jette, par exemple Wikipédia exemple de code), pas de trucs ou extra variables temporaires sont nécessaires à tous. La mise en œuvre quel que soit le compilateur niveau d'optimisation, y compris à l' -O3 et -Ofast.)
C99 stipule de façon explicite (p. ex. dans 5.4.2.2) que le casting et l'affectation à la fois de supprimer toutes supplément de portée et de précision. Cela signifie que vous pouvez utiliser long double de l'arithmétique par la définition de vos variables temporaires utilisés pendant le calcul d' long double, aussi le moulage de votre entrée de variables de ce type, chaque fois qu'une norme IEEE-754 binary64 est nécessaire, il suffit de lancer pour un double.
Sur '387, le cast génère une affectation et une charge sur les deux compilateurs; ce n'correctement rond de 80 bits de la valeur à la norme IEEE-754 binary64. Ce coût est très raisonnable à mon avis. L'heure exacte dépend de l'architecture et du code environnant; en général, il est et peut être entrelacé avec d'autres code de réduire les coûts pour négligeable les niveaux. Lorsque MMX, SSE ou AVX sont disponibles, leurs registres sont séparés de l'80 bits 80387 registres, et le casting est généralement assuré par le déplacement de la valeur vers le MMX/SSE/AVX registre.
(Je préfère le code de production en vue d'utiliser une variable de type point, disons tempdouble ou, pour les variables temporaires, de sorte qu'il peut être défini soit double ou long double en fonction de l'architecture et de la vitesse/précision compromis désiré.)
En un mot:
Ne pas assumer (expression) est double de précision, tout simplement parce que toutes les variables et les constantes littérales sont. L'écrire comme (double)(expression) si vous souhaitez que le résultat en double de précision.
Cela s'applique à des expressions composées, aussi, et peut parfois conduire à manier les expressions avec de nombreux niveaux de jette.
Si vous avez expr1 et expr2 que vous souhaitez calculer à 80 bits de précision, mais aussi besoin que le produit de chaque arrondi à 64 bits, utilisez
long double expr1;
long double expr2;
double product = (double)(expr1) * (double)(expr2);
Remarque, product est calculée comme le produit de deux valeurs de 64 bits; pas calculée à 80 bits de précision, puis arrondis à la baisse. Calculer le produit à 80 bits de précision, puis en arrondissant, serait
double other = expr1 * expr2;
ou, en ajoutant descriptif jette que vous dire exactement ce qui se passe,
double other = (double)((long double)(expr1) * (long double)(expr2));
Il devrait être évident que l' product et other diffèrent souvent.
Le C99 règles de conversion sont juste un autre outil, vous devez apprendre à manier, si vous travaillez avec un mélange de 32-bit/64-bit/80 bits/128 bits des valeurs à virgule flottante. Vraiment, vous rencontre exactement les mêmes problèmes si vous mélangez binary32 et binary64 flotteurs (float et double sur la plupart des architectures)!
Peut-être la réécriture de Pascal Cuoq du code relatif à l'exploration, à appliquer correctement les règles de casting, rend cela plus clair?
#include <stdio.h>
#define TEST(eq) printf("%-56s%s\n", "" # eq ":", (eq) ? "true" : "false")
int main(void)
{
double d = 1.0 / 10.0;
long double ld = 1.0L / 10.0L;
printf("sizeof (double) = %d\n", (int)sizeof (double));
printf("sizeof (long double) == %d\n", (int)sizeof (long double));
printf("\nExpect true:\n");
TEST(d == (double)(0.1));
TEST(ld == (long double)(0.1L));
TEST(d == (double)(1.0 / 10.0));
TEST(ld == (long double)(1.0L / 10.0L));
TEST(d == (double)(ld));
TEST((double)(1.0L/10.0L) == (double)(0.1));
TEST((long double)(1.0L/10.0L) == (long double)(0.1L));
printf("\nExpect false:\n");
TEST(d == ld);
TEST((long double)(d) == ld);
TEST(d == 0.1L);
TEST(ld == 0.1);
TEST(d == (long double)(1.0L / 10.0L));
TEST(ld == (double)(1.0L / 10.0));
return 0;
}
La sortie, à la fois avec GCC et clang, est
sizeof (double) = 8
sizeof (long double) == 12
Expect true:
d == (double)(0.1): true
ld == (long double)(0.1L): true
d == (double)(1.0 / 10.0): true
ld == (long double)(1.0L / 10.0L): true
d == (double)(ld): true
(double)(1.0L/10.0L) == (double)(0.1): true
(long double)(1.0L/10.0L) == (long double)(0.1L): true
Expect false:
d == ld: false
(long double)(d) == ld: false
d == 0.1L: false
ld == 0.1: false
d == (long double)(1.0L / 10.0L): false
ld == (double)(1.0L / 10.0): false
sauf que les versions récentes de GCC promouvoir le côté droit de l' ld == 0.1 à long double première (c'est à dire à l' ld == 0.1L), rendement true, et que la SSE/AVX, long double est de 128 bits.
Pour le pur '387 tests, j'ai utilisé
gcc -W -Wall -m32 -mfpmath=387 -mno-sse ... test.c -o test
clang -W -Wall -m32 -mfpmath=387 -mno-sse ... test.c -o test
avec l'optimisation des différents drapeau combinaisons ..., y compris l' -fomit-frame-pointer, -O0, -O1, -O2, -O3, et -Os.
À l'aide de tous les autres drapeaux ou C99 compilateurs aboutissent au même résultat, à l'exception de long double de la taille (et ld == 1.0 pour les actuelles versions de GCC). Si vous rencontrez des différences, je serais très reconnaissant à entendre parler d'eux; j'ai peut-être besoin d'avertir mes utilisateurs de ces compilateurs (des versions de compilateur). Notez que Microsoft ne prend pas en charge le C99, de sorte qu'ils sont complètement inintéressant pour moi.
Pascal Cuoq apporte un intéressant problème dans le commentaire de la chaîne ci-dessous, que je n'ai pas reconnaître immédiatement.
Lors de l'évaluation d'une expression, à la fois GCC et clang avec -mfpmath=387 spécifier que toutes les expressions sont évaluées à l'aide de 80 bits de précision. Cela conduit par exemple
7491907632491941888 = 0x1.9fe2693112e14p+62 = 110011111111000100110100100110001000100101110000101000000000000
5698883734965350400 = 0x1.3c5a02407b71cp+62 = 100111100010110100000001001000000011110110111000111000000000000
7491907632491941888 * 5698883734965350400 = 42695510550671093541385598890357555200 = 100000000111101101101100110001101000010100100001011110111111111111110011000111000001011101010101100011000000000000000000000000
produisant des résultats incorrects, parce que la chaîne de ceux dans le milieu de la binaire résultat est tout simplement à la différence entre 53 et 64 bits mantissas (64 et 80 bits des nombres à virgule flottante, respectivement). Ainsi, alors que le résultat attendu est
42695510550671088819251326462451515392 = 0x1.00f6d98d0a42fp+125 = 100000000111101101101100110001101000010100100001011110000000000000000000000000000000000000000000000000000000000000000000000000
le résultat obtenu avec seulement -std=c99 -m32 -mno-sse -mfpmath=387 est
42695510550671098263984292201741942784 = 0x1.00f6d98d0a43p+125 = 100000000111101101101100110001101000010100100001100000000000000000000000000000000000000000000000000000000000000000000000000000
En théorie, vous devriez être en mesure de dire à gcc et clang pour faire respecter le bon C99 règles d'arrondi à l'aide d'options
-std=c99 -m32 -mno-sse -mfpmath=387 -ffloat-store -fexcess-precision=standard
Cependant, ceci n'affecte que les expressions le compilateur optimise, et ne semble pas résoudre le 387 manipulation. Si vous utilisez par exemple clang -O1 -std=c99 -m32 -mno-sse -mfpmath=387 -ffloat-store -fexcess-precision=standard test.c -o test && ./test avec test.c cours de Pascal Cuoq l'exemple du programme, vous obtiendrez le résultat correct conformément à la norme IEEE-754 règles, mais seulement parce que le compilateur optimise loin de l'expression, et non pas à l'aide de l'387.
Tout simplement, au lieu de l'informatique
(double)d1 * (double)d2
gcc et clang dire le '387 à calculer
(double)((long double)d1 * (long double)d2)
C'est, en effet, je crois que c'est un compilateur bug affectant à la fois la gcc 4.6.3 et clang-llvm-3.0, et facile à reproduire. (Pascal Cuoq points qu' FLT_EVAL_METHOD=2 moyen des opérations sur le double de la précision des arguments est toujours fait à précision étendue, mais je ne vois pas de toute saine raison, en plus d'avoir à réécrire des parties d' libm sur le '387 -- faire qu'en C99 et compte tenu de la norme IEEE-754 règles sont réalisables par le matériel! Après tout, le correct fonctionnement est facilement réalisable par le compilateur, en modifiant les 387 mot de contrôle pour correspondre à la précision de l'expression. Et, étant donné les options du compilateur qui doit forcer ce comportement -- -std=c99 -ffloat-store -fexcess-precision=standard -- aucun sens s' FLT_EVAL_METHOD=2 comportement est réellement désiré, il n'y a pas de rétro-compatibilité des questions, soit.) Il est important de noter que, compte tenu de la bonne drapeaux du compilateur, les expressions évaluées au moment de la compilation ne évaluée correctement, et que seules les expressions évaluées au moment de l'exécution obtenir des résultats incorrects.
La solution la plus simple, et le portable, est d'utiliser fesetround(FE_TOWARDZERO) (à partir de fenv.h) à l'année, tous les résultats vers zéro.
Dans certains cas, l'arrondi vers zéro peut aider avec la prévisibilité et les cas pathologiques. En particulier, pour des intervalles comme x = [0,1), arrondi vers zéro signifie que la limite supérieure n'est jamais atteint par le biais de l'arrondissement; important si vous évaluer par exemple par morceaux splines.
Pour les autres modes d'arrondi, vous avez besoin de contrôler l'387 directement au matériel.
Vous pouvez utiliser __FPU_SETCW() de #include <fpu_control.h>, ou ouvrez-code. Par exemple, precision.c:
#include <stdlib.h>
#include <stdio.h>
#include <limits.h>
#define FP387_NEAREST 0x0000
#define FP387_ZERO 0x0C00
#define FP387_UP 0x0800
#define FP387_DOWN 0x0400
#define FP387_SINGLE 0x0000
#define FP387_DOUBLE 0x0200
#define FP387_EXTENDED 0x0300
static inline void fp387(const unsigned short control)
{
unsigned short cw = (control & 0x0F00) | 0x007f;
__asm__ volatile ("fldcw %0" : : "m" (*&cw));
}
const char *bits(const double value)
{
const unsigned char *const data = (const unsigned char *)&value;
static char buffer[CHAR_BIT * sizeof value + 1];
char *p = buffer;
size_t i = CHAR_BIT * sizeof value;
while (i-->0)
*(p++) = '0' + !!(data[i / CHAR_BIT] & (1U << (i % CHAR_BIT)));
*p = '\0';
return (const char *)buffer;
}
int main(int argc, char *argv[])
{
double d1, d2;
char dummy;
if (argc != 3) {
fprintf(stderr, "\nUsage: %s 7491907632491941888 5698883734965350400\n\n", argv[0]);
return EXIT_FAILURE;
}
if (sscanf(argv[1], " %lf %c", &d1, &dummy) != 1) {
fprintf(stderr, "%s: Not a number.\n", argv[1]);
return EXIT_FAILURE;
}
if (sscanf(argv[2], " %lf %c", &d2, &dummy) != 1) {
fprintf(stderr, "%s: Not a number.\n", argv[2]);
return EXIT_FAILURE;
}
printf("%s:\td1 = %.0f\n\t %s in binary\n", argv[1], d1, bits(d1));
printf("%s:\td2 = %.0f\n\t %s in binary\n", argv[2], d2, bits(d2));
printf("\nDefaults:\n");
printf("Product = %.0f\n\t %s in binary\n", d1 * d2, bits(d1 * d2));
printf("\nExtended precision, rounding to nearest integer:\n");
fp387(FP387_EXTENDED | FP387_NEAREST);
printf("Product = %.0f\n\t %s in binary\n", d1 * d2, bits(d1 * d2));
printf("\nDouble precision, rounding to nearest integer:\n");
fp387(FP387_DOUBLE | FP387_NEAREST);
printf("Product = %.0f\n\t %s in binary\n", d1 * d2, bits(d1 * d2));
printf("\nExtended precision, rounding to zero:\n");
fp387(FP387_EXTENDED | FP387_ZERO);
printf("Product = %.0f\n\t %s in binary\n", d1 * d2, bits(d1 * d2));
printf("\nDouble precision, rounding to zero:\n");
fp387(FP387_DOUBLE | FP387_ZERO);
printf("Product = %.0f\n\t %s in binary\n", d1 * d2, bits(d1 * d2));
return 0;
}
À l'aide de clang-llvm-3.0 pour compiler et exécuter, j'obtiens des résultats corrects,
clang -std=c99 -m32 -mno-sse -mfpmath=387 -O3 -W -Wall precision.c -o precision
./precision 7491907632491941888 5698883734965350400
7491907632491941888: d1 = 7491907632491941888
0100001111011001111111100010011010010011000100010010111000010100 in binary
5698883734965350400: d2 = 5698883734965350400
0100001111010011110001011010000000100100000001111011011100011100 in binary
Defaults:
Product = 42695510550671098263984292201741942784
0100011111000000000011110110110110011000110100001010010000110000 in binary
Extended precision, rounding to nearest integer:
Product = 42695510550671098263984292201741942784
0100011111000000000011110110110110011000110100001010010000110000 in binary
Double precision, rounding to nearest integer:
Product = 42695510550671088819251326462451515392
0100011111000000000011110110110110011000110100001010010000101111 in binary
Extended precision, rounding to zero:
Product = 42695510550671088819251326462451515392
0100011111000000000011110110110110011000110100001010010000101111 in binary
Double precision, rounding to zero:
Product = 42695510550671088819251326462451515392
0100011111000000000011110110110110011000110100001010010000101111 in binary
En d'autres termes, vous pouvez contourner les problèmes de compilateur en utilisant fp387() définir la précision et le mode d'arrondi.
L'inconvénient est que certaines bibliothèques de mathématiques (libm.a, libm.so) peut être écrit avec l'hypothèse que les résultats intermédiaires sont toujours calculés à 80 bits de précision. Au moins la bibliothèque C de GNU fpu_control.h sur x86_64 a le commentaire "libm exige de la précision étendue". Heureusement, vous pouvez prendre l '387 implémentations par exemple, d'une bibliothèque C de GNU, et les mettre en œuvre dans un fichier d'en-tête ou d'écrire un travail, libm, si vous avez besoin de l' math.h fonctionnalités; en fait, je pense que je pourrais être en mesure de les y aider.
